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.
