Hi, all While using gmp for my own work gmp 4.1.2 for djgpp I have found small bugs in the interface gmp.pas which cause compiling errors when compiling for gmp 4 (it compiles fine with gmp 3) The diff file is attached: gmp.pas.minimal.diff While working with it I have found an other error in the provided arctan(x) function. It works only for |x| <= 1, diverges above. For |x|>1 the simplest way is to use arctan(x)+arctan(1/x)=pi/2 but the convergence of the corresponding series is exceedingly slow at (near) |x|=1 I have programmed the same trick as in the QLib CEPHES library (which uses a home made 100 digits format instead of gmp), namely - use the previous formula to fold x above x=tan(3pi/8) = sqrt(2)+1 - use the formula arctan(x)=pi/4 +arctan((x-1)/(x+1)) to fold x comprised between tan(pi/8)=sqrt(2)-1 and tan(3pi/8)=sqrt(2)+1 The series has to be computed only for |x|<=0.414 where the speed of convergence is reasonable. In the CEPHES library Moshier uses a continuous fraction development instead of the series. I suppose this is a tradeoff speed/accuracy (not sure): the continuous fraction has only positive instead of alternating terms, which minimises rounding errors, but it is programmed with a fixed number of terms, of the order of the maximum number of terms in the series in the tests I have done, forbidding to benefit of the faster convergence for small |x|. I thus kept the series. Notice that pi included in these formulas is computed by a double call to arctan ! This is not a closed loop because pi needs only x=1/5 and x=1/239 below the lowest threshold. (Moshier uses a constant for pi because he has a fixed number of digits). I needed also the sin and cos functions. I have programmed the usual series development for sin, after reduction to the range [0,pi/2[, and cos(x)=sin(pi/2-x). All this is included in gmp.pas.full.diff, which includes also the diff contained in the minimal.diff (do not apply both). With this we have for gmp all the basic transcendental functions required for pascal (in the real type !) Hope this helps Maurice -- Maurice Lombardi Laboratoire de Spectrometrie Physique, Universite Joseph Fourier de Grenoble, BP87 38402 Saint Martin d'Heres Cedex FRANCE Tel: 33 (0)4 76 51 47 51 Fax: 33 (0)4 76 63 54 95 mailto:Maurice.Lombardi@ujf-grenoble.fr --- gmp.pas.orig 2004-02-07 18:27:00.000000000 +0000 +++ gmp.pas 2004-04-10 16:06:24.000000000 +0000 @@ -330,7 +330,7 @@ {$endif} { New declarations in GMP 3.x. @@ Mostly untested! } -{$ifdef HAVE_GMP3} +{$if not defined (HAVE_GMP2)} { Available random number generation algorithms. } type @@ -359,8 +359,8 @@ end; procedure gmp_randinit (var State: gmp_randstate_t; Alg: gmp_randalg_t; ...); external name '__gmp_randinit'; -procedure gmp_randinit_lc (var State: gmp_randstate_t; a: {$ifdef HAVE_GMP4} protected var {$endif} mpz_t; c: MedCard; m: {$ifdef HAVE_GMP4} protected var {$endif} mpz_t); external name '__gmp_randinit_lc'; -procedure gmp_randinit_lc_2exp (var State: gmp_randstate_t; a: {$ifdef HAVE_GMP4} protected var {$endif} mpz_t; c: MedCard; M2Exp: MedCard); external name '__gmp_randinit_lc_2exp'; +procedure gmp_randinit_lc (var State: gmp_randstate_t; {$ifdef HAVE_GMP4} protected var {$endif} a: mpz_t; c: MedCard; {$ifdef HAVE_GMP4} protected var {$endif} m: mpz_t); external name '__gmp_randinit_lc'; +procedure gmp_randinit_lc_2exp (var State: gmp_randstate_t; {$ifdef HAVE_GMP4} protected var {$endif} a: mpz_t; c: MedCard; M2Exp: MedCard); external name '__gmp_randinit_lc_2exp'; procedure gmp_randseed (var State: gmp_randstate_t; Seed: mpz_t); external name '__gmp_randseed'; procedure gmp_randseed_ui (var State: gmp_randstate_t; Seed: MedCard); external name '__gmp_randseed_ui'; procedure gmp_randclear (var State: gmp_randstate_t); external name '__gmp_randclear'; --- gmp.pas.orig 2004-02-07 18:27:00.000000000 +0000 +++ gmp.pas 2004-04-10 17:17:52.000000000 +0000 @@ -330,7 +330,7 @@ {$endif} { New declarations in GMP 3.x. @@ Mostly untested! } -{$ifdef HAVE_GMP3} +{$if not defined (HAVE_GMP2)} { Available random number generation algorithms. } type @@ -359,8 +359,8 @@ end; procedure gmp_randinit (var State: gmp_randstate_t; Alg: gmp_randalg_t; ...); external name '__gmp_randinit'; -procedure gmp_randinit_lc (var State: gmp_randstate_t; a: {$ifdef HAVE_GMP4} protected var {$endif} mpz_t; c: MedCard; m: {$ifdef HAVE_GMP4} protected var {$endif} mpz_t); external name '__gmp_randinit_lc'; -procedure gmp_randinit_lc_2exp (var State: gmp_randstate_t; a: {$ifdef HAVE_GMP4} protected var {$endif} mpz_t; c: MedCard; M2Exp: MedCard); external name '__gmp_randinit_lc_2exp'; +procedure gmp_randinit_lc (var State: gmp_randstate_t; {$ifdef HAVE_GMP4} protected var {$endif} a: mpz_t; c: MedCard; {$ifdef HAVE_GMP4} protected var {$endif} m: mpz_t); external name '__gmp_randinit_lc'; +procedure gmp_randinit_lc_2exp (var State: gmp_randstate_t; {$ifdef HAVE_GMP4} protected var {$endif} a: mpz_t; c: MedCard; M2Exp: MedCard); external name '__gmp_randinit_lc_2exp'; procedure gmp_randseed (var State: gmp_randstate_t; Seed: mpz_t); external name '__gmp_randseed'; procedure gmp_randseed_ui (var State: gmp_randstate_t; Seed: MedCard); external name '__gmp_randseed_ui'; procedure gmp_randclear (var State: gmp_randstate_t); external name '__gmp_randclear'; @@ -395,6 +395,8 @@ procedure mpf_ceil (var Dest: mpf_t; protected var Src: mpf_t); external name '__gmpf_ceil'; procedure mpf_floor (var Dest: mpf_t; protected var Src: mpf_t); external name '__gmpf_floor'; +function mpf_get_si (protected var Src: mpf_t): MedInt; external name '__gmpf_get_si'; +function mpf_get_ui (protected var Src: mpf_t): MedCard; external name '__gmpf_get_ui'; procedure mpf_pow_ui (var Dest: mpf_t; protected var Src1: mpf_t; Src2: MedCard); external name '__gmpf_pow_ui'; procedure mpf_trunc (var Dest: mpf_t; protected var Src: mpf_t); external name '__gmpf_trunc'; procedure mpf_urandomb (ROP: mpf_t; var State: gmp_randstate_t; n: MedCard); external name '__gmpf_urandomb'; @@ -419,6 +421,8 @@ procedure mpf_pow (var Dest: mpf_t; protected var Src1, Src2: mpf_t); procedure mpf_arctan (var c: mpf_t; protected var x: mpf_t); procedure mpf_pi (var c: mpf_t); +procedure mpf_sin (var c: mpf_t; protected var x: mpf_t); +procedure mpf_cos (var Dest: mpf_t; protected var Src: mpf_t); implementation @@ -588,32 +592,57 @@ end; procedure mpf_arctan (var c: mpf_t; protected var x: mpf_t); -{ $$\arctan x = \sum_{n=0}^{\infty} (-1)^n \frac{x^{2n+1}}{2n+1}$$ } -var - p, mx2, c0, s: mpf_t; - Precision, n: MedCard; -begin - Precision := mpf_get_prec (c); - mpf_init2 (p, Precision); - mpf_set (p, x); - mpf_init2 (mx2, Precision); - mpf_mul (mx2, x, x); - mpf_neg (mx2, mx2); - mpf_init2 (c0, Precision); - mpf_init2 (s, Precision); - mpf_set (c, x); - n := 1; - repeat - mpf_mul (p, p, mx2); - mpf_div_ui (s, p, 2 * n + 1); - mpf_set (c0, c); - mpf_add (c, c, s); - Inc (n) - until mpf_eq (c0, c, Precision) <> 0; - mpf_clear (s); - mpf_clear (c0); - mpf_clear (mx2); - mpf_clear (p) +{ $$\arctan x = \sum_{n=0}^{\infty} (-1)^n \frac{x^{2n+1}}{2n+1}$$ + after double range reduction + around $\tan(3\pi/8)=\sqrt(2)+1$ with $\arctan x = \pi/2+\arctan(-1/x)$ + around $\tan(\pi/8)=\sqrt(2)-1$ with $\arctan x = \pi/4+\arctan((x-1)/(x+1))$ } +var + Precision, n: MedCard; + xx, mx2, a, b, aa, bb: mpf_t; +begin + Precision := mpf_get_prec (c); + mpf_init2 (xx, Precision); + mpf_init2 (mx2, Precision); + mpf_init2 (a, Precision); + mpf_init2 (b, Precision); + mpf_abs (xx, x); + mpf_sqrt_ui (a, 2); + mpf_sub_ui (b, a, 1); + mpf_add_ui (a, a, 1); + if mpf_cmp (xx, a) > 0 then begin + mpf_pi (c); + mpf_div_2exp (c, c, 1); + mpf_ui_div (xx, 1, xx); + mpf_neg (xx, xx) + end else if mpf_cmp (xx, b) > 0 then begin + mpf_pi (c); + mpf_div_2exp (c, c, 2); + mpf_init2 (aa, Precision); + mpf_init2 (bb, Precision); + mpf_sub_ui (aa, xx, 1); + mpf_add_ui (bb, xx, 1); + mpf_div (xx, aa, bb); + mpf_clear (aa); + mpf_clear (bb); + end else + mpf_set_ui (c, 0); + mpf_mul (mx2, xx, xx); + mpf_neg (mx2, mx2); + mpf_add (c, c, xx); + n := 1; + repeat + mpf_mul (xx, xx, mx2); + mpf_div_ui (a, xx, 2*n+1); + mpf_set (b, c); + mpf_add (c, c, a); + Inc (n) + until mpf_eq (b, c, Precision) <> 0; + if mpf_sgn (x) < 0 then + mpf_neg (c, c); + mpf_clear (xx); + mpf_clear (mx2); + mpf_clear (a); + mpf_clear (b); end; procedure mpf_pi (var c: mpf_t); @@ -633,4 +662,77 @@ mpf_clear (b) end; +procedure mpf_sin (var c: mpf_t; protected var x: mpf_t); +{ $$\sin x = \sum_{n=0}^{\infty} (-1)^n \frac{x^{2n+1}}{(2n+1)!}$$ + after reduction to range $[0,\pi/2[$ } +var + Precision, quadrant: MedCard; + sign: integer; + a, b, z, xx, c0: mpf_t; +begin + Precision := mpf_get_prec (c); + mpf_init2 (a, Precision); + mpf_init2 (b, Precision); + mpf_init2 (z, Precision); + mpf_init2 (xx, Precision); + mpf_init2 (c0, Precision); + sign := mpf_sgn(x); + mpf_abs (xx, x); + mpf_pi (z); + mpf_div_2exp (z, z, 1); + mpf_div (a, xx, z); + mpf_floor (xx, a); + if mpf_cmp_ui (xx, 4) >= 0 then begin + mpf_div_2exp (b, xx, 2); + mpf_floor (b, b); + mpf_mul_2exp (b, b, 2); + end else + mpf_set_ui (b, 0); + mpf_sub (b, xx, b); +{$if not defined (HAVE_GMP2)} + quadrant := mpf_get_ui (b); +{$else} + quadrant := round (mpf_get_d (b)); +{$endif} + mpf_sub (b, a, xx); + mpf_mul (xx, z, b); + if quadrant > 1 then + sign := -sign; + if (quadrant and 1) <> 0 then + mpf_sub (xx, z, xx); + mpf_mul (z, xx, xx); + mpf_neg (z, z); + mpf_set_ui (a, 1); + mpf_set_ui (b, 1); + mpf_set_ui (c, 1); + repeat + mpf_add_ui (a, a, 1); + mpf_div (b, b, a); + mpf_add_ui (a, a, 1); + mpf_div (b, b, a); + mpf_mul (b, b, z); + mpf_set (c0, c); + mpf_add (c, c, b); + until mpf_eq (c0, c, Precision) <> 0; + mpf_mul (c, c, xx); + if sign < 0 then + mpf_neg (c, c); + mpf_clear (a); + mpf_clear (b); + mpf_clear (z); + mpf_clear (xx); + mpf_clear (c0); +end; + +procedure mpf_cos (var Dest: mpf_t; protected var Src: mpf_t); +var Temp: mpf_t; +begin + mpf_init2 (Temp, mpf_get_prec (Dest)); + mpf_pi (Temp); + mpf_div_2exp (Temp, Temp, 1); + mpf_sub (Temp, Temp, Src); + mpf_sin (Dest, Temp); + mpf_clear (Temp); +end; + end.