contrib/gdtoa/dmisc.c

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