1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
|
#include <stdio.h>
#include <math.h>
#include "gmp.h"
#include "gmp-impl.h"
#include "mpfr.h"
/* #define DEBUG */
void print_double(x) mpfr_t x;
{
int e, i; double d;
d = mpfr_get_d(x);
e = (int) ceil(log(d)/log(2.0));
/* d <= 2^e */
e -= 53;
if (e>0) for (i=0;i<e;i++) d /= 2.0;
else for (i=0;i<-e;i++) d *= 2.0;
printf("%1.0f*2^(%d)",d,e);
}
void mpfr_sqrt(mpfr_ptr x, mpfr_srcptr a, unsigned char rnd_mode)
{
int p, q, err, i, e, n; mpfr_t t, u;
#ifdef DEBUG
printf("enter mpfr_sqrt, a=%1.20e, rnd=%d\n",mpfr_get_d(a), rnd_mode);
printf("a="); print_double(a); putchar('\n');
#endif
/* use Newton's iteration x[n+1] = 1/2*(x[n]+a/x[n]),
the error e[n] = x[n]-sqrt(a) satisfies e[n+1] <= e[n]/2/sqrt(a) */
if (SIGN(a)<0) { SET_NAN(x); return; }
e = EXP(a)/2; if (2*e<EXP(a)) e++;
#ifdef DEBUG
printf("e=%d\n",e);
#endif
/* now 2^(2*e-2) <= a <= 2^(2*e) i.e. 1/4 <= a/2^(2e) <= 1 */
q = p = PREC(x);
for (i=0; i<3; i++)
q = p + (int) ceil(log(4.0*ceil(log((double)q)/log(2.0))+2.0)/log(2.0));
err = q-p; /* the error is at most 2^err ulp */
q = ((q-1)/mp_bits_per_limb)*mp_bits_per_limb; /* adjust to entire limb */
mpfr_init(t); mpfr_init(u);
do {
q += mp_bits_per_limb;
#ifdef DEBUG
printf("prec=%d q=%d err=%d\n",p,q,err);
#endif
mpfr_set_prec(t, q, GMP_RNDU);
mpfr_set_prec(x, q, GMP_RNDU);
mpfr_set_prec(u, q, GMP_RNDU);
mpfr_set_ui(x, 1, GMP_RNDU);
EXP(x) += e;
n = (int) ceil(log((double) q)/log(2.0));
for (i=0; i<n; i++) {
mpfr_div(t, a, x, GMP_RNDU);
mpfr_add(u, x, t, GMP_RNDU);
mpfr_div_2exp(x, u, 1, GMP_RNDU);
#ifdef DEBUG
printf("t="); print_double(t); putchar('\n');
printf("u="); print_double(u); putchar('\n');
printf("x="); print_double(x); putchar('\n');
printf("i=%d t=%1.20e u=%1.20e x=%1.20e\n",i,mpfr_get_d(t),mpfr_get_d(u),
mpfr_get_d(x));
printf("t="); mpfr_print_raw(t); putchar('\n');
printf("u="); mpfr_print_raw(u); putchar('\n');
printf("x="); mpfr_print_raw(x); putchar('\n');
#endif
}
} while (mpfr_can_round(x, q-err, GMP_RNDU, rnd_mode, p)==-1);
mpfr_round(x, rnd_mode, p);
mpfr_clear(t); mpfr_clear(u);
}
|
