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) 95219557Sdas#ifdef FREE 96219557Sdas FREE((void*)v); 97219557Sdas#else 98196916Sattilio free((void*)v); 99219557Sdas#endif 100196916Sattilio else { 101196916Sattilio ACQUIRE_DTOA_LOCK(0); 102196916Sattilio v->next = freelist[v->k]; 103196916Sattilio freelist[v->k] = v; 104196916Sattilio FREE_DTOA_LOCK(0); 105196916Sattilio } 106112158Sdas } 107112158Sdas } 108112158Sdas 109112158Sdas int 110112158Sdaslo0bits 111112158Sdas#ifdef KR_headers 112112158Sdas (y) ULong *y; 113112158Sdas#else 114112158Sdas (ULong *y) 115112158Sdas#endif 116112158Sdas{ 117219557Sdas int k; 118219557Sdas ULong x = *y; 119112158Sdas 120112158Sdas if (x & 7) { 121112158Sdas if (x & 1) 122112158Sdas return 0; 123112158Sdas if (x & 2) { 124112158Sdas *y = x >> 1; 125112158Sdas return 1; 126112158Sdas } 127112158Sdas *y = x >> 2; 128112158Sdas return 2; 129112158Sdas } 130112158Sdas k = 0; 131112158Sdas if (!(x & 0xffff)) { 132112158Sdas k = 16; 133112158Sdas x >>= 16; 134112158Sdas } 135112158Sdas if (!(x & 0xff)) { 136112158Sdas k += 8; 137112158Sdas x >>= 8; 138112158Sdas } 139112158Sdas if (!(x & 0xf)) { 140112158Sdas k += 4; 141112158Sdas x >>= 4; 142112158Sdas } 143112158Sdas if (!(x & 0x3)) { 144112158Sdas k += 2; 145112158Sdas x >>= 2; 146112158Sdas } 147112158Sdas if (!(x & 1)) { 148112158Sdas k++; 149112158Sdas x >>= 1; 150165743Sdas if (!x) 151112158Sdas return 32; 152112158Sdas } 153112158Sdas *y = x; 154112158Sdas return k; 155112158Sdas } 156112158Sdas 157112158Sdas Bigint * 158112158Sdasmultadd 159112158Sdas#ifdef KR_headers 160112158Sdas (b, m, a) Bigint *b; int m, a; 161112158Sdas#else 162112158Sdas (Bigint *b, int m, int a) /* multiply by m and add a */ 163112158Sdas#endif 164112158Sdas{ 165112158Sdas int i, wds; 166112158Sdas#ifdef ULLong 167112158Sdas ULong *x; 168112158Sdas ULLong carry, y; 169112158Sdas#else 170112158Sdas ULong carry, *x, y; 171112158Sdas#ifdef Pack_32 172112158Sdas ULong xi, z; 173112158Sdas#endif 174112158Sdas#endif 175112158Sdas Bigint *b1; 176112158Sdas 177112158Sdas wds = b->wds; 178112158Sdas x = b->x; 179112158Sdas i = 0; 180112158Sdas carry = a; 181112158Sdas do { 182112158Sdas#ifdef ULLong 183112158Sdas y = *x * (ULLong)m + carry; 184112158Sdas carry = y >> 32; 185112158Sdas *x++ = y & 0xffffffffUL; 186112158Sdas#else 187112158Sdas#ifdef Pack_32 188112158Sdas xi = *x; 189112158Sdas y = (xi & 0xffff) * m + carry; 190112158Sdas z = (xi >> 16) * m + (y >> 16); 191112158Sdas carry = z >> 16; 192112158Sdas *x++ = (z << 16) + (y & 0xffff); 193112158Sdas#else 194112158Sdas y = *x * m + carry; 195112158Sdas carry = y >> 16; 196112158Sdas *x++ = y & 0xffff; 197112158Sdas#endif 198112158Sdas#endif 199112158Sdas } 200112158Sdas while(++i < wds); 201112158Sdas if (carry) { 202112158Sdas if (wds >= b->maxwds) { 203112158Sdas b1 = Balloc(b->k+1); 204112158Sdas Bcopy(b1, b); 205112158Sdas Bfree(b); 206112158Sdas b = b1; 207112158Sdas } 208112158Sdas b->x[wds++] = carry; 209112158Sdas b->wds = wds; 210112158Sdas } 211112158Sdas return b; 212112158Sdas } 213112158Sdas 214112158Sdas int 215165743Sdashi0bits_D2A 216112158Sdas#ifdef KR_headers 217219557Sdas (x) ULong x; 218112158Sdas#else 219219557Sdas (ULong x) 220112158Sdas#endif 221112158Sdas{ 222219557Sdas int k = 0; 223112158Sdas 224112158Sdas if (!(x & 0xffff0000)) { 225112158Sdas k = 16; 226112158Sdas x <<= 16; 227112158Sdas } 228112158Sdas if (!(x & 0xff000000)) { 229112158Sdas k += 8; 230112158Sdas x <<= 8; 231112158Sdas } 232112158Sdas if (!(x & 0xf0000000)) { 233112158Sdas k += 4; 234112158Sdas x <<= 4; 235112158Sdas } 236112158Sdas if (!(x & 0xc0000000)) { 237112158Sdas k += 2; 238112158Sdas x <<= 2; 239112158Sdas } 240112158Sdas if (!(x & 0x80000000)) { 241112158Sdas k++; 242112158Sdas if (!(x & 0x40000000)) 243112158Sdas return 32; 244112158Sdas } 245112158Sdas return k; 246112158Sdas } 247112158Sdas 248112158Sdas Bigint * 249112158Sdasi2b 250112158Sdas#ifdef KR_headers 251112158Sdas (i) int i; 252112158Sdas#else 253112158Sdas (int i) 254112158Sdas#endif 255112158Sdas{ 256112158Sdas Bigint *b; 257112158Sdas 258112158Sdas b = Balloc(1); 259112158Sdas b->x[0] = i; 260112158Sdas b->wds = 1; 261112158Sdas return b; 262112158Sdas } 263112158Sdas 264112158Sdas Bigint * 265112158Sdasmult 266112158Sdas#ifdef KR_headers 267112158Sdas (a, b) Bigint *a, *b; 268112158Sdas#else 269112158Sdas (Bigint *a, Bigint *b) 270112158Sdas#endif 271112158Sdas{ 272112158Sdas Bigint *c; 273112158Sdas int k, wa, wb, wc; 274112158Sdas ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; 275112158Sdas ULong y; 276112158Sdas#ifdef ULLong 277112158Sdas ULLong carry, z; 278112158Sdas#else 279112158Sdas ULong carry, z; 280112158Sdas#ifdef Pack_32 281112158Sdas ULong z2; 282112158Sdas#endif 283112158Sdas#endif 284112158Sdas 285112158Sdas if (a->wds < b->wds) { 286112158Sdas c = a; 287112158Sdas a = b; 288112158Sdas b = c; 289112158Sdas } 290112158Sdas k = a->k; 291112158Sdas wa = a->wds; 292112158Sdas wb = b->wds; 293112158Sdas wc = wa + wb; 294112158Sdas if (wc > a->maxwds) 295112158Sdas k++; 296112158Sdas c = Balloc(k); 297112158Sdas for(x = c->x, xa = x + wc; x < xa; x++) 298112158Sdas *x = 0; 299112158Sdas xa = a->x; 300112158Sdas xae = xa + wa; 301112158Sdas xb = b->x; 302112158Sdas xbe = xb + wb; 303112158Sdas xc0 = c->x; 304112158Sdas#ifdef ULLong 305112158Sdas for(; xb < xbe; xc0++) { 306112158Sdas if ( (y = *xb++) !=0) { 307112158Sdas x = xa; 308112158Sdas xc = xc0; 309112158Sdas carry = 0; 310112158Sdas do { 311112158Sdas z = *x++ * (ULLong)y + *xc + carry; 312112158Sdas carry = z >> 32; 313112158Sdas *xc++ = z & 0xffffffffUL; 314112158Sdas } 315112158Sdas while(x < xae); 316112158Sdas *xc = carry; 317112158Sdas } 318112158Sdas } 319112158Sdas#else 320112158Sdas#ifdef Pack_32 321112158Sdas for(; xb < xbe; xb++, xc0++) { 322112158Sdas if ( (y = *xb & 0xffff) !=0) { 323112158Sdas x = xa; 324112158Sdas xc = xc0; 325112158Sdas carry = 0; 326112158Sdas do { 327112158Sdas z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; 328112158Sdas carry = z >> 16; 329112158Sdas z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; 330112158Sdas carry = z2 >> 16; 331112158Sdas Storeinc(xc, z2, z); 332112158Sdas } 333112158Sdas while(x < xae); 334112158Sdas *xc = carry; 335112158Sdas } 336112158Sdas if ( (y = *xb >> 16) !=0) { 337112158Sdas x = xa; 338112158Sdas xc = xc0; 339112158Sdas carry = 0; 340112158Sdas z2 = *xc; 341112158Sdas do { 342112158Sdas z = (*x & 0xffff) * y + (*xc >> 16) + carry; 343112158Sdas carry = z >> 16; 344112158Sdas Storeinc(xc, z, z2); 345112158Sdas z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; 346112158Sdas carry = z2 >> 16; 347112158Sdas } 348112158Sdas while(x < xae); 349112158Sdas *xc = z2; 350112158Sdas } 351112158Sdas } 352112158Sdas#else 353112158Sdas for(; xb < xbe; xc0++) { 354112158Sdas if ( (y = *xb++) !=0) { 355112158Sdas x = xa; 356112158Sdas xc = xc0; 357112158Sdas carry = 0; 358112158Sdas do { 359112158Sdas z = *x++ * y + *xc + carry; 360112158Sdas carry = z >> 16; 361112158Sdas *xc++ = z & 0xffff; 362112158Sdas } 363112158Sdas while(x < xae); 364112158Sdas *xc = carry; 365112158Sdas } 366112158Sdas } 367112158Sdas#endif 368112158Sdas#endif 369112158Sdas for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ; 370112158Sdas c->wds = wc; 371112158Sdas return c; 372112158Sdas } 373112158Sdas 374112158Sdas static Bigint *p5s; 375112158Sdas 376112158Sdas Bigint * 377112158Sdaspow5mult 378112158Sdas#ifdef KR_headers 379112158Sdas (b, k) Bigint *b; int k; 380112158Sdas#else 381112158Sdas (Bigint *b, int k) 382112158Sdas#endif 383112158Sdas{ 384112158Sdas Bigint *b1, *p5, *p51; 385112158Sdas int i; 386112158Sdas static int p05[3] = { 5, 25, 125 }; 387112158Sdas 388112158Sdas if ( (i = k & 3) !=0) 389112158Sdas b = multadd(b, p05[i-1], 0); 390112158Sdas 391112158Sdas if (!(k >>= 2)) 392112158Sdas return b; 393112158Sdas if ((p5 = p5s) == 0) { 394112158Sdas /* first time */ 395112158Sdas#ifdef MULTIPLE_THREADS 396112158Sdas ACQUIRE_DTOA_LOCK(1); 397112158Sdas if (!(p5 = p5s)) { 398112158Sdas p5 = p5s = i2b(625); 399112158Sdas p5->next = 0; 400112158Sdas } 401112158Sdas FREE_DTOA_LOCK(1); 402112158Sdas#else 403112158Sdas p5 = p5s = i2b(625); 404112158Sdas p5->next = 0; 405112158Sdas#endif 406112158Sdas } 407112158Sdas for(;;) { 408112158Sdas if (k & 1) { 409112158Sdas b1 = mult(b, p5); 410112158Sdas Bfree(b); 411112158Sdas b = b1; 412112158Sdas } 413112158Sdas if (!(k >>= 1)) 414112158Sdas break; 415112158Sdas if ((p51 = p5->next) == 0) { 416112158Sdas#ifdef MULTIPLE_THREADS 417112158Sdas ACQUIRE_DTOA_LOCK(1); 418112158Sdas if (!(p51 = p5->next)) { 419112158Sdas p51 = p5->next = mult(p5,p5); 420112158Sdas p51->next = 0; 421112158Sdas } 422112158Sdas FREE_DTOA_LOCK(1); 423112158Sdas#else 424112158Sdas p51 = p5->next = mult(p5,p5); 425112158Sdas p51->next = 0; 426112158Sdas#endif 427112158Sdas } 428112158Sdas p5 = p51; 429112158Sdas } 430112158Sdas return b; 431112158Sdas } 432112158Sdas 433112158Sdas Bigint * 434112158Sdaslshift 435112158Sdas#ifdef KR_headers 436112158Sdas (b, k) Bigint *b; int k; 437112158Sdas#else 438112158Sdas (Bigint *b, int k) 439112158Sdas#endif 440112158Sdas{ 441112158Sdas int i, k1, n, n1; 442112158Sdas Bigint *b1; 443112158Sdas ULong *x, *x1, *xe, z; 444112158Sdas 445112158Sdas n = k >> kshift; 446112158Sdas k1 = b->k; 447112158Sdas n1 = n + b->wds + 1; 448112158Sdas for(i = b->maxwds; n1 > i; i <<= 1) 449112158Sdas k1++; 450112158Sdas b1 = Balloc(k1); 451112158Sdas x1 = b1->x; 452112158Sdas for(i = 0; i < n; i++) 453112158Sdas *x1++ = 0; 454112158Sdas x = b->x; 455112158Sdas xe = x + b->wds; 456112158Sdas if (k &= kmask) { 457112158Sdas#ifdef Pack_32 458112158Sdas k1 = 32 - k; 459112158Sdas z = 0; 460112158Sdas do { 461112158Sdas *x1++ = *x << k | z; 462112158Sdas z = *x++ >> k1; 463112158Sdas } 464112158Sdas while(x < xe); 465112158Sdas if ((*x1 = z) !=0) 466112158Sdas ++n1; 467112158Sdas#else 468112158Sdas k1 = 16 - k; 469112158Sdas z = 0; 470112158Sdas do { 471112158Sdas *x1++ = *x << k & 0xffff | z; 472112158Sdas z = *x++ >> k1; 473112158Sdas } 474112158Sdas while(x < xe); 475112158Sdas if (*x1 = z) 476112158Sdas ++n1; 477112158Sdas#endif 478112158Sdas } 479112158Sdas else do 480112158Sdas *x1++ = *x++; 481112158Sdas while(x < xe); 482112158Sdas b1->wds = n1 - 1; 483112158Sdas Bfree(b); 484112158Sdas return b1; 485112158Sdas } 486112158Sdas 487112158Sdas int 488112158Sdascmp 489112158Sdas#ifdef KR_headers 490112158Sdas (a, b) Bigint *a, *b; 491112158Sdas#else 492112158Sdas (Bigint *a, Bigint *b) 493112158Sdas#endif 494112158Sdas{ 495112158Sdas ULong *xa, *xa0, *xb, *xb0; 496112158Sdas int i, j; 497112158Sdas 498112158Sdas i = a->wds; 499112158Sdas j = b->wds; 500112158Sdas#ifdef DEBUG 501112158Sdas if (i > 1 && !a->x[i-1]) 502112158Sdas Bug("cmp called with a->x[a->wds-1] == 0"); 503112158Sdas if (j > 1 && !b->x[j-1]) 504112158Sdas Bug("cmp called with b->x[b->wds-1] == 0"); 505112158Sdas#endif 506112158Sdas if (i -= j) 507112158Sdas return i; 508112158Sdas xa0 = a->x; 509112158Sdas xa = xa0 + j; 510112158Sdas xb0 = b->x; 511112158Sdas xb = xb0 + j; 512112158Sdas for(;;) { 513112158Sdas if (*--xa != *--xb) 514112158Sdas return *xa < *xb ? -1 : 1; 515112158Sdas if (xa <= xa0) 516112158Sdas break; 517112158Sdas } 518112158Sdas return 0; 519112158Sdas } 520112158Sdas 521112158Sdas Bigint * 522112158Sdasdiff 523112158Sdas#ifdef KR_headers 524112158Sdas (a, b) Bigint *a, *b; 525112158Sdas#else 526112158Sdas (Bigint *a, Bigint *b) 527112158Sdas#endif 528112158Sdas{ 529112158Sdas Bigint *c; 530112158Sdas int i, wa, wb; 531112158Sdas ULong *xa, *xae, *xb, *xbe, *xc; 532112158Sdas#ifdef ULLong 533112158Sdas ULLong borrow, y; 534112158Sdas#else 535112158Sdas ULong borrow, y; 536112158Sdas#ifdef Pack_32 537112158Sdas ULong z; 538112158Sdas#endif 539112158Sdas#endif 540112158Sdas 541112158Sdas i = cmp(a,b); 542112158Sdas if (!i) { 543112158Sdas c = Balloc(0); 544112158Sdas c->wds = 1; 545112158Sdas c->x[0] = 0; 546112158Sdas return c; 547112158Sdas } 548112158Sdas if (i < 0) { 549112158Sdas c = a; 550112158Sdas a = b; 551112158Sdas b = c; 552112158Sdas i = 1; 553112158Sdas } 554112158Sdas else 555112158Sdas i = 0; 556112158Sdas c = Balloc(a->k); 557112158Sdas c->sign = i; 558112158Sdas wa = a->wds; 559112158Sdas xa = a->x; 560112158Sdas xae = xa + wa; 561112158Sdas wb = b->wds; 562112158Sdas xb = b->x; 563112158Sdas xbe = xb + wb; 564112158Sdas xc = c->x; 565112158Sdas borrow = 0; 566112158Sdas#ifdef ULLong 567112158Sdas do { 568112158Sdas y = (ULLong)*xa++ - *xb++ - borrow; 569112158Sdas borrow = y >> 32 & 1UL; 570112158Sdas *xc++ = y & 0xffffffffUL; 571112158Sdas } 572112158Sdas while(xb < xbe); 573112158Sdas while(xa < xae) { 574112158Sdas y = *xa++ - borrow; 575112158Sdas borrow = y >> 32 & 1UL; 576112158Sdas *xc++ = y & 0xffffffffUL; 577112158Sdas } 578112158Sdas#else 579112158Sdas#ifdef Pack_32 580112158Sdas do { 581112158Sdas y = (*xa & 0xffff) - (*xb & 0xffff) - borrow; 582112158Sdas borrow = (y & 0x10000) >> 16; 583112158Sdas z = (*xa++ >> 16) - (*xb++ >> 16) - borrow; 584112158Sdas borrow = (z & 0x10000) >> 16; 585112158Sdas Storeinc(xc, z, y); 586112158Sdas } 587112158Sdas while(xb < xbe); 588112158Sdas while(xa < xae) { 589112158Sdas y = (*xa & 0xffff) - borrow; 590112158Sdas borrow = (y & 0x10000) >> 16; 591112158Sdas z = (*xa++ >> 16) - borrow; 592112158Sdas borrow = (z & 0x10000) >> 16; 593112158Sdas Storeinc(xc, z, y); 594112158Sdas } 595112158Sdas#else 596112158Sdas do { 597112158Sdas y = *xa++ - *xb++ - borrow; 598112158Sdas borrow = (y & 0x10000) >> 16; 599112158Sdas *xc++ = y & 0xffff; 600112158Sdas } 601112158Sdas while(xb < xbe); 602112158Sdas while(xa < xae) { 603112158Sdas y = *xa++ - borrow; 604112158Sdas borrow = (y & 0x10000) >> 16; 605112158Sdas *xc++ = y & 0xffff; 606112158Sdas } 607112158Sdas#endif 608112158Sdas#endif 609112158Sdas while(!*--xc) 610112158Sdas wa--; 611112158Sdas c->wds = wa; 612112158Sdas return c; 613112158Sdas } 614112158Sdas 615112158Sdas double 616112158Sdasb2d 617112158Sdas#ifdef KR_headers 618112158Sdas (a, e) Bigint *a; int *e; 619112158Sdas#else 620112158Sdas (Bigint *a, int *e) 621112158Sdas#endif 622112158Sdas{ 623112158Sdas ULong *xa, *xa0, w, y, z; 624112158Sdas int k; 625219557Sdas U d; 626112158Sdas#ifdef VAX 627112158Sdas ULong d0, d1; 628112158Sdas#else 629219557Sdas#define d0 word0(&d) 630219557Sdas#define d1 word1(&d) 631112158Sdas#endif 632112158Sdas 633112158Sdas xa0 = a->x; 634112158Sdas xa = xa0 + a->wds; 635112158Sdas y = *--xa; 636112158Sdas#ifdef DEBUG 637112158Sdas if (!y) Bug("zero y in b2d"); 638112158Sdas#endif 639112158Sdas k = hi0bits(y); 640112158Sdas *e = 32 - k; 641112158Sdas#ifdef Pack_32 642112158Sdas if (k < Ebits) { 643219557Sdas d0 = Exp_1 | y >> (Ebits - k); 644112158Sdas w = xa > xa0 ? *--xa : 0; 645219557Sdas d1 = y << ((32-Ebits) + k) | w >> (Ebits - k); 646112158Sdas goto ret_d; 647112158Sdas } 648112158Sdas z = xa > xa0 ? *--xa : 0; 649112158Sdas if (k -= Ebits) { 650219557Sdas d0 = Exp_1 | y << k | z >> (32 - k); 651112158Sdas y = xa > xa0 ? *--xa : 0; 652219557Sdas d1 = z << k | y >> (32 - k); 653112158Sdas } 654112158Sdas else { 655112158Sdas d0 = Exp_1 | y; 656112158Sdas d1 = z; 657112158Sdas } 658112158Sdas#else 659112158Sdas if (k < Ebits + 16) { 660112158Sdas z = xa > xa0 ? *--xa : 0; 661112158Sdas d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; 662112158Sdas w = xa > xa0 ? *--xa : 0; 663112158Sdas y = xa > xa0 ? *--xa : 0; 664112158Sdas d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; 665112158Sdas goto ret_d; 666112158Sdas } 667112158Sdas z = xa > xa0 ? *--xa : 0; 668112158Sdas w = xa > xa0 ? *--xa : 0; 669112158Sdas k -= Ebits + 16; 670112158Sdas d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; 671112158Sdas y = xa > xa0 ? *--xa : 0; 672112158Sdas d1 = w << k + 16 | y << k; 673112158Sdas#endif 674112158Sdas ret_d: 675112158Sdas#ifdef VAX 676219557Sdas word0(&d) = d0 >> 16 | d0 << 16; 677219557Sdas word1(&d) = d1 >> 16 | d1 << 16; 678112158Sdas#endif 679219557Sdas return dval(&d); 680112158Sdas } 681112158Sdas#undef d0 682112158Sdas#undef d1 683112158Sdas 684112158Sdas Bigint * 685112158Sdasd2b 686112158Sdas#ifdef KR_headers 687219557Sdas (dd, e, bits) double dd; int *e, *bits; 688112158Sdas#else 689219557Sdas (double dd, int *e, int *bits) 690112158Sdas#endif 691112158Sdas{ 692112158Sdas Bigint *b; 693219557Sdas U d; 694165743Sdas#ifndef Sudden_Underflow 695165743Sdas int i; 696165743Sdas#endif 697165743Sdas int de, k; 698112158Sdas ULong *x, y, z; 699112158Sdas#ifdef VAX 700112158Sdas ULong d0, d1; 701112158Sdas#else 702219557Sdas#define d0 word0(&d) 703219557Sdas#define d1 word1(&d) 704112158Sdas#endif 705219557Sdas d.d = dd; 706219557Sdas#ifdef VAX 707219557Sdas d0 = word0(&d) >> 16 | word0(&d) << 16; 708219557Sdas d1 = word1(&d) >> 16 | word1(&d) << 16; 709219557Sdas#endif 710112158Sdas 711112158Sdas#ifdef Pack_32 712112158Sdas b = Balloc(1); 713112158Sdas#else 714112158Sdas b = Balloc(2); 715112158Sdas#endif 716112158Sdas x = b->x; 717112158Sdas 718112158Sdas z = d0 & Frac_mask; 719112158Sdas d0 &= 0x7fffffff; /* clear sign bit, which we ignore */ 720112158Sdas#ifdef Sudden_Underflow 721112158Sdas de = (int)(d0 >> Exp_shift); 722112158Sdas#ifndef IBM 723112158Sdas z |= Exp_msk11; 724112158Sdas#endif 725112158Sdas#else 726112158Sdas if ( (de = (int)(d0 >> Exp_shift)) !=0) 727112158Sdas z |= Exp_msk1; 728112158Sdas#endif 729112158Sdas#ifdef Pack_32 730112158Sdas if ( (y = d1) !=0) { 731112158Sdas if ( (k = lo0bits(&y)) !=0) { 732219557Sdas x[0] = y | z << (32 - k); 733112158Sdas z >>= k; 734112158Sdas } 735112158Sdas else 736112158Sdas x[0] = y; 737165743Sdas#ifndef Sudden_Underflow 738165743Sdas i = 739165743Sdas#endif 740165743Sdas b->wds = (x[1] = z) !=0 ? 2 : 1; 741112158Sdas } 742112158Sdas else { 743112158Sdas k = lo0bits(&z); 744112158Sdas x[0] = z; 745165743Sdas#ifndef Sudden_Underflow 746165743Sdas i = 747165743Sdas#endif 748165743Sdas b->wds = 1; 749112158Sdas k += 32; 750112158Sdas } 751112158Sdas#else 752112158Sdas if ( (y = d1) !=0) { 753112158Sdas if ( (k = lo0bits(&y)) !=0) 754112158Sdas if (k >= 16) { 755112158Sdas x[0] = y | z << 32 - k & 0xffff; 756112158Sdas x[1] = z >> k - 16 & 0xffff; 757112158Sdas x[2] = z >> k; 758112158Sdas i = 2; 759112158Sdas } 760112158Sdas else { 761112158Sdas x[0] = y & 0xffff; 762112158Sdas x[1] = y >> 16 | z << 16 - k & 0xffff; 763112158Sdas x[2] = z >> k & 0xffff; 764112158Sdas x[3] = z >> k+16; 765112158Sdas i = 3; 766112158Sdas } 767112158Sdas else { 768112158Sdas x[0] = y & 0xffff; 769112158Sdas x[1] = y >> 16; 770112158Sdas x[2] = z & 0xffff; 771112158Sdas x[3] = z >> 16; 772112158Sdas i = 3; 773112158Sdas } 774112158Sdas } 775112158Sdas else { 776112158Sdas#ifdef DEBUG 777112158Sdas if (!z) 778112158Sdas Bug("Zero passed to d2b"); 779112158Sdas#endif 780112158Sdas k = lo0bits(&z); 781112158Sdas if (k >= 16) { 782112158Sdas x[0] = z; 783112158Sdas i = 0; 784112158Sdas } 785112158Sdas else { 786112158Sdas x[0] = z & 0xffff; 787112158Sdas x[1] = z >> 16; 788112158Sdas i = 1; 789112158Sdas } 790112158Sdas k += 32; 791112158Sdas } 792112158Sdas while(!x[i]) 793112158Sdas --i; 794112158Sdas b->wds = i + 1; 795112158Sdas#endif 796112158Sdas#ifndef Sudden_Underflow 797112158Sdas if (de) { 798112158Sdas#endif 799112158Sdas#ifdef IBM 800112158Sdas *e = (de - Bias - (P-1) << 2) + k; 801219557Sdas *bits = 4*P + 8 - k - hi0bits(word0(&d) & Frac_mask); 802112158Sdas#else 803112158Sdas *e = de - Bias - (P-1) + k; 804112158Sdas *bits = P - k; 805112158Sdas#endif 806112158Sdas#ifndef Sudden_Underflow 807112158Sdas } 808112158Sdas else { 809112158Sdas *e = de - Bias - (P-1) + 1 + k; 810112158Sdas#ifdef Pack_32 811112158Sdas *bits = 32*i - hi0bits(x[i-1]); 812112158Sdas#else 813112158Sdas *bits = (i+2)*16 - hi0bits(x[i]); 814112158Sdas#endif 815112158Sdas } 816112158Sdas#endif 817112158Sdas return b; 818112158Sdas } 819112158Sdas#undef d0 820112158Sdas#undef d1 821112158Sdas 822112158Sdas CONST double 823112158Sdas#ifdef IEEE_Arith 824112158Sdasbigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; 825112158SdasCONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 826112158Sdas }; 827112158Sdas#else 828112158Sdas#ifdef IBM 829112158Sdasbigtens[] = { 1e16, 1e32, 1e64 }; 830112158SdasCONST double tinytens[] = { 1e-16, 1e-32, 1e-64 }; 831112158Sdas#else 832112158Sdasbigtens[] = { 1e16, 1e32 }; 833112158SdasCONST double tinytens[] = { 1e-16, 1e-32 }; 834112158Sdas#endif 835112158Sdas#endif 836112158Sdas 837112158Sdas CONST double 838112158Sdastens[] = { 839112158Sdas 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 840112158Sdas 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 841112158Sdas 1e20, 1e21, 1e22 842112158Sdas#ifdef VAX 843112158Sdas , 1e23, 1e24 844112158Sdas#endif 845112158Sdas }; 846112158Sdas 847112158Sdas char * 848112158Sdas#ifdef KR_headers 849112158Sdasstrcp_D2A(a, b) char *a; char *b; 850112158Sdas#else 851112158Sdasstrcp_D2A(char *a, CONST char *b) 852112158Sdas#endif 853112158Sdas{ 854219557Sdas while((*a = *b++)) 855112158Sdas a++; 856112158Sdas return a; 857112158Sdas } 858112158Sdas 859112158Sdas#ifdef NO_STRING_H 860112158Sdas 861112158Sdas Char * 862112158Sdas#ifdef KR_headers 863112158Sdasmemcpy_D2A(a, b, len) Char *a; Char *b; size_t len; 864112158Sdas#else 865112158Sdasmemcpy_D2A(void *a1, void *b1, size_t len) 866112158Sdas#endif 867112158Sdas{ 868219557Sdas char *a = (char*)a1, *ae = a + len; 869219557Sdas char *b = (char*)b1, *a0 = a; 870112158Sdas while(a < ae) 871112158Sdas *a++ = *b++; 872112158Sdas return a0; 873112158Sdas } 874112158Sdas 875112158Sdas#endif /* NO_STRING_H */ 876