libmpfr library interface and tests updated in kForth-32

Post by 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.

...

Post by Krishna Myneni
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

2023-03-20 13:08:32 UTC

Post by 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.

...

Post by Krishna Myneni
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

Impressive.

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

Krishna Myneni

2023-03-21 13:26:36 UTC

Post by 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.

...

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

2023-03-24 18:19:42 UTC

Post by 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.

...

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

"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

2023-03-24 19:21:12 UTC

Post by 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.

...

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

senile me .. exponent of course

Krishna Myneni

2023-03-27 12:39:25 UTC

...

Post by Krishna Myneni
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

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

2023-03-30 07:43:34 UTC

Post by Krishna Myneni
...

senile me .. exponent of course

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 ;-)

minforth

2024-05-27 13:46:26 UTC

Post by 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

I am completing my bigfloat library with gmp/mpfr. Compiled with -static
and
debugging information stripped off, the size of the executable grew
4-fold (!) from
192k (with libbf) to 783k (with mpfr) with otherwise identical
functionality.

Admittedly mpfr already includes some higher functions like
Bessels00000, Gamma and such,
but when you don't need them, mpfr is a heavy load.

I'll do some more tests incl. some speed tests. I was surprised to find
that the
mpfr example given in
https://www.mpfr.org/sample.html
seems to be slightly wrong, after comparing the Euler number with
Wolfram Alpha.
Accumulated rounding error.

I use a libbf (or mpfr) wrapper in the style of the standard
floating-point number
wordset. So my equivalent Forth program is much more compact than the
original:

\ EULER BigFloat Test

62 SET-BFPRECISION

: EULER ( n -- B: e ; compute euler number by series expansion )
B1 ( sum )
B1 ( inverse factorial )
1 DO
i bs/ b+ bup ( undrop )
LOOP bdrop ;

100 EULER
( Euler expansion 62 digits:)
CR B.
CR .( 2.7182818284590452353602874713526624977572470936999595749669676)

Run:

MinForth 3.6 (64 bit) (bf libbf)
# fload e.mf Euler expansion 62 digits:
2.7182818284590452353602874713526624977572470936999595749669676
2.7182818284590452353602874713526624977572470936999595749669676 ok
#

Krishna Myneni

2024-05-27 18:38:56 UTC

Is libbf a C library? Is your executable a test program?

Post by minforth
Admittedly mpfr already includes some higher functions like
Bessels00000, Gamma and such,
but when you don't need them, mpfr is a heavy load.

Special functions occur frequently in many of my programs, but once you
have basic arithmetic, power, and common trig functions (sin and cos),
many of the higher order functions can be implemented. There's even
Forth source within the FSL for it.

Post by minforth
I'll do some more tests incl. some speed tests. I was surprised to find
that the
mpfr example given in
https://www.mpfr.org/sample.html
seems to be slightly wrong, after comparing the Euler number with
Wolfram Alpha.
Accumulated rounding error.

Interesting.

Post by minforth
I use a libbf (or mpfr) wrapper in the style of the standard
floating-point number
wordset. So my equivalent Forth program is much more compact than the
\ EULER BigFloat Test
62 SET-BFPRECISION

So you set precision not in binary digits but in decimal digits.

Post by minforth
: EULER ( n -- B: e ; compute euler number by series expansion )
B1 ( sum )
B1 ( inverse factorial )
1 DO
i bs/ b+ bup ( undrop )
LOOP bdrop ;
100 EULER ( Euler expansion 62 digits:)
CR B.
CR .( 2.7182818284590452353602874713526624977572470936999595749669676)
MinForth 3.6 (64 bit) (bf libbf)
2.7182818284590452353602874713526624977572470936999595749669676
2.7182818284590452353602874713526624977572470936999595749669676 ok
#

Congrats!

--
Krishna

minforth

2024-05-27 20:12:15 UTC

Post by minforth
I am completing my bigfloat library with gmp/mpfr. Compiled with -static
and
debugging information stripped off, the size of the executable grew
4-fold (!) from
192k (with libbf) to 783k (with mpfr) with otherwise identical
functionality.

Is libbf a C library? Is your executable a test program?

Libbf is a C library, actually kind of "mpfr light".
See: https://www.bellard.org/libbf/

My 192k executable is a full Forth system (CLI) with core,
exceptions, float64 and bigfloat wordsets, plus some few utilities.

Post by minforth
62 SET-BFPRECISION

So you set precision not in binary digits but in decimal digits.

Of course you can also set the internal bit-width mpfr-prec (and
rounding mode) directly. But more often than that I am interested
in precise decimal digits. SET-BFPRECISION calculates and sets a
suitable bid-width for a given decimal precision. A convenience.

Post by minforth
: EULER ( n -- B: e ; compute euler number by series expansion )
B1 ( sum )
B1 ( inverse factorial )
1 DO
i bs/ b+ bup ( undrop )
LOOP bdrop ;

Congrats!

Thanks. Actually it is a tad shorter:
: EULER ( n -- B: e ) 1 s>b bdup FOR n bs/ b+ bup NEXT bdrop ;

minforth

2024-05-27 21:20:46 UTC