lib/libc/gdtoa/dmisc.c

   1 /* $NetBSD: dmisc.c,v 1.3 2007/02/03 16:44:02 christos Exp $ */
   2
   3 /****************************************************************
   4
   5 The author of this software is David M. Gay.
   6
   7 Copyright (C) 1998 by Lucent Technologies
   8 All Rights Reserved
   9
  10 Permission to use, copy, modify, and distribute this software and
  11 its documentation for any purpose and without fee is hereby
  12 granted, provided that the above copyright notice appear in all
  13 copies and that both that the copyright notice and this
  14 permission notice and warranty disclaimer appear in supporting
  15 documentation, and that the name of Lucent or any of its entities
  16 not be used in advertising or publicity pertaining to
  17 distribution of the software without specific, written prior
  18 permission.
  19
  20 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
  21 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
  22 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
  23 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  24 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
  25 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
  26 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
  27 THIS SOFTWARE.
  28
  29 ****************************************************************/
  30
  31 /* Please send bug reports to David M. Gay (dmg at acm dot org,
  32  * with " at " changed at "@" and " dot " changed to ".").      */
  33
  34 #include "gdtoaimp.h"
  35
  36 #ifndef MULTIPLE_THREADS
  37  char *dtoa_result;
  38 #endif
  39
  40  char *
  41 #ifdef KR_headers
  42 rv_alloc(i) size_t i;
  43 #else
  44 rv_alloc(size_t i)
  45 #endif
  46 {
  47         size_t j;
  48         int k, *r;
  49
  50         j = sizeof(ULong);
  51         for(k = 0;
  52                 sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i;
  53                 j <<= 1)
  54                         k++;
  55         r = (int*)(void*)Balloc(k);
  56         if (r == NULL)
  57                 return NULL;
  58         *r = k;
  59         return
  60 #ifndef MULTIPLE_THREADS
  61         dtoa_result =
  62 #endif
  63                 (char *)(void *)(r+1);
  64         }
  65
  66  char *
  67 #ifdef KR_headers
  68 nrv_alloc(s, rve, n) CONST char *s; char **rve; size_t n;
  69 #else
  70 nrv_alloc(CONST char *s, char **rve, size_t n)
  71 #endif
  72 {
  73         char *rv, *t;
  74
  75         t = rv = rv_alloc(n);
  76         if (t == NULL)
  77                 return NULL;
  78         while((*t = *s++) !=0)
  79                 t++;
  80         if (rve)
  81                 *rve = t;
  82         return rv;
  83         }
  84
  85 /* freedtoa(s) must be used to free values s returned by dtoa
  86  * when MULTIPLE_THREADS is #defined.  It should be used in all cases,
  87  * but for consistency with earlier versions of dtoa, it is optional
  88  * when MULTIPLE_THREADS is not defined.
  89  */
  90
  91  void
  92 #ifdef KR_headers
  93 freedtoa(s) char *s;
  94 #else
  95 freedtoa(char *s)
  96 #endif
  97 {
  98         Bigint *b = (Bigint *)(void *)((int *)(void *)s - 1);
  99         b->maxwds = 1 << (b->k = *(int*)(void*)b);
 100         Bfree(b);
 101 #ifndef MULTIPLE_THREADS
 102         if (s == dtoa_result)
 103                 dtoa_result = 0;
 104 #endif
 105         }
 106
 107  int
 108 quorem
 109 #ifdef KR_headers
 110         (b, S) Bigint *b, *S;
 111 #else
 112         (Bigint *b, Bigint *S)
 113 #endif
 114 {
 115         int n;
 116         ULong *bx, *bxe, q, *sx, *sxe;
 117 #ifdef ULLong
 118         ULLong borrow, carry, y, ys;
 119 #else
 120         ULong borrow, carry, y, ys;
 121 #ifdef Pack_32
 122         ULong si, z, zs;
 123 #endif
 124 #endif
 125
 126         n = S->wds;
 127 #ifdef DEBUG
 128         /*debug*/ if (b->wds > n)
 129         /*debug*/       Bug("oversize b in quorem");
 130 #endif
 131         if (b->wds < n)
 132                 return 0;
 133         sx = S->x;
 134         sxe = sx + --n;
 135         bx = b->x;
 136         bxe = bx + n;
 137         q = *bxe / (*sxe + 1);  /* ensure q <= true quotient */
 138 #ifdef DEBUG
 139         /*debug*/ if (q > 9)
 140         /*debug*/       Bug("oversized quotient in quorem");
 141 #endif
 142         if (q) {
 143                 borrow = 0;
 144                 carry = 0;
 145                 do {
 146 #ifdef ULLong
 147                         ys = *sx++ * (ULLong)q + carry;
 148                         carry = ys >> 32;
 149                         /* LINTED conversion */
 150                         y = *bx - (ys & 0xffffffffUL) - borrow;
 151                         borrow = y >> 32 & 1UL;
 152                         /* LINTED conversion */
 153                         *bx++ = y & 0xffffffffUL;
 154 #else
 155 #ifdef Pack_32
 156                         si = *sx++;
 157                         ys = (si & 0xffff) * q + carry;
 158                         zs = (si >> 16) * q + (ys >> 16);
 159                         carry = zs >> 16;
 160                         y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
 161                         borrow = (y & 0x10000) >> 16;
 162                         z = (*bx >> 16) - (zs & 0xffff) - borrow;
 163                         borrow = (z & 0x10000) >> 16;
 164                         Storeinc(bx, z, y);
 165 #else
 166                         ys = *sx++ * q + carry;
 167                         carry = ys >> 16;
 168                         y = *bx - (ys & 0xffff) - borrow;
 169                         borrow = (y & 0x10000) >> 16;
 170                         *bx++ = y & 0xffff;
 171 #endif
 172 #endif
 173                         }
 174                         while(sx <= sxe);
 175                 if (!*bxe) {
 176                         bx = b->x;
 177                         while(--bxe > bx && !*bxe)
 178                                 --n;
 179                         b->wds = n;
 180                         }
 181                 }
 182         if (cmp(b, S) >= 0) {
 183                 q++;
 184                 borrow = 0;
 185                 carry = 0;
 186                 bx = b->x;
 187                 sx = S->x;
 188                 do {
 189 #ifdef ULLong
 190                         ys = *sx++ + carry;
 191                         carry = ys >> 32;
 192                         /* LINTED conversion */
 193                         y = *bx - (ys & 0xffffffffUL) - borrow;
 194                         borrow = y >> 32 & 1UL;
 195                         /* LINTED conversion */
 196                         *bx++ = y & 0xffffffffUL;
 197 #else
 198 #ifdef Pack_32
 199                         si = *sx++;
 200                         ys = (si & 0xffff) + carry;
 201                         zs = (si >> 16) + (ys >> 16);
 202                         carry = zs >> 16;
 203                         y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
 204                         borrow = (y & 0x10000) >> 16;
 205                         z = (*bx >> 16) - (zs & 0xffff) - borrow;
 206                         borrow = (z & 0x10000) >> 16;
 207                         Storeinc(bx, z, y);
 208 #else
 209                         ys = *sx++ + carry;
 210                         carry = ys >> 16;
 211                         y = *bx - (ys & 0xffff) - borrow;
 212                         borrow = (y & 0x10000) >> 16;
 213                         *bx++ = y & 0xffff;
 214 #endif
 215 #endif
 216                         }
 217                         while(sx <= sxe);
 218                 bx = b->x;
 219                 bxe = bx + n;
 220                 if (!*bxe) {
 221                         while(--bxe > bx && !*bxe)
 222                                 --n;
 223                         b->wds = n;
 224                         }
 225                 }
 226         return q;
 227         }