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