1 /* arith07.c Copyright (C) 1990-2008 Codemist Ltd */ 2 3 /* 4 * Arithmetic functions. negation plus a load of Common Lisp things 5 * for support of complex numbers. 6 * 7 */ 8 9 /************************************************************************** 10 * Copyright (C) 2008, Codemist Ltd. A C Norman * 11 * * 12 * Redistribution and use in source and binary forms, with or without * 13 * modification, are permitted provided that the following conditions are * 14 * met: * 15 * * 16 * * Redistributions of source code must retain the relevant * 17 * copyright notice, this list of conditions and the following * 18 * disclaimer. * 19 * * Redistributions in binary form must reproduce the above * 20 * copyright notice, this list of conditions and the following * 21 * disclaimer in the documentation and/or other materials provided * 22 * with the distribution. * 23 * * 24 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS * 25 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT * 26 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS * 27 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE * 28 * COPYRIGHT OWNERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, * 29 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, * 30 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS * 31 * OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND * 32 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR * 33 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF * 34 * THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH * 35 * DAMAGE. * 36 *************************************************************************/ 37 38 39 40 /* Signature: 55b237df 24-May-2008 */ 41 42 #include "headers.h" 43 44 45 Lisp_Object copyb(Lisp_Object a) 46 /* 47 * copy a bignum. 48 */ 49 { 50 Lisp_Object b, nil; 51 int32_t len = bignum_length(a), i; 52 push(a); 53 b = getvector(TAG_NUMBERS, TYPE_BIGNUM, len); 54 pop(a); 55 errexit(); 56 len = (len-CELL)/4; 57 for (i=0; i<len; i++) 58 bignum_digits(b)[i] = bignum_digits(a)[i]; 59 return b; 60 } 61 62 Lisp_Object negateb(Lisp_Object a) 63 /* 64 * Negate a bignum. Note that negating the 1-word bignum 65 * value of 0x08000000 will produce a fixnum as a result, 66 * which might confuse the caller... in a similar way negating 67 * the value -0x40000000 will need to promote from a one-word 68 * bignum to a 2-word bignum. How messy just for negation! 69 */ 70 { 71 Lisp_Object b, nil; 72 int32_t len = bignum_length(a), i, carry; 73 if (len == CELL+4) /* one-word bignum - do specially */ 74 { len = -(int32_t)bignum_digits(a)[0]; 75 if (len == fix_mask) return fixnum_of_int(len); 76 else if (len == 0x40000000) return make_two_word_bignum(0, len); 77 else return make_one_word_bignum(len); 78 } 79 push(a); 80 b = getvector(TAG_NUMBERS, TYPE_BIGNUM, len); 81 pop(a); 82 errexit(); 83 len = (len-CELL)/4-1; 84 carry = -1; 85 for (i=0; i<len; i++) 86 { carry = clear_top_bit(~bignum_digits(a)[i]) + top_bit(carry); 87 bignum_digits(b)[i] = clear_top_bit(carry); 88 } 89 /* 90 * Handle the top digit separately since it is signed. 91 */ 92 carry = ~bignum_digits(a)[i] + top_bit(carry); 93 if (!signed_overflow(carry)) 94 { 95 /* 96 * If the most significant word ends up as -1 then I just might 97 * have 0x40000000 in the next word down and so I may need to shrink 98 * the number. Since I handled 1-word bignums specially I have at 99 * least two words to deal with here. 100 */ 101 if (carry == -1 && (bignum_digits(b)[i-1] & 0x40000000) != 0) 102 { bignum_digits(b)[i-1] |= ~0x7fffffff; 103 numhdr(b) -= pack_hdrlength(1); 104 if (SIXTY_FOUR_BIT) 105 { if ((i & 1) != 0) bignum_digits(b)[i] = 0; 106 else bignum_digits(b)[i] = make_bighdr(2); 107 } 108 else 109 { if ((i & 1) == 0) bignum_digits(b)[i] = 0; 110 else bignum_digits(b)[i] = make_bighdr(2); 111 } 112 } 113 else bignum_digits(b)[i] = carry; /* no shrinking needed */ 114 return b; 115 } 116 /* 117 * Here I have overflow: this can only happen when I negate a number 118 * that started off with 0xc0000000 in the most significant digit, 119 * and I have to pad a zero word onto the front. 120 */ 121 bignum_digits(b)[i] = clear_top_bit(carry); 122 return lengthen_by_one_bit(b, carry); 123 } 124 125 /* 126 * generic negation 127 */ 128 129 Lisp_Object negate(Lisp_Object a) 130 { 131 #ifdef COMMON 132 Lisp_Object nil; /* needed for errexit() */ 133 #endif 134 switch ((int)a & TAG_BITS) 135 { 136 case TAG_FIXNUM: 137 { int32_t aa = -int_of_fixnum(a); 138 /* 139 * negating the number -#x8000000 (which is a fixnum) yields a value 140 * which just fails to be a fixnum. 141 */ 142 if (aa != 0x08000000) return fixnum_of_int(aa); 143 else return make_one_word_bignum(aa); 144 } 145 #ifdef COMMON 146 case TAG_SFLOAT: 147 { Float_union aa; 148 aa.i = a - TAG_SFLOAT; 149 aa.f = (float) (-aa.f); 150 return (aa.i & ~(int32_t)0xf) + TAG_SFLOAT; 151 } 152 #endif 153 case TAG_NUMBERS: 154 { int32_t ha = type_of_header(numhdr(a)); 155 switch (ha) 156 { 157 case TYPE_BIGNUM: 158 return negateb(a); 159 #ifdef COMMON 160 case TYPE_RATNUM: 161 { Lisp_Object n = numerator(a), 162 d = denominator(a); 163 push(d); 164 n = negate(n); 165 pop(d); 166 errexit(); 167 return make_ratio(n, d); 168 } 169 case TYPE_COMPLEX_NUM: 170 { Lisp_Object r = real_part(a), 171 i = imag_part(a); 172 push(i); 173 r = negate(r); 174 pop(i); 175 errexit(); 176 push(r); 177 i = negate(i); 178 pop(r); 179 errexit(); 180 return make_complex(r, i); 181 } 182 #endif 183 default: 184 return aerror1("bad arg for minus", a); 185 } 186 } 187 case TAG_BOXFLOAT: 188 { double d = float_of_number(a); 189 return make_boxfloat(-d, type_of_header(flthdr(a))); 190 } 191 default: 192 return aerror1("bad arg for minus", a); 193 } 194 } 195 196 197 /*****************************************************************************/ 198 /*** Transcendental functions etcetera. ***/ 199 /*****************************************************************************/ 200 201 202 /* 203 * Much of the code here is extracted from the portable Fortran library 204 * used by Codemist with its Fortran compiler. 205 */ 206 207 /* 208 * The object of the following macro is to adjust the floating point 209 * variables concerned so that the more significant one can be squared 210 * with NO LOSS OF PRECISION. It is only used when there is no danger 211 * of over- or under-flow. 212 * 213 * This code is NOT PORTABLE but can be modified for use elsewhere 214 * It should, however, serve for IEEE and IBM FP formats. 215 */ 216 217 #define _fp_normalize(high, low) \ 218 { double temp; /* access to representation */ \ 219 temp = high; /* take original number */ \ 220 ((int32_t *)&temp)[current_fp_rep & FP_WORD_ORDER] = 0; \ 221 /* make low part of mantissa 0 */ \ 222 low += (high - temp); /* add into low-order result */ \ 223 high = temp; \ 224 } 225 226 #ifdef COMMON 227 228 double MS_CDECL Cabs(Complex z) 229 { 230 /* 231 * Obtain the absolute value of a complex number - note that the main 232 * agony here is in ensuring that neither overflow nor underflow can 233 * wreck the calculation. Given ideal arithmetic the sum could be carried 234 * through as just sqrt(x^2 + y^2). 235 */ 236 double x = z.real, y = z.imag; 237 double scale; 238 int n1, n2; 239 if (x==0.0) return fabs(y); 240 else if (y==0.0) return fabs(x); 241 (void)frexp(x, &n1); 242 (void)frexp(y, &n2); 243 /* The exact range of values returned by frexp does not matter here */ 244 if (n2>n1) n1 = n2; 245 /* n1 is now the exponent of the larger (in absolute value) of x, y */ 246 scale = ldexp(1.0, n1); /* can not be 0.0 */ 247 x /= scale; 248 y /= scale; 249 /* The above scaling operation introduces no rounding error (since the */ 250 /* scale factor is exactly a power of 2). It reduces the larger of x, y */ 251 /* to be somewhere near 1.0 so overflow in x*x+y*y is impossible. It is */ 252 /* still possible that one of x*x and y*y will underflow (but not both) */ 253 /* but this is harmless. */ 254 return scale * sqrt(x*x + y*y); 255 } 256 257 Complex MS_CDECL Ccos(Complex z) 258 { 259 double x = z.real, y = z.imag; 260 /* 261 * cos(x + iy) = cos(x)*cosh(y) - i sin(x)*sinh(y) 262 * For smallish y this can be used directly. For |y| > 50 I will 263 * compute sinh and cosh as just +/- exp(|y|)/2 264 */ 265 double s = sin(x), c = cos(x); 266 double absy = fabs(y); 267 if (absy <= 50.0) 268 { double sh = sinh(y), ch = cosh(y); 269 z.real = c*ch; 270 z.imag = - s*sh; 271 return z; 272 } 273 else 274 { 275 double w; 276 int n = _reduced_exp(absy, &w) - 1; 277 z.real = ldexp(c*w, n); 278 if (y < 0.0) z.imag = ldexp(s*w, n); 279 else z.imag = ldexp(-s*w, n); 280 return z; 281 } 282 } 283 284 static double reduced_power(double a, int n) 285 { 286 /* 287 * Compute (1 + a)^n - 1 avoiding undue roundoff error. 288 * Assumes n >= 1 on entry and that a is small. 289 */ 290 if (n == 1) return a; 291 { double d = reduced_power(a, n/2); 292 d = (2.0 + d)*d; 293 if (n & 1) d += (1.0 + d)*a; 294 return d; 295 } 296 } 297 298 /* 299 * The following value is included for documentation purposes - it 300 * give the largest args that can be given to exp() without leading to 301 * overflow on IEEE-arithmetic machines. 302 * #define _exp_arg_limit 709.78271289338397 303 * Note that in any case exp(50.0) will not overflow (it is only 5.2e21), 304 * so it can be evaluated the simple direct way. 305 */ 306 307 int _reduced_exp(double x, double *r) 308 { 309 /* 310 * (*r) = exp(x)/2^n; return n; 311 * where n will be selected so that *r gets set to a fairly small value 312 * (precise range of r unimportant provided it will be WELL away from 313 * chances of overflow, even when exp(x) would actually overflow). This 314 * function may only be called with argument x that is positive and 315 * large enough that n will end up satisfying n>=1. The coding here 316 * will ensure that if x>4.0, and in general the use of this function 317 * will only be for x > 50. 318 * For IBM hardware it would be good to be able to control the value 319 * of n mod 4, maybe, to help counter wobbling precision. This is not 320 * done here. 321 */ 322 int n; 323 double f; 324 n = (int)(x / 7.625 + 0.5); 325 /* 326 * 7.625 = 61/8 and is expected to have an exact floating point 327 * representation here, so f is computed without any rounding error. 328 * (do I need something like the (x - 0.5) - 0.5 trick here?) 329 */ 330 f = exp(x - 7.625*(double)n); 331 /* 332 * the magic constant is ((exp(61/8) / 2048) - 1) and it arises because 333 * 61/88 is a decent rational approximation to log(2), hence exp(61/8) 334 * is almost 2^11. Thus I compute exp(x) as 335 * 2^(11*n) * (exp(61/8)/2^11)^n * exp(f) 336 * The first factor is exact, the second is (1+e)^n where e is small and 337 * n is an integer, so can be computer accurately, and the residue f at the 338 * end is small enough not to give over-bad trouble. 339 * The numeric constant given here was calculated with the REDUCE 3.3 340 * bigfloat package. 341 */ 342 #define _e61q8 3.81086435594567676751e-4 343 *r = reduced_power(_e61q8, n)*f + f; 344 #undef _e61q8 345 return 11*n; 346 } 347 348 Complex MS_CDECL Cexp(Complex z) 349 { 350 double x = z.real, y = z.imag; 351 /* 352 * value is exp(x)*(cos(y) + i sin(y)) but have care with overflow 353 * Here (and throughout the complex library) there is an opportunity 354 * to save time by computing sin(y) and cos(y) together. Since this 355 * code is (to begin with) to sit on top of an arbitrary C library, 356 * perhaps with hardware support for the calculation of real-valued 357 * trig functions I am not going to try to realise this saving. 358 */ 359 double s = sin(y), c = cos(y); 360 /* 361 * if x > 50 I will use a cautious sceme which computes exp(x) with 362 * its (binary) exponent separated. Note that 50.0 is chosen as a 363 * number noticably smaller than _exp_arg_limit (exp(50) = 5.18e21), 364 * but is not a critical very special number. 365 */ 366 if (x <= 50.0) /* includes x < 0.0, of course */ 367 { double w = exp(x); 368 z.real = w*c; 369 z.imag = w*s; 370 return z; 371 } 372 else 373 { double w; 374 int n = _reduced_exp(x, &w); 375 z.real = ldexp(w*c, n); 376 z.imag = ldexp(w*s, n); 377 return z; 378 } 379 } 380 381 Complex MS_CDECL Cln(Complex z) 382 { 383 double x = z.real, y = z.imag; 384 /* 385 * if x and y are both very large then cabs(z) may be out of range 386 * even though log or if is OK. Thus it is necessary to perform an 387 * elaborate scaled calculation here, and not just 388 * z.real = log(cabs(z)); 389 */ 390 double scale, r; 391 int n1, n2; 392 if (x==0.0) r = log(fabs(y)); 393 else if (y==0.0) r = log(fabs(x)); 394 else 395 { 396 (void)frexp(x, &n1); 397 (void)frexp(y, &n2); 398 /* The exact range of values returned by frexp does not matter here */ 399 if (n2>n1) n1 = n2; 400 scale = ldexp(1.0, n1); 401 x /= scale; 402 y /= scale; 403 r = log(scale) + 0.5*log(x*x + y*y); 404 } 405 z.real = r; 406 /* 407 * The C standard is not very explicit about the behaviour of atan2(0.0, -n) 408 * while for Fortran it is necessary that this returns +pi not -pi. Hence 409 * with extreme caution I put a special test here. 410 */ 411 if (y == 0.0) 412 if (x < 0.0) z.imag = _pi; 413 else z.imag = 0.0; 414 else z.imag = atan2(y, x); 415 return z; 416 } 417 418 419 /* 420 * Complex raising to a power. This seems to be pretty nasty 421 * to get right, and the code includes extra high precision variants 422 * on atan() and log(). Further refinements wrt efficiency may be 423 * possible later on. 424 * This code has been partially tested, and seems to be uniformly 425 * better than using just a**b = exp(b*log(a)), but much more careful 426 * study is needed before it can possibly be claimed that it is 427 * right in the sense of not throwing away accuracy when it does not 428 * have to. I also need to make careful checks to verify that the 429 * correct (principal) value is computed. 430 */ 431 432 /* 433 * The next function is used after arithmetic has been done on extra- 434 * precision numbers so that the relationship between high and low parts 435 * is no longer known. Re-instate it. 436 */ 437 438 #define _two_minus_25 2.98023223876953125e-8 /* 2^(-25) */ 439 440 static double fp_add(double a, double b, double *lowres) 441 { 442 /* Result is the high part of a+b, with the low part assigned to *lowres */ 443 double absa, absb; 444 if (a >= 0.0) absa = a; else absa = -a; 445 if (b >= 0.0) absb = b; else absb = -b; 446 if (absa < absb) 447 { double t = a; a = b; b = t; 448 } 449 /* Now a is the larger (in absolute value) of the two numbers */ 450 if (absb > absa * _two_minus_25) 451 { double al = 0.0, bl = 0.0; 452 /* 453 * If the exponent difference beweeen a and b is no more than 25 then 454 * I can add the top part (20 or 24 bits) of a to the top part of b 455 * without going beyond the 52 or 56 bits that a full mantissa can hold. 456 */ 457 _fp_normalize(a, al); 458 _fp_normalize(b, bl); 459 a = a + b; /* No rounding needed here */ 460 b = al + bl; 461 if (a == 0.0) 462 { a = b; 463 b = 0.0; 464 } 465 } 466 /* 467 * The above step leaves b small wrt the value in a (unless a+b led 468 * to substantial cancellation of leading digits), but leaves the high 469 * part a with bits everywhere. Force low part of a to zero. 470 */ 471 { double al = 0.0; 472 _fp_normalize(a, al); 473 b = b + al; 474 } 475 if (a >= 0.0) absa = a; else absa = -a; 476 if (b >= 0.0) absb = b; else absb = -b; 477 if (absb > absa * _two_minus_25) 478 /* 479 * If on input a is close to -b, then a+b is close to zero. In this 480 * case the exponents of a abd b matched, and so earlier calculations 481 * have all been done exactly. Go around again to split residue into 482 * high and low parts. 483 */ 484 { double al = 0.0, bl = 0.0; 485 _fp_normalize(b, bl); 486 a = a + b; 487 _fp_normalize(a, al); 488 b = bl + al; 489 } 490 *lowres = b; 491 return a; 492 } 493 494 #undef _two_minus_25 495 496 static void extended_atan2(double b, double a, double *thetah, double *thetal) 497 { 498 int octant; 499 double rh, rl, thh, thl; 500 /* 501 * First reduce the argument to the first octant (i.e. a, b both +ve, 502 * and b <= a). 503 */ 504 if (b < 0.0) 505 { octant = 4; 506 a = -a; 507 b = -b; 508 } 509 else octant = 0; 510 if (a < 0.0) 511 { double t = a; 512 octant += 2; 513 a = b; 514 b = -t; 515 } 516 if (b > a) 517 { double t = a; 518 octant += 1; 519 a = b; 520 b = t; 521 } 522 523 { static struct { double h; double l;} _atan[] = { 524 /* 525 * The table here gives atan(n/16) in 1.5-precision for n=0..16 526 * Note that all the magic numbers used in this file were calculated 527 * using the REDUCE bigfloat package, and all the 'exact' parts have 528 * a denominator of at worst 2^26 and so are expected to have lots 529 * of trailing zero bits in their floating point representation. 530 */ 531 { 0.0, 0.0 }, 532 { 0.06241881847381591796875, -0.84778585694947708870e-8 }, 533 { 0.124355018138885498046875, -2.35921240630155201508e-8 }, 534 { 0.185347974300384521484375, -2.43046897565983490387e-8 }, 535 { 0.2449786663055419921875, -0.31786778380154175187e-8 }, 536 { 0.302884876728057861328125, -0.83530864557675689054e-8 }, 537 { 0.358770668506622314453125, 0.17639499059427950639e-8 }, 538 { 0.412410438060760498046875, 0.35366268088529162896e-8 }, 539 { 0.4636476039886474609375, 0.50121586552767562314e-8 }, 540 { 0.512389481067657470703125, -2.07569197640365239794e-8 }, 541 { 0.558599293231964111328125, 2.21115983246433832164e-8 }, 542 { 0.602287352085113525390625, -0.59501493437085023057e-8 }, 543 { 0.643501102924346923828125, 0.58689374629746842287e-8 }, 544 { 0.6823165416717529296875, 1.32029951485689299817e-8 }, 545 { 0.71882998943328857421875, 1.01883359311982641515e-8 }, 546 { 0.75315129756927490234375, -1.66070805128190106297e-8 }, 547 { 0.785398185253143310546875, -2.18556950009312141541e-8 } 548 }; 549 int k = (int)(16.0*(b/a + 0.03125)); /* 0 to 16 */ 550 double kd = (double)k/16.0; 551 double ah = a, al = 0.0, 552 bh = b, bl = 0.0, 553 ch, cl, q, q2; 554 _fp_normalize(ah, al); 555 _fp_normalize(bh, bl); 556 ch = bh - ah*kd; cl = bl - al*kd; 557 ah = ah + bh*kd; al = al + bl*kd; 558 bh = ch; bl = cl; 559 /* Now |(a/b)| <= 1/32 */ 560 ah = fp_add(ah, al, &al); /* Re-normalise */ 561 bh = fp_add(bh, bl, &bl); 562 /* Compute approximation to b/a */ 563 rh = (bh + bl)/(ah + al); rl = 0.0; 564 _fp_normalize(rh, rl); 565 bh -= ah*rh; bl -= al*rh; 566 rl = (bh + bl)/(ah + al); /* Quotient now formed */ 567 /* 568 * Now it is necessary to compute atan(q) to one-and-a-half precision. 569 * Since |q| < 1/32 I will leave just q as the high order word of 570 * the result and compute atan(q)-q as a single precision value. This 571 * gives about 16 bits accuracy beyond regular single precision work. 572 */ 573 q = rh + rl; 574 q2 = q*q; 575 /* The expansion the follows could be done better using a minimax poly */ 576 rl -= q*q2*(0.33333333333333333333 - 577 q2*(0.20000000000000000000 - 578 q2*(0.14285714285714285714 - 579 q2*(0.11111111111111111111 - 580 q2* 0.09090909090909090909)))); 581 /* OK - now (rh, rl) is atan(reduced b/a). Need to add on atan(kd) */ 582 rh += _atan[k].h; 583 rl += _atan[k].l; 584 } 585 /* 586 * The following constants give high precision versions of pi and pi/2, 587 * and the high partwords (p2h and pih) have lots of low order zero bits 588 * in their binary representation. Expect (=require) that the arithmetic 589 * that computes thh is done without introduced rounding error. 590 */ 591 #define _p2h 1.57079632580280303955078125 592 #define _p2l 9.92093579680540441639751e-10 593 594 #define _pih 3.14159265160560607910156250 595 #define _pil 1.984187159361080883279502e-9 596 597 switch (octant) 598 { 599 default: 600 case 0: thh = rh; thl = rl; break; 601 case 1: thh = _p2h - rh; thl = _p2l - rl; break; 602 case 2: thh = _p2h + rh; thl = _p2l + rl; break; 603 case 3: thh = _pih - rh; thl = _pil - rl; break; 604 case 4: thh = -_pih + rh; thl = -_pil + rl; break; 605 case 5: thh = -_p2h - rh; thl = -_p2l - rl; break; 606 case 6: thh = -_p2h + rh; thl = -_p2l + rl; break; 607 case 7: thh = -rh; thl = -rl; break; 608 } 609 610 #undef _p2h 611 #undef _p2l 612 #undef _pih 613 #undef _pil 614 615 *thetah = fp_add(thh, thl, thetal); 616 } 617 618 static void extended_log(int k, double a, double b, 619 double *logrh, double *logrl) 620 { 621 /* 622 * If we had exact arithmetic this procedure could be: 623 * k*log(2) + 0.5*log(a^2 + b^2) 624 */ 625 double al = 0.0, bl = 0.0, all = 0.0, bll = 0.0, c, ch, cl, cll; 626 double w, q, qh, ql, rh, rl; 627 int n; 628 /* 629 * First (a^2 + b^2) is calculated, using extra precision. 630 * Because any rounding at this stage can lead to bad errors in 631 * the power that I eventually want to compute, I use 3-word 632 * arithmetic here, and with the version of _fp_normalize given 633 * above and IEEE or IBM370 arithmetic this part of the 634 * computation is exact. 635 */ 636 _fp_normalize(a, al); _fp_normalize(al, all); 637 _fp_normalize(b, bl); _fp_normalize(bl, bll); 638 ch = a*a + b*b; 639 cl = 2.0*(a*al + b*bl); 640 cll = (al*al + bl*bl) + 641 all*(2.0*(a + al) + all) + 642 bll*(2.0*(b + bl) + bll); 643 _fp_normalize(ch, cl); 644 _fp_normalize(cl, cll); 645 c = ch + (cl + cll); /* single precision approximation */ 646 /* 647 * At this stage the scaling of the input value will mean that we 648 * have 0.25 <= c <= 2.0 649 * 650 * Now rewrite things as 651 * (2*k + n)*log(s) + 0.5*log((a^2 + b^2)/2^n)) 652 * where s = sqrt(2) 653 * and where the arg to the log is in sqrt(0.5), sqrt(2) 654 */ 655 656 #define _sqrt_half 0.70710678118654752440 657 #define _sqrt_two 1.41421356237309504880 658 k = 2*k; 659 while (c < _sqrt_half) 660 { k -= 1; 661 ch *= 2.0; 662 cl *= 2.0; 663 cll *= 2.0; 664 c *= 2.0; 665 } 666 while (c > _sqrt_two) 667 { k += 1; 668 ch *= 0.5; 669 cl *= 0.5; 670 cll *= 0.5; 671 c *= 0.5; 672 } 673 #undef _sqrt_half 674 #undef _sqrt_two 675 n = (int)(16.0/c + 0.5); 676 w = (double)n / 16.0; 677 ch *= w; 678 cl *= w; 679 cll *= w; /* Now |c-1| < 0.04317 */ 680 ch = (ch - 0.5) - 0.5; 681 ch = fp_add(ch, cl, &cl); 682 cl = cl + cll; 683 /* 684 * (ch, cl) is now the reduced argument ready for calculating log(1+c), 685 * and now that the reduction is over I can drop back to the use of just 686 * two doubleprecision values to represent c. 687 */ 688 c = ch + cl; 689 qh = c / (2.0 + c); 690 ql = 0.0; 691 _fp_normalize(qh, ql); 692 ql = ((ch - qh*(2.0 + ch)) + cl - qh*cl) / (2.0 + c); 693 /* (qh, ql) is now c/(2.0 + c) */ 694 rh = qh; /* 18 bits bigger than low part will end up */ 695 q = qh + ql; 696 w = q*q; 697 rl = ql + q*w*(0.33333333333333333333 + 698 w*(0.20000000000000000000 + 699 w*(0.14285714285714285714 + 700 w*(0.11111111111111111111 + 701 w*(0.09090909090909090909))))); 702 /* 703 * (rh, rl) is now atan(c) correct to double precision plus about 18 bits. 704 */ 705 706 { double temp; 707 static struct { double h; double l; } _log_table[] = { 708 /* 709 * The following values are (in extra precision) -log(n/16)/2 for n 710 * in the range 11 to 23 (i.e. roughly over sqrt(2)/2 to sqrt(2)) 711 */ 712 { 0.1873466968536376953125, 2.786706765149099245e-8 }, 713 { 0.1438410282135009765625, 0.801238948715710950e-8 }, 714 { 0.103819668292999267578125, 1.409612298322959552e-8 }, 715 { 0.066765725612640380859375, -2.930037906928620319e-8 }, 716 { 0.0322692394256591796875, 2.114312640614896196e-8 }, 717 { 0.0, 0.0 }, 718 { -0.0303122997283935546875, -1.117982386660280307e-8 }, 719 { -0.0588915348052978515625, 1.697710612429310295e-8 }, 720 { -0.08592510223388671875, -2.622944289242004947e-8 }, 721 { -0.111571788787841796875, 1.313073691899185245e-8 }, 722 { -0.135966837406158447265625, -2.033566243215020975e-8 }, 723 { -0.159226894378662109375, 2.881939480146987639e-8 }, 724 { -0.18145275115966796875, 0.431498374218108783e-8 }}; 725 726 rh = fp_add(rh, _log_table[n-11].h, &temp); 727 rl = temp + _log_table[n-11].l + rl; 728 } 729 730 #define _exact_part_logroot2 0.3465735912322998046875 731 #define _approx_part_logroot2 (-9.5232714997888393927e-10) 732 /* Multiply this by k and add it in */ 733 { double temp, kd = (double)k; 734 rh = fp_add(rh, kd*_exact_part_logroot2, &temp); 735 rl = rl + temp + kd*_approx_part_logroot2; 736 } 737 #undef _exact_part_logroot2 738 #undef _approx_part_logroot2 739 740 *logrh = rh; 741 *logrl = rl; 742 return; 743 } 744 745 746 Complex MS_CDECL Cpow(Complex z1, Complex z2) 747 { 748 double a = z1.real, b = z1.imag, 749 c = z2.real, d = z2.imag; 750 int k, m, n; 751 double scale; 752 double logrh, logrl, thetah, thetal; 753 double r, i, rh, rl, ih, il, clow, dlow, q; 754 double cw, sw, cost, sint; 755 756 if (b == 0.0 && d == 0.0 && a >= 0.0)/* Simple case if both args are real */ 757 { z1.real = pow(a, c); 758 z1.imag = 0.0; 759 return z1; 760 } 761 /* 762 * Start by scaling z1 by dividing out a power of 2 so that |z1| is 763 * fairly close to 1 764 */ 765 if (a == 0.0) 766 { if (b == 0.0) return z1; /* 0.0**anything is really an error */ 767 /* The exact values returned by frexp do not matter here */ 768 (void)frexp(b, &k); 769 } 770 else 771 { (void)frexp(a, &k); 772 if (b != 0.0) 773 { int n; 774 (void)frexp(b, &n); 775 if (n > k) k = n; 776 } 777 } 778 scale = ldexp(1.0, k); 779 a /= scale; 780 b /= scale; 781 /* 782 * The idea next is to express z1 as r*exp(i theta), then z1**z2 783 * is exp(z2*log(z1)) and the exponent simplifies to 784 * (c*log r - d*theta) + i(c theta + d log r). 785 * Note r = scale*sqrt(a*a + b*b) now. 786 * The argument for exp() must be computed to noticably more than 787 * regular precision, since otherwise exp() greatly magnifies the 788 * error in the real part (if it is large and positive) and range 789 * reduction in the imaginary part can equally lead to accuracy 790 * problems. Neither c nor d can usefully be anywhere near the 791 * extreme range of floating point values, so I will not even try 792 * to scale them, thus I will end up happy(ish) if I can compute 793 * atan2(b, a) and log(|z1|) in one-and-a-half precision form. 794 * It would be best to compute theta in units of pi/4 rather than in 795 * raidians, since then the case of z^n (integer n) with arg(z) an 796 * exact multiple of pi/4 would work out better. However at present I 797 * just hope that the 1.5-precision working is good enough. 798 */ 799 extended_atan2(b, a, &thetah, &thetal); 800 extended_log(k, a, b, &logrh, &logrl); 801 802 /* Normalise all numbers */ 803 clow = 0.0; 804 dlow = 0.0; 805 _fp_normalize(c, clow); 806 _fp_normalize(d, dlow); 807 808 /* (rh, rl) = c*logr - d*theta; */ 809 rh = c*logrh - d*thetah; /* No rounding in this computation */ 810 rl = c*logrl + clow*(logrh + logrl) - d*thetal - dlow*(thetah + thetal); 811 812 /* (ih, il) = c*theta + d*logr; */ 813 ih = c*thetah + d*logrh; /* No rounding in this computation */ 814 il = c*thetal + clow*(thetah + thetal) + d*logrl + dlow*(logrh + logrl); 815 816 /* 817 * Now it remains to take the exponential of the extended precision 818 * value in ((rh, rl), (ih, il)). 819 * Range reduce the real part by multiples of log(2), and the imaginary 820 * part by multiples of pi/2. 821 */ 822 823 rh = fp_add(rh, rl, &rl); /* renormalise */ 824 r = rh + rl; /* Approximate value */ 825 826 #define _recip_log_2 1.4426950408889634074 827 828 q = r * _recip_log_2; 829 m = (q < 0.0) ? (int)(q - 0.5) : (int)(q + 0.5); 830 q = (double)m; 831 #undef _recip_log_2 832 /* 833 * log 2 = 11629080/2^24 - 0.000 00000 19046 54299 95776 78785 834 * to reasonable accuracy. It is vital that the exact part be 835 * read in as a number with lots of low zero bits, so when it gets 836 * multiplied by the integer q there is NO rounding needed. 837 */ 838 839 #define _exact_part_log2 0.693147182464599609375 840 #define _approx_part_log2 (-1.9046542999577678785418e-9) 841 rh = rh - q*_exact_part_log2; 842 rl = rl - q*_approx_part_log2; 843 #undef _exact_part_log2 844 #undef _approx_part_log2 845 r = exp(rh + rl); /* This should now be accurate enough */ 846 847 ih = fp_add(ih, il, &il); 848 i = ih + il; /* Approximate value */ 849 850 #define _recip_pi_by_2 0.6366197723675813431 851 852 q = i * _recip_pi_by_2; 853 n = (q < 0.0) ? (int)(q - 0.5) : (int)(q + 0.5); 854 q = (double)n; 855 /* 856 * pi/2 = 105414357/2^26 + 0.000 00000 09920 93579 68054 04416 39751 857 * to reasonable accuracy. It is vital that the exact part be 858 * read in as a number with lots of low zero bits, so when it gets 859 * multiplied by the integer q there is NO rounding needed. 860 */ 861 862 #define _exact_part_pi2 1.57079632580280303955078125 863 #define _approx_part_pi2 9.92093579680540441639751e-10 864 ih = ih - q*_exact_part_pi2; 865 il = il - q*_approx_part_pi2; 866 #undef _recip_pi_by_2 867 #undef _exact_part_pi2 868 #undef _approx_part_pi2 869 i = ih + il; /* Now accurate enough */ 870 /* 871 * Having done super-careful range reduction I can call the regular 872 * sin/cos routines here. If sin/cos could both be computed together 873 * that could speed things up a little bit, but by this stage they have 874 * not much in common. 875 */ 876 cw = cos(i); 877 sw = sin(i); 878 switch (n & 3) /* quadrant control */ 879 { 880 default: 881 case 0: cost = cw; sint = sw; break; 882 case 1: cost = -sw; sint = cw; break; 883 case 2: cost = -cw; sint = -sw; break; 884 case 3: cost = sw; sint = -cw; break; 885 } 886 887 /* 888 * Now, at long last, I can assemble the results and return. 889 */ 890 891 z1.real = ldexp(r*cost, m); 892 z1.imag = ldexp(r*sint, m); 893 return z1; 894 } 895 896 /* 897 * End of complex-to-complex-power code. 898 */ 899 900 #endif 901 902 /* end of arith07.c */ 903