misc.c revision 112158
1112158Sdas/**************************************************************** 2112158Sdas 3112158SdasThe author of this software is David M. Gay. 4112158Sdas 5112158SdasCopyright (C) 1998, 1999 by Lucent Technologies 6112158SdasAll Rights Reserved 7112158Sdas 8112158SdasPermission to use, copy, modify, and distribute this software and 9112158Sdasits documentation for any purpose and without fee is hereby 10112158Sdasgranted, provided that the above copyright notice appear in all 11112158Sdascopies and that both that the copyright notice and this 12112158Sdaspermission notice and warranty disclaimer appear in supporting 13112158Sdasdocumentation, and that the name of Lucent or any of its entities 14112158Sdasnot be used in advertising or publicity pertaining to 15112158Sdasdistribution of the software without specific, written prior 16112158Sdaspermission. 17112158Sdas 18112158SdasLUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, 19112158SdasINCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. 20112158SdasIN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY 21112158SdasSPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 22112158SdasWHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER 23112158SdasIN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, 24112158SdasARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF 25112158SdasTHIS SOFTWARE. 26112158Sdas 27112158Sdas****************************************************************/ 28112158Sdas 29112158Sdas/* Please send bug reports to 30112158Sdas David M. Gay 31112158Sdas Bell Laboratories, Room 2C-463 32112158Sdas 600 Mountain Avenue 33112158Sdas Murray Hill, NJ 07974-0636 34112158Sdas U.S.A. 35112158Sdas dmg@bell-labs.com 36112158Sdas */ 37112158Sdas 38112158Sdas#include "gdtoaimp.h" 39112158Sdas 40112158Sdas static Bigint *freelist[Kmax+1]; 41112158Sdas#ifndef Omit_Private_Memory 42112158Sdas#ifndef PRIVATE_MEM 43112158Sdas#define PRIVATE_MEM 2304 44112158Sdas#endif 45112158Sdas#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double)) 46112158Sdasstatic double private_mem[PRIVATE_mem], *pmem_next = private_mem; 47112158Sdas#endif 48112158Sdas 49112158Sdas Bigint * 50112158SdasBalloc 51112158Sdas#ifdef KR_headers 52112158Sdas (k) int k; 53112158Sdas#else 54112158Sdas (int k) 55112158Sdas#endif 56112158Sdas{ 57112158Sdas int x; 58112158Sdas Bigint *rv; 59112158Sdas#ifndef Omit_Private_Memory 60112158Sdas unsigned int len; 61112158Sdas#endif 62112158Sdas 63112158Sdas ACQUIRE_DTOA_LOCK(0); 64112158Sdas if ( (rv = freelist[k]) !=0) { 65112158Sdas freelist[k] = rv->next; 66112158Sdas } 67112158Sdas else { 68112158Sdas x = 1 << k; 69112158Sdas#ifdef Omit_Private_Memory 70112158Sdas rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong)); 71112158Sdas#else 72112158Sdas len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1) 73112158Sdas /sizeof(double); 74112158Sdas if (pmem_next - private_mem + len <= PRIVATE_mem) { 75112158Sdas rv = (Bigint*)pmem_next; 76112158Sdas pmem_next += len; 77112158Sdas } 78112158Sdas else 79112158Sdas rv = (Bigint*)MALLOC(len*sizeof(double)); 80112158Sdas#endif 81112158Sdas rv->k = k; 82112158Sdas rv->maxwds = x; 83112158Sdas } 84112158Sdas FREE_DTOA_LOCK(0); 85112158Sdas rv->sign = rv->wds = 0; 86112158Sdas return rv; 87112158Sdas } 88112158Sdas 89112158Sdas void 90112158SdasBfree 91112158Sdas#ifdef KR_headers 92112158Sdas (v) Bigint *v; 93112158Sdas#else 94112158Sdas (Bigint *v) 95112158Sdas#endif 96112158Sdas{ 97112158Sdas if (v) { 98112158Sdas ACQUIRE_DTOA_LOCK(0); 99112158Sdas v->next = freelist[v->k]; 100112158Sdas freelist[v->k] = v; 101112158Sdas FREE_DTOA_LOCK(0); 102112158Sdas } 103112158Sdas } 104112158Sdas 105112158Sdas int 106112158Sdaslo0bits 107112158Sdas#ifdef KR_headers 108112158Sdas (y) ULong *y; 109112158Sdas#else 110112158Sdas (ULong *y) 111112158Sdas#endif 112112158Sdas{ 113112158Sdas register int k; 114112158Sdas register ULong x = *y; 115112158Sdas 116112158Sdas if (x & 7) { 117112158Sdas if (x & 1) 118112158Sdas return 0; 119112158Sdas if (x & 2) { 120112158Sdas *y = x >> 1; 121112158Sdas return 1; 122112158Sdas } 123112158Sdas *y = x >> 2; 124112158Sdas return 2; 125112158Sdas } 126112158Sdas k = 0; 127112158Sdas if (!(x & 0xffff)) { 128112158Sdas k = 16; 129112158Sdas x >>= 16; 130112158Sdas } 131112158Sdas if (!(x & 0xff)) { 132112158Sdas k += 8; 133112158Sdas x >>= 8; 134112158Sdas } 135112158Sdas if (!(x & 0xf)) { 136112158Sdas k += 4; 137112158Sdas x >>= 4; 138112158Sdas } 139112158Sdas if (!(x & 0x3)) { 140112158Sdas k += 2; 141112158Sdas x >>= 2; 142112158Sdas } 143112158Sdas if (!(x & 1)) { 144112158Sdas k++; 145112158Sdas x >>= 1; 146112158Sdas if (!x & 1) 147112158Sdas return 32; 148112158Sdas } 149112158Sdas *y = x; 150112158Sdas return k; 151112158Sdas } 152112158Sdas 153112158Sdas Bigint * 154112158Sdasmultadd 155112158Sdas#ifdef KR_headers 156112158Sdas (b, m, a) Bigint *b; int m, a; 157112158Sdas#else 158112158Sdas (Bigint *b, int m, int a) /* multiply by m and add a */ 159112158Sdas#endif 160112158Sdas{ 161112158Sdas int i, wds; 162112158Sdas#ifdef ULLong 163112158Sdas ULong *x; 164112158Sdas ULLong carry, y; 165112158Sdas#else 166112158Sdas ULong carry, *x, y; 167112158Sdas#ifdef Pack_32 168112158Sdas ULong xi, z; 169112158Sdas#endif 170112158Sdas#endif 171112158Sdas Bigint *b1; 172112158Sdas 173112158Sdas wds = b->wds; 174112158Sdas x = b->x; 175112158Sdas i = 0; 176112158Sdas carry = a; 177112158Sdas do { 178112158Sdas#ifdef ULLong 179112158Sdas y = *x * (ULLong)m + carry; 180112158Sdas carry = y >> 32; 181112158Sdas *x++ = y & 0xffffffffUL; 182112158Sdas#else 183112158Sdas#ifdef Pack_32 184112158Sdas xi = *x; 185112158Sdas y = (xi & 0xffff) * m + carry; 186112158Sdas z = (xi >> 16) * m + (y >> 16); 187112158Sdas carry = z >> 16; 188112158Sdas *x++ = (z << 16) + (y & 0xffff); 189112158Sdas#else 190112158Sdas y = *x * m + carry; 191112158Sdas carry = y >> 16; 192112158Sdas *x++ = y & 0xffff; 193112158Sdas#endif 194112158Sdas#endif 195112158Sdas } 196112158Sdas while(++i < wds); 197112158Sdas if (carry) { 198112158Sdas if (wds >= b->maxwds) { 199112158Sdas b1 = Balloc(b->k+1); 200112158Sdas Bcopy(b1, b); 201112158Sdas Bfree(b); 202112158Sdas b = b1; 203112158Sdas } 204112158Sdas b->x[wds++] = carry; 205112158Sdas b->wds = wds; 206112158Sdas } 207112158Sdas return b; 208112158Sdas } 209112158Sdas 210112158Sdas int 211112158Sdashi0bits 212112158Sdas#ifdef KR_headers 213112158Sdas (x) register ULong x; 214112158Sdas#else 215112158Sdas (register ULong x) 216112158Sdas#endif 217112158Sdas{ 218112158Sdas register int k = 0; 219112158Sdas 220112158Sdas if (!(x & 0xffff0000)) { 221112158Sdas k = 16; 222112158Sdas x <<= 16; 223112158Sdas } 224112158Sdas if (!(x & 0xff000000)) { 225112158Sdas k += 8; 226112158Sdas x <<= 8; 227112158Sdas } 228112158Sdas if (!(x & 0xf0000000)) { 229112158Sdas k += 4; 230112158Sdas x <<= 4; 231112158Sdas } 232112158Sdas if (!(x & 0xc0000000)) { 233112158Sdas k += 2; 234112158Sdas x <<= 2; 235112158Sdas } 236112158Sdas if (!(x & 0x80000000)) { 237112158Sdas k++; 238112158Sdas if (!(x & 0x40000000)) 239112158Sdas return 32; 240112158Sdas } 241112158Sdas return k; 242112158Sdas } 243112158Sdas 244112158Sdas Bigint * 245112158Sdasi2b 246112158Sdas#ifdef KR_headers 247112158Sdas (i) int i; 248112158Sdas#else 249112158Sdas (int i) 250112158Sdas#endif 251112158Sdas{ 252112158Sdas Bigint *b; 253112158Sdas 254112158Sdas b = Balloc(1); 255112158Sdas b->x[0] = i; 256112158Sdas b->wds = 1; 257112158Sdas return b; 258112158Sdas } 259112158Sdas 260112158Sdas Bigint * 261112158Sdasmult 262112158Sdas#ifdef KR_headers 263112158Sdas (a, b) Bigint *a, *b; 264112158Sdas#else 265112158Sdas (Bigint *a, Bigint *b) 266112158Sdas#endif 267112158Sdas{ 268112158Sdas Bigint *c; 269112158Sdas int k, wa, wb, wc; 270112158Sdas ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; 271112158Sdas ULong y; 272112158Sdas#ifdef ULLong 273112158Sdas ULLong carry, z; 274112158Sdas#else 275112158Sdas ULong carry, z; 276112158Sdas#ifdef Pack_32 277112158Sdas ULong z2; 278112158Sdas#endif 279112158Sdas#endif 280112158Sdas 281112158Sdas if (a->wds < b->wds) { 282112158Sdas c = a; 283112158Sdas a = b; 284112158Sdas b = c; 285112158Sdas } 286112158Sdas k = a->k; 287112158Sdas wa = a->wds; 288112158Sdas wb = b->wds; 289112158Sdas wc = wa + wb; 290112158Sdas if (wc > a->maxwds) 291112158Sdas k++; 292112158Sdas c = Balloc(k); 293112158Sdas for(x = c->x, xa = x + wc; x < xa; x++) 294112158Sdas *x = 0; 295112158Sdas xa = a->x; 296112158Sdas xae = xa + wa; 297112158Sdas xb = b->x; 298112158Sdas xbe = xb + wb; 299112158Sdas xc0 = c->x; 300112158Sdas#ifdef ULLong 301112158Sdas for(; xb < xbe; xc0++) { 302112158Sdas if ( (y = *xb++) !=0) { 303112158Sdas x = xa; 304112158Sdas xc = xc0; 305112158Sdas carry = 0; 306112158Sdas do { 307112158Sdas z = *x++ * (ULLong)y + *xc + carry; 308112158Sdas carry = z >> 32; 309112158Sdas *xc++ = z & 0xffffffffUL; 310112158Sdas } 311112158Sdas while(x < xae); 312112158Sdas *xc = carry; 313112158Sdas } 314112158Sdas } 315112158Sdas#else 316112158Sdas#ifdef Pack_32 317112158Sdas for(; xb < xbe; xb++, xc0++) { 318112158Sdas if ( (y = *xb & 0xffff) !=0) { 319112158Sdas x = xa; 320112158Sdas xc = xc0; 321112158Sdas carry = 0; 322112158Sdas do { 323112158Sdas z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; 324112158Sdas carry = z >> 16; 325112158Sdas z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; 326112158Sdas carry = z2 >> 16; 327112158Sdas Storeinc(xc, z2, z); 328112158Sdas } 329112158Sdas while(x < xae); 330112158Sdas *xc = carry; 331112158Sdas } 332112158Sdas if ( (y = *xb >> 16) !=0) { 333112158Sdas x = xa; 334112158Sdas xc = xc0; 335112158Sdas carry = 0; 336112158Sdas z2 = *xc; 337112158Sdas do { 338112158Sdas z = (*x & 0xffff) * y + (*xc >> 16) + carry; 339112158Sdas carry = z >> 16; 340112158Sdas Storeinc(xc, z, z2); 341112158Sdas z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; 342112158Sdas carry = z2 >> 16; 343112158Sdas } 344112158Sdas while(x < xae); 345112158Sdas *xc = z2; 346112158Sdas } 347112158Sdas } 348112158Sdas#else 349112158Sdas for(; xb < xbe; xc0++) { 350112158Sdas if ( (y = *xb++) !=0) { 351112158Sdas x = xa; 352112158Sdas xc = xc0; 353112158Sdas carry = 0; 354112158Sdas do { 355112158Sdas z = *x++ * y + *xc + carry; 356112158Sdas carry = z >> 16; 357112158Sdas *xc++ = z & 0xffff; 358112158Sdas } 359112158Sdas while(x < xae); 360112158Sdas *xc = carry; 361112158Sdas } 362112158Sdas } 363112158Sdas#endif 364112158Sdas#endif 365112158Sdas for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ; 366112158Sdas c->wds = wc; 367112158Sdas return c; 368112158Sdas } 369112158Sdas 370112158Sdas static Bigint *p5s; 371112158Sdas 372112158Sdas Bigint * 373112158Sdaspow5mult 374112158Sdas#ifdef KR_headers 375112158Sdas (b, k) Bigint *b; int k; 376112158Sdas#else 377112158Sdas (Bigint *b, int k) 378112158Sdas#endif 379112158Sdas{ 380112158Sdas Bigint *b1, *p5, *p51; 381112158Sdas int i; 382112158Sdas static int p05[3] = { 5, 25, 125 }; 383112158Sdas 384112158Sdas if ( (i = k & 3) !=0) 385112158Sdas b = multadd(b, p05[i-1], 0); 386112158Sdas 387112158Sdas if (!(k >>= 2)) 388112158Sdas return b; 389112158Sdas if ((p5 = p5s) == 0) { 390112158Sdas /* first time */ 391112158Sdas#ifdef MULTIPLE_THREADS 392112158Sdas ACQUIRE_DTOA_LOCK(1); 393112158Sdas if (!(p5 = p5s)) { 394112158Sdas p5 = p5s = i2b(625); 395112158Sdas p5->next = 0; 396112158Sdas } 397112158Sdas FREE_DTOA_LOCK(1); 398112158Sdas#else 399112158Sdas p5 = p5s = i2b(625); 400112158Sdas p5->next = 0; 401112158Sdas#endif 402112158Sdas } 403112158Sdas for(;;) { 404112158Sdas if (k & 1) { 405112158Sdas b1 = mult(b, p5); 406112158Sdas Bfree(b); 407112158Sdas b = b1; 408112158Sdas } 409112158Sdas if (!(k >>= 1)) 410112158Sdas break; 411112158Sdas if ((p51 = p5->next) == 0) { 412112158Sdas#ifdef MULTIPLE_THREADS 413112158Sdas ACQUIRE_DTOA_LOCK(1); 414112158Sdas if (!(p51 = p5->next)) { 415112158Sdas p51 = p5->next = mult(p5,p5); 416112158Sdas p51->next = 0; 417112158Sdas } 418112158Sdas FREE_DTOA_LOCK(1); 419112158Sdas#else 420112158Sdas p51 = p5->next = mult(p5,p5); 421112158Sdas p51->next = 0; 422112158Sdas#endif 423112158Sdas } 424112158Sdas p5 = p51; 425112158Sdas } 426112158Sdas return b; 427112158Sdas } 428112158Sdas 429112158Sdas Bigint * 430112158Sdaslshift 431112158Sdas#ifdef KR_headers 432112158Sdas (b, k) Bigint *b; int k; 433112158Sdas#else 434112158Sdas (Bigint *b, int k) 435112158Sdas#endif 436112158Sdas{ 437112158Sdas int i, k1, n, n1; 438112158Sdas Bigint *b1; 439112158Sdas ULong *x, *x1, *xe, z; 440112158Sdas 441112158Sdas n = k >> kshift; 442112158Sdas k1 = b->k; 443112158Sdas n1 = n + b->wds + 1; 444112158Sdas for(i = b->maxwds; n1 > i; i <<= 1) 445112158Sdas k1++; 446112158Sdas b1 = Balloc(k1); 447112158Sdas x1 = b1->x; 448112158Sdas for(i = 0; i < n; i++) 449112158Sdas *x1++ = 0; 450112158Sdas x = b->x; 451112158Sdas xe = x + b->wds; 452112158Sdas if (k &= kmask) { 453112158Sdas#ifdef Pack_32 454112158Sdas k1 = 32 - k; 455112158Sdas z = 0; 456112158Sdas do { 457112158Sdas *x1++ = *x << k | z; 458112158Sdas z = *x++ >> k1; 459112158Sdas } 460112158Sdas while(x < xe); 461112158Sdas if ((*x1 = z) !=0) 462112158Sdas ++n1; 463112158Sdas#else 464112158Sdas k1 = 16 - k; 465112158Sdas z = 0; 466112158Sdas do { 467112158Sdas *x1++ = *x << k & 0xffff | z; 468112158Sdas z = *x++ >> k1; 469112158Sdas } 470112158Sdas while(x < xe); 471112158Sdas if (*x1 = z) 472112158Sdas ++n1; 473112158Sdas#endif 474112158Sdas } 475112158Sdas else do 476112158Sdas *x1++ = *x++; 477112158Sdas while(x < xe); 478112158Sdas b1->wds = n1 - 1; 479112158Sdas Bfree(b); 480112158Sdas return b1; 481112158Sdas } 482112158Sdas 483112158Sdas int 484112158Sdascmp 485112158Sdas#ifdef KR_headers 486112158Sdas (a, b) Bigint *a, *b; 487112158Sdas#else 488112158Sdas (Bigint *a, Bigint *b) 489112158Sdas#endif 490112158Sdas{ 491112158Sdas ULong *xa, *xa0, *xb, *xb0; 492112158Sdas int i, j; 493112158Sdas 494112158Sdas i = a->wds; 495112158Sdas j = b->wds; 496112158Sdas#ifdef DEBUG 497112158Sdas if (i > 1 && !a->x[i-1]) 498112158Sdas Bug("cmp called with a->x[a->wds-1] == 0"); 499112158Sdas if (j > 1 && !b->x[j-1]) 500112158Sdas Bug("cmp called with b->x[b->wds-1] == 0"); 501112158Sdas#endif 502112158Sdas if (i -= j) 503112158Sdas return i; 504112158Sdas xa0 = a->x; 505112158Sdas xa = xa0 + j; 506112158Sdas xb0 = b->x; 507112158Sdas xb = xb0 + j; 508112158Sdas for(;;) { 509112158Sdas if (*--xa != *--xb) 510112158Sdas return *xa < *xb ? -1 : 1; 511112158Sdas if (xa <= xa0) 512112158Sdas break; 513112158Sdas } 514112158Sdas return 0; 515112158Sdas } 516112158Sdas 517112158Sdas Bigint * 518112158Sdasdiff 519112158Sdas#ifdef KR_headers 520112158Sdas (a, b) Bigint *a, *b; 521112158Sdas#else 522112158Sdas (Bigint *a, Bigint *b) 523112158Sdas#endif 524112158Sdas{ 525112158Sdas Bigint *c; 526112158Sdas int i, wa, wb; 527112158Sdas ULong *xa, *xae, *xb, *xbe, *xc; 528112158Sdas#ifdef ULLong 529112158Sdas ULLong borrow, y; 530112158Sdas#else 531112158Sdas ULong borrow, y; 532112158Sdas#ifdef Pack_32 533112158Sdas ULong z; 534112158Sdas#endif 535112158Sdas#endif 536112158Sdas 537112158Sdas i = cmp(a,b); 538112158Sdas if (!i) { 539112158Sdas c = Balloc(0); 540112158Sdas c->wds = 1; 541112158Sdas c->x[0] = 0; 542112158Sdas return c; 543112158Sdas } 544112158Sdas if (i < 0) { 545112158Sdas c = a; 546112158Sdas a = b; 547112158Sdas b = c; 548112158Sdas i = 1; 549112158Sdas } 550112158Sdas else 551112158Sdas i = 0; 552112158Sdas c = Balloc(a->k); 553112158Sdas c->sign = i; 554112158Sdas wa = a->wds; 555112158Sdas xa = a->x; 556112158Sdas xae = xa + wa; 557112158Sdas wb = b->wds; 558112158Sdas xb = b->x; 559112158Sdas xbe = xb + wb; 560112158Sdas xc = c->x; 561112158Sdas borrow = 0; 562112158Sdas#ifdef ULLong 563112158Sdas do { 564112158Sdas y = (ULLong)*xa++ - *xb++ - borrow; 565112158Sdas borrow = y >> 32 & 1UL; 566112158Sdas *xc++ = y & 0xffffffffUL; 567112158Sdas } 568112158Sdas while(xb < xbe); 569112158Sdas while(xa < xae) { 570112158Sdas y = *xa++ - borrow; 571112158Sdas borrow = y >> 32 & 1UL; 572112158Sdas *xc++ = y & 0xffffffffUL; 573112158Sdas } 574112158Sdas#else 575112158Sdas#ifdef Pack_32 576112158Sdas do { 577112158Sdas y = (*xa & 0xffff) - (*xb & 0xffff) - borrow; 578112158Sdas borrow = (y & 0x10000) >> 16; 579112158Sdas z = (*xa++ >> 16) - (*xb++ >> 16) - borrow; 580112158Sdas borrow = (z & 0x10000) >> 16; 581112158Sdas Storeinc(xc, z, y); 582112158Sdas } 583112158Sdas while(xb < xbe); 584112158Sdas while(xa < xae) { 585112158Sdas y = (*xa & 0xffff) - borrow; 586112158Sdas borrow = (y & 0x10000) >> 16; 587112158Sdas z = (*xa++ >> 16) - borrow; 588112158Sdas borrow = (z & 0x10000) >> 16; 589112158Sdas Storeinc(xc, z, y); 590112158Sdas } 591112158Sdas#else 592112158Sdas do { 593112158Sdas y = *xa++ - *xb++ - borrow; 594112158Sdas borrow = (y & 0x10000) >> 16; 595112158Sdas *xc++ = y & 0xffff; 596112158Sdas } 597112158Sdas while(xb < xbe); 598112158Sdas while(xa < xae) { 599112158Sdas y = *xa++ - borrow; 600112158Sdas borrow = (y & 0x10000) >> 16; 601112158Sdas *xc++ = y & 0xffff; 602112158Sdas } 603112158Sdas#endif 604112158Sdas#endif 605112158Sdas while(!*--xc) 606112158Sdas wa--; 607112158Sdas c->wds = wa; 608112158Sdas return c; 609112158Sdas } 610112158Sdas 611112158Sdas double 612112158Sdasb2d 613112158Sdas#ifdef KR_headers 614112158Sdas (a, e) Bigint *a; int *e; 615112158Sdas#else 616112158Sdas (Bigint *a, int *e) 617112158Sdas#endif 618112158Sdas{ 619112158Sdas ULong *xa, *xa0, w, y, z; 620112158Sdas int k; 621112158Sdas double d; 622112158Sdas#ifdef VAX 623112158Sdas ULong d0, d1; 624112158Sdas#else 625112158Sdas#define d0 word0(d) 626112158Sdas#define d1 word1(d) 627112158Sdas#endif 628112158Sdas 629112158Sdas xa0 = a->x; 630112158Sdas xa = xa0 + a->wds; 631112158Sdas y = *--xa; 632112158Sdas#ifdef DEBUG 633112158Sdas if (!y) Bug("zero y in b2d"); 634112158Sdas#endif 635112158Sdas k = hi0bits(y); 636112158Sdas *e = 32 - k; 637112158Sdas#ifdef Pack_32 638112158Sdas if (k < Ebits) { 639112158Sdas d0 = Exp_1 | y >> Ebits - k; 640112158Sdas w = xa > xa0 ? *--xa : 0; 641112158Sdas d1 = y << (32-Ebits) + k | w >> Ebits - k; 642112158Sdas goto ret_d; 643112158Sdas } 644112158Sdas z = xa > xa0 ? *--xa : 0; 645112158Sdas if (k -= Ebits) { 646112158Sdas d0 = Exp_1 | y << k | z >> 32 - k; 647112158Sdas y = xa > xa0 ? *--xa : 0; 648112158Sdas d1 = z << k | y >> 32 - k; 649112158Sdas } 650112158Sdas else { 651112158Sdas d0 = Exp_1 | y; 652112158Sdas d1 = z; 653112158Sdas } 654112158Sdas#else 655112158Sdas if (k < Ebits + 16) { 656112158Sdas z = xa > xa0 ? *--xa : 0; 657112158Sdas d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; 658112158Sdas w = xa > xa0 ? *--xa : 0; 659112158Sdas y = xa > xa0 ? *--xa : 0; 660112158Sdas d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; 661112158Sdas goto ret_d; 662112158Sdas } 663112158Sdas z = xa > xa0 ? *--xa : 0; 664112158Sdas w = xa > xa0 ? *--xa : 0; 665112158Sdas k -= Ebits + 16; 666112158Sdas d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; 667112158Sdas y = xa > xa0 ? *--xa : 0; 668112158Sdas d1 = w << k + 16 | y << k; 669112158Sdas#endif 670112158Sdas ret_d: 671112158Sdas#ifdef VAX 672112158Sdas word0(d) = d0 >> 16 | d0 << 16; 673112158Sdas word1(d) = d1 >> 16 | d1 << 16; 674112158Sdas#endif 675112158Sdas return dval(d); 676112158Sdas } 677112158Sdas#undef d0 678112158Sdas#undef d1 679112158Sdas 680112158Sdas Bigint * 681112158Sdasd2b 682112158Sdas#ifdef KR_headers 683112158Sdas (d, e, bits) double d; int *e, *bits; 684112158Sdas#else 685112158Sdas (double d, int *e, int *bits) 686112158Sdas#endif 687112158Sdas{ 688112158Sdas Bigint *b; 689112158Sdas int de, i, k; 690112158Sdas ULong *x, y, z; 691112158Sdas#ifdef VAX 692112158Sdas ULong d0, d1; 693112158Sdas d0 = word0(d) >> 16 | word0(d) << 16; 694112158Sdas d1 = word1(d) >> 16 | word1(d) << 16; 695112158Sdas#else 696112158Sdas#define d0 word0(d) 697112158Sdas#define d1 word1(d) 698112158Sdas#endif 699112158Sdas 700112158Sdas#ifdef Pack_32 701112158Sdas b = Balloc(1); 702112158Sdas#else 703112158Sdas b = Balloc(2); 704112158Sdas#endif 705112158Sdas x = b->x; 706112158Sdas 707112158Sdas z = d0 & Frac_mask; 708112158Sdas d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ 709112158Sdas#ifdef Sudden_Underflow 710112158Sdas de = (int)(d0 >> Exp_shift); 711112158Sdas#ifndef IBM 712112158Sdas z |= Exp_msk11; 713112158Sdas#endif 714112158Sdas#else 715112158Sdas if ( (de = (int)(d0 >> Exp_shift)) !=0) 716112158Sdas z |= Exp_msk1; 717112158Sdas#endif 718112158Sdas#ifdef Pack_32 719112158Sdas if ( (y = d1) !=0) { 720112158Sdas if ( (k = lo0bits(&y)) !=0) { 721112158Sdas x[0] = y | z << 32 - k; 722112158Sdas z >>= k; 723112158Sdas } 724112158Sdas else 725112158Sdas x[0] = y; 726112158Sdas i = b->wds = (x[1] = z) !=0 ? 2 : 1; 727112158Sdas } 728112158Sdas else { 729112158Sdas#ifdef DEBUG 730112158Sdas if (!z) 731112158Sdas Bug("Zero passed to d2b"); 732112158Sdas#endif 733112158Sdas k = lo0bits(&z); 734112158Sdas x[0] = z; 735112158Sdas i = b->wds = 1; 736112158Sdas k += 32; 737112158Sdas } 738112158Sdas#else 739112158Sdas if ( (y = d1) !=0) { 740112158Sdas if ( (k = lo0bits(&y)) !=0) 741112158Sdas if (k >= 16) { 742112158Sdas x[0] = y | z << 32 - k & 0xffff; 743112158Sdas x[1] = z >> k - 16 & 0xffff; 744112158Sdas x[2] = z >> k; 745112158Sdas i = 2; 746112158Sdas } 747112158Sdas else { 748112158Sdas x[0] = y & 0xffff; 749112158Sdas x[1] = y >> 16 | z << 16 - k & 0xffff; 750112158Sdas x[2] = z >> k & 0xffff; 751112158Sdas x[3] = z >> k+16; 752112158Sdas i = 3; 753112158Sdas } 754112158Sdas else { 755112158Sdas x[0] = y & 0xffff; 756112158Sdas x[1] = y >> 16; 757112158Sdas x[2] = z & 0xffff; 758112158Sdas x[3] = z >> 16; 759112158Sdas i = 3; 760112158Sdas } 761112158Sdas } 762112158Sdas else { 763112158Sdas#ifdef DEBUG 764112158Sdas if (!z) 765112158Sdas Bug("Zero passed to d2b"); 766112158Sdas#endif 767112158Sdas k = lo0bits(&z); 768112158Sdas if (k >= 16) { 769112158Sdas x[0] = z; 770112158Sdas i = 0; 771112158Sdas } 772112158Sdas else { 773112158Sdas x[0] = z & 0xffff; 774112158Sdas x[1] = z >> 16; 775112158Sdas i = 1; 776112158Sdas } 777112158Sdas k += 32; 778112158Sdas } 779112158Sdas while(!x[i]) 780112158Sdas --i; 781112158Sdas b->wds = i + 1; 782112158Sdas#endif 783112158Sdas#ifndef Sudden_Underflow 784112158Sdas if (de) { 785112158Sdas#endif 786112158Sdas#ifdef IBM 787112158Sdas *e = (de - Bias - (P-1) << 2) + k; 788112158Sdas *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask); 789112158Sdas#else 790112158Sdas *e = de - Bias - (P-1) + k; 791112158Sdas *bits = P - k; 792112158Sdas#endif 793112158Sdas#ifndef Sudden_Underflow 794112158Sdas } 795112158Sdas else { 796112158Sdas *e = de - Bias - (P-1) + 1 + k; 797112158Sdas#ifdef Pack_32 798112158Sdas *bits = 32*i - hi0bits(x[i-1]); 799112158Sdas#else 800112158Sdas *bits = (i+2)*16 - hi0bits(x[i]); 801112158Sdas#endif 802112158Sdas } 803112158Sdas#endif 804112158Sdas return b; 805112158Sdas } 806112158Sdas#undef d0 807112158Sdas#undef d1 808112158Sdas 809112158Sdas CONST double 810112158Sdas#ifdef IEEE_Arith 811112158Sdasbigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; 812112158SdasCONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 813112158Sdas }; 814112158Sdas#else 815112158Sdas#ifdef IBM 816112158Sdasbigtens[] = { 1e16, 1e32, 1e64 }; 817112158SdasCONST double tinytens[] = { 1e-16, 1e-32, 1e-64 }; 818112158Sdas#else 819112158Sdasbigtens[] = { 1e16, 1e32 }; 820112158SdasCONST double tinytens[] = { 1e-16, 1e-32 }; 821112158Sdas#endif 822112158Sdas#endif 823112158Sdas 824112158Sdas CONST double 825112158Sdastens[] = { 826112158Sdas 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 827112158Sdas 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 828112158Sdas 1e20, 1e21, 1e22 829112158Sdas#ifdef VAX 830112158Sdas , 1e23, 1e24 831112158Sdas#endif 832112158Sdas }; 833112158Sdas 834112158Sdas char * 835112158Sdas#ifdef KR_headers 836112158Sdasstrcp_D2A(a, b) char *a; char *b; 837112158Sdas#else 838112158Sdasstrcp_D2A(char *a, CONST char *b) 839112158Sdas#endif 840112158Sdas{ 841112158Sdas while(*a = *b++) 842112158Sdas a++; 843112158Sdas return a; 844112158Sdas } 845112158Sdas 846112158Sdas#ifdef NO_STRING_H 847112158Sdas 848112158Sdas Char * 849112158Sdas#ifdef KR_headers 850112158Sdasmemcpy_D2A(a, b, len) Char *a; Char *b; size_t len; 851112158Sdas#else 852112158Sdasmemcpy_D2A(void *a1, void *b1, size_t len) 853112158Sdas#endif 854112158Sdas{ 855112158Sdas register char *a = (char*)a1, *ae = a + len; 856112158Sdas register char *b = (char*)b1, *a0 = a; 857112158Sdas while(a < ae) 858112158Sdas *a++ = *b++; 859112158Sdas return a0; 860112158Sdas } 861112158Sdas 862112158Sdas#endif /* NO_STRING_H */ 863