Frank Heckenbach a écrit:
Maurice Lombardi wrote:
OK. However, it might be worthwhile to cache pi (as long as
precision suffices), like I do with LnHalf in mpf_ln. (Quite a
similar situation, BTW: Also not a closed loop, though even a direct
recursion there.) (Same for pi in mpf_sin and \sqrt 2 in
mpf_arctan.)
done (inside mpf_pi which is called from several places)
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).
OK. I wonder if it would be better to store the value of `a' in the
final loop in an integer variable instead. The range will certainly
be enough (since it counts by 2 each turn). I haven't done timing
tests, but I guess that dividing by an integer should at least not
be slower than by an mpf_t that happens to represent an integer
value, should it?
done
in addition gmp 3 gives a simpler way do do the internal function GetExp
All is included in the attached gmp.pas.diff2.diff to be applied
to the original gmp.pas
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-14 19:42:02.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,9 @@
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';
+function mpf_get_d_2exp (var exp: MedInt; protected var Src: mpf_t): Real; external name '__gmpf_get_d_2exp';
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';
@@ -417,8 +420,10 @@
procedure mpf_exp (var Dest: mpf_t; protected var Src: mpf_t);
procedure mpf_ln (var Dest: mpf_t; protected var Src: mpf_t);
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_arctan (var Dest: mpf_t; protected var Src: mpf_t);
+procedure mpf_pi (var Dest: mpf_t);
+procedure mpf_sin (var Dest: mpf_t; protected var Src: mpf_t);
+procedure mpf_cos (var Dest: mpf_t; protected var Src: mpf_t);
implementation
@@ -474,11 +479,19 @@
mpf_clear (Temp)
end;
+{$ifdef HAVE_GMP2}
function GetExp (protected var x: mpf_t) = Exp: mp_exp_t; attribute (inline);
{ @@ This is a kludge, but how to get the exponent (of base 2) in a better way? }
begin
Dispose (mpf_get_str (nil, Exp, 2, 0, x))
end;
+{$else}
+function GetExp (protected var x: mpf_t) = Exp: MedInt; attribute (inline);
+ var Dummy: Real;
+begin
+ Dummy := mpf_get_d_2exp (Exp, x);
+end;
+{$endif}
procedure mpf_exp (var Dest: mpf_t; protected var Src: mpf_t);
{ $$ \exp x = \sum_{n = 0}^{\infty} \frac{x^n}{n!} $$
@@ -546,7 +559,8 @@
mpf_init2 (Half, Precision);
mpf_set_d (Half, 0.5);
mpf_ln (LnHalf, Half);
- mpf_clear (Half)
+ mpf_clear (Half);
+ LnHalfInited := True;
end;
mpf_set (Dest, LnHalf);
mpf_mul_ui (Dest, Dest, Abs (Exp));
@@ -587,50 +601,166 @@
mpf_clear (Temp)
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)
+procedure mpf_arctan (var Dest: mpf_t; protected var Src: mpf_t);
+{ $$\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;
+ SqrtTwo: mpf_t; attribute (static);
+ SqrtTwoInited: boolean = False; attribute (static);
+begin
+ Precision := mpf_get_prec (Dest);
+ mpf_init2 (xx, Precision);
+ mpf_init2 (mx2, Precision);
+ mpf_init2 (a, Precision);
+ mpf_init2 (b, Precision);
+ mpf_abs (xx, Src);
+ if not SqrtTwoInited or (mpf_get_prec (SqrtTwo) < Precision) then begin
+ if SqrtTwoInited then mpf_clear (SqrtTwo);
+ mpf_init2 (SqrtTwo, Precision);
+ mpf_sqrt_ui (SqrtTwo, 2);
+ SqrtTwoInited := true;
+ end;
+ mpf_add_ui (a, SqrtTwo, 1);
+ mpf_sub_ui (b, SqrtTwo, 1);
+ if mpf_cmp (xx, a) > 0 then begin
+ mpf_pi (Dest);
+ mpf_div_2exp (Dest, Dest, 1);
+ mpf_ui_div (xx, 1, xx);
+ mpf_neg (xx, xx)
+ end else if mpf_cmp (xx, b) > 0 then begin
+ mpf_pi (Dest);
+ mpf_div_2exp (Dest, Dest, 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 (Dest, 0);
+ mpf_mul (mx2, xx, xx);
+ mpf_neg (mx2, mx2);
+ mpf_add (Dest, Dest, xx);
+ n := 1;
+ repeat
+ mpf_mul (xx, xx, mx2);
+ mpf_div_ui (a, xx, 2*n+1);
+ mpf_set (b, Dest);
+ mpf_add (Dest, Dest, a);
+ Inc (n)
+ until mpf_eq (b, Dest, Precision) <> 0;
+ if mpf_sgn (Src) < 0 then
+ mpf_neg (Dest, Dest);
+ mpf_clear (xx);
+ mpf_clear (mx2);
+ mpf_clear (a);
+ mpf_clear (b);
end;
-procedure mpf_pi (var c: mpf_t);
+procedure mpf_pi (var Dest: mpf_t);
{ 4 arctan 1/5 - arctan 1/239 = pi/4 }
-var b: mpf_t;
+var
+ b: mpf_t;
+ Precision: MedCard;
+ Pi: mpf_t; attribute (static);
+ PiInited: boolean = False; attribute (static);
+begin
+ Precision := mpf_get_prec (Dest);
+ if not PiInited or (mpf_get_prec (Pi) < Precision) then begin
+ if PiInited then mpf_clear (Pi);
+ mpf_init2 (Pi, Precision);
+ mpf_set_ui (Pi, 1);
+ mpf_div_ui (Pi, Pi, 5);
+ mpf_arctan (Pi, Pi);
+ mpf_mul_ui (Pi, Pi, 4);
+ mpf_init2 (b, Precision);
+ mpf_set_ui (b, 1);
+ mpf_div_ui (b, b, 239);
+ mpf_arctan (b, b);
+ mpf_sub (Pi, Pi, b);
+ mpf_mul_ui (Pi, Pi, 4);
+ mpf_clear (b);
+ PiInited := true;
+ end;
+ mpf_set (Dest, Pi);
+end;
+
+procedure mpf_sin (var Dest: mpf_t; protected var Src: 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, n: MedCard;
+ sign: integer;
+ a, b, z, xx, c0: mpf_t;
+begin
+ Precision := mpf_get_prec (Dest);
+ mpf_init2 (a, Precision);
+ mpf_init2 (b, Precision);
+ mpf_init2 (z, Precision);
+ mpf_init2 (xx, Precision);
+ mpf_init2 (c0, Precision);
+ sign := mpf_sgn(Src);
+ mpf_abs (xx, Src);
+ 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);
+{$ifdef HAVE_GMP2}
+ quadrant := round (mpf_get_d (b));
+{$else}
+ quadrant := mpf_get_ui (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);
+ n := 1;
+ mpf_set_ui (b, 1);
+ mpf_set_ui (Dest, 1);
+ repeat
+ Inc (n);
+ mpf_div_ui (b, b, n);
+ Inc (n);
+ mpf_div_ui (b, b, n);
+ mpf_mul (b, b, z);
+ mpf_set (c0, Dest);
+ mpf_add (Dest, Dest, b);
+ until mpf_eq (c0, Dest, Precision) <> 0;
+ mpf_mul (Dest, Dest, xx);
+ if sign < 0 then
+ mpf_neg (Dest, Dest);
+ 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_set_ui (c, 1);
- mpf_div_ui (c, c, 5);
- mpf_arctan (c, c);
- mpf_mul_ui (c, c, 4);
- mpf_init2 (b, mpf_get_prec (c));
- mpf_set_ui (b, 1);
- mpf_div_ui (b, b, 239);
- mpf_arctan (b, b);
- mpf_sub (c, c, b);
- mpf_mul_ui (c, c, 4);
- mpf_clear (b)
+ 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.
Frank Heckenbach
Maurice Lombardi