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