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. */ 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 */ 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 */ 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. */ 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 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))) */ 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))) */ 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 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 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 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 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 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 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 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 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 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