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
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
|
// file kernel/x/c/sqrt.c: square root of extensible integers
/*-----------------------------------------------------------------------+
| Copyright 2005-2006, Michel Quercia (michel.quercia@prepas.org) |
| |
| This file is part of Numerix. Numerix is free software; you can |
| redistribute it and/or modify it under the terms of the GNU Lesser |
| General Public License as published by the Free Software Foundation; |
| either version 2.1 of the License, or (at your option) any later |
| version. |
| |
| The Numerix Library is distributed in the hope that it will be |
| useful, but WITHOUT ANY WARRANTY; without even the implied warranty |
| of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| Lesser General Public License for more details. |
| |
| You should have received a copy of the GNU Lesser General Public |
| License along with the GNU MP Library; see the file COPYING. If not, |
| write to the Free Software Foundation, Inc., 59 Temple Place - |
| Suite 330, Boston, MA 02111-1307, USA. |
+-----------------------------------------------------------------------+
| |
| Racine carre |
| |
+-----------------------------------------------------------------------*/
/*
entre :
a = entier extensible
_b = NULL ou pointeur sur un entier extensible
mode = long
sortie :
si mode & 3 = 0: b <- floor(sqrt(a))
si mode & 3 = 1: b <- floor(sqrt(a)+1/2)
si mode & 3 = 2: b <- ceil(sqrt(a))
si mode & 3 = 3: b <- ceil(sqrt(a)-1/2) (= floor(sqrt(a)+1/2))
si _b != NULL, *_b <- b
retourne b
erreur :
NEGATIVE_BASE si a < 0
*/
xint xx(private_sqrt)(xint *_b, xint a, long mode) {
long la = xx_lg(a), lb, n,p;
chiffre *x,r;
int inc_b;
xint b;
xx_push_roots_2(a,_b);
#ifdef caml_api
#define a __lr.a
#define _b __lr._b
#endif
if (xx_sgn(a)) xx(failwith)(NEGATIVE_BASE);
if (la == 0) {
b = xx(enlarge)(_b,0);
b->hd = 0;
xx_update_and_return(_b,b);
}
/* dcale a pour avoir une longueur paire et un chiffre de poids
fort compris entre BASE/16 et BASE/4 */
r = a->val[la-1];
if (r & (3*(BASE_2/2))) {
n = HW/2-1;
p = 1 - (la&1);
lb = (la+p+1)/2;
x = xn(alloc)(2*lb);
x[0] = 0;
x[la+p] = xn(shift_up)(a->val,la,x+p,2*n);
}
else {
for (n=0; (r & (3*(BASE_2/8))) == 0; n++, r<<=2);
p = la&1;
lb = (la+p)/2;
x = xn(alloc)(2*lb);
x[0] = 0;
xn(shift_up)(a->val,la,x+p,2*n);
}
/* b <- 2*floor(sqrt(a*2^(n+p*HW))) */
b = xx(enlarge)(_b,lb);
(2*lb > zimsqrt_lim) ? xn(modsqrt)(x,2*lb,b->val)
: xn(sqrt_n2)(x,2*lb,b->val);
/* arrondi selon le mode */
n += p*(HW/2)+1;
switch(mode & 3) {
case 0: inc_b = 0; break;
case 2: inc_b = (xn(cmp)(x,lb,x,0)); break;
default:inc_b = (n>1) ? (b->val[0] >> (n-1)) & 1
: (xn(cmp2)(x,lb,b->val,lb) > 0);
break;
}
/* libre la mmoire temporaire */
xn(free)(x);
/* dcale b et rectifie la longueur */
xn(shift_down)(b->val,lb,b->val,n);
if (inc_b) xn(inc1)(b->val,lb);
xx(make_head)(b,lb,0);
xx_update_and_return(_b,b);
#undef a
#undef _b
}
#if defined(caml_api) || defined(ocaml_api)
xint xx(sqrt) ( xint *_b, xint a) {return xx(private_sqrt)(_b, a,0);}
xint xx(gsqrt) (value mode, xint *_b, xint a) {return xx(private_sqrt)(_b, a,Round_val(mode));}
xint xx(f_sqrt) ( xint a) {return xx(private_sqrt)(xx_null,a,0);}
xint xx(f_gsqrt)(value mode, xint a) {return xx(private_sqrt)(xx_null,a,Round_val(mode));}
#endif /* api */
|
