gonzui

    1: /****************************************************************
    2: 
    3: The author of this software is David M. Gay.
    4: 
    5: Copyright (C) 1998, 1999 by Lucent Technologies
    6: All Rights Reserved
    7: 
    8: Permission to use, copy, modify, and distribute this software and
    9: its documentation for any purpose and without fee is hereby
   10: granted, provided that the above copyright notice appear in all
   11: copies and that both that the copyright notice and this
   12: permission notice and warranty disclaimer appear in supporting
   13: documentation, and that the name of Lucent or any of its entities
   14: not be used in advertising or publicity pertaining to
   15: distribution of the software without specific, written prior
   16: permission.
   17: 
   18: LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
   19: INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
   20: IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
   21: SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
   22: WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
   23: IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
   24: ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
   25: THIS SOFTWARE.
   26: 
   27: ****************************************************************/
   28: 
   29: /* Please send bug reports to David M. Gay (dmg at acm dot org,
   30:  * with " at " changed at "@" and " dot " changed to ".").      */
   31: 
   32: #include "gdtoaimp.h"
   33: 
   34:  static Bigint *freelist[Kmax+1];
   35: #ifndef Omit_Private_Memory
   36: #ifndef PRIVATE_MEM
   37: //#define PRIVATE_MEM   2304
   38: #define PRIVATE_MEM     (32 + 32 + 40)      /* k=0,1,2 */
   39: #endif
   40: #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
   41: static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
   42: #endif
   43: 
   44:  Bigint *
   45: Balloc
   46: #ifdef KR_headers
   47:         (k) int k;
   48: #else
   49:         (int k)
   50: #endif
   51: {
   52:         int x;
   53:         Bigint *rv;
   54: #ifndef Omit_Private_Memory
   55:         unsigned int len;
   56: #endif
   57: 
   58:         ACQUIRE_DTOA_LOCK(0);
   59:         /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
   60:         /* but this case seems very unlikely. */
   61:         if (k <= Kmax && (rv = freelist[k]) !=0) {
   62:                 freelist[k] = rv->next;
   63:                 }
   64:         else {
   65:                 x = 1 << k;
   66: #ifdef Omit_Private_Memory
   67:                 rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
   68:                 if (rv == NULL)
   69:                         return (NULL);
   70: #else
   71:                 len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
   72:                         /sizeof(double);
   73:                 if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
   74:                         rv = (Bigint*)pmem_next;
   75:                         pmem_next += len;
   76:                         }
   77:                 else {
   78:                         rv = (Bigint*)MALLOC(len*sizeof(double));
   79:                         if (rv == NULL)
   80:                                 return (NULL);
   81:                 }
   82: #endif
   83:                 rv->k = k;
   84:                 rv->maxwds = x;
   85:                 }
   86:         FREE_DTOA_LOCK(0);
   87:         rv->sign = rv->wds = 0;
   88:         return rv;
   89:         }
   90: 
   91:  void
   92: Bfree
   93: #ifdef KR_headers
   94:         (v) Bigint *v;
   95: #else
   96:         (Bigint *v)
   97: #endif
   98: {
   99:         if (v) {
  100:                 if (v->k > Kmax)
  101: #ifdef FREE
  102:                         FREE((void*)v);
  103: #else
  104:                         free((void*)v);
  105: #endif
  106:                 else {
  107:                         ACQUIRE_DTOA_LOCK(0);
  108:                         v->next = freelist[v->k];
  109:                         freelist[v->k] = v;
  110:                         FREE_DTOA_LOCK(0);
  111:                         }
  112:                 }
  113:         }
  114: 
  115:  int
  116: lo0bits
  117: #ifdef KR_headers
  118:         (y) ULong *y;
  119: #else
  120:         (ULong *y)
  121: #endif
  122: {
  123:         int k;
  124:         ULong x = *y;
  125: 
  126:         if (x & 7) {
  127:                 if (x & 1)
  128:                         return 0;
  129:                 if (x & 2) {
  130:                         *y = x >> 1;
  131:                         return 1;
  132:                         }
  133:                 *y = x >> 2;
  134:                 return 2;
  135:                 }
  136:         k = 0;
  137:         if (!(x & 0xffff)) {
  138:                 k = 16;
  139:                 x >>= 16;
  140:                 }
  141:         if (!(x & 0xff)) {
  142:                 k += 8;
  143:                 x >>= 8;
  144:                 }
  145:         if (!(x & 0xf)) {
  146:                 k += 4;
  147:                 x >>= 4;
  148:                 }
  149:         if (!(x & 0x3)) {
  150:                 k += 2;
  151:                 x >>= 2;
  152:                 }
  153:         if (!(x & 1)) {
  154:                 k++;
  155:                 x >>= 1;
  156:                 if (!x)
  157:                         return 32;
  158:                 }
  159:         *y = x;
  160:         return k;
  161:         }
  162: 
  163:  Bigint *
  164: multadd
  165: #ifdef KR_headers
  166:         (b, m, a) Bigint *b; int m, a;
  167: #else
  168:         (Bigint *b, int m, int a)      /* multiply by m and add a */
  169: #endif
  170: {
  171:         int i, wds;
  172: #ifdef ULLong
  173:         ULong *x;
  174:         ULLong carry, y;
  175: #else
  176:         ULong carry, *x, y;
  177: #ifdef Pack_32
  178:         ULong xi, z;
  179: #endif
  180: #endif
  181:         Bigint *b1;
  182: 
  183:         wds = b->wds;
  184:         x = b->x;
  185:         i = 0;
  186:         carry = a;
  187:         do {
  188: #ifdef ULLong
  189:                 y = *x * (ULLong)m + carry;
  190:                 carry = y >> 32;
  191:                 *x++ = y & 0xffffffffUL;
  192: #else
  193: #ifdef Pack_32
  194:                 xi = *x;
  195:                 y = (xi & 0xffff) * m + carry;
  196:                 z = (xi >> 16) * m + (y >> 16);
  197:                 carry = z >> 16;
  198:                 *x++ = (z << 16) + (y & 0xffff);
  199: #else
  200:                 y = *x * m + carry;
  201:                 carry = y >> 16;
  202:                 *x++ = y & 0xffff;
  203: #endif
  204: #endif
  205:                 }
  206:                 while(++i < wds);
  207:         if (carry) {
  208:                 if (wds >= b->maxwds) {
  209:                         b1 = Balloc(b->k+1);
  210:                         if (b1 == NULL)
  211:                                 return (NULL);
  212:                         Bcopy(b1, b);
  213:                         Bfree(b);
  214:                         b = b1;
  215:                         }
  216:                 b->x[wds++] = carry;
  217:                 b->wds = wds;
  218:                 }
  219:         return b;
  220:         }
  221: 
  222:  int
  223: hi0bits_D2A
  224: #ifdef KR_headers
  225:         (x) ULong x;
  226: #else
  227:         (ULong x)
  228: #endif
  229: {
  230:         int k = 0;
  231: 
  232:         if (!(x & 0xffff0000)) {
  233:                 k = 16;
  234:                 x <<= 16;
  235:                 }
  236:         if (!(x & 0xff000000)) {
  237:                 k += 8;
  238:                 x <<= 8;
  239:                 }
  240:         if (!(x & 0xf0000000)) {
  241:                 k += 4;
  242:                 x <<= 4;
  243:                 }
  244:         if (!(x & 0xc0000000)) {
  245:                 k += 2;
  246:                 x <<= 2;
  247:                 }
  248:         if (!(x & 0x80000000)) {
  249:                 k++;
  250:                 if (!(x & 0x40000000))
  251:                         return 32;
  252:                 }
  253:         return k;
  254:         }
  255: 
  256:  Bigint *
  257: i2b
  258: #ifdef KR_headers
  259:         (i) int i;
  260: #else
  261:         (int i)
  262: #endif
  263: {
  264:         Bigint *b;
  265: 
  266:         b = Balloc(1);
  267:         if (b == NULL)
  268:                 return (NULL);
  269:         b->x[0] = i;
  270:         b->wds = 1;
  271:         return b;
  272:         }
  273: 
  274:  Bigint *
  275: mult
  276: #ifdef KR_headers
  277:         (a, b) Bigint *a, *b;
  278: #else
  279:         (Bigint *a, Bigint *b)
  280: #endif
  281: {
  282:         Bigint *c;
  283:         int k, wa, wb, wc;
  284:         ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
  285:         ULong y;
  286: #ifdef ULLong
  287:         ULLong carry, z;
  288: #else
  289:         ULong carry, z;
  290: #ifdef Pack_32
  291:         ULong z2;
  292: #endif
  293: #endif
  294: 
  295:         if (a->wds < b->wds) {
  296:                 c = a;
  297:                 a = b;
  298:                 b = c;
  299:                 }
  300:         k = a->k;
  301:         wa = a->wds;
  302:         wb = b->wds;
  303:         wc = wa + wb;
  304:         if (wc > a->maxwds)
  305:                 k++;
  306:         c = Balloc(k);
  307:         if (c == NULL)
  308:                 return (NULL);
  309:         for(x = c->x, xa = x + wc; x < xa; x++)
  310:                 *x = 0;
  311:         xa = a->x;
  312:         xae = xa + wa;
  313:         xb = b->x;
  314:         xbe = xb + wb;
  315:         xc0 = c->x;
  316: #ifdef ULLong
  317:         for(; xb < xbe; xc0++) {
  318:                 if ( (y = *xb++) !=0) {
  319:                         x = xa;
  320:                         xc = xc0;
  321:                         carry = 0;
  322:                         do {
  323:                                 z = *x++ * (ULLong)y + *xc + carry;
  324:                                 carry = z >> 32;
  325:                                 *xc++ = z & 0xffffffffUL;
  326:                                 }
  327:                                 while(x < xae);
  328:                         *xc = carry;
  329:                         }
  330:                 }
  331: #else
  332: #ifdef Pack_32
  333:         for(; xb < xbe; xb++, xc0++) {
  334:                 if ( (y = *xb & 0xffff) !=0) {
  335:                         x = xa;
  336:                         xc = xc0;
  337:                         carry = 0;
  338:                         do {
  339:                                 z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
  340:                                 carry = z >> 16;
  341:                                 z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
  342:                                 carry = z2 >> 16;
  343:                                 Storeinc(xc, z2, z);
  344:                                 }
  345:                                 while(x < xae);
  346:                         *xc = carry;
  347:                         }
  348:                 if ( (y = *xb >> 16) !=0) {
  349:                         x = xa;
  350:                         xc = xc0;
  351:                         carry = 0;
  352:                         z2 = *xc;
  353:                         do {
  354:                                 z = (*x & 0xffff) * y + (*xc >> 16) + carry;
  355:                                 carry = z >> 16;
  356:                                 Storeinc(xc, z, z2);
  357:                                 z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
  358:                                 carry = z2 >> 16;
  359:                                 }
  360:                                 while(x < xae);
  361:                         *xc = z2;
  362:                         }
  363:                 }
  364: #else
  365:         for(; xb < xbe; xc0++) {
  366:                 if ( (y = *xb++) !=0) {
  367:                         x = xa;
  368:                         xc = xc0;
  369:                         carry = 0;
  370:                         do {
  371:                                 z = *x++ * y + *xc + carry;
  372:                                 carry = z >> 16;
  373:                                 *xc++ = z & 0xffff;
  374:                                 }
  375:                                 while(x < xae);
  376:                         *xc = carry;
  377:                         }
  378:                 }
  379: #endif
  380: #endif
  381:         for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
  382:         c->wds = wc;
  383:         return c;
  384:         }
  385: 
  386:  static Bigint *p5s;
  387: 
  388:  Bigint *
  389: pow5mult
  390: #ifdef KR_headers
  391:         (b, k) Bigint *b; int k;
  392: #else
  393:         (Bigint *b, int k)
  394: #endif
  395: {
  396:         Bigint *b1, *p5, *p51;
  397:         int i;
  398:         static int p05[3] = { 5, 25, 125 };
  399: 
  400:         if ( (i = k & 3) !=0) {
  401:                 b = multadd(b, p05[i-1], 0);
  402:                 if (b == NULL)
  403:                         return (NULL);
  404:                 }
  405: 
  406:         if (!(k >>= 2))
  407:                 return b;
  408:         if ((p5 = p5s) == 0) {
  409:                 /* first time */
  410: #ifdef MULTIPLE_THREADS
  411:                 ACQUIRE_DTOA_LOCK(1);
  412:                 if (!(p5 = p5s)) {
  413:                         p5 = p5s = i2b(625);
  414:                         if (p5 == NULL)
  415:                                 return (NULL);
  416:                         p5->next = 0;
  417:                         }
  418:                 FREE_DTOA_LOCK(1);
  419: #else
  420:                 p5 = p5s = i2b(625);
  421:                 if (p5 == NULL)
  422:                         return (NULL);
  423:                 p5->next = 0;
  424: #endif
  425:                 }
  426:         for(;;) {
  427:                 if (k & 1) {
  428:                         b1 = mult(b, p5);
  429:                         if (b1 == NULL)
  430:                                 return (NULL);
  431:                         Bfree(b);
  432:                         b = b1;
  433:                         }
  434:                 if (!(k >>= 1))
  435:                         break;
  436:                 if ((p51 = p5->next) == 0) {
  437: #ifdef MULTIPLE_THREADS
  438:                         ACQUIRE_DTOA_LOCK(1);
  439:                         if (!(p51 = p5->next)) {
  440:                                 p51 = p5->next = mult(p5,p5);
  441:                                 if (p51 == NULL)
  442:                                         return (NULL);
  443:                                 p51->next = 0;
  444:                                 }
  445:                         FREE_DTOA_LOCK(1);
  446: #else
  447:                         p51 = p5->next = mult(p5,p5);
  448:                         if (p51 == NULL)
  449:                                 return (NULL);
  450:                         p51->next = 0;
  451: #endif
  452:                         }
  453:                 p5 = p51;
  454:                 }
  455:         return b;
  456:         }
  457: 
  458:  Bigint *
  459: lshift
  460: #ifdef KR_headers
  461:         (b, k) Bigint *b; int k;
  462: #else
  463:         (Bigint *b, int k)
  464: #endif
  465: {
  466:         int i, k1, n, n1;
  467:         Bigint *b1;
  468:         ULong *x, *x1, *xe, z;
  469: 
  470:         n = k >> kshift;
  471:         k1 = b->k;
  472:         n1 = n + b->wds + 1;
  473:         for(i = b->maxwds; n1 > i; i <<= 1)
  474:                 k1++;
  475:         b1 = Balloc(k1);
  476:         if (b1 == NULL)
  477:                 return (NULL);
  478:         x1 = b1->x;
  479:         for(i = 0; i < n; i++)
  480:                 *x1++ = 0;
  481:         x = b->x;
  482:         xe = x + b->wds;
  483:         if (k &= kmask) {
  484: #ifdef Pack_32
  485:                 k1 = 32 - k;
  486:                 z = 0;
  487:                 do {
  488:                         *x1++ = *x << k | z;
  489:                         z = *x++ >> k1;
  490:                         }
  491:                         while(x < xe);
  492:                 if ((*x1 = z) !=0)
  493:                         ++n1;
  494: #else
  495:                 k1 = 16 - k;
  496:                 z = 0;
  497:                 do {
  498:                         *x1++ = *x << k  & 0xffff | z;
  499:                         z = *x++ >> k1;
  500:                         }
  501:                         while(x < xe);
  502:                 if (*x1 = z)
  503:                         ++n1;
  504: #endif
  505:                 }
  506:         else do
  507:                 *x1++ = *x++;
  508:                 while(x < xe);
  509:         b1->wds = n1 - 1;
  510:         Bfree(b);
  511:         return b1;
  512:         }
  513: 
  514:  int
  515: cmp
  516: #ifdef KR_headers
  517:         (a, b) Bigint *a, *b;
  518: #else
  519:         (Bigint *a, Bigint *b)
  520: #endif
  521: {
  522:         ULong *xa, *xa0, *xb, *xb0;
  523:         int i, j;
  524: 
  525:         i = a->wds;
  526:         j = b->wds;
  527: #ifdef DEBUG
  528:         if (i > 1 && !a->x[i-1])
  529:                 Bug("cmp called with a->x[a->wds-1] == 0");
  530:         if (j > 1 && !b->x[j-1])
  531:                 Bug("cmp called with b->x[b->wds-1] == 0");
  532: #endif
  533:         if (i -= j)
  534:                 return i;
  535:         xa0 = a->x;
  536:         xa = xa0 + j;
  537:         xb0 = b->x;
  538:         xb = xb0 + j;
  539:         for(;;) {
  540:                 if (*--xa != *--xb)
  541:                         return *xa < *xb ? -1 : 1;
  542:                 if (xa <= xa0)
  543:                         break;
  544:                 }
  545:         return 0;
  546:         }
  547: 
  548:  Bigint *
  549: diff
  550: #ifdef KR_headers
  551:         (a, b) Bigint *a, *b;
  552: #else
  553:         (Bigint *a, Bigint *b)
  554: #endif
  555: {
  556:         Bigint *c;
  557:         int i, wa, wb;
  558:         ULong *xa, *xae, *xb, *xbe, *xc;
  559: #ifdef ULLong
  560:         ULLong borrow, y;
  561: #else
  562:         ULong borrow, y;
  563: #ifdef Pack_32
  564:         ULong z;
  565: #endif
  566: #endif
  567: 
  568:         i = cmp(a,b);
  569:         if (!i) {
  570:                 c = Balloc(0);
  571:                 if (c == NULL)
  572:                         return (NULL);
  573:                 c->wds = 1;
  574:                 c->x[0] = 0;
  575:                 return c;
  576:                 }
  577:         if (i < 0) {
  578:                 c = a;
  579:                 a = b;
  580:                 b = c;
  581:                 i = 1;
  582:                 }
  583:         else
  584:                 i = 0;
  585:         c = Balloc(a->k);
  586:         if (c == NULL)
  587:                 return (NULL);
  588:         c->sign = i;
  589:         wa = a->wds;
  590:         xa = a->x;
  591:         xae = xa + wa;
  592:         wb = b->wds;
  593:         xb = b->x;
  594:         xbe = xb + wb;
  595:         xc = c->x;
  596:         borrow = 0;
  597: #ifdef ULLong
  598:         do {
  599:                 y = (ULLong)*xa++ - *xb++ - borrow;
  600:                 borrow = y >> 32 & 1UL;
  601:                 *xc++ = y & 0xffffffffUL;
  602:                 }
  603:                 while(xb < xbe);
  604:         while(xa < xae) {
  605:                 y = *xa++ - borrow;
  606:                 borrow = y >> 32 & 1UL;
  607:                 *xc++ = y & 0xffffffffUL;
  608:                 }
  609: #else
  610: #ifdef Pack_32
  611:         do {
  612:                 y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
  613:                 borrow = (y & 0x10000) >> 16;
  614:                 z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
  615:                 borrow = (z & 0x10000) >> 16;
  616:                 Storeinc(xc, z, y);
  617:                 }
  618:                 while(xb < xbe);
  619:         while(xa < xae) {
  620:                 y = (*xa & 0xffff) - borrow;
  621:                 borrow = (y & 0x10000) >> 16;
  622:                 z = (*xa++ >> 16) - borrow;
  623:                 borrow = (z & 0x10000) >> 16;
  624:                 Storeinc(xc, z, y);
  625:                 }
  626: #else
  627:         do {
  628:                 y = *xa++ - *xb++ - borrow;
  629:                 borrow = (y & 0x10000) >> 16;
  630:                 *xc++ = y & 0xffff;
  631:                 }
  632:                 while(xb < xbe);
  633:         while(xa < xae) {
  634:                 y = *xa++ - borrow;
  635:                 borrow = (y & 0x10000) >> 16;
  636:                 *xc++ = y & 0xffff;
  637:                 }
  638: #endif
  639: #endif
  640:         while(!*--xc)
  641:                 wa--;
  642:         c->wds = wa;
  643:         return c;
  644:         }
  645: 
  646:  double
  647: b2d
  648: #ifdef KR_headers
  649:         (a, e) Bigint *a; int *e;
  650: #else
  651:         (Bigint *a, int *e)
  652: #endif
  653: {
  654:         ULong *xa, *xa0, w, y, z;
  655:         int k;
  656:         U d;
  657: #ifdef VAX
  658:         ULong d0, d1;
  659: #else
  660: #define d0 word0(&d)
  661: #define d1 word1(&d)
  662: #endif
  663: 
  664:         xa0 = a->x;
  665:         xa = xa0 + a->wds;
  666:         y = *--xa;
  667: #ifdef DEBUG
  668:         if (!y) Bug("zero y in b2d");
  669: #endif
  670:         k = hi0bits(y);
  671:         *e = 32 - k;
  672: #ifdef Pack_32
  673:         if (k < Ebits) {
  674:                 d0 = Exp_1 | y >> (Ebits - k);
  675:                 w = xa > xa0 ? *--xa : 0;
  676:                 d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
  677:                 goto ret_d;
  678:                 }
  679:         z = xa > xa0 ? *--xa : 0;
  680:         if (k -= Ebits) {
  681:                 d0 = Exp_1 | y << k | z >> (32 - k);
  682:                 y = xa > xa0 ? *--xa : 0;
  683:                 d1 = z << k | y >> (32 - k);
  684:                 }
  685:         else {
  686:                 d0 = Exp_1 | y;
  687:                 d1 = z;
  688:                 }
  689: #else
  690:         if (k < Ebits + 16) {
  691:                 z = xa > xa0 ? *--xa : 0;
  692:                 d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
  693:                 w = xa > xa0 ? *--xa : 0;
  694:                 y = xa > xa0 ? *--xa : 0;
  695:                 d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
  696:                 goto ret_d;
  697:                 }
  698:         z = xa > xa0 ? *--xa : 0;
  699:         w = xa > xa0 ? *--xa : 0;
  700:         k -= Ebits + 16;
  701:         d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
  702:         y = xa > xa0 ? *--xa : 0;
  703:         d1 = w << k + 16 | y << k;
  704: #endif
  705:  ret_d:
  706: #ifdef VAX
  707:         word0(&d) = d0 >> 16 | d0 << 16;
  708:         word1(&d) = d1 >> 16 | d1 << 16;
  709: #endif
  710:         return dval(&d);
  711:         }
  712: #undef d0
  713: #undef d1
  714: 
  715:  Bigint *
  716: d2b
  717: #ifdef KR_headers
  718:         (dd, e, bits) double dd; int *e, *bits;
  719: #else
  720:         (double dd, int *e, int *bits)
  721: #endif
  722: {
  723:         Bigint *b;
  724:         U d;
  725: #ifndef Sudden_Underflow
  726:         int i;
  727: #endif
  728:         int de, k;
  729:         ULong *x, y, z;
  730: #ifdef VAX
  731:         ULong d0, d1;
  732: #else
  733: #define d0 word0(&d)
  734: #define d1 word1(&d)
  735: #endif
  736:         d.d = dd;
  737: #ifdef VAX
  738:         d0 = word0(&d) >> 16 | word0(&d) << 16;
  739:         d1 = word1(&d) >> 16 | word1(&d) << 16;
  740: #endif
  741: 
  742: #ifdef Pack_32
  743:         b = Balloc(1);
  744: #else
  745:         b = Balloc(2);
  746: #endif
  747:         if (b == NULL)
  748:                 return (NULL);
  749:         x = b->x;
  750: 
  751:         z = d0 & Frac_mask;
  752:         d0 &= 0x7fffffff;      /* clear sign bit, which we ignore */
  753: #ifdef Sudden_Underflow
  754:         de = (int)(d0 >> Exp_shift);
  755: #ifndef IBM
  756:         z |= Exp_msk11;
  757: #endif
  758: #else
  759:         if ( (de = (int)(d0 >> Exp_shift)) !=0)
  760:                 z |= Exp_msk1;
  761: #endif
  762: #ifdef Pack_32
  763:         if ( (y = d1) !=0) {
  764:                 if ( (k = lo0bits(&y)) !=0) {
  765:                         x[0] = y | z << (32 - k);
  766:                         z >>= k;
  767:                         }
  768:                 else
  769:                         x[0] = y;
  770: #ifndef Sudden_Underflow
  771:                 i =
  772: #endif
  773:                      b->wds = (x[1] = z) !=0 ? 2 : 1;
  774:                 }
  775:         else {
  776:                 k = lo0bits(&z);
  777:                 x[0] = z;
  778: #ifndef Sudden_Underflow
  779:                 i =
  780: #endif
  781:                     b->wds = 1;
  782:                 k += 32;
  783:                 }
  784: #else
  785:         if ( (y = d1) !=0) {
  786:                 if ( (k = lo0bits(&y)) !=0)
  787:                         if (k >= 16) {
  788:                                 x[0] = y | z << 32 - k & 0xffff;
  789:                                 x[1] = z >> k - 16 & 0xffff;
  790:                                 x[2] = z >> k;
  791:                                 i = 2;
  792:                                 }
  793:                         else {
  794:                                 x[0] = y & 0xffff;
  795:                                 x[1] = y >> 16 | z << 16 - k & 0xffff;
  796:                                 x[2] = z >> k & 0xffff;
  797:                                 x[3] = z >> k+16;
  798:                                 i = 3;
  799:                                 }
  800:                 else {
  801:                         x[0] = y & 0xffff;
  802:                         x[1] = y >> 16;
  803:                         x[2] = z & 0xffff;
  804:                         x[3] = z >> 16;
  805:                         i = 3;
  806:                         }
  807:                 }
  808:         else {
  809: #ifdef DEBUG
  810:                 if (!z)
  811:                         Bug("Zero passed to d2b");
  812: #endif
  813:                 k = lo0bits(&z);
  814:                 if (k >= 16) {
  815:                         x[0] = z;
  816:                         i = 0;
  817:                         }
  818:                 else {
  819:                         x[0] = z & 0xffff;
  820:                         x[1] = z >> 16;
  821:                         i = 1;
  822:                         }
  823:                 k += 32;
  824:                 }
  825:         while(!x[i])
  826:                 --i;
  827:         b->wds = i + 1;
  828: #endif
  829: #ifndef Sudden_Underflow
  830:         if (de) {
  831: #endif
  832: #ifdef IBM
  833:                 *e = (de - Bias - (P-1) << 2) + k;
  834:                 *bits = 4*P + 8 - k - hi0bits(word0(&d) & Frac_mask);
  835: #else
  836:                 *e = de - Bias - (P-1) + k;
  837:                 *bits = P - k;
  838: #endif
  839: #ifndef Sudden_Underflow
  840:                 }
  841:         else {
  842:                 *e = de - Bias - (P-1) + 1 + k;
  843: #ifdef Pack_32
  844:                 *bits = 32*i - hi0bits(x[i-1]);
  845: #else
  846:                 *bits = (i+2)*16 - hi0bits(x[i]);
  847: #endif
  848:                 }
  849: #endif
  850:         return b;
  851:         }
  852: #undef d0
  853: #undef d1
  854: 
  855:  CONST double
  856: #ifdef IEEE_Arith
  857: bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
  858: CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256
  859:                 };
  860: #else
  861: #ifdef IBM
  862: bigtens[] = { 1e16, 1e32, 1e64 };
  863: CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
  864: #else
  865: bigtens[] = { 1e16, 1e32 };
  866: CONST double tinytens[] = { 1e-16, 1e-32 };
  867: #endif
  868: #endif
  869: 
  870:  CONST double
  871: tens[] = {
  872:                 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
  873:                 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
  874:                 1e20, 1e21, 1e22
  875: #ifdef VAX
  876:                 , 1e23, 1e24
  877: #endif
  878:                 };
  879: 
  880:  char *
  881: #ifdef KR_headers
  882: strcp_D2A(a, b) char *a; char *b;
  883: #else
  884: strcp_D2A(char *a, CONST char *b)
  885: #endif
  886: {
  887:         while((*a = *b++))
  888:                 a++;
  889:         return a;
  890:         }
  891: 
  892: #ifdef NO_STRING_H
  893: 
  894:  Char *
  895: #ifdef KR_headers
  896: memcpy_D2A(a, b, len) Char *a; Char *b; size_t len;
  897: #else
  898: memcpy_D2A(void *a1, void *b1, size_t len)
  899: #endif
  900: {
  901:         char *a = (char*)a1, *ae = a + len;
  902:         char *b = (char*)b1, *a0 = a;
  903:         while(a < ae)
  904:                 *a++ = *b++;
  905:         return a0;
  906:         }
  907: 
  908: #endif /* NO_STRING_H */