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