1
2 /*
3 * IBM Accurate Mathematical Library
4 * written by International Business Machines Corp.
5 * Copyright (C) 2001 Free Software Foundation
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
18 * along with this program; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
20 */
21 /************************************************************************/
22 /* MODULE_NAME: mpa.c */
23 /* */
24 /* FUNCTIONS: */
25 /* mcr */
26 /* acr */
27 /* cr */
28 /* cpy */
29 /* cpymn */
30 /* norm */
31 /* denorm */
32 /* mp_dbl */
33 /* dbl_mp */
34 /* add_magnitudes */
35 /* sub_magnitudes */
36 /* add */
37 /* sub */
38 /* mul */
39 /* inv */
40 /* dvd */
41 /* */
42 /* Arithmetic functions for multiple precision numbers. */
43 /* Relative errors are bounded */
44 /************************************************************************/
45
46
47 #include "endian.h"
48 #include "mpa.h"
49 #include "mpa2.h"
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. */
mcr(const mp_no * x,const mp_no * y,int p)54 static int mcr(const mp_no *x, const mp_no *y, int p) {
55 int i;
56 for (i=1; i<=p; i++) {
57 if (X[i] == Y[i]) continue;
58 else if (X[i] > Y[i]) return 1;
59 else return -1; }
60 return 0;
61 }
62
63
64
65 /* acr() compares the absolute values of two multiple precision numbers */
__acr(const mp_no * x,const mp_no * y,int p)66 int __acr(const mp_no *x, const mp_no *y, int p) {
67 int i;
68
69 if (X[0] == ZERO) {
70 if (Y[0] == ZERO) i= 0;
71 else i=-1;
72 }
73 else if (Y[0] == ZERO) i= 1;
74 else {
75 if (EX > EY) i= 1;
76 else if (EX < EY) i=-1;
77 else i= mcr(x,y,p);
78 }
79
80 return i;
81 }
82
83
84 /* cr90 compares the values of two multiple precision numbers */
__cr(const mp_no * x,const mp_no * y,int p)85 int __cr(const mp_no *x, const mp_no *y, int p) {
86 int i;
87
88 if (X[0] > Y[0]) i= 1;
89 else if (X[0] < Y[0]) i=-1;
90 else if (X[0] < ZERO ) i= __acr(y,x,p);
91 else i= __acr(x,y,p);
92
93 return i;
94 }
95
96
97 /* Copy a multiple precision number. Set *y=*x. x=y is permissible. */
__cpy(const mp_no * x,mp_no * y,int p)98 void __cpy(const mp_no *x, mp_no *y, int p) {
99 int i;
100
101 EY = EX;
102 for (i=0; i <= p; i++) Y[i] = X[i];
103
104 return;
105 }
106
107
108 /* Copy a multiple precision number x of precision m into a */
109 /* multiple precision number y of precision n. In case n>m, */
110 /* the digits of y beyond the m'th are set to zero. In case */
111 /* n<m, the digits of x beyond the n'th are ignored. */
112 /* x=y is permissible. */
113
__cpymn(const mp_no * x,int m,mp_no * y,int n)114 void __cpymn(const mp_no *x, int m, mp_no *y, int n) {
115
116 int i,k;
117
118 EY = EX; k=MIN(m,n);
119 for (i=0; i <= k; i++) Y[i] = X[i];
120 for ( ; i <= n; i++) Y[i] = ZERO;
121
122 return;
123 }
124
125 /* Convert a multiple precision number *x into a double precision */
126 /* number *y, normalized case (|x| >= 2**(-1022))) */
norm(const mp_no * x,double * y,int p)127 static void norm(const mp_no *x, double *y, int p)
128 {
129 #define R radixi.d
130 int i;
131 #if 0
132 int k;
133 #endif
134 double a,c,u,v,z[5];
135 if (p<5) {
136 if (p==1) c = X[1];
137 else if (p==2) c = X[1] + R* X[2];
138 else if (p==3) c = X[1] + R*(X[2] + R* X[3]);
139 else if (p==4) c =(X[1] + R* X[2]) + R*R*(X[3] + R*X[4]);
140 }
141 else {
142 for (a=ONE, z[1]=X[1]; z[1] < TWO23; )
143 {a *= TWO; z[1] *= TWO; }
144
145 for (i=2; i<5; i++) {
146 z[i] = X[i]*a;
147 u = (z[i] + CUTTER)-CUTTER;
148 if (u > z[i]) u -= RADIX;
149 z[i] -= u;
150 z[i-1] += u*RADIXI;
151 }
152
153 u = (z[3] + TWO71) - TWO71;
154 if (u > z[3]) u -= TWO19;
155 v = z[3]-u;
156
157 if (v == TWO18) {
158 if (z[4] == ZERO) {
159 for (i=5; i <= p; i++) {
160 if (X[i] == ZERO) continue;
161 else {z[3] += ONE; break; }
162 }
163 }
164 else z[3] += ONE;
165 }
166
167 c = (z[1] + R *(z[2] + R * z[3]))/a;
168 }
169
170 c *= X[0];
171
172 for (i=1; i<EX; i++) c *= RADIX;
173 for (i=1; i>EX; i--) c *= RADIXI;
174
175 *y = c;
176 return;
177 #undef R
178 }
179
180 /* Convert a multiple precision number *x into a double precision */
181 /* number *y, denormalized case (|x| < 2**(-1022))) */
denorm(const mp_no * x,double * y,int p)182 static void denorm(const mp_no *x, double *y, int p)
183 {
184 int i,k;
185 double c,u,z[5];
186 #if 0
187 double a,v;
188 #endif
189
190 #define R radixi.d
191 if (EX<-44 || (EX==-44 && X[1]<TWO5))
192 { *y=ZERO; return; }
193
194 if (p==1) {
195 if (EX==-42) {z[1]=X[1]+TWO10; z[2]=ZERO; z[3]=ZERO; k=3;}
196 else if (EX==-43) {z[1]= TWO10; z[2]=X[1]; z[3]=ZERO; k=2;}
197 else {z[1]= TWO10; z[2]=ZERO; z[3]=X[1]; k=1;}
198 }
199 else if (p==2) {
200 if (EX==-42) {z[1]=X[1]+TWO10; z[2]=X[2]; z[3]=ZERO; k=3;}
201 else if (EX==-43) {z[1]= TWO10; z[2]=X[1]; z[3]=X[2]; k=2;}
202 else {z[1]= TWO10; z[2]=ZERO; z[3]=X[1]; k=1;}
203 }
204 else {
205 if (EX==-42) {z[1]=X[1]+TWO10; z[2]=X[2]; k=3;}
206 else if (EX==-43) {z[1]= TWO10; z[2]=X[1]; k=2;}
207 else {z[1]= TWO10; z[2]=ZERO; k=1;}
208 z[3] = X[k];
209 }
210
211 u = (z[3] + TWO57) - TWO57;
212 if (u > z[3]) u -= TWO5;
213
214 if (u==z[3]) {
215 for (i=k+1; i <= p; i++) {
216 if (X[i] == ZERO) continue;
217 else {z[3] += ONE; break; }
218 }
219 }
220
221 c = X[0]*((z[1] + R*(z[2] + R*z[3])) - TWO10);
222
223 *y = c*TWOM1032;
224 return;
225
226 #undef R
227 }
228
229 /* Convert a multiple precision number *x into a double precision number *y. */
230 /* The result is correctly rounded to the nearest/even. *x is left unchanged */
231
__mp_dbl(const mp_no * x,double * y,int p)232 void __mp_dbl(const mp_no *x, double *y, int p) {
233 #if 0
234 int i,k;
235 double a,c,u,v,z[5];
236 #endif
237
238 if (X[0] == ZERO) {*y = ZERO; return; }
239
240 if (EX> -42) norm(x,y,p);
241 else if (EX==-42 && X[1]>=TWO10) norm(x,y,p);
242 else denorm(x,y,p);
243 }
244
245
246 /* dbl_mp() converts a double precision number x into a multiple precision */
247 /* number *y. If the precision p is too small the result is truncated. x is */
248 /* left unchanged. */
249
__dbl_mp(double x,mp_no * y,int p)250 void __dbl_mp(double x, mp_no *y, int p) {
251
252 int i,n;
253 double u;
254
255 /* Sign */
256 if (x == ZERO) {Y[0] = ZERO; return; }
257 else if (x > ZERO) Y[0] = ONE;
258 else {Y[0] = MONE; x=-x; }
259
260 /* Exponent */
261 for (EY=ONE; x >= RADIX; EY += ONE) x *= RADIXI;
262 for ( ; x < ONE; EY -= ONE) x *= RADIX;
263
264 /* Digits */
265 n=MIN(p,4);
266 for (i=1; i<=n; i++) {
267 u = (x + TWO52) - TWO52;
268 if (u>x) u -= ONE;
269 Y[i] = u; x -= u; x *= RADIX; }
270 for ( ; i<=p; i++) Y[i] = ZERO;
271 return;
272 }
273
274
275 /* add_magnitudes() adds the magnitudes of *x & *y assuming that */
276 /* abs(*x) >= abs(*y) > 0. */
277 /* The sign of the sum *z is undefined. x&y may overlap but not x&z or y&z. */
278 /* No guard digit is used. The result equals the exact sum, truncated. */
279 /* *x & *y are left unchanged. */
280
add_magnitudes(const mp_no * x,const mp_no * y,mp_no * z,int p)281 static void add_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
282
283 int i,j,k;
284
285 EZ = EX;
286
287 i=p; j=p+ EY - EX; k=p+1;
288
289 if (j<1)
290 {__cpy(x,z,p); return; }
291 else Z[k] = ZERO;
292
293 for (; j>0; i--,j--) {
294 Z[k] += X[i] + Y[j];
295 if (Z[k] >= RADIX) {
296 Z[k] -= RADIX;
297 Z[--k] = ONE; }
298 else
299 Z[--k] = ZERO;
300 }
301
302 for (; i>0; i--) {
303 Z[k] += X[i];
304 if (Z[k] >= RADIX) {
305 Z[k] -= RADIX;
306 Z[--k] = ONE; }
307 else
308 Z[--k] = ZERO;
309 }
310
311 if (Z[1] == ZERO) {
312 for (i=1; i<=p; i++) Z[i] = Z[i+1]; }
313 else EZ += ONE;
314 }
315
316
317 /* sub_magnitudes() subtracts the magnitudes of *x & *y assuming that */
318 /* abs(*x) > abs(*y) > 0. */
319 /* The sign of the difference *z is undefined. x&y may overlap but not x&z */
320 /* or y&z. One guard digit is used. The error is less than one ulp. */
321 /* *x & *y are left unchanged. */
322
sub_magnitudes(const mp_no * x,const mp_no * y,mp_no * z,int p)323 static void sub_magnitudes(const mp_no *x, const mp_no *y, mp_no *z, int p) {
324
325 int i,j,k;
326
327 EZ = EX;
328
329 if (EX == EY) {
330 i=j=k=p;
331 Z[k] = Z[k+1] = ZERO; }
332 else {
333 j= EX - EY;
334 if (j > p) {__cpy(x,z,p); return; }
335 else {
336 i=p; j=p+1-j; k=p;
337 if (Y[j] > ZERO) {
338 Z[k+1] = RADIX - Y[j--];
339 Z[k] = MONE; }
340 else {
341 Z[k+1] = ZERO;
342 Z[k] = ZERO; j--;}
343 }
344 }
345
346 for (; j>0; i--,j--) {
347 Z[k] += (X[i] - Y[j]);
348 if (Z[k] < ZERO) {
349 Z[k] += RADIX;
350 Z[--k] = MONE; }
351 else
352 Z[--k] = ZERO;
353 }
354
355 for (; i>0; i--) {
356 Z[k] += X[i];
357 if (Z[k] < ZERO) {
358 Z[k] += RADIX;
359 Z[--k] = MONE; }
360 else
361 Z[--k] = ZERO;
362 }
363
364 for (i=1; Z[i] == ZERO; i++) ;
365 EZ = EZ - i + 1;
366 for (k=1; i <= p+1; )
367 Z[k++] = Z[i++];
368 for (; k <= p; )
369 Z[k++] = ZERO;
370
371 return;
372 }
373
374
375 /* Add two multiple precision numbers. Set *z = *x + *y. x&y may overlap */
376 /* but not x&z or y&z. One guard digit is used. The error is less than */
377 /* one ulp. *x & *y are left unchanged. */
378
__add(const mp_no * x,const mp_no * y,mp_no * z,int p)379 void __add(const mp_no *x, const mp_no *y, mp_no *z, int p) {
380
381 int n;
382
383 if (X[0] == ZERO) {__cpy(y,z,p); return; }
384 else if (Y[0] == ZERO) {__cpy(x,z,p); return; }
385
386 if (X[0] == Y[0]) {
387 if (__acr(x,y,p) > 0) {add_magnitudes(x,y,z,p); Z[0] = X[0]; }
388 else {add_magnitudes(y,x,z,p); Z[0] = Y[0]; }
389 }
390 else {
391 if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p); Z[0] = X[0]; }
392 else if (n == -1) {sub_magnitudes(y,x,z,p); Z[0] = Y[0]; }
393 else Z[0] = ZERO;
394 }
395 return;
396 }
397
398
399 /* Subtract two multiple precision numbers. *z is set to *x - *y. x&y may */
400 /* overlap but not x&z or y&z. One guard digit is used. The error is */
401 /* less than one ulp. *x & *y are left unchanged. */
402
__sub(const mp_no * x,const mp_no * y,mp_no * z,int p)403 void __sub(const mp_no *x, const mp_no *y, mp_no *z, int p) {
404
405 int n;
406
407 if (X[0] == ZERO) {__cpy(y,z,p); Z[0] = -Z[0]; return; }
408 else if (Y[0] == ZERO) {__cpy(x,z,p); return; }
409
410 if (X[0] != Y[0]) {
411 if (__acr(x,y,p) > 0) {add_magnitudes(x,y,z,p); Z[0] = X[0]; }
412 else {add_magnitudes(y,x,z,p); Z[0] = -Y[0]; }
413 }
414 else {
415 if ((n=__acr(x,y,p)) == 1) {sub_magnitudes(x,y,z,p); Z[0] = X[0]; }
416 else if (n == -1) {sub_magnitudes(y,x,z,p); Z[0] = -Y[0]; }
417 else Z[0] = ZERO;
418 }
419 return;
420 }
421
422
423 /* Multiply two multiple precision numbers. *z is set to *x * *y. x&y */
424 /* may overlap but not x&z or y&z. In case p=1,2,3 the exact result is */
425 /* truncated to p digits. In case p>3 the error is bounded by 1.001 ulp. */
426 /* *x & *y are left unchanged. */
427
__mul(const mp_no * x,const mp_no * y,mp_no * z,int p)428 void __mul(const mp_no *x, const mp_no *y, mp_no *z, int p) {
429
430 int i, i1, i2, j, k, k2;
431 double u;
432
433 /* Is z=0? */
434 if (X[0]*Y[0]==ZERO)
435 { Z[0]=ZERO; return; }
436
437 /* Multiply, add and carry */
438 k2 = (p<3) ? p+p : p+3;
439 Z[k2]=ZERO;
440 for (k=k2; k>1; ) {
441 if (k > p) {i1=k-p; i2=p+1; }
442 else {i1=1; i2=k; }
443 for (i=i1,j=i2-1; i<i2; i++,j--) Z[k] += X[i]*Y[j];
444
445 u = (Z[k] + CUTTER)-CUTTER;
446 if (u > Z[k]) u -= RADIX;
447 Z[k] -= u;
448 Z[--k] = u*RADIXI;
449 }
450
451 /* Is there a carry beyond the most significant digit? */
452 if (Z[1] == ZERO) {
453 for (i=1; i<=p; i++) Z[i]=Z[i+1];
454 EZ = EX + EY - 1; }
455 else
456 EZ = EX + EY;
457
458 Z[0] = X[0] * Y[0];
459 return;
460 }
461
462
463 /* Invert a multiple precision number. Set *y = 1 / *x. */
464 /* Relative error bound = 1.001*r**(1-p) for p=2, 1.063*r**(1-p) for p=3, */
465 /* 2.001*r**(1-p) for p>3. */
466 /* *x=0 is not permissible. *x is left unchanged. */
467
__inv(const mp_no * x,mp_no * y,int p)468 void __inv(const mp_no *x, mp_no *y, int p) {
469 int i;
470 #if 0
471 int l;
472 #endif
473 double t;
474 mp_no z,w;
475 static const int np1[] = {0,0,0,0,1,2,2,2,2,3,3,3,3,3,3,3,3,3,
476 4,4,4,4,4,4,4,4,4,4,4,4,4,4,4};
477 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,
478 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
479 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
480 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}};
481
482 __cpy(x,&z,p); z.e=0; __mp_dbl(&z,&t,p);
483 t=ONE/t; __dbl_mp(t,y,p); EY -= EX;
484
485 for (i=0; i<np1[p]; i++) {
486 __cpy(y,&w,p);
487 __mul(x,&w,y,p);
488 __sub(&mptwo,y,&z,p);
489 __mul(&w,&z,y,p);
490 }
491 return;
492 }
493
494
495 /* Divide one multiple precision number by another.Set *z = *x / *y. *x & *y */
496 /* are left unchanged. x&y may overlap but not x&z or y&z. */
497 /* Relative error bound = 2.001*r**(1-p) for p=2, 2.063*r**(1-p) for p=3 */
498 /* and 3.001*r**(1-p) for p>3. *y=0 is not permissible. */
499
__dvd(const mp_no * x,const mp_no * y,mp_no * z,int p)500 void __dvd(const mp_no *x, const mp_no *y, mp_no *z, int p) {
501
502 mp_no w;
503
504 if (X[0] == ZERO) Z[0] = ZERO;
505 else {__inv(y,&w,p); __mul(x,&w,z,p);}
506 return;
507 }
508