gonzui

    1: /****************************************************************
    2: 
    3: The author of this software is David M. Gay.
    4: 
    5: Copyright (C) 1998 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: #ifndef MULTIPLE_THREADS
   35:  char *dtoa_result;
   36: #endif
   37: 
   38:  char *
   39: #ifdef KR_headers
   40: rv_alloc(i) int i;
   41: #else
   42: rv_alloc(int i)
   43: #endif
   44: {
   45:         int j, k, *r;
   46: 
   47:         j = sizeof(ULong);
   48:         for(k = 0;
   49:                 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i;
   50:                 j <<= 1)
   51:                         k++;
   52:         r = (int*)Balloc(k);
   53:         if (r == NULL)
   54:                 return (
   55: #ifndef MULTIPLE_THREADS
   56:                 dtoa_result =
   57: #endif
   58:                         NULL);
   59:         *r = k;
   60:         return
   61: #ifndef MULTIPLE_THREADS
   62:         dtoa_result =
   63: #endif
   64:                 (char *)(r+1);
   65:         }
   66: 
   67:  char *
   68: #ifdef KR_headers
   69: nrv_alloc(s, rve, n) char *s, **rve; int n;
   70: #else
   71: nrv_alloc(char *s, char **rve, int n)
   72: #endif
   73: {
   74:         char *rv, *t;
   75: 
   76:         t = rv = rv_alloc(n);
   77:         if (t == NULL)
   78:                 return (NULL);
   79:         while((*t = *s++) !=0)
   80:                 t++;
   81:         if (rve)
   82:                 *rve = t;
   83:         return rv;
   84:         }
   85: 
   86: /* freedtoa(s) must be used to free values s returned by dtoa
   87:  * when MULTIPLE_THREADS is #defined.  It should be used in all cases,
   88:  * but for consistency with earlier versions of dtoa, it is optional
   89:  * when MULTIPLE_THREADS is not defined.
   90:  */
   91: 
   92:  void
   93: #ifdef KR_headers
   94: freedtoa(s) char *s;
   95: #else
   96: freedtoa(char *s)
   97: #endif
   98: {
   99:         Bigint *b = (Bigint *)((int *)s - 1);
  100:         b->maxwds = 1 << (b->k = *(int*)b);
  101:         Bfree(b);
  102: #ifndef MULTIPLE_THREADS
  103:         if (s == dtoa_result)
  104:                 dtoa_result = 0;
  105: #endif
  106:         }
  107: 
  108:  int
  109: quorem
  110: #ifdef KR_headers
  111:         (b, S) Bigint *b, *S;
  112: #else
  113:         (Bigint *b, Bigint *S)
  114: #endif
  115: {
  116:         int n;
  117:         ULong *bx, *bxe, q, *sx, *sxe;
  118: #ifdef ULLong
  119:         ULLong borrow, carry, y, ys;
  120: #else
  121:         ULong borrow, carry, y, ys;
  122: #ifdef Pack_32
  123:         ULong si, z, zs;
  124: #endif
  125: #endif
  126: 
  127:         n = S->wds;
  128: #ifdef DEBUG
  129:         /*debug*/ if (b->wds > n)
  130:         /*debug*/      Bug("oversize b in quorem");
  131: #endif
  132:         if (b->wds < n)
  133:                 return 0;
  134:         sx = S->x;
  135:         sxe = sx + --n;
  136:         bx = b->x;
  137:         bxe = bx + n;
  138:         q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
  139: #ifdef DEBUG
  140:         /*debug*/ if (q > 9)
  141:         /*debug*/      Bug("oversized quotient in quorem");
  142: #endif
  143:         if (q) {
  144:                 borrow = 0;
  145:                 carry = 0;
  146:                 do {
  147: #ifdef ULLong
  148:                         ys = *sx++ * (ULLong)q + carry;
  149:                         carry = ys >> 32;
  150:                         y = *bx - (ys & 0xffffffffUL) - borrow;
  151:                         borrow = y >> 32 & 1UL;
  152:                         *bx++ = y & 0xffffffffUL;
  153: #else
  154: #ifdef Pack_32
  155:                         si = *sx++;
  156:                         ys = (si & 0xffff) * q + carry;
  157:                         zs = (si >> 16) * q + (ys >> 16);
  158:                         carry = zs >> 16;
  159:                         y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
  160:                         borrow = (y & 0x10000) >> 16;
  161:                         z = (*bx >> 16) - (zs & 0xffff) - borrow;
  162:                         borrow = (z & 0x10000) >> 16;
  163:                         Storeinc(bx, z, y);
  164: #else
  165:                         ys = *sx++ * q + carry;
  166:                         carry = ys >> 16;
  167:                         y = *bx - (ys & 0xffff) - borrow;
  168:                         borrow = (y & 0x10000) >> 16;
  169:                         *bx++ = y & 0xffff;
  170: #endif
  171: #endif
  172:                         }
  173:                         while(sx <= sxe);
  174:                 if (!*bxe) {
  175:                         bx = b->x;
  176:                         while(--bxe > bx && !*bxe)
  177:                                 --n;
  178:                         b->wds = n;
  179:                         }
  180:                 }
  181:         if (cmp(b, S) >= 0) {
  182:                 q++;
  183:                 borrow = 0;
  184:                 carry = 0;
  185:                 bx = b->x;
  186:                 sx = S->x;
  187:                 do {
  188: #ifdef ULLong
  189:                         ys = *sx++ + carry;
  190:                         carry = ys >> 32;
  191:                         y = *bx - (ys & 0xffffffffUL) - borrow;
  192:                         borrow = y >> 32 & 1UL;
  193:                         *bx++ = y & 0xffffffffUL;
  194: #else
  195: #ifdef Pack_32
  196:                         si = *sx++;
  197:                         ys = (si & 0xffff) + carry;
  198:                         zs = (si >> 16) + (ys >> 16);
  199:                         carry = zs >> 16;
  200:                         y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
  201:                         borrow = (y & 0x10000) >> 16;
  202:                         z = (*bx >> 16) - (zs & 0xffff) - borrow;
  203:                         borrow = (z & 0x10000) >> 16;
  204:                         Storeinc(bx, z, y);
  205: #else
  206:                         ys = *sx++ + carry;
  207:                         carry = ys >> 16;
  208:                         y = *bx - (ys & 0xffff) - borrow;
  209:                         borrow = (y & 0x10000) >> 16;
  210:                         *bx++ = y & 0xffff;
  211: #endif
  212: #endif
  213:                         }
  214:                         while(sx <= sxe);
  215:                 bx = b->x;
  216:                 bxe = bx + n;
  217:                 if (!*bxe) {
  218:                         while(--bxe > bx && !*bxe)
  219:                                 --n;
  220:                         b->wds = n;
  221:                         }
  222:                 }
  223:         return q;
  224:         }