libmpfr library interface and tests updated in kForth-32

[Krishna Myneni]

Having need again of the GNU MPFR library, I made some updates to the
kForth-32 interface to the library and its tests. Links to the files are
given below.

The kForth interface is libmpfr.4th, and the tests are in
libmpfr-test.4th. The interface has only changed via substituting older
structures with Forth 200x structures, and appears to work for versions
of the MPFR library >= 3.0.0.

Tests of the mpfr_xx words were incomplete and buggy for the
transcendental functions, and I have fixed these problems. The tests
only scratch the surface of the functions provided by MPFR, but they
provide a minimum set to start working with it.

https://github.com/mynenik/kForth-32/blob/master/forth-src/libs/gmp/libmpfr.4th

https://github.com/mynenik/kForth-32/blob/master/forth-src/libs/gmp/libmpfr-test.4th

--
Krishna Myneni

[Krishna Myneni]

On 3/18/23 16:02, Krishna Myneni wrote:
> Having need again of the GNU MPFR library, I made some updates to the
> kForth-32 interface to the library and its tests. Links to the files are
> given below.

...
>
> https://github.com/mynenik/kForth-32/blob/master/forth-src/libs/gmp/libmpfr.4th
>
> https://github.com/mynenik/kForth-32/blob/master/forth-src/libs/gmp/libmpfr-test.4th
>

I had to update the definition of MPFR. and various other files in the
forth-src/libs/gmp folder as well, so the entire set of files should be
refreshed if you want to try out the GMP and MPFR interfaces under
kForth-32.

--
KM

[minforth]

Impressive.

For smaller requirements Fabrice Bellard's libbf is a worthy alternative
https://bellard.org/libbf/

[Krishna Myneni]

On 3/20/23 08:08, minforth wrote:
> Krishna Myneni schrieb am Samstag, 18. März 2023 um 23:33:20 UTC+1:
>> On 3/18/23 16:02, Krishna Myneni wrote:
>>> Having need again of the GNU MPFR library, I made some updates to the
>>> kForth-32 interface to the library and its tests. Links to the files are
>>> given below.
>> ...

>

> For smaller requirements Fabrice Bellard's libbf is a worthy alternative
> https://bellard.org/libbf/

Thanks. I looked at the info on libbf. The benchmarks suggest that it
has some advantages over mpfr in execution speed when an extremely large
number of digits are required. For < 100 digits precision, their graph
suggest that mpfr performs quite a bit better.

I ran into a problem with double precision fp limits when I needed to
compute the spherical Bessel and Neumann functions, j_l(x) and n_l(x),
for arguments x >= 2000, for l of a few hundred. The functions
themselves have values in the 1e-4 to 1e-5 range; however, the
intermediate results of the calculation overflow the 11-bit exponent
(around 1e308). There's also a tremendous loss of accuracy with the
double precision routine for high l.

Using the double precision sph_bes_neu.4th, with MAX-L set to 5000,
j_l=340 (x = 2100) and n_l=340 (x = 2100) become -INF and NAN,
respectively. Yes, such large arguments for l and x are needed in
certain calculations. I considered using the double-double Forth
arithmetic library of Julian Noble's until I realized that the exponent
overflow in the intermediate quantities was the cause of the problem --
double double's have the same exponent range as ordinary double
precision fp numbers. The MPFR version of the spherical Bessel function
code (mpfr_sph_bes_neu.4th) has a 32-bit exponent range, and provides
full double precision accuracy even for l and x arguments up to several
thousand.

--
Krishna

https://github.com/mynenik/kForth-32/blob/master/forth-src/fsl/extras/sph_bes_neu.4th

https://github.com/mynenik/kForth-32/blob/master/forth-src/libs/gmp/mpfr_sph_bes_neu.4th

[minforth]

"Somewhere over the rainbow" there are 128-bit quadruple floating-point numbers
with 15-bit mantissa as specified by IEEE754. But apart from gcc libquadmath
there does not seem to be wide usage. I have only played with it in the past, so I
couldn't tell whether it could present an alternative to linking with such "huge hog libraries"
as gnu gmp and mpfr.

[minforth]

senile me .. exponent of course

[Krishna Myneni]

On 3/24/23 14:21, minforth wrote:
> minforth schrieb am Freitag, 24. März 2023 um 19:19:43 UTC+1:

...

>>> The MPFR version of the spherical Bessel function
>>> code (mpfr_sph_bes_neu.4th) has a 32-bit exponent range, and provides
>>> full double precision accuracy even for l and x arguments up to several
>>> thousand.
>>>
>>> --
>>> Krishna
>>>
>>> https://github.com/mynenik/kForth-32/blob/master/forth-src/fsl/extras/sph_bes_neu.4th
>>>
>>> https://github.com/mynenik/kForth-32/blob/master/forth-src/libs/gmp/mpfr_sph_bes_neu.4th
>> "Somewhere over the rainbow" there are 128-bit quadruple floating-point numbers
>> with 15-bit mantissa as specified by IEEE754. But apart from gcc libquadmath
>> there does not seem to be wide usage. I have only played with it in the past, so I
>> couldn't tell whether it could present an alternative to linking with such "huge hog libraries"
>> as gnu gmp and mpfr.
>
> senile me .. exponent of course

The 15-bit exponent of IEEE quad fp is likely to be sufficient, but to
achieve full double-precision over the range of x ~ few thousand and l ~
few thousand, my tests indicate that quad precision will not be
sufficient. In the MPFR version, I am currently using a 128-bit mantissa
(and a 32-bit exponent). Now, full double precision in j_l(x) and
n_l(x), over the desired range may not be needed for my application, but
I can think of situations where it may be important to start out with
full double precision for those values. libquadmath is important for
applications, but not a replacement for MPFR.

--
Krishna

[minforth]

I can understand. Long time ago I was bitten by a similar problem while
solving numeric differential equations for electric fields around high-voltage cables.
Luckily I could revive the old battle-horse: FORTRAN ;-)