| 1 | /**************************************************************** |
| 2 | |
| 3 | The author of this software is David M. Gay. |
| 4 | |
| 5 | Copyright (C) 1998, 1999 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 | |
| 33 | #if defined(__MINGW32__) || defined(__MINGW64__) |
| 34 | /* we have to include windows.h before gdtoa |
| 35 | headers, otherwise defines cause conflicts. */ |
| 36 | #ifndef WIN32_LEAN_AND_MEAN |
| 37 | #define WIN32_LEAN_AND_MEAN |
| 38 | #endif |
| 39 | #include <windows.h> |
| 40 | |
| 41 | #define NLOCKS 2 |
| 42 | |
| 43 | #ifdef USE_WIN32_SL |
| 44 | /* Use spin locks. */ |
| 45 | static long dtoa_sl[NLOCKS]; |
| 46 | |
| 47 | #define ACQUIRE_DTOA_LOCK(n) \ |
| 48 | while (InterlockedCompareExchange (&dtoa_sl[n], 1, 0) != 0) \ |
| 49 | Sleep (0); |
| 50 | #define FREE_DTOA_LOCK(n) InterlockedExchange (&dtoa_sl[n], 0); |
| 51 | |
| 52 | #else	/* USE_WIN32_SL */ |
| 53 | |
| 54 | #include <stdlib.h> |
| 55 | static CRITICAL_SECTION dtoa_CritSec[NLOCKS]; |
| 56 | static long dtoa_CS_init = 0; |
| 57 | /* |
| 58 | 1 = initializing |
| 59 | 2 = initialized |
| 60 | 3 = deleted |
| 61 | */ |
| 62 | static void dtoa_lock_cleanup (void) |
| 63 | { |
| 64 | 	long last_CS_init = InterlockedExchange (&dtoa_CS_init,3); |
| 65 | 	if (2 == last_CS_init) { |
| 66 | 		int i; |
| 67 | 		for (i = 0; i < NLOCKS; i++) |
| 68 | 			DeleteCriticalSection (&dtoa_CritSec[i]); |
| 69 | 	} |
| 70 | } |
| 71 | |
| 72 | static void dtoa_lock (unsigned int n) |
| 73 | { |
| 74 | 	if (2 == dtoa_CS_init) { |
| 75 | 		EnterCriticalSection (&dtoa_CritSec[n]); |
| 76 | 		return; |
| 77 | 	} |
| 78 | 	else if (0 == dtoa_CS_init) { |
| 79 | 		long last_CS_init = InterlockedExchange (&dtoa_CS_init, 1); |
| 80 | 		if (0 == last_CS_init) { |
| 81 | 			int i; |
| 82 | 			for (i = 0; i < NLOCKS; i++) |
| 83 | 				InitializeCriticalSection (&dtoa_CritSec[i]); |
| 84 | 			atexit (dtoa_lock_cleanup); |
| 85 | 			dtoa_CS_init = 2; |
| 86 | 		} |
| 87 | 		else if (2 == last_CS_init) |
| 88 | 			dtoa_CS_init = 2; |
| 89 | 	} |
| 90 | 	/* Another thread is initializing. Wait. */ |
| 91 | 	while (1 == dtoa_CS_init) |
| 92 | 		Sleep (1); |
| 93 | |
| 94 | 	/* It had better be initialized now. */ |
| 95 | 	if (2 == dtoa_CS_init) |
| 96 | 		EnterCriticalSection(&dtoa_CritSec[n]); |
| 97 | } |
| 98 | |
| 99 | static void dtoa_unlock (unsigned int n) |
| 100 | { |
| 101 | 	if (2 == dtoa_CS_init) |
| 102 | 		LeaveCriticalSection (&dtoa_CritSec[n]); |
| 103 | } |
| 104 | |
| 105 | #define ACQUIRE_DTOA_LOCK(n) dtoa_lock(n) |
| 106 | #define FREE_DTOA_LOCK(n) dtoa_unlock(n) |
| 107 | #endif	/* USE_WIN32_SL */ |
| 108 | |
| 109 | #endif	/* __MINGW32__ / __MINGW64__ */ |
| 110 | |
| 111 | #include "gdtoaimp.h" |
| 112 | |
| 113 | static Bigint *freelist[Kmax+1]; |
| 114 | #ifndef Omit_Private_Memory |
| 115 | #ifndef PRIVATE_MEM |
| 116 | #define PRIVATE_MEM 2304 |
| 117 | #endif |
| 118 | #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double)) |
| 119 | static double private_mem[PRIVATE_mem], *pmem_next = private_mem; |
| 120 | #endif |
| 121 | |
| 122 | Bigint *Balloc (int k) |
| 123 | { |
| 124 | 	int x; |
| 125 | 	Bigint *rv; |
| 126 | #ifndef Omit_Private_Memory |
| 127 | 	unsigned int len; |
| 128 | #endif |
| 129 | |
| 130 | 	ACQUIRE_DTOA_LOCK(0); |
| 131 | 	/* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */ |
| 132 | 	/* but this case seems very unlikely. */ |
| 133 | 	if (k <= Kmax && (rv = freelist[k]) !=0) { |
| 134 | 		freelist[k] = rv->next; |
| 135 | 	} |
| 136 | 	else { |
| 137 | 		x = 1 << k; |
| 138 | #ifdef Omit_Private_Memory |
| 139 | 		rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong)); |
| 140 | 		if (rv == NULL) |
| 141 | 			return NULL; |
| 142 | #else |
| 143 | 		len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1) |
| 144 | 			/sizeof(double); |
| 145 | 		if (k <= Kmax |
| 146 | 		 && (size_t) (pmem_next - private_mem + len) <= PRIVATE_mem) |
| 147 | 		{ |
| 148 | 			rv = (Bigint*)pmem_next; |
| 149 | 			pmem_next += len; |
| 150 | 		} |
| 151 | 		else |
| 152 | 		{ |
| 153 | 			rv = (Bigint*)MALLOC(len*sizeof(double)); |
| 154 | 			if (rv == NULL) |
| 155 | 				return NULL; |
| 156 | 		} |
| 157 | #endif |
| 158 | 		rv->k = k; |
| 159 | 		rv->maxwds = x; |
| 160 | 	} |
| 161 | 	FREE_DTOA_LOCK(0); |
| 162 | 	rv->sign = rv->wds = 0; |
| 163 | 	return rv; |
| 164 | } |
| 165 | |
| 166 | void Bfree (Bigint *v) |
| 167 | { |
| 168 | 	if (v) { |
| 169 | 		if (v->k > Kmax) |
| 170 | 			free((void*)v); |
| 171 | 		else { |
| 172 | 			ACQUIRE_DTOA_LOCK(0); |
| 173 | 			v->next = freelist[v->k]; |
| 174 | 			freelist[v->k] = v; |
| 175 | 			FREE_DTOA_LOCK(0); |
| 176 | 		} |
| 177 | 	} |
| 178 | } |
| 179 | |
| 180 | /* lo0bits(): Shift y so lowest bit is 1 and return the |
| 181 | *		 number of bits y was shifted. |
| 182 | * With GCC, we use an inline wrapper for __builtin_clz() |
| 183 | */ |
| 184 | #ifndef __GNUC__ |
| 185 | int lo0bits (ULong *y) |
| 186 | { |
| 187 | 	int k; |
| 188 | 	ULong x = *y; |
| 189 | |
| 190 | 	if (x & 7) { |
| 191 | 		if (x & 1) |
| 192 | 			return 0; |
| 193 | 		if (x & 2) { |
| 194 | 			*y = x >> 1; |
| 195 | 			return 1; |
| 196 | 		} |
| 197 | 		*y = x >> 2; |
| 198 | 		return 2; |
| 199 | 	} |
| 200 | 	k = 0; |
| 201 | 	if (!(x & 0xffff)) { |
| 202 | 		k = 16; |
| 203 | 		x >>= 16; |
| 204 | 	} |
| 205 | 	if (!(x & 0xff)) { |
| 206 | 		k += 8; |
| 207 | 		x >>= 8; |
| 208 | 	} |
| 209 | 	if (!(x & 0xf)) { |
| 210 | 		k += 4; |
| 211 | 		x >>= 4; |
| 212 | 	} |
| 213 | 	if (!(x & 0x3)) { |
| 214 | 		k += 2; |
| 215 | 		x >>= 2; |
| 216 | 	} |
| 217 | 	if (!(x & 1)) { |
| 218 | 		k++; |
| 219 | 		x >>= 1; |
| 220 | 		if (!x) |
| 221 | 			return 32; |
| 222 | 	} |
| 223 | 	*y = x; |
| 224 | 	return k; |
| 225 | } |
| 226 | #endif	/* __GNUC__ */ |
| 227 | |
| 228 | Bigint *multadd (Bigint *b, int m, int a)	/* multiply by m and add a */ |
| 229 | { |
| 230 | 	int i, wds; |
| 231 | #ifdef ULLong |
| 232 | 	ULong *x; |
| 233 | 	ULLong carry, y; |
| 234 | #else |
| 235 | 	ULong carry, *x, y; |
| 236 | #ifdef Pack_32 |
| 237 | 	ULong xi, z; |
| 238 | #endif |
| 239 | #endif |
| 240 | 	Bigint *b1; |
| 241 | |
| 242 | 	wds = b->wds; |
| 243 | 	x = b->x; |
| 244 | 	i = 0; |
| 245 | 	carry = a; |
| 246 | 	do { |
| 247 | #ifdef ULLong |
| 248 | 		y = *x * (ULLong)m + carry; |
| 249 | 		carry = y >> 32; |
| 250 | 		*x++ = y & 0xffffffffUL; |
| 251 | #else |
| 252 | #ifdef Pack_32 |
| 253 | 		xi = *x; |
| 254 | 		y = (xi & 0xffff) * m + carry; |
| 255 | 		z = (xi >> 16) * m + (y >> 16); |
| 256 | 		carry = z >> 16; |
| 257 | 		*x++ = (z << 16) + (y & 0xffff); |
| 258 | #else |
| 259 | 		y = *x * m + carry; |
| 260 | 		carry = y >> 16; |
| 261 | 		*x++ = y & 0xffff; |
| 262 | #endif |
| 263 | #endif |
| 264 | 	} while(++i < wds); |
| 265 | 	if (carry) { |
| 266 | 		if (wds >= b->maxwds) { |
| 267 | 			b1 = Balloc(b->k+1); |
| 268 | 			if (b1 == NULL) |
| 269 | 				return NULL; |
| 270 | 			Bcopy(b1, b); |
| 271 | 			Bfree(b); |
| 272 | 			b = b1; |
| 273 | 		} |
| 274 | 		b->x[wds++] = carry; |
| 275 | 		b->wds = wds; |
| 276 | 	} |
| 277 | 	return b; |
| 278 | } |
| 279 | |
| 280 | /* hi0bits(); |
| 281 | * With GCC, we use an inline wrapper for __builtin_clz() |
| 282 | */ |
| 283 | #ifndef __GNUC__ |
| 284 | int hi0bits_D2A (ULong x) |
| 285 | { |
| 286 | 	int k = 0; |
| 287 | |
| 288 | 	if (!(x & 0xffff0000)) { |
| 289 | 		k = 16; |
| 290 | 		x <<= 16; |
| 291 | 	} |
| 292 | 	if (!(x & 0xff000000)) { |
| 293 | 		k += 8; |
| 294 | 		x <<= 8; |
| 295 | 	} |
| 296 | 	if (!(x & 0xf0000000)) { |
| 297 | 		k += 4; |
| 298 | 		x <<= 4; |
| 299 | 	} |
| 300 | 	if (!(x & 0xc0000000)) { |
| 301 | 		k += 2; |
| 302 | 		x <<= 2; |
| 303 | 	} |
| 304 | 	if (!(x & 0x80000000)) { |
| 305 | 		k++; |
| 306 | 		if (!(x & 0x40000000)) |
| 307 | 			return 32; |
| 308 | 	} |
| 309 | 	return k; |
| 310 | } |
| 311 | #endif	/* __GNUC__ */ |
| 312 | |
| 313 | Bigint *i2b (int i) |
| 314 | { |
| 315 | 	Bigint *b; |
| 316 | |
| 317 | 	b = Balloc(1); |
| 318 | 	if (b == NULL) |
| 319 | 		return NULL; |
| 320 | 	b->x[0] = i; |
| 321 | 	b->wds = 1; |
| 322 | 	return b; |
| 323 | } |
| 324 | |
| 325 | Bigint *mult (Bigint *a, Bigint *b) |
| 326 | { |
| 327 | 	Bigint *c; |
| 328 | 	int k, wa, wb, wc; |
| 329 | 	ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0; |
| 330 | 	ULong y; |
| 331 | #ifdef ULLong |
| 332 | 	ULLong carry, z; |
| 333 | #else |
| 334 | 	ULong carry, z; |
| 335 | #ifdef Pack_32 |
| 336 | 	ULong z2; |
| 337 | #endif |
| 338 | #endif |
| 339 | |
| 340 | 	if (a->wds < b->wds) { |
| 341 | 		c = a; |
| 342 | 		a = b; |
| 343 | 		b = c; |
| 344 | 	} |
| 345 | 	k = a->k; |
| 346 | 	wa = a->wds; |
| 347 | 	wb = b->wds; |
| 348 | 	wc = wa + wb; |
| 349 | 	if (wc > a->maxwds) |
| 350 | 		k++; |
| 351 | 	c = Balloc(k); |
| 352 | 	if (c == NULL) |
| 353 | 		return NULL; |
| 354 | 	for(x = c->x, xa = x + wc; x < xa; x++) |
| 355 | 		*x = 0; |
| 356 | 	xa = a->x; |
| 357 | 	xae = xa + wa; |
| 358 | 	xb = b->x; |
| 359 | 	xbe = xb + wb; |
| 360 | 	xc0 = c->x; |
| 361 | #ifdef ULLong |
| 362 | 	for(; xb < xbe; xc0++) { |
| 363 | 		if ( (y = *xb++) !=0) { |
| 364 | 			x = xa; |
| 365 | 			xc = xc0; |
| 366 | 			carry = 0; |
| 367 | 			do { |
| 368 | 				z = *x++ * (ULLong)y + *xc + carry; |
| 369 | 				carry = z >> 32; |
| 370 | 				*xc++ = z & 0xffffffffUL; |
| 371 | 			} while(x < xae); |
| 372 | 			*xc = carry; |
| 373 | 		} |
| 374 | 	} |
| 375 | #else |
| 376 | #ifdef Pack_32 |
| 377 | 	for(; xb < xbe; xb++, xc0++) { |
| 378 | 		if ( (y = *xb & 0xffff) !=0) { |
| 379 | 			x = xa; |
| 380 | 			xc = xc0; |
| 381 | 			carry = 0; |
| 382 | 			do { |
| 383 | 				z = (*x & 0xffff) * y + (*xc & 0xffff) + carry; |
| 384 | 				carry = z >> 16; |
| 385 | 				z2 = (*x++ >> 16) * y + (*xc >> 16) + carry; |
| 386 | 				carry = z2 >> 16; |
| 387 | 				Storeinc(xc, z2, z); |
| 388 | 			} while(x < xae); |
| 389 | 			*xc = carry; |
| 390 | 		} |
| 391 | 		if ( (y = *xb >> 16) !=0) { |
| 392 | 			x = xa; |
| 393 | 			xc = xc0; |
| 394 | 			carry = 0; |
| 395 | 			z2 = *xc; |
| 396 | 			do { |
| 397 | 				z = (*x & 0xffff) * y + (*xc >> 16) + carry; |
| 398 | 				carry = z >> 16; |
| 399 | 				Storeinc(xc, z, z2); |
| 400 | 				z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry; |
| 401 | 				carry = z2 >> 16; |
| 402 | 			} while(x < xae); |
| 403 | 			*xc = z2; |
| 404 | 		} |
| 405 | 	} |
| 406 | #else |
| 407 | 	for(; xb < xbe; xc0++) { |
| 408 | 		if ( (y = *xb++) !=0) { |
| 409 | 			x = xa; |
| 410 | 			xc = xc0; |
| 411 | 			carry = 0; |
| 412 | 			do { |
| 413 | 				z = *x++ * y + *xc + carry; |
| 414 | 				carry = z >> 16; |
| 415 | 				*xc++ = z & 0xffff; |
| 416 | 			} while(x < xae); |
| 417 | 			*xc = carry; |
| 418 | 		} |
| 419 | 	} |
| 420 | #endif |
| 421 | #endif |
| 422 | 	for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ; |
| 423 | 	c->wds = wc; |
| 424 | 	return c; |
| 425 | } |
| 426 | |
| 427 | static Bigint *p5s; |
| 428 | |
| 429 | Bigint *pow5mult (Bigint *b, int k) |
| 430 | { |
| 431 | 	Bigint *b1, *p5, *p51; |
| 432 | 	int i; |
| 433 | 	static int p05[3] = { 5, 25, 125 }; |
| 434 | |
| 435 | 	if ( (i = k & 3) !=0) |
| 436 | 	{ |
| 437 | 		b = multadd(b, p05[i-1], 0); |
| 438 | 		if (b == NULL) |
| 439 | 		return NULL; |
| 440 | 	} |
| 441 | 	if (!(k >>= 2)) |
| 442 | 		return b; |
| 443 | 	if ((p5 = p5s) == 0) { |
| 444 | 		/* first time */ |
| 445 | #ifdef MULTIPLE_THREADS |
| 446 | 		ACQUIRE_DTOA_LOCK(1); |
| 447 | 		if (!(p5 = p5s)) { |
| 448 | 			p5 = p5s = i2b(625); |
| 449 | 			if (p5 == NULL) |
| 450 | 				return NULL; |
| 451 | 			p5->next = 0; |
| 452 | 		} |
| 453 | 		FREE_DTOA_LOCK(1); |
| 454 | #else |
| 455 | 		p5 = p5s = i2b(625); |
| 456 | 		if (p5 == NULL) |
| 457 | 			return NULL; |
| 458 | 		p5->next = 0; |
| 459 | #endif |
| 460 | 	} |
| 461 | 	for(;;) { |
| 462 | 		if (k & 1) { |
| 463 | 			b1 = mult(b, p5); |
| 464 | 			if (b1 == NULL) |
| 465 | 				return NULL; |
| 466 | 			Bfree(b); |
| 467 | 			b = b1; |
| 468 | 		} |
| 469 | 		if (!(k >>= 1)) |
| 470 | 			break; |
| 471 | 		if ((p51 = p5->next) == 0) { |
| 472 | #ifdef MULTIPLE_THREADS |
| 473 | 			ACQUIRE_DTOA_LOCK(1); |
| 474 | 			if (!(p51 = p5->next)) { |
| 475 | 				p51 = p5->next = mult(p5,p5); |
| 476 | 				if (p51 == NULL) |
| 477 | 					return NULL; |
| 478 | 				p51->next = 0; |
| 479 | 			} |
| 480 | 			FREE_DTOA_LOCK(1); |
| 481 | #else |
| 482 | 			p51 = p5->next = mult(p5,p5); |
| 483 | 			if (p51 == NULL) |
| 484 | 				return NULL; |
| 485 | 			p51->next = 0; |
| 486 | #endif |
| 487 | 		} |
| 488 | 		p5 = p51; |
| 489 | 	} |
| 490 | 	return b; |
| 491 | } |
| 492 | |
| 493 | Bigint *lshift (Bigint *b, int k) |
| 494 | { |
| 495 | 	int i, k1, n, n1; |
| 496 | 	Bigint *b1; |
| 497 | 	ULong *x, *x1, *xe, z; |
| 498 | |
| 499 | 	n = k >> kshift; |
| 500 | 	k1 = b->k; |
| 501 | 	n1 = n + b->wds + 1; |
| 502 | 	for(i = b->maxwds; n1 > i; i <<= 1) |
| 503 | 		k1++; |
| 504 | 	b1 = Balloc(k1); |
| 505 | 	if (b1 == NULL) |
| 506 | 		return NULL; |
| 507 | 	x1 = b1->x; |
| 508 | 	for(i = 0; i < n; i++) |
| 509 | 		*x1++ = 0; |
| 510 | 	x = b->x; |
| 511 | 	xe = x + b->wds; |
| 512 | 	if (k &= kmask) { |
| 513 | #ifdef Pack_32 |
| 514 | 		k1 = 32 - k; |
| 515 | 		z = 0; |
| 516 | 		do { |
| 517 | 			*x1++ = *x << k | z; |
| 518 | 			z = *x++ >> k1; |
| 519 | 		} while(x < xe); |
| 520 | 		if ((*x1 = z) !=0) |
| 521 | 			++n1; |
| 522 | #else |
| 523 | 		k1 = 16 - k; |
| 524 | 		z = 0; |
| 525 | 		do { |
| 526 | 			*x1++ = *x << k & 0xffff | z; |
| 527 | 			z = *x++ >> k1; |
| 528 | 		} while(x < xe); |
| 529 | 		if (*x1 = z) |
| 530 | 			++n1; |
| 531 | #endif |
| 532 | 	} |
| 533 | 	else do |
| 534 | 		*x1++ = *x++; |
| 535 | 		while(x < xe); |
| 536 | 	b1->wds = n1 - 1; |
| 537 | 	Bfree(b); |
| 538 | 	return b1; |
| 539 | } |
| 540 | |
| 541 | int cmp (Bigint *a, Bigint *b) |
| 542 | { |
| 543 | 	ULong *xa, *xa0, *xb, *xb0; |
| 544 | 	int i, j; |
| 545 | |
| 546 | 	i = a->wds; |
| 547 | 	j = b->wds; |
| 548 | #ifdef DEBUG |
| 549 | 	if (i > 1 && !a->x[i-1]) |
| 550 | 		Bug("cmp called with a->x[a->wds-1] == 0"); |
| 551 | 	if (j > 1 && !b->x[j-1]) |
| 552 | 		Bug("cmp called with b->x[b->wds-1] == 0"); |
| 553 | #endif |
| 554 | 	if (i -= j) |
| 555 | 		return i; |
| 556 | 	xa0 = a->x; |
| 557 | 	xa = xa0 + j; |
| 558 | 	xb0 = b->x; |
| 559 | 	xb = xb0 + j; |
| 560 | 	for(;;) { |
| 561 | 		if (*--xa != *--xb) |
| 562 | 			return *xa < *xb ? -1 : 1; |
| 563 | 		if (xa <= xa0) |
| 564 | 			break; |
| 565 | 	} |
| 566 | 	return 0; |
| 567 | } |
| 568 | |
| 569 | Bigint *diff (Bigint *a, Bigint *b) |
| 570 | { |
| 571 | 	Bigint *c; |
| 572 | 	int i, wa, wb; |
| 573 | 	ULong *xa, *xae, *xb, *xbe, *xc; |
| 574 | #ifdef ULLong |
| 575 | 	ULLong borrow, y; |
| 576 | #else |
| 577 | 	ULong borrow, y; |
| 578 | #ifdef Pack_32 |
| 579 | 	ULong z; |
| 580 | #endif |
| 581 | #endif |
| 582 | |
| 583 | 	i = cmp(a,b); |
| 584 | 	if (!i) { |
| 585 | 		c = Balloc(0); |
| 586 | 		if (c == NULL) |
| 587 | 			return NULL; |
| 588 | 		c->wds = 1; |
| 589 | 		c->x[0] = 0; |
| 590 | 		return c; |
| 591 | 	} |
| 592 | 	if (i < 0) { |
| 593 | 		c = a; |
| 594 | 		a = b; |
| 595 | 		b = c; |
| 596 | 		i = 1; |
| 597 | 	} |
| 598 | 	else |
| 599 | 		i = 0; |
| 600 | 	c = Balloc(a->k); |
| 601 | 	if (c == NULL) |
| 602 | 		return NULL; |
| 603 | 	c->sign = i; |
| 604 | 	wa = a->wds; |
| 605 | 	xa = a->x; |
| 606 | 	xae = xa + wa; |
| 607 | 	wb = b->wds; |
| 608 | 	xb = b->x; |
| 609 | 	xbe = xb + wb; |
| 610 | 	xc = c->x; |
| 611 | 	borrow = 0; |
| 612 | #ifdef ULLong |
| 613 | 	do { |
| 614 | 		y = (ULLong)*xa++ - *xb++ - borrow; |
| 615 | 		borrow = y >> 32 & 1UL; |
| 616 | 		*xc++ = y & 0xffffffffUL; |
| 617 | 	} while(xb < xbe); |
| 618 | 	while(xa < xae) { |
| 619 | 		y = *xa++ - borrow; |
| 620 | 		borrow = y >> 32 & 1UL; |
| 621 | 		*xc++ = y & 0xffffffffUL; |
| 622 | 	} |
| 623 | #else |
| 624 | #ifdef Pack_32 |
| 625 | 	do { |
| 626 | 		y = (*xa & 0xffff) - (*xb & 0xffff) - borrow; |
| 627 | 		borrow = (y & 0x10000) >> 16; |
| 628 | 		z = (*xa++ >> 16) - (*xb++ >> 16) - borrow; |
| 629 | 		borrow = (z & 0x10000) >> 16; |
| 630 | 		Storeinc(xc, z, y); |
| 631 | 	} while(xb < xbe); |
| 632 | 	while(xa < xae) { |
| 633 | 		y = (*xa & 0xffff) - borrow; |
| 634 | 		borrow = (y & 0x10000) >> 16; |
| 635 | 		z = (*xa++ >> 16) - borrow; |
| 636 | 		borrow = (z & 0x10000) >> 16; |
| 637 | 		Storeinc(xc, z, y); |
| 638 | 	} |
| 639 | #else |
| 640 | 	do { |
| 641 | 		y = *xa++ - *xb++ - borrow; |
| 642 | 		borrow = (y & 0x10000) >> 16; |
| 643 | 		*xc++ = y & 0xffff; |
| 644 | 	} while(xb < xbe); |
| 645 | 	while(xa < xae) { |
| 646 | 		y = *xa++ - borrow; |
| 647 | 		borrow = (y & 0x10000) >> 16; |
| 648 | 		*xc++ = y & 0xffff; |
| 649 | 	} |
| 650 | #endif |
| 651 | #endif |
| 652 | 	while(!*--xc) |
| 653 | 		wa--; |
| 654 | 	c->wds = wa; |
| 655 | 	return c; |
| 656 | } |
| 657 | |
| 658 | double b2d (Bigint *a, int *e) |
| 659 | { |
| 660 | 	ULong *xa, *xa0, w, y, z; |
| 661 | 	int k; |
| 662 | 	union _dbl_union d; |
| 663 | #define d0 word0(&d) |
| 664 | #define d1 word1(&d) |
| 665 | |
| 666 | 	xa0 = a->x; |
| 667 | 	xa = xa0 + a->wds; |
| 668 | 	y = *--xa; |
| 669 | #ifdef DEBUG |
| 670 | 	if (!y) Bug("zero y in b2d"); |
| 671 | #endif |
| 672 | 	k = hi0bits(y); |
| 673 | 	*e = 32 - k; |
| 674 | #ifdef Pack_32 |
| 675 | 	if (k < Ebits) { |
| 676 | 		d0 = Exp_1 | y >> (Ebits - k); |
| 677 | 		w = xa > xa0 ? *--xa : 0; |
| 678 | 		d1 = y << ((32-Ebits) + k) | w >> (Ebits - k); |
| 679 | 		goto ret_d; |
| 680 | 	} |
| 681 | 	z = xa > xa0 ? *--xa : 0; |
| 682 | 	if (k -= Ebits) { |
| 683 | 		d0 = Exp_1 | y << k | z >> (32 - k); |
| 684 | 		y = xa > xa0 ? *--xa : 0; |
| 685 | 		d1 = z << k | y >> (32 - k); |
| 686 | 	} |
| 687 | 	else { |
| 688 | 		d0 = Exp_1 | y; |
| 689 | 		d1 = z; |
| 690 | 	} |
| 691 | #else |
| 692 | 	if (k < Ebits + 16) { |
| 693 | 		z = xa > xa0 ? *--xa : 0; |
| 694 | 		d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k; |
| 695 | 		w = xa > xa0 ? *--xa : 0; |
| 696 | 		y = xa > xa0 ? *--xa : 0; |
| 697 | 		d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k; |
| 698 | 		goto ret_d; |
| 699 | 	} |
| 700 | 	z = xa > xa0 ? *--xa : 0; |
| 701 | 	w = xa > xa0 ? *--xa : 0; |
| 702 | 	k -= Ebits + 16; |
| 703 | 	d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k; |
| 704 | 	y = xa > xa0 ? *--xa : 0; |
| 705 | 	d1 = w << k + 16 | y << k; |
| 706 | #endif |
| 707 | ret_d: |
| 708 | 	return dval(&d); |
| 709 | #undef d0 |
| 710 | #undef d1 |
| 711 | } |
| 712 | |
| 713 | Bigint *d2b (double dd, int *e, int *bits) |
| 714 | { |
| 715 | 	Bigint *b; |
| 716 | 	union _dbl_union d; |
| 717 | #ifndef Sudden_Underflow |
| 718 | 	int i; |
| 719 | #endif |
| 720 | 	int de, k; |
| 721 | 	ULong *x, y, z; |
| 722 | #define d0 word0(&d) |
| 723 | #define d1 word1(&d) |
| 724 | 	d.d = dd; |
| 725 | |
| 726 | #ifdef Pack_32 |
| 727 | 	b = Balloc(1); |
| 728 | #else |
| 729 | 	b = Balloc(2); |
| 730 | #endif |
| 731 | 	if (b == NULL) |
| 732 | 		return NULL; |
| 733 | 	x = b->x; |
| 734 | |
| 735 | 	z = d0 & Frac_mask; |
| 736 | 	d0 &= 0x7fffffff;	/* clear sign bit, which we ignore */ |
| 737 | #ifdef Sudden_Underflow |
| 738 | 	de = (int)(d0 >> Exp_shift); |
| 739 | 	z |= Exp_msk11; |
| 740 | #else |
| 741 | 	if ( (de = (int)(d0 >> Exp_shift)) !=0) |
| 742 | 		z |= Exp_msk1; |
| 743 | #endif |
| 744 | #ifdef Pack_32 |
| 745 | 	if ( (y = d1) !=0) { |
| 746 | 		if ( (k = lo0bits(&y)) !=0) { |
| 747 | 			x[0] = y | z << (32 - k); |
| 748 | 			z >>= k; |
| 749 | 		} |
| 750 | 		else |
| 751 | 			x[0] = y; |
| 752 | #ifndef Sudden_Underflow |
| 753 | 		i = |
| 754 | #endif |
| 755 | 		 b->wds = (x[1] = z) !=0 ? 2 : 1; |
| 756 | 	} |
| 757 | 	else { |
| 758 | 		k = lo0bits(&z); |
| 759 | 		x[0] = z; |
| 760 | #ifndef Sudden_Underflow |
| 761 | 		i = |
| 762 | #endif |
| 763 | 		 b->wds = 1; |
| 764 | 		k += 32; |
| 765 | 	} |
| 766 | #else |
| 767 | 	if ( (y = d1) !=0) { |
| 768 | 		if ( (k = lo0bits(&y)) !=0) |
| 769 | 			if (k >= 16) { |
| 770 | 				x[0] = y | z << 32 - k & 0xffff; |
| 771 | 				x[1] = z >> k - 16 & 0xffff; |
| 772 | 				x[2] = z >> k; |
| 773 | 				i = 2; |
| 774 | 			} |
| 775 | 			else { |
| 776 | 				x[0] = y & 0xffff; |
| 777 | 				x[1] = y >> 16 | z << 16 - k & 0xffff; |
| 778 | 				x[2] = z >> k & 0xffff; |
| 779 | 				x[3] = z >> k+16; |
| 780 | 				i = 3; |
| 781 | 			} |
| 782 | 		else { |
| 783 | 			x[0] = y & 0xffff; |
| 784 | 			x[1] = y >> 16; |
| 785 | 			x[2] = z & 0xffff; |
| 786 | 			x[3] = z >> 16; |
| 787 | 			i = 3; |
| 788 | 		} |
| 789 | 	} |
| 790 | 	else { |
| 791 | #ifdef DEBUG |
| 792 | 		if (!z) |
| 793 | 			Bug("Zero passed to d2b"); |
| 794 | #endif |
| 795 | 		k = lo0bits(&z); |
| 796 | 		if (k >= 16) { |
| 797 | 			x[0] = z; |
| 798 | 			i = 0; |
| 799 | 		} |
| 800 | 		else { |
| 801 | 			x[0] = z & 0xffff; |
| 802 | 			x[1] = z >> 16; |
| 803 | 			i = 1; |
| 804 | 		} |
| 805 | 		k += 32; |
| 806 | 	} |
| 807 | 	while(!x[i]) |
| 808 | 		--i; |
| 809 | 	b->wds = i + 1; |
| 810 | #endif |
| 811 | #ifndef Sudden_Underflow |
| 812 | 	if (de) { |
| 813 | #endif |
| 814 | 		*e = de - Bias - (P-1) + k; |
| 815 | 		*bits = P - k; |
| 816 | #ifndef Sudden_Underflow |
| 817 | 	} |
| 818 | 	else { |
| 819 | 		*e = de - Bias - (P-1) + 1 + k; |
| 820 | #ifdef Pack_32 |
| 821 | 		*bits = 32*i - hi0bits(x[i-1]); |
| 822 | #else |
| 823 | 		*bits = (i+2)*16 - hi0bits(x[i]); |
| 824 | #endif |
| 825 | 	} |
| 826 | #endif |
| 827 | 	return b; |
| 828 | #undef d0 |
| 829 | #undef d1 |
| 830 | } |
| 831 | |
| 832 | const double |
| 833 | bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 }; |
| 834 | const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 }; |
| 835 | |
| 836 | const double |
| 837 | tens[] = { |
| 838 | 		1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, |
| 839 | 		1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, |
| 840 | 		1e20, 1e21, 1e22 |
| 841 | }; |
| 842 | |
| 843 | char *strcp_D2A (char *a, const char *b) |
| 844 | { |
| 845 | 	while((*a = *b++)) |
| 846 | 		a++; |
| 847 | 	return a; |
| 848 | } |
| 849 | |
| 850 | #ifdef NO_STRING_H |
| 851 | void *memcpy_D2A (void *a1, void *b1, size_t len) |
| 852 | { |
| 853 | 	char *a = (char*)a1, *ae = a + len; |
| 854 | 	char *b = (char*)b1, *a0 = a; |
| 855 | 	while(a < ae) |
| 856 | 		*a++ = *b++; |
| 857 | 	return a0; |
| 858 | } |
| 859 | #endif /* NO_STRING_H */ |
| 860 | |