misc.c revision 196916
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); 58196916Sattilio /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */ 59196916Sattilio /* but this case seems very unlikely. */ 60196916Sattilio if (k <= Kmax && (rv = freelist[k]) !=0) { 61112158Sdas freelist[k] = rv->next; 62112158Sdas } 63112158Sdas else { 64112158Sdas x = 1 << k; 65112158Sdas#ifdef Omit_Private_Memory 66112158Sdas rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong)); 67112158Sdas#else 68112158Sdas len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1) 69112158Sdas /sizeof(double); 70196916Sattilio if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) { 71112158Sdas rv = (Bigint*)pmem_next; 72112158Sdas pmem_next += len; 73112158Sdas } 74112158Sdas else 75112158Sdas rv = (Bigint*)MALLOC(len*sizeof(double)); 76112158Sdas#endif 77112158Sdas rv->k = k; 78112158Sdas rv->maxwds = x; 79112158Sdas } 80112158Sdas FREE_DTOA_LOCK(0); 81112158Sdas rv->sign = rv->wds = 0; 82112158Sdas return rv; 83112158Sdas } 84112158Sdas 85112158Sdas void 86112158SdasBfree 87112158Sdas#ifdef KR_headers 88112158Sdas (v) Bigint *v; 89112158Sdas#else 90112158Sdas (Bigint *v) 91112158Sdas#endif 92112158Sdas{ 93112158Sdas if (v) { 94196916Sattilio if (v->k > Kmax) 95196916Sattilio free((void*)v); 96196916Sattilio else { 97196916Sattilio ACQUIRE_DTOA_LOCK(0); 98196916Sattilio v->next = freelist[v->k]; 99196916Sattilio freelist[v->k] = v; 100196916Sattilio FREE_DTOA_LOCK(0); 101196916Sattilio } 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; 146165743Sdas if (!x) 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 211165743Sdashi0bits_D2A 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; 689165743Sdas#ifndef Sudden_Underflow 690165743Sdas int i; 691165743Sdas#endif 692165743Sdas int de, k; 693112158Sdas ULong *x, y, z; 694112158Sdas#ifdef VAX 695112158Sdas ULong d0, d1; 696112158Sdas d0 = word0(d) >> 16 | word0(d) << 16; 697112158Sdas d1 = word1(d) >> 16 | word1(d) << 16; 698112158Sdas#else 699112158Sdas#define d0 word0(d) 700112158Sdas#define d1 word1(d) 701112158Sdas#endif 702112158Sdas 703112158Sdas#ifdef Pack_32 704112158Sdas b = Balloc(1); 705112158Sdas#else 706112158Sdas b = Balloc(2); 707112158Sdas#endif 708112158Sdas x = b->x; 709112158Sdas 710112158Sdas z = d0 & Frac_mask; 711112158Sdas d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ 712112158Sdas#ifdef Sudden_Underflow 713112158Sdas de = (int)(d0 >> Exp_shift); 714112158Sdas#ifndef IBM 715112158Sdas z |= Exp_msk11; 716112158Sdas#endif 717112158Sdas#else 718112158Sdas if ( (de = (int)(d0 >> Exp_shift)) !=0) 719112158Sdas z |= Exp_msk1; 720112158Sdas#endif 721112158Sdas#ifdef Pack_32 722112158Sdas if ( (y = d1) !=0) { 723112158Sdas if ( (k = lo0bits(&y)) !=0) { 724112158Sdas x[0] = y | z << 32 - k; 725112158Sdas z >>= k; 726112158Sdas } 727112158Sdas else 728112158Sdas x[0] = y; 729165743Sdas#ifndef Sudden_Underflow 730165743Sdas i = 731165743Sdas#endif 732165743Sdas b->wds = (x[1] = z) !=0 ? 2 : 1; 733112158Sdas } 734112158Sdas else { 735112158Sdas#ifdef DEBUG 736112158Sdas if (!z) 737112158Sdas Bug("Zero passed to d2b"); 738112158Sdas#endif 739112158Sdas k = lo0bits(&z); 740112158Sdas x[0] = z; 741165743Sdas#ifndef Sudden_Underflow 742165743Sdas i = 743165743Sdas#endif 744165743Sdas b->wds = 1; 745112158Sdas k += 32; 746112158Sdas } 747112158Sdas#else 748112158Sdas if ( (y = d1) !=0) { 749112158Sdas if ( (k = lo0bits(&y)) !=0) 750112158Sdas if (k >= 16) { 751112158Sdas x[0] = y | z << 32 - k & 0xffff; 752112158Sdas x[1] = z >> k - 16 & 0xffff; 753112158Sdas x[2] = z >> k; 754112158Sdas i = 2; 755112158Sdas } 756112158Sdas else { 757112158Sdas x[0] = y & 0xffff; 758112158Sdas x[1] = y >> 16 | z << 16 - k & 0xffff; 759112158Sdas x[2] = z >> k & 0xffff; 760112158Sdas x[3] = z >> k+16; 761112158Sdas i = 3; 762112158Sdas } 763112158Sdas else { 764112158Sdas x[0] = y & 0xffff; 765112158Sdas x[1] = y >> 16; 766112158Sdas x[2] = z & 0xffff; 767112158Sdas x[3] = z >> 16; 768112158Sdas i = 3; 769112158Sdas } 770112158Sdas } 771112158Sdas else { 772112158Sdas#ifdef DEBUG 773112158Sdas if (!z) 774112158Sdas Bug("Zero passed to d2b"); 775112158Sdas#endif 776112158Sdas k = lo0bits(&z); 777112158Sdas if (k >= 16) { 778112158Sdas x[0] = z; 779112158Sdas i = 0; 780112158Sdas } 781112158Sdas else { 782112158Sdas x[0] = z & 0xffff; 783112158Sdas x[1] = z >> 16; 784112158Sdas i = 1; 785112158Sdas } 786112158Sdas k += 32; 787112158Sdas } 788112158Sdas while(!x[i]) 789112158Sdas --i; 790112158Sdas b->wds = i + 1; 791112158Sdas#endif 792112158Sdas#ifndef Sudden_Underflow 793112158Sdas if (de) { 794112158Sdas#endif 795112158Sdas#ifdef IBM 796112158Sdas *e = (de - Bias - (P-1) << 2) + k; 797112158Sdas *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask); 798112158Sdas#else 799112158Sdas *e = de - Bias - (P-1) + k; 800112158Sdas *bits = P - k; 801112158Sdas#endif 802112158Sdas#ifndef Sudden_Underflow 803112158Sdas } 804112158Sdas else { 805112158Sdas *e = de - Bias - (P-1) + 1 + k; 806112158Sdas#ifdef Pack_32 807112158Sdas *bits = 32*i - hi0bits(x[i-1]); 808112158Sdas#else 809112158Sdas *bits = (i+2)*16 - hi0bits(x[i]); 810112158Sdas#endif 811112158Sdas } 812112158Sdas#endif 813112158Sdas return b; 814112158Sdas } 815112158Sdas#undef d0 816112158Sdas#undef d1 817112158Sdas 818112158Sdas CONST double 819112158Sdas#ifdef IEEE_Arith 820112158Sdasbigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; 821112158SdasCONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 822112158Sdas }; 823112158Sdas#else 824112158Sdas#ifdef IBM 825112158Sdasbigtens[] = { 1e16, 1e32, 1e64 }; 826112158SdasCONST double tinytens[] = { 1e-16, 1e-32, 1e-64 }; 827112158Sdas#else 828112158Sdasbigtens[] = { 1e16, 1e32 }; 829112158SdasCONST double tinytens[] = { 1e-16, 1e-32 }; 830112158Sdas#endif 831112158Sdas#endif 832112158Sdas 833112158Sdas CONST double 834112158Sdastens[] = { 835112158Sdas 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 836112158Sdas 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 837112158Sdas 1e20, 1e21, 1e22 838112158Sdas#ifdef VAX 839112158Sdas , 1e23, 1e24 840112158Sdas#endif 841112158Sdas }; 842112158Sdas 843112158Sdas char * 844112158Sdas#ifdef KR_headers 845112158Sdasstrcp_D2A(a, b) char *a; char *b; 846112158Sdas#else 847112158Sdasstrcp_D2A(char *a, CONST char *b) 848112158Sdas#endif 849112158Sdas{ 850112158Sdas while(*a = *b++) 851112158Sdas a++; 852112158Sdas return a; 853112158Sdas } 854112158Sdas 855112158Sdas#ifdef NO_STRING_H 856112158Sdas 857112158Sdas Char * 858112158Sdas#ifdef KR_headers 859112158Sdasmemcpy_D2A(a, b, len) Char *a; Char *b; size_t len; 860112158Sdas#else 861112158Sdasmemcpy_D2A(void *a1, void *b1, size_t len) 862112158Sdas#endif 863112158Sdas{ 864112158Sdas register char *a = (char*)a1, *ae = a + len; 865112158Sdas register char *b = (char*)b1, *a0 = a; 866112158Sdas while(a < ae) 867112158Sdas *a++ = *b++; 868112158Sdas return a0; 869112158Sdas } 870112158Sdas 871112158Sdas#endif /* NO_STRING_H */ 872