]>
Commit | Line | Data |
---|---|---|
04067002 UD |
1 | |
2 | /* | |
3 | * IBM Accurate Mathematical Library | |
4 | * written by International Business Machines Corp. | |
f4cf5f2d | 5 | * Copyright (C) 2001-2012 Free Software Foundation, Inc. |
04067002 UD |
6 | * |
7 | * This program is free software; you can redistribute it and/or modify | |
8 | * it under the terms of the GNU Lesser General Public License as published by | |
9 | * the Free Software Foundation; either version 2.1 of the License, or | |
10 | * (at your option) any later version. | |
11 | * | |
12 | * This program is distributed in the hope that it will be useful, | |
13 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | |
14 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
15 | * GNU Lesser General Public License for more details. | |
16 | * | |
17 | * You should have received a copy of the GNU Lesser General Public License | |
59ba27a6 | 18 | * along with this program; if not, see <http://www.gnu.org/licenses/>. |
04067002 UD |
19 | */ |
20 | /************************************************************************/ | |
21 | /* MODULE_NAME: mpa.c */ | |
22 | /* */ | |
23 | /* FUNCTIONS: */ | |
24 | /* mcr */ | |
25 | /* acr */ | |
26 | /* cr */ | |
27 | /* cpy */ | |
28 | /* cpymn */ | |
29 | /* norm */ | |
30 | /* denorm */ | |
31 | /* mp_dbl */ | |
32 | /* dbl_mp */ | |
33 | /* add_magnitudes */ | |
34 | /* sub_magnitudes */ | |
35 | /* add */ | |
36 | /* sub */ | |
37 | /* mul */ | |
38 | /* inv */ | |
39 | /* dvd */ | |
40 | /* */ | |
41 | /* Arithmetic functions for multiple precision numbers. */ | |
42 | /* Relative errors are bounded */ | |
43 | /************************************************************************/ | |
44 | ||
45 | ||
46 | #include "endian.h" | |
47 | #include "mpa.h" | |
48 | #include "mpa2.h" | |
49 | #include <sys/param.h> /* For MIN() */ | |
50 | /* mcr() compares the sizes of the mantissas of two multiple precision */ | |
51 | /* numbers. Mantissas are compared regardless of the signs of the */ | |
52 | /* numbers, even if x->d[0] or y->d[0] are zero. Exponents are also */ | |
53 | /* disregarded. */ | |
54 | static int mcr(const mp_no *x, const mp_no *y, int p) { | |
55 | long i; | |
56 | long p2 = p; | |
57 | for (i=1; i<=p2; i++) { | |
58 | if (X[i] == Y[i]) continue; | |
59 | else if (X[i] > Y[i]) return 1; | |
60 | else return -1; } | |
61 | return 0; | |
62 | } | |
63 | ||
64 | ||
65 | ||
66 | /* acr() compares the absolute values of two multiple precision numbers */ | |
67 | int __acr(const mp_no *x, const mp_no *y, int p) { | |
68 | long i; | |
69 | ||
70 | if (X[0] == ZERO) { | |
71 | if (Y[0] == ZERO) i= 0; | |
72 | else i=-1; | |
73 | } | |
74 | else if (Y[0] == ZERO) i= 1; | |
75 | else { | |
76 | if (EX > EY) i= 1; | |
77 | else if (EX < EY) i=-1; | |
78 | else i= mcr(x,y,p); | |
79 | } | |
80 | ||
81 | return i; | |
82 | } | |
83 | ||
84 | ||
85 | /* cr90 compares the values of two multiple precision numbers */ | |
86 | int __cr(const mp_no *x, const mp_no *y, int p) { | |
87 | int i; | |
88 | ||
89 | if (X[0] > Y[0]) i= 1; | |
90 | else if (X[0] < Y[0]) i=-1; | |
91 | else if (X[0] < ZERO ) i= __acr(y,x,p); | |
92 | else i= __acr(x,y,p); | |
93 | ||
94 | return i; | |
95 | } | |
96 | ||
97 | ||
98 | /* Copy a multiple precision number. Set *y=*x. x=y is permissible. */ | |
99 | void __cpy(const mp_no *x, mp_no *y, int p) { | |
100 | long i; | |
101 | ||
102 | EY = EX; | |
103 | for (i=0; i <= p; i++) Y[i] = X[i]; | |
104 | ||
105 | return; | |
106 | } | |
107 | ||
108 | ||
109 | /* Copy a multiple precision number x of precision m into a */ | |
110 | /* multiple precision number y of precision n. In case n>m, */ | |
111 | /* the digits of y beyond the m'th are set to zero. In case */ | |
112 | /* n<m, the digits of x beyond the n'th are ignored. */ | |
113 | /* x=y is permissible. */ | |
114 | ||
115 | void __cpymn(const mp_no *x, int m, mp_no *y, int n) { | |
116 | ||
117 | long i,k; | |
118 | long n2 = n; | |
119 | long m2 = m; | |
120 | ||
121 | EY = EX; k=MIN(m2,n2); | |
122 | for (i=0; i <= k; i++) Y[i] = X[i]; | |
123 | for ( ; i <= n2; i++) Y[i] = ZERO; | |
124 | ||
125 | return; | |
126 | } | |
127 | ||
128 | /* Convert a multiple precision number *x into a double precision */ | |
129 | /* number *y, normalized case (|x| >= 2**(-1022))) */ | |
130 | static void norm(const mp_no *x, double *y, int p) | |
131 | { | |
132 | #define R radixi.d | |
133 | long i; | |
134 | #if 0 | |
135 | int k; | |
136 | #endif | |
137 | double a,c,u,v,z[5]; | |
138 | if (p<5) { | |
139 | if (p==1) c = X[1]; | |
140 | else if (p==2) c = X[1] + R* X[2]; | |
141 | else if (p==3) c = X[1] + R*(X[2] + R* X[3]); | |
142 | else if (p==4) c =(X[1] + R* X[2]) + R*R*(X[3] + R*X[4]); | |
143 | } | |
144 | else { | |
145 | for (a=ONE, z[1]=X[1]; z[1] < TWO23; ) | |
146 | {a *= TWO; z[1] *= TWO; } | |
147 | ||
148 | for (i=2; i<5; i++) { | |
149 | z[i] = X[i]*a; | |
150 | u = (z[i] + CUTTER)-CUTTER; | |
151 | if (u > z[i]) u -= RADIX; | |
152 | z[i] -= u; | |
153 | z[i-1] += u*RADIXI; | |
154 | } | |
155 | ||
156 | u = (z[3] + TWO71) - TWO71; | |
157 | if (u > z[3]) u -= TWO19; | |
158 | v = z[3]-u; | |
159 | ||
160 | if (v == TWO18) { | |
161 | if (z[4] == ZERO) { | |
162 | for (i=5; i <= p; i++) { | |
163 | if (X[i] == ZERO) continue; | |
164 | else {z[3] += ONE; break; } | |
165 | } | |
166 | } | |
167 | else z[3] += ONE; | |
168 | } | |
169 | ||
170 | c = (z[1] + R *(z[2] + R * z[3]))/a; | |
171 | } | |
172 | ||
173 | c *= X[0]; | |
174 | ||
175 | for (i=1; i<EX; i++) c *= RADIX; | |
176 | for (i=1; i>EX; i--) c *= RADIXI; | |
177 | ||
178 | *y = c; | |
179 | return; | |
180 | #undef R | |
181 | } | |
182 | ||
183 | /* Convert a multiple precision number *x into a double precision */ | |
184 | /* number *y, denormalized case (|x| < 2**(-1022))) */ | |
185 | static void denorm(const mp_no *x, double *y, int p) | |
186 | { | |
187 | long i,k; | |
188 | long p2 = p; | |
189 | double c,u,z[5]; | |
190 | #if 0 | |
191 | double a,v; | |
192 | #endif | |
193 | ||
194 | #define R radixi.d | |
195 | if (EX<-44 || (EX==-44 && X[1]<TWO5)) | |
196 | { *y=ZERO; return; } | |
197 | ||
198 | if (p2==1) { | |
199 | if (EX==-42) {z[1]=X[1]+TWO10; z[2]=ZERO; z[3]=ZERO; k=3;} | |
200 | else if (EX==-43) {z[1]= TWO10; z[2]=X[1]; z[3]=ZERO; k=2;} | |
201 | else {z[1]= TWO10; z[2]=ZERO; z[3]=X[1]; k=1;} | |
202 | } | |
203 | else if (p2==2) { | |
204 | if (EX==-42) {z[1]=X[1]+TWO10; z[2]=X[2]; z[3]=ZERO; k=3;} | |
205 | else if (EX==-43) {z[1]= TWO10; z[2]=X[1]; z[3]=X[2]; k=2;} | |
206 | else {z[1]= TWO10; z[2]=ZERO; z[3]=X[1]; k=1;} | |
207 | } | |
208 | else { | |
209 | if (EX==-42) {z[1]=X[1]+TWO10; z[2]=X[2]; k=3;} | |
210 | else if (EX==-43) {z[1]= TWO10; z[2]=X[1]; k=2;} | |
211 | else {z[1]= TWO10; z[2]=ZERO; k=1;} | |
212 | z[3] = X[k]; | |
213 | } | |
214 | ||
215 | u = (z[3] + TWO57) - TWO57; | |
216 | if (u > z[3]) u -= TWO5; | |
217 | ||
218 | if (u==z[3]) { | |
219 | for (i=k+1; i <= p2; i++) { | |
220 | if (X[i] == ZERO) continue; | |
221 | else {z[3] += ONE; break; } | |
222 | } | |
223 | } | |
224 | ||
225 | c = X[0]*((z[1] + R*(z[2] + R*z[3])) - TWO10); | |
226 | ||
227 | *y = c*TWOM1032; | |
228 | return; | |
229 | ||
230 | #undef R | |
231 | } | |
232 | ||
233 | /* Convert a multiple precision number *x into a double precision number *y. */ | |
234 | /* The result is correctly rounded to the nearest/even. *x is left unchanged */ | |
235 | ||
236 | void __mp_dbl(const mp_no *x, double *y, int p) { | |
237 | #if 0 | |
238 | int i,k; | |
239 | double a,c,u,v,z[5]; | |
240 | #endif | |
241 | ||
242 | if (X[0] == ZERO) {*y = ZERO; return; } | |
243 | ||
244 | if (EX> -42) norm(x,y,p); | |
245 | else if (EX==-42 && X[1]>=TWO10) norm(x,y,p); | |
246 | else denorm(x,y,p); | |
247 | } | |
248 | ||
249 | ||
250 | /* dbl_mp() converts a double precision number x into a multiple precision */ | |
251 | /* number *y. If the precision p is too small the result is truncated. x is */ | |
252 | /* left unchanged. */ | |
253 | ||
254 | void __dbl_mp(double x, mp_no *y, int p) { | |
255 | ||
256 | long i,n; | |
257 | long p2 = p; | |
258 | double u; | |
259 | ||
260 | /* Sign */ | |
261 | if (x == ZERO) {Y[0] = ZERO; return; } | |
262 | else if (x > ZERO) Y[0] = ONE; | |
263 | else {Y[0] = MONE; x=-x; } | |
264 | ||
265 | /* Exponent */ | |
266 | for (EY=ONE; x >= RADIX; EY += ONE) x *= RADIXI; | |
267 | for ( ; x < ONE; EY -= ONE) x *= RADIX; | |
268 | ||
269 | /* Digits */ | |
270 | n=MIN(p2,4); | |
271 | for (i=1; i<=n; i++) { | |
272 | u = (x + TWO52) - TWO52; | |
273 | if (u>x) u -= ONE; | |
274 | Y[i] = u; x -= u; x *= RADIX; } | |
275 | for ( ; i<=p2; i++) Y[i] = ZERO; | |
276 | return; | |
277 | } | |
278 | ||
279 | ||
280 | /* add_magnitudes() adds the magnitudes of *x & *y assuming that */ | |
281 | /* abs(*x) >= abs(*y) > 0. */ | |
282 | /* The sign of the sum *z is undefined. x&y may overlap but not x&z or y&z. */ | |
283 | /* No guard digit is used. The result equals the exact sum, truncated. */ | |
284 | /* *x & *y are left unchanged. */ | |
285 | ||
286 | static void add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { | |
287 | ||
288 | long i,j,k; | |
289 | long p2 = p; | |
290 | ||
291 | EZ = EX; | |
292 | ||
293 | i=p2; j=p2+ EY - EX; k=p2+1; | |
294 | ||
295 | if (j<1) | |
296 | {__cpy(x,z,p); return; } | |
297 | else Z[k] = ZERO; | |
298 | ||
299 | for (; j>0; i--,j--) { | |
300 | Z[k] += X[i] + Y[j]; | |
301 | if (Z[k] >= RADIX) { | |
302 | Z[k] -= RADIX; | |
303 | Z[--k] = ONE; } | |
304 | else | |
305 | Z[--k] = ZERO; | |
306 | } | |
307 | ||
308 | for (; i>0; i--) { | |
309 | Z[k] += X[i]; | |
310 | if (Z[k] >= RADIX) { | |
311 | Z[k] -= RADIX; | |
312 | Z[--k] = ONE; } | |
313 | else | |
314 | Z[--k] = ZERO; | |
315 | } | |
316 | ||
317 | if (Z[1] == ZERO) { | |
318 | for (i=1; i<=p2; i++) Z[i] = Z[i+1]; } | |
319 | else EZ += ONE; | |
320 | } | |
321 | ||
322 | ||
323 | /* sub_magnitudes() subtracts the magnitudes of *x & *y assuming that */ | |
324 | /* abs(*x) > abs(*y) > 0. */ | |
325 | /* The sign of the difference *z is undefined. x&y may overlap but not x&z */ | |
326 | /* or y&z. One guard digit is used. The error is less than one ulp. */ | |
327 | /* *x & *y are left unchanged. */ | |
328 | ||
329 | static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) { | |
330 | ||
331 | long i,j,k; | |
332 | long p2 = p; | |
333 | ||
334 | EZ = EX; | |
335 | ||
336 | if (EX == EY) { | |
337 | i=j=k=p2; | |
338 | Z[k] = Z[k+1] = ZERO; } | |
339 | else { | |
340 | j= EX - EY; | |
341 | if (j > p2) {__cpy(x,z,p); return; } | |
342 | else { | |
343 | i=p2; j=p2+1-j; k=p2; | |
344 | if (Y[j] > ZERO) { | |
345 | Z[k+1] = RADIX - Y[j--]; | |
346 | Z[k] = MONE; } | |
347 | else { | |
348 | Z[k+1] = ZERO; | |
349 | Z[k] = ZERO; j--;} | |
350 | } | |
351 | } | |
352 | ||
353 | for (; j>0; i--,j--) { | |
354 | Z[k] += (X[i] - Y[j]); | |
355 | if (Z[k] < ZERO) { | |
356 | Z[k] += RADIX; | |
357 | Z[--k] = MONE; } | |
358 | else | |
359 | Z[--k] = ZERO; | |
360 | } | |
361 | ||
362 | for (; i>0; i--) { | |
363 | Z[k] += X[i]; | |
364 | if (Z[k] < ZERO) { | |
365 | Z[k] += RADIX; | |
366 | Z[--k] = MONE; } | |
367 | else | |
368 | Z[--k] = ZERO; | |
369 | } | |
370 | ||
371 | for (i=1; Z[i] == ZERO; i++) ; | |
372 | EZ = EZ - i + 1; | |
373 | for (k=1; i <= p2+1; ) | |
374 | Z[k++] = Z[i++]; | |
375 | for (; k <= p2; ) | |
376 | Z[k++] = ZERO; | |
377 | ||
378 | return; | |
379 | } | |
380 | ||
381 | ||
382 | /* Add two multiple precision numbers. Set *z = *x + *y. x&y may overlap */ | |
383 | /* but not x&z or y&z. One guard digit is used. The error is less than */ | |
384 | /* one ulp. *x & *y are left unchanged. */ | |
385 | ||
386 | void __add(const mp_no *x, const mp_no *y, mp_no *z, int p) { | |
387 | ||
388 | int n; | |
389 | ||
390 | if (X[0] == ZERO) {__cpy(y,z,p); return; } | |
391 | else if (Y[0] == ZERO) {__cpy(x,z,p); return; } | |
392 | ||
393 | if (X[0] == Y[0]) { | |
394 | if (__acr(x,y,p) > 0) {add_magnitudes(x,y,z,p); Z[0] = X[0]; } | |
395 | else {add_magnitudes(y,x,z,p); Z[0] = Y[0]; } | |
396 | } | |
397 | else { | |
398 | if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p); Z[0] = X[0]; } | |
399 | else if (n == -1) {sub_magnitudes(y,x,z,p); Z[0] = Y[0]; } | |
400 | else Z[0] = ZERO; | |
401 | } | |
402 | return; | |
403 | } | |
404 | ||
405 | ||
406 | /* Subtract two multiple precision numbers. *z is set to *x - *y. x&y may */ | |
407 | /* overlap but not x&z or y&z. One guard digit is used. The error is */ | |
408 | /* less than one ulp. *x & *y are left unchanged. */ | |
409 | ||
410 | void __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) { | |
411 | ||
412 | int n; | |
413 | ||
414 | if (X[0] == ZERO) {__cpy(y,z,p); Z[0] = -Z[0]; return; } | |
415 | else if (Y[0] == ZERO) {__cpy(x,z,p); return; } | |
416 | ||
417 | if (X[0] != Y[0]) { | |
418 | if (__acr(x,y,p) > 0) {add_magnitudes(x,y,z,p); Z[0] = X[0]; } | |
419 | else {add_magnitudes(y,x,z,p); Z[0] = -Y[0]; } | |
420 | } | |
421 | else { | |
422 | if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p); Z[0] = X[0]; } | |
423 | else if (n == -1) {sub_magnitudes(y,x,z,p); Z[0] = -Y[0]; } | |
424 | else Z[0] = ZERO; | |
425 | } | |
426 | return; | |
427 | } | |
428 | ||
429 | ||
430 | /* Multiply two multiple precision numbers. *z is set to *x * *y. x&y */ | |
431 | /* may overlap but not x&z or y&z. In case p=1,2,3 the exact result is */ | |
432 | /* truncated to p digits. In case p>3 the error is bounded by 1.001 ulp. */ | |
433 | /* *x & *y are left unchanged. */ | |
434 | ||
435 | void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) { | |
436 | ||
437 | long i, i1, i2, j, k, k2; | |
438 | long p2 = p; | |
439 | double u, zk, zk2; | |
440 | ||
441 | /* Is z=0? */ | |
442 | if (X[0]*Y[0]==ZERO) | |
443 | { Z[0]=ZERO; return; } | |
444 | ||
445 | /* Multiply, add and carry */ | |
446 | k2 = (p2<3) ? p2+p2 : p2+3; | |
447 | zk = Z[k2]=ZERO; | |
448 | for (k=k2; k>1; ) { | |
449 | if (k > p2) {i1=k-p2; i2=p2+1; } | |
450 | else {i1=1; i2=k; } | |
451 | #if 1 | |
452 | /* rearange this inner loop to allow the fmadd instructions to be | |
453 | independent and execute in parallel on processors that have | |
454 | dual symetrical FP pipelines. */ | |
455 | if (i1 < (i2-1)) | |
456 | { | |
457 | /* make sure we have at least 2 iterations */ | |
458 | if (((i2 - i1) & 1L) == 1L) | |
459 | { | |
460 | /* Handle the odd iterations case. */ | |
461 | zk2 = x->d[i2-1]*y->d[i1]; | |
462 | } | |
463 | else | |
464 | zk2 = zero.d; | |
465 | /* Do two multiply/adds per loop iteration, using independent | |
466 | accumulators; zk and zk2. */ | |
467 | for (i=i1,j=i2-1; i<i2-1; i+=2,j-=2) | |
468 | { | |
469 | zk += x->d[i]*y->d[j]; | |
470 | zk2 += x->d[i+1]*y->d[j-1]; | |
471 | } | |
472 | zk += zk2; /* final sum. */ | |
473 | } | |
474 | else | |
475 | { | |
476 | /* Special case when iterations is 1. */ | |
477 | zk += x->d[i1]*y->d[i1]; | |
478 | } | |
479 | #else | |
480 | /* The orginal code. */ | |
481 | for (i=i1,j=i2-1; i<i2; i++,j--) zk += X[i]*Y[j]; | |
482 | #endif | |
483 | ||
484 | u = (zk + CUTTER)-CUTTER; | |
485 | if (u > zk) u -= RADIX; | |
486 | Z[k] = zk - u; | |
487 | zk = u*RADIXI; | |
488 | --k; | |
489 | } | |
490 | Z[k] = zk; | |
491 | ||
492 | /* Is there a carry beyond the most significant digit? */ | |
493 | if (Z[1] == ZERO) { | |
494 | for (i=1; i<=p2; i++) Z[i]=Z[i+1]; | |
495 | EZ = EX + EY - 1; } | |
496 | else | |
497 | EZ = EX + EY; | |
498 | ||
499 | Z[0] = X[0] * Y[0]; | |
500 | return; | |
501 | } | |
502 | ||
503 | ||
504 | /* Invert a multiple precision number. Set *y = 1 / *x. */ | |
505 | /* Relative error bound = 1.001*r**(1-p) for p=2, 1.063*r**(1-p) for p=3, */ | |
506 | /* 2.001*r**(1-p) for p>3. */ | |
507 | /* *x=0 is not permissible. *x is left unchanged. */ | |
508 | ||
509 | void __inv(const mp_no *x, mp_no *y, int p) { | |
510 | long i; | |
511 | #if 0 | |
512 | int l; | |
513 | #endif | |
514 | double t; | |
515 | mp_no z,w; | |
516 | static const int np1[] = {0,0,0,0,1,2,2,2,2,3,3,3,3,3,3,3,3,3, | |
517 | 4,4,4,4,4,4,4,4,4,4,4,4,4,4,4}; | |
518 | const mp_no mptwo = {1,{1.0,2.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, | |
519 | 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, | |
520 | 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, | |
521 | 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; | |
522 | ||
523 | __cpy(x,&z,p); z.e=0; __mp_dbl(&z,&t,p); | |
524 | t=ONE/t; __dbl_mp(t,y,p); EY -= EX; | |
525 | ||
526 | for (i=0; i<np1[p]; i++) { | |
527 | __cpy(y,&w,p); | |
528 | __mul(x,&w,y,p); | |
529 | __sub(&mptwo,y,&z,p); | |
530 | __mul(&w,&z,y,p); | |
531 | } | |
532 | return; | |
533 | } | |
534 | ||
535 | ||
536 | /* Divide one multiple precision number by another.Set *z = *x / *y. *x & *y */ | |
537 | /* are left unchanged. x&y may overlap but not x&z or y&z. */ | |
538 | /* Relative error bound = 2.001*r**(1-p) for p=2, 2.063*r**(1-p) for p=3 */ | |
539 | /* and 3.001*r**(1-p) for p>3. *y=0 is not permissible. */ | |
540 | ||
541 | void __dvd(const mp_no *x, const mp_no *y, mp_no *z, int p) { | |
542 | ||
543 | mp_no w; | |
544 | ||
545 | if (X[0] == ZERO) Z[0] = ZERO; | |
546 | else {__inv(y,&w,p); __mul(x,&w,z,p);} | |
547 | return; | |
548 | } |