1 /*
2 * bignum.c - multiple precision exact integer arithmetic
3 *
4 * Copyright (c) 2000-2020 Shiro Kawai <shiro@acm.org>
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 *
10 * 1. Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 *
13 * 2. Redistributions in binary form must reproduce the above copyright
14 * notice, this list of conditions and the following disclaimer in the
15 * documentation and/or other materials provided with the distribution.
16 *
17 * 3. Neither the name of the authors nor the names of its contributors
18 * may be used to endorse or promote products derived from this
19 * software without specific prior written permission.
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
24 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
25 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
26 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
27 * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
28 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
29 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
30 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
31 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 */
33
34 /* Bignum library. Not optimized well yet---I think bignum performance
35 * is not very critical for Gauche, except a few special cases (like
36 * the cases used in the numeric I/O routine). So the implementation
37 * emphasizes robustness rather than performance.
38 *
39 * Bignum is represented by ScmBignum structure. There are "normalized"
40 * and "denormalized" bignums. Scheme part only sees the normalized
41 * bignums. Normalized bignum uses the minimum words to represent the
42 * given number, and no normalized bignums for the numbers that can be
43 * representable in fixnum. Most external bignum API expects normalized
44 * bignums, and return normalized bignums. Normalized bignums should
45 * be seen as read-only construct.
46 *
47 * Denormalized bignums are used to keep intermediate results.
48 * Denormalized forms shouldn't "leak out" to the Scheme world; but
49 * can be useful to write efficient routine.
50 */
51 /* Cf: Knuth: The Art of Computer Programming, sectin 4.3 */
52
53 /* They say AIX requires this to be the first thing in the file, so
54 I include gauche/config.h before the real "gauche.h" */
55 #include <gauche/config.h>
56 #ifndef __GNUC__
57 # if HAVE_ALLOCA_H
58 # include <alloca.h>
59 # else
60 # ifdef _AIX
61 #pragma alloca
62 # else
63 # ifndef alloca /* predefined by HP cc +Olibcalls */
64 char *alloca ();
65 # endif
66 # endif
67 # endif
68 #else
69 # if HAVE_ALLOCA_H
70 # include <alloca.h>
71 # endif
72 # if HAVE_MALLOC_H
73 /* MinGW helds alloca() in "malloc.h" instead of "alloca.h" */
74 # include <malloc.h>
75 # endif
76 #endif
77
78 #include <stdlib.h>
79 #include <math.h>
80 #include <limits.h>
81 #define LIBGAUCHE_BODY
82 #include "gauche.h"
83 #include "gauche/priv/arith.h"
84 #include "gauche/bits.h"
85 #include "gauche/bits_inline.h"
86 #include "gauche/bignum.h"
87
88 #undef min
89 #define min(x, y) (((x) < (y))? (x) : (y))
90 #undef max
91 #define max(x, y) (((x) > (y))? (x) : (y))
92
93 static ScmBignum *bignum_rshift(ScmBignum *br, const ScmBignum *bx, int amount);
94 static ScmBignum *bignum_lshift(ScmBignum *br, const ScmBignum *bx, int amount);
95 static int bignum_safe_size_for_add(const ScmBignum *x, const ScmBignum *y);
96 static ScmBignum *bignum_add_int(ScmBignum *br, const ScmBignum *bx, const ScmBignum *by);
97 static ScmBignum *bignum_2scmpl(ScmBignum *br);
98
99 /*---------------------------------------------------------------------
100 * Constructor
101 *
102 * Scm_MakeBignum* always returns bignum, possibly denormalized.
103 */
bignum_clear(ScmBignum * b)104 static ScmBignum *bignum_clear(ScmBignum *b)
105 {
106 for (u_int i=0; i<b->size; i++) b->values[i] = 0;
107 return b;
108 }
109
110 #define BIGNUM_SIZE(size) (sizeof(ScmBignum)+((size)-1)*sizeof(long))
111
make_bignum(int size)112 static ScmBignum *make_bignum(int size)
113 {
114 if (size < 0) Scm_Error("invalid bignum size (internal error): %d", size);
115 if (size > (int)SCM_BIGNUM_MAX_DIGITS) Scm_Error("too large bignum");
116 ScmBignum *b = SCM_NEW_ATOMIC2(ScmBignum*, BIGNUM_SIZE(size));
117 SCM_SET_CLASS(b, SCM_CLASS_INTEGER);
118 b->size = size;
119 b->sign = 1;
120 return bignum_clear(b);
121 }
122
123
124 /* Allocate temporary bignum in the current function's stack frame
125 if alloca() is available. */
126 #ifdef HAVE_ALLOCA
127 #define ALLOC_TEMP_BIGNUM(var_, size_) \
128 (var_) = SCM_BIGNUM(alloca(BIGNUM_SIZE(size_))); \
129 SCM_SET_CLASS(var_, SCM_CLASS_INTEGER); \
130 (var_)->size = (size_); \
131 (var_)->sign = 1; \
132 bignum_clear(var_)
133 #else /*!HAVE_ALLOCA*/
134 #define ALLOC_TEMP_BIGNUM(var_, size_) (var_) = make_bignum(size_);
135 #endif /*!HAVE_ALLOCA*/
136
Scm_MakeBignumFromSI(long val)137 ScmObj Scm_MakeBignumFromSI(long val)
138 {
139 ScmBignum *b;
140 if (val == LONG_MIN) {
141 b = make_bignum(1);
142 b->sign = -1;
143 b->values[0] = (unsigned long)LONG_MAX+1;
144 } else if (val < 0) {
145 b = make_bignum(1);
146 b->sign = -1;
147 b->values[0] = -val;
148 } else {
149 b = make_bignum(1);
150 b->sign = 1;
151 b->values[0] = val;
152 }
153 return SCM_OBJ(b);
154 }
155
Scm_MakeBignumFromUI(u_long val)156 ScmObj Scm_MakeBignumFromUI(u_long val)
157 {
158 ScmBignum *b = make_bignum(1);
159 b->sign = 1;
160 b->values[0] = val;
161 return SCM_OBJ(b);
162 }
163
164 /* If sign > 0 or sign < 0, values[] has absolute value.
165 If sign == 0, values[] has 2's complement signed representation */
Scm_MakeBignumFromUIArray(int sign,const u_long * values,int size)166 ScmObj Scm_MakeBignumFromUIArray(int sign, const u_long *values, int size)
167 {
168 ScmBignum *b = make_bignum(size);
169 if (sign != 0) {
170 b->sign = (sign > 0)? 1 : -1;
171 for (int i=0; i<size; i++) b->values[i] = values[i];
172 } else {
173 int nonzerop = FALSE;
174 for (int i=0; i<size; i++) {
175 if ((b->values[i] = values[i]) != 0) nonzerop = TRUE;
176 }
177 if (nonzerop) {
178 if (values[size-1] <= LONG_MAX) b->sign = 1;
179 else { b->sign = -1; bignum_2scmpl(b); }
180 } else {
181 b->sign = 0;
182 }
183 }
184 return SCM_OBJ(b);
185 }
186
Scm_MakeBignumFromDouble(double val)187 ScmObj Scm_MakeBignumFromDouble(double val)
188 {
189 if (LONG_MIN <= val
190 #if SIZEOF_LONG == 4
191 && val <= LONG_MAX
192 #else
193 && val <= nextafter((double)LONG_MAX, 0.0)
194 #endif
195 )
196 return Scm_MakeBignumFromSI((long)val);
197
198 int exponent, sign;
199 ScmObj mantissa = Scm_DecodeFlonum(val, &exponent, &sign);
200 if (!SCM_NUMBERP(mantissa)) {
201 Scm_Error("can't convert %lf to an integer", val);
202 }
203 ScmObj b = Scm_Ash(mantissa, exponent);
204 if (sign < 0) b = Scm_Negate(b);
205 /* always returns bignum */
206 if (SCM_INTP(b)) {
207 return Scm_MakeBignumFromSI(SCM_INT_VALUE(b));
208 } else {
209 return b;
210 }
211 }
212
Scm_BignumCopy(const ScmBignum * b)213 ScmObj Scm_BignumCopy(const ScmBignum *b)
214 {
215 ScmBignum *c = make_bignum(b->size);
216 c->sign = b->sign;
217 for (u_int i=0; i<b->size; i++) c->values[i] = b->values[i];
218 return SCM_OBJ(c);
219 }
220
221 /*-----------------------------------------------------------------------
222 * Conversion
223 */
224
225 /* Modifies B and return it. */
Scm_NormalizeBignum(ScmBignum * b)226 ScmObj Scm_NormalizeBignum(ScmBignum *b)
227 {
228 int size = b->size;
229 int i;
230 for (i=size-1; i>0; i--) {
231 if (b->values[i] == 0) size--;
232 else break;
233 }
234 if (i==0) {
235 if (b->sign == 0) {
236 return SCM_MAKE_INT(0);
237 }
238 if (b->sign > 0 && b->values[0] <= (u_long)SCM_SMALL_INT_MAX) {
239 return SCM_MAKE_INT(b->values[0]);
240 }
241 if (b->sign < 0 && b->values[0] <= (u_long)-SCM_SMALL_INT_MIN) {
242 return SCM_MAKE_INT(-((signed long)b->values[0]));
243 }
244 }
245 b->size = size;
246 return SCM_OBJ(b);
247 }
248
249 /* b must be normalized. */
Scm_BignumToSI(const ScmBignum * b,int clamp,int * oor)250 long Scm_BignumToSI(const ScmBignum *b, int clamp, int *oor)
251 {
252 if (clamp == SCM_CLAMP_NONE && oor != NULL) *oor = FALSE;
253 if (b->sign >= 0) {
254 if (b->values[0] > LONG_MAX || b->size >= 2) {
255 if (clamp & SCM_CLAMP_HI) return LONG_MAX;
256 else goto err;
257 } else {
258 return (long)b->values[0];
259 }
260 } else {
261 if (b->values[0] > (u_long)LONG_MAX+1 || b->size >= 2) {
262 if (clamp & SCM_CLAMP_LO) return LONG_MIN;
263 else goto err;
264 } else {
265 return -(long)b->values[0];
266 }
267 }
268 err:
269 if (clamp == SCM_CLAMP_NONE && oor != NULL) {
270 *oor = TRUE;
271 } else {
272 Scm_Error("argument out of range: %S", SCM_OBJ(b));
273 }
274 return 0;
275 }
276
277 /* b must be normalized. */
Scm_BignumToUI(const ScmBignum * b,int clamp,int * oor)278 u_long Scm_BignumToUI(const ScmBignum *b, int clamp, int *oor)
279 {
280 if (clamp == SCM_CLAMP_NONE && oor != NULL) *oor = FALSE;
281 if (b->sign >= 0) {
282 if (b->size >= 2) {
283 if (clamp & SCM_CLAMP_HI) return SCM_ULONG_MAX;
284 else goto err;
285 } else {
286 return b->values[0];
287 }
288 } else {
289 if (clamp & SCM_CLAMP_LO) return 0;
290 else goto err;
291 }
292 err:
293 if (clamp == SCM_CLAMP_NONE && oor != NULL) {
294 *oor = TRUE;
295 } else {
296 Scm_Error("argument out of range: %S", SCM_OBJ(b));
297 }
298 return 0;
299 }
300
301 #if SIZEOF_LONG == 4
302 /* we need special routines for int64 */
Scm_BignumToSI64(const ScmBignum * b,int clamp,int * oor)303 int64_t Scm_BignumToSI64(const ScmBignum *b, int clamp, int *oor)
304 {
305 int64_t r = 0;
306 if (clamp == SCM_CLAMP_NONE && oor != NULL) *oor = FALSE;
307 if (b->sign > 0) {
308 if (b->size == 1) {
309 r = b->values[0];
310 } else if (b->size > 2 || b->values[1] > LONG_MAX) {
311 if (!(clamp & SCM_CLAMP_HI)) goto err;
312 SCM_SET_INT64_MAX(r);
313 } else {
314 r = ((int64_t)b->values[1] << 32) + (uint64_t)b->values[0];
315 }
316 } else { /* b->sign < 0 */
317 if (b->size == 1) {
318 r = -(int64_t)b->values[0];
319 } else if (b->size > 2 || (b->values[1] > LONG_MAX && b->values[0] > 0)) {
320 if (!(clamp&SCM_CLAMP_LO)) goto err;
321 SCM_SET_INT64_MIN(r);
322 } else {
323 r = -(int64_t)(((int64_t)b->values[1] << 32) + (uint64_t)b->values[0]);
324 }
325 }
326 return r;
327 err:
328 if (clamp == SCM_CLAMP_NONE && oor != NULL) {
329 *oor = TRUE;
330 } else {
331 Scm_Error("argument out of range: %S", SCM_OBJ(b));
332 }
333 return r;
334 }
335
Scm_BignumToUI64(const ScmBignum * b,int clamp,int * oor)336 uint64_t Scm_BignumToUI64(const ScmBignum *b, int clamp, int *oor)
337 {
338 uint64_t r = 0;
339 if (clamp == SCM_CLAMP_NONE && oor != NULL) *oor = FALSE;
340 if (b->sign > 0) {
341 if (b->size > 2) {
342 if (!(clamp&SCM_CLAMP_HI)) goto err;
343 SCM_SET_UINT64_MAX(r);
344 } else if (b->size == 2) {
345 r = (((uint64_t)b->values[1]) << 32) + (uint64_t)b->values[0];
346 } else {
347 r = (uint64_t)b->values[0];
348 }
349 } else { /* b->sign < 0 */
350 if (!(clamp&SCM_CLAMP_LO)) goto err;
351 }
352 return r;
353 err:
354 if (clamp == SCM_CLAMP_NONE && oor != NULL) {
355 *oor = TRUE;
356 } else {
357 Scm_Error("argument out of range: %S", SCM_OBJ(b));
358 }
359 return r;
360 }
361 #endif /* SIZEOF_LONG == 4 */
362
363 extern double Scm__EncodeDouble(u_long, u_long, int, int);
364
365 /* Converts a bignum b to a double. b must be normalized.
366 We don't rely on double arithmetic, for it may result
367 an error in the LSB from multiple roundings. Instead we
368 extract bits directly from the bignum values array.
369 (We use ScmBits API by casting b->values to ScmBits*).
370 */
Scm_BignumToDouble(const ScmBignum * b)371 double Scm_BignumToDouble(const ScmBignum *b)
372 {
373 ScmBits *bits = (ScmBits*)b->values;
374 ScmBits dst[2];
375
376 /* first, filter out a special case. */
377 if (b->size == 0) return 0.0;
378
379 int maxbit = Scm_BitsHighest1(bits, 0, b->size*WORD_BITS);
380 int exponent = maxbit+1023;
381 #if SIZEOF_LONG >= 8
382 SCM_ASSERT(maxbit >= 54); /* because b is normalized */
383 dst[0] = 0;
384 Scm_BitsCopyX(dst, 0, bits, maxbit-52, maxbit);
385 /* Rounding. We have to round up if maxbit-53 == 1, EXCEPT
386 the special case where maxbit-52==0, maxbit-53==1, and all
387 bits below are 0 (round-to-even rule) */
388 if (SCM_BITS_TEST(bits, maxbit-53)
389 && ((dst[0]&1) == 1
390 || (Scm_BitsCount1(bits, 0, maxbit-53) > 0))) {
391 dst[0]++;
392 if (dst[0] >= (1UL<<52)) {
393 /* Overflow. We mask the hidden bit, then shift. */
394 dst[0] &= ~(1UL<<52);
395 dst[0] >>= 1;
396 exponent++;
397 }
398 }
399 if (exponent > 2046) return Scm__EncodeDouble(0, 0, 2047, (b->sign < 0));
400 else return Scm__EncodeDouble(dst[0], 0, exponent, (b->sign < 0));
401 #else /*SIZEOF_LONG == 4 */
402 dst[0] = dst[1] = 0;
403 if (maxbit < 53) {
404 Scm_BitsCopyX(dst, 52-maxbit, bits, 0, maxbit);
405 } else {
406 Scm_BitsCopyX(dst, 0, bits, maxbit-52, maxbit);
407 /* Rounding. See the above comment. */
408 if (SCM_BITS_TEST(bits, maxbit-53)
409 && ((dst[0]&1) == 1
410 || (maxbit > 53 && Scm_BitsCount1(bits, 0, maxbit-53) > 0))) {
411 u_long inc = dst[0] + 1;
412 if (inc < dst[0]) dst[1]++;
413 dst[0] = inc;
414 if (dst[1] >= (1UL<<(52-32))) {
415 /* Overflow. We mask the hidden bit, then shift. */
416 dst[1] &= ~(1UL<<(52-32));
417 dst[0] = (dst[0] >> 1) | ((u_long)(dst[1]&1) << 31);
418 dst[1] >>= 1;
419 exponent++;
420 }
421 }
422 }
423 if (exponent > 2046) return Scm__EncodeDouble(0, 0, 2047, (b->sign < 0));
424 else return Scm__EncodeDouble(dst[0], dst[1], exponent, (b->sign < 0));
425 #endif /*SIZEOF_LONG==4*/
426 }
427
428 /* return -b, normalized */
Scm_BignumNegate(const ScmBignum * b)429 ScmObj Scm_BignumNegate(const ScmBignum *b)
430 {
431 ScmObj c = Scm_BignumCopy(b);
432 SCM_BIGNUM_SIGN(c) = -SCM_BIGNUM_SIGN(c);
433 return Scm_NormalizeBignum(SCM_BIGNUM(c));
434 }
435
436 /*-----------------------------------------------------------------------
437 * Compare
438 */
439
440 /* bx and by must be normalized */
Scm_BignumCmp(const ScmBignum * bx,const ScmBignum * by)441 int Scm_BignumCmp(const ScmBignum *bx, const ScmBignum *by)
442 {
443 if (bx->sign < by->sign) return -1;
444 if (bx->sign > by->sign) return 1;
445 if (bx->size < by->size) return (bx->sign > 0) ? -1 : 1;
446 if (bx->size > by->size) return (bx->sign > 0) ? 1 : -1;
447
448 for (int i=(int)bx->size-1; i>=0; i--) {
449 if (bx->values[i] < by->values[i]) return (bx->sign > 0) ? -1 : 1;
450 if (bx->values[i] > by->values[i]) return (bx->sign > 0) ? 1 : -1;
451 }
452 return 0;
453 }
454
455 /* compare absolute values. assume bx and by are nomalized. */
Scm_BignumAbsCmp(const ScmBignum * bx,const ScmBignum * by)456 int Scm_BignumAbsCmp(const ScmBignum *bx, const ScmBignum *by)
457 {
458 if (bx->size < by->size) return -1;
459 if (bx->size > by->size) return 1;
460 for (int i=(int)bx->size-1; i>=0; i--) {
461 if (bx->values[i] < by->values[i]) return -1;
462 if (bx->values[i] > by->values[i]) return 1;
463 }
464 return 0;
465 }
466
467 /* Compare bx + off and by. All arguments must be positive. bx and by
468 must be normalized. off may be denormalized if it is created directly
469 by Scm_MakeBignumFromUI call.
470 Expect bx >> off for most cases.
471 Screen out the obvious case without actually calculating bx+off.
472 Experimentary, the following set of conditions avoid 93% of cases from
473 doing actual bignum addition. */
Scm_BignumCmp3U(const ScmBignum * bx,const ScmBignum * off,const ScmBignum * by)474 int Scm_BignumCmp3U(const ScmBignum *bx,
475 const ScmBignum *off,
476 const ScmBignum *by)
477 {
478 u_int xsize = SCM_BIGNUM_SIZE(bx), ysize = SCM_BIGNUM_SIZE(by);
479 u_int osize = SCM_BIGNUM_SIZE(off);
480
481 if (xsize > ysize) return 1;
482 if (xsize < ysize) {
483 if (osize < ysize && by->values[ysize-1] > 1) {
484 return -1;
485 }
486 if (osize == ysize) {
487 if (off->values[osize-1] > by->values[ysize-1]) return 1;
488 if (off->values[osize-1] < by->values[ysize-1]-1) return -1;
489 }
490 /* fallthrough */
491 } else {
492 /* xsize == ysize */
493 u_long w, c = 0;
494 if (osize > ysize) return 1;
495 if (bx->values[xsize-1] > by->values[ysize-1]) return 1;
496 if (osize < xsize) {
497 if (bx->values[xsize-1] < by->values[ysize-1]-1) return -1;
498 } else {
499 /* osize == xsize */
500 u_long xx = bx->values[xsize-1], oo = off->values[osize-1];
501 UADD(w, c, xx, oo);
502 if (c > 0 || w > by->values[ysize-1]) return 1;
503 if (w < by->values[ysize-1] - 1) return -1;
504 }
505 /* fallthrough */
506 }
507 u_int tsize = bignum_safe_size_for_add(bx, off);
508 ScmBignum *br;
509 ALLOC_TEMP_BIGNUM(br, tsize);
510 bignum_add_int(br, bx, off);
511
512 if (br->size < by->size) return -1;
513 for (int i=(int)br->size-1; i>=0; i--) {
514 if (i >= (int)by->size) {
515 if (br->values[i]) return 1;
516 continue;
517 }
518 if (br->values[i] < by->values[i]) return -1;
519 if (br->values[i] > by->values[i]) return 1;
520 }
521 return 0;
522 }
523
524 /*-----------------------------------------------------------------------
525 * Add & subtract
526 */
bignum_safe_size_for_add(const ScmBignum * x,const ScmBignum * y)527 static int bignum_safe_size_for_add(const ScmBignum *x, const ScmBignum *y)
528 {
529 int xsize = SCM_BIGNUM_SIZE(x);
530 int ysize = SCM_BIGNUM_SIZE(y);
531 if (xsize > ysize) {
532 if (x->values[xsize-1] == SCM_ULONG_MAX) return xsize+1;
533 else return xsize;
534 } else if (ysize > xsize) {
535 if (y->values[ysize-1] == SCM_ULONG_MAX) return ysize+1;
536 else return ysize;
537 } else {
538 return xsize+1;
539 }
540 }
541
542 /* take 2's complement */
bignum_2scmpl(ScmBignum * br)543 static ScmBignum *bignum_2scmpl(ScmBignum *br)
544 {
545 int rsize = SCM_BIGNUM_SIZE(br);
546 u_long c = 1;
547 for (int i=0; i<rsize; i++) {
548 unsigned long x = ~br->values[i];
549 UADD(br->values[i], c, x, 0);
550 }
551 return br;
552 }
553
Scm_BignumComplement(const ScmBignum * bx)554 ScmObj Scm_BignumComplement(const ScmBignum *bx)
555 {
556 ScmBignum *r = SCM_BIGNUM(Scm_BignumCopy(bx));
557 return SCM_OBJ(bignum_2scmpl(r));
558 }
559
560 /* br = abs(bx) + abs(by), assuming br has enough size. br and bx can be
561 the same object. */
bignum_add_int(ScmBignum * br,const ScmBignum * bx,const ScmBignum * by)562 static ScmBignum *bignum_add_int(ScmBignum *br,
563 const ScmBignum *bx, const ScmBignum *by)
564 {
565 int rsize = SCM_BIGNUM_SIZE(br);
566 int xsize = SCM_BIGNUM_SIZE(bx);
567 int ysize = SCM_BIGNUM_SIZE(by);
568 u_long c = 0;
569
570 for (int i=0; i<rsize; i++, xsize--, ysize--) {
571 if (xsize <= 0) {
572 if (ysize <= 0) {
573 UADD(br->values[i], c, 0, 0);
574 continue;
575 }
576 u_long y = by->values[i];
577 UADD(br->values[i], c, 0, y);
578 continue;
579 }
580 if (ysize <= 0) {
581 u_long x = bx->values[i];
582 UADD(br->values[i], c, x, 0);
583 continue;
584 }
585 u_long x = bx->values[i];
586 u_long y = by->values[i];
587 UADD(br->values[i], c, x, y);
588 }
589 return br;
590 }
591
592 /* br = abs(bx) - abs(by), assuming br has enough size. br and bx can be
593 the same object. */
bignum_sub_int(ScmBignum * br,const ScmBignum * bx,const ScmBignum * by)594 static ScmBignum *bignum_sub_int(ScmBignum *br,
595 const ScmBignum *bx, const ScmBignum *by)
596 {
597 int rsize = SCM_BIGNUM_SIZE(br);
598 int xsize = SCM_BIGNUM_SIZE(bx);
599 int ysize = SCM_BIGNUM_SIZE(by);
600 u_long c = 0;
601
602 for (int i=0; i<rsize; i++, xsize--, ysize--) {
603 if (xsize <= 0) {
604 if (ysize <= 0) {
605 USUB(br->values[i], c, 0, 0);
606 continue;
607 }
608 u_long y = by->values[i];
609 USUB(br->values[i], c, 0, y);
610 continue;
611 }
612 if (ysize <= 0) {
613 u_long x = bx->values[i];
614 USUB(br->values[i], c, x, 0);
615 continue;
616 }
617 u_long x = bx->values[i];
618 u_long y = by->values[i];
619 USUB(br->values[i], c, x, y);
620 }
621 if (c != 0) {
622 bignum_2scmpl(br);
623 br->sign = 0 - br->sign; /* flip sign */
624 }
625 return br;
626 }
627
628 /* returns bx + by, not normalized */
bignum_add(const ScmBignum * bx,const ScmBignum * by)629 static ScmBignum *bignum_add(const ScmBignum *bx, const ScmBignum *by)
630 {
631 int rsize = bignum_safe_size_for_add(bx, by);
632 ScmBignum *br = make_bignum(rsize);
633 br->sign = SCM_BIGNUM_SIGN(bx);
634 if (SCM_BIGNUM_SIGN(bx) == SCM_BIGNUM_SIGN(by)) {
635 bignum_add_int(br, bx, by);
636 } else {
637 bignum_sub_int(br, bx, by);
638 }
639 return br;
640 }
641
642 /* returns bx - by, not normalized */
bignum_sub(const ScmBignum * bx,const ScmBignum * by)643 static ScmBignum *bignum_sub(const ScmBignum *bx, const ScmBignum *by)
644 {
645 int rsize = bignum_safe_size_for_add(bx, by);
646 ScmBignum *br = make_bignum(rsize);
647 br->sign = SCM_BIGNUM_SIGN(bx);
648 if (SCM_BIGNUM_SIGN(bx) == SCM_BIGNUM_SIGN(by)) {
649 bignum_sub_int(br, bx, by);
650 } else {
651 bignum_add_int(br, bx, by);
652 }
653 return br;
654 }
655
656 /* returns bx + y, not nomalized */
bignum_add_si(const ScmBignum * bx,long y)657 static ScmBignum *bignum_add_si(const ScmBignum *bx, long y)
658 {
659 long c = 0;
660 u_int rsize = bx->size+1;
661 u_long yabs = ((y < 0)? -y : y);
662 int ysign = ((y < 0)? -1 : 1);
663 ScmBignum *br = make_bignum(rsize);
664 br->sign = bx->sign;
665 if (SCM_BIGNUM_SIGN(bx) == ysign) {
666 for (u_int i=0; i<bx->size; i++) {
667 UADD(br->values[i], c, bx->values[i], yabs);
668 yabs = 0;
669 }
670 } else {
671 for (u_int i=0; i<bx->size; i++) {
672 USUB(br->values[i], c, bx->values[i], yabs);
673 yabs = 0;
674 }
675 }
676 br->values[rsize-1] = c;
677 return br;
678 }
679
Scm_BignumAdd(const ScmBignum * bx,const ScmBignum * by)680 ScmObj Scm_BignumAdd(const ScmBignum *bx, const ScmBignum *by)
681 {
682 return Scm_NormalizeBignum(bignum_add(bx, by));
683 }
684
Scm_BignumSub(const ScmBignum * bx,const ScmBignum * by)685 ScmObj Scm_BignumSub(const ScmBignum *bx, const ScmBignum *by)
686 {
687 return Scm_NormalizeBignum(bignum_sub(bx, by));
688 }
689
Scm_BignumAddSI(const ScmBignum * bx,long y)690 ScmObj Scm_BignumAddSI(const ScmBignum *bx, long y)
691 {
692 if (y == 0) return SCM_OBJ(bx);
693 else return Scm_NormalizeBignum(bignum_add_si(bx, y));
694 }
695
Scm_BignumSubSI(const ScmBignum * bx,long y)696 ScmObj Scm_BignumSubSI(const ScmBignum *bx, long y)
697 {
698 if (y == 0) return SCM_OBJ(bx);
699 else return Scm_NormalizeBignum(bignum_add_si(bx, -y));
700 }
701
702 /*-----------------------------------------------------------------------
703 * Shifter
704 */
705
706 /* br = bx >> amount. amount >= 0. no normalization. assumes br
707 has enough size to hold the result. br and bx can be the same object. */
bignum_rshift(ScmBignum * br,const ScmBignum * bx,int amount)708 static ScmBignum *bignum_rshift(ScmBignum *br, const ScmBignum *bx, int amount)
709 {
710 u_int nwords = amount / WORD_BITS;
711 u_int nbits = amount % WORD_BITS;
712 int i;
713
714 if (bx->size <= nwords) {
715 br->size = 0; br->values[0] = 0;
716 } else if (nbits == 0) {
717 for (i = (int)nwords; i < (int)bx->size; i++) {
718 br->values[i-nwords] = bx->values[i];
719 }
720 br->size = bx->size - nwords;
721 br->sign = bx->sign;
722 } else {
723 u_long x;
724 for (i = (int)nwords; i < (int)bx->size-1; i++) {
725 x = (bx->values[i+1]<<(WORD_BITS-nbits))|(bx->values[i]>>nbits);
726 br->values[i-nwords] = x;
727 }
728 br->values[i-nwords] = bx->values[i]>>nbits;
729 br->size = bx->size - nwords;
730 br->sign = bx->sign;
731 }
732 return br;
733 }
734
735 /* br = bx << amount, amount > 0. no normalization. assumes br
736 has enough size. br and bx can be the same object. */
bignum_lshift(ScmBignum * br,const ScmBignum * bx,int amount)737 static ScmBignum *bignum_lshift(ScmBignum *br, const ScmBignum *bx, int amount)
738 {
739 int nwords = amount / WORD_BITS;
740 int nbits = amount % WORD_BITS;
741
742 if (nbits == 0) {
743 /* short path */
744 for (int i = (int)bx->size-1; i>=0; i--) {
745 if ((int)br->size > i+nwords) br->values[i+nwords] = bx->values[i];
746 }
747 for (int i = nwords-1; i>=0; i--) br->values[i] = 0;
748 } else {
749 if (br->size > bx->size + nwords) {
750 br->values[bx->size+nwords] =
751 bx->values[bx->size-1]>>(WORD_BITS-nbits);
752 }
753 int i;
754 for (i = (int)bx->size-1; i > 0; i--) {
755 u_long x = (bx->values[i]<<nbits)|(bx->values[i-1]>>(WORD_BITS-nbits));
756 if ((int)br->size > i+nwords) br->values[i+nwords] = x;
757 }
758 br->values[nwords] = bx->values[0]<<nbits;
759 for (i = nwords-1; i>=0; i--) br->values[i] = 0;
760 }
761 if (br != bx) {
762 br->sign = bx->sign;
763 }
764 return br;
765 }
766
767 /*-----------------------------------------------------------------------
768 * Multiplication
769 */
770
771 /* br += bx * (y << off*WORD_BITS). br must have enough size. */
bignum_mul_word(ScmBignum * br,const ScmBignum * bx,u_long y,int off)772 static ScmBignum *bignum_mul_word(ScmBignum *br, const ScmBignum *bx,
773 u_long y, int off)
774 {
775 for (int i=0; i<bx->size; i++) {
776 u_long hi, lo, r1;
777 u_long x = bx->values[i];
778 UMUL(hi, lo, x, y);
779
780 u_long c = 0;
781 u_long r0 = br->values[i+off];
782 UADD(r1, c, r0, lo);
783 br->values[i+off] = r1;
784
785 r0 = br->values[i+off+1];
786 UADD(r1, c, r0, hi);
787 br->values[i+off+1] = r1;
788
789 for (int j=i+off+2; c && j<br->size; j++) {
790 r0 = br->values[j];
791 UADD(r1, c, r0, 0);
792 br->values[j] = r1;
793 }
794 }
795 return br;
796 }
797
798 /* returns bx * by. not normalized */
bignum_mul(const ScmBignum * bx,const ScmBignum * by)799 static ScmBignum *bignum_mul(const ScmBignum *bx, const ScmBignum *by)
800 {
801 ScmBignum *br = make_bignum(bx->size + by->size);
802 for (u_int i=0; i<by->size; i++) {
803 bignum_mul_word(br, bx, by->values[i], i);
804 }
805 br->sign = bx->sign * by->sign;
806 return br;
807 }
808
809 /* return bx * y, y != 0 and y != 1 */
bignum_mul_si(const ScmBignum * bx,long y)810 static ScmBignum *bignum_mul_si(const ScmBignum *bx, long y)
811 {
812 if (y == -1) {
813 ScmBignum *br = SCM_BIGNUM(Scm_BignumCopy(bx));
814 br->sign = -br->sign;
815 return br;
816 }
817 /* TODO: optimize for 2^n case !*/
818 ScmBignum *br;
819 br = make_bignum(bx->size + 1); /* TODO: more accurate estimation */
820 u_long yabs = (y<0)? -y:y;
821 br->sign = bx->sign;
822 bignum_mul_word(br, bx, yabs, 0);
823 if (y<0) br->sign = -br->sign;
824 return br;
825 }
826
Scm_BignumMul(const ScmBignum * bx,const ScmBignum * by)827 ScmObj Scm_BignumMul(const ScmBignum *bx, const ScmBignum *by)
828 {
829 ScmBignum *br = bignum_mul(bx, by);
830 return Scm_NormalizeBignum(br);
831 }
832
Scm_BignumMulSI(const ScmBignum * bx,long y)833 ScmObj Scm_BignumMulSI(const ScmBignum *bx, long y)
834 {
835 if (y == 1) return SCM_OBJ(bx);
836 else if (y == 0) return SCM_MAKE_INT(0);
837 else return Scm_NormalizeBignum(bignum_mul_si(bx, y));
838 }
839
840 /*-----------------------------------------------------------------------
841 * Division
842 */
843
844 /* returns # of bits in the leftmost '1' in the word, counting from MSB. */
div_normalization_factor(u_long w)845 static inline int div_normalization_factor(u_long w)
846 {
847 u_long b = (1L<<(WORD_BITS-1)), c = 0;
848 for (; b > 0; b>>=1, c++) {
849 if (w & b) return c;
850 }
851 /* something got wrong here */
852 Scm_Panic("bignum.c: div_normalization_factor: can't be here");
853 return 0; /* dummy */
854 }
855
856 /* General case of division. We use each half word as a digit.
857 Assumes digitsof(dividend) >= digitsof(divisor) > 1.
858 Assumes enough digits are allocated to quotient.
859 Remainder is returned (not normalized) */
bignum_gdiv(const ScmBignum * dividend,const ScmBignum * divisor,ScmBignum * quotient)860 static ScmBignum *bignum_gdiv(const ScmBignum *dividend,
861 const ScmBignum *divisor,
862 ScmBignum *quotient)
863 {
864 int d = div_normalization_factor(divisor->values[divisor->size-1]);
865 int n, m;
866 u_long vv, uj, uj2, cy;
867
868 #define DIGIT(num, n) (((n)%2)? HI((num)->values[(n)/2]) : LO((num)->values[(n)/2]))
869 #define DIGIT2(num, n) \
870 (((n)%2)? \
871 ((LO((num)->values[(n)/2+1])<<HALF_BITS)|HI((num)->values[(n)/2])): \
872 (num)->values[(n)/2])
873 #define SETDIGIT(num, n, v) \
874 (((n)%2)? \
875 (num->values[(n)/2] = (num->values[(n)/2] & LOMASK)|((v) << HALF_BITS)) :\
876 (num->values[(n)/2] = (num->values[(n)/2] & HIMASK)|((v) & LOMASK)))
877 #define SETDIGIT2(num, n, v) \
878 (((n)%2)? \
879 ((num->values[(n)/2] = LO(num->values[(n)/2])|((v)<<HALF_BITS)), \
880 (num->values[(n)/2+1] = (num->values[(n)/2+1] & HIMASK)|HI(v))) : \
881 (num->values[(n)/2] = (v)))
882
883 /* normalize */
884 ScmBignum *u, *v;
885 u = make_bignum(dividend->size + 1); /* will be returned as a remainder */
886 ALLOC_TEMP_BIGNUM(v, divisor->size);
887 if (d >= HALF_BITS) {
888 d -= HALF_BITS;
889 n = divisor->size*2 - 1;
890 m = dividend->size*2 - n;
891 } else {
892 n = divisor->size*2;
893 m = dividend->size*2 - n;
894 }
895 bignum_lshift(u, dividend, d);
896 bignum_lshift(v, divisor, d);
897 u_long vn_1 = DIGIT(v, n-1);
898 u_long vn_2 = DIGIT(v, n-2);
899 #undef DIV_DEBUG
900 #ifdef DIV_DEBUG
901 Scm_Printf(SCM_CUROUT, "shift=%d, n=%d, m=%d\n", d, n, m);
902 Scm_Printf(SCM_CUROUT, "u="); Scm_DumpBignum(u, SCM_CUROUT);
903 Scm_Printf(SCM_CUROUT, "\nv="); Scm_DumpBignum(v, SCM_CUROUT);
904 Scm_Printf(SCM_CUROUT, "\nvn_1=%08lx, vn_2=%08lx\n", vn_1, vn_2);
905 #endif
906
907 for (int j = m; j >= 0; j--) {
908 u_long uu = (DIGIT(u, j+n) << HALF_BITS) + DIGIT(u, j+n-1);
909 u_long qq = uu/vn_1;
910 u_long rr = uu%vn_1;
911 #ifdef DIV_DEBUG
912 Scm_Printf(SCM_CUROUT, "loop on j=%d, uu=%08lx, qq=%08lx, rr=%08lx\n",
913 j, uu, qq, rr);
914 #endif
915 while (qq >= HALF_WORD) { qq--; rr += vn_1; }
916 while ((qq*vn_2 > (rr<<HALF_BITS)+DIGIT(u, j+n-2)) && (rr < HALF_WORD)) {
917 qq--; rr += vn_1;
918 }
919 #ifdef DIV_DEBUG
920 Scm_Printf(SCM_CUROUT, "estimate uu=%08lx, qq=%08lx, rr=%08lx\n",
921 uu, qq, rr);
922 #endif
923 cy = 0;
924 for (int k = 0; k < n; k++) {
925 vv = qq * DIGIT(v, k);
926 uj = DIGIT2(u, j+k);
927 uj2 = uj - vv - cy;
928 cy = (uj2 > uj)? HALF_WORD : 0;
929 SETDIGIT2(u, j+k, uj2);
930 }
931 #ifdef DIV_DEBUG
932 Scm_Printf(SCM_CUROUT, "subtract cy = %d, ", cy);
933 Scm_Printf(SCM_CUROUT, "u="); Scm_DumpBignum(u, SCM_CUROUT);
934 Scm_Printf(SCM_CUROUT, "\n");
935 #endif
936 if (cy) {
937 qq--;
938 cy = 0;
939 for (int k = 0; k < n; k++) {
940 vv = DIGIT(v, k);
941 uj = DIGIT(u, j+k) + vv + cy;
942 cy = (uj >= HALF_WORD)? 1 : 0;
943 SETDIGIT(u, j+k, uj);
944 }
945 uj = DIGIT(u, j+n) + cy;
946 SETDIGIT(u, j+n, uj);
947 }
948 SETDIGIT(quotient, j, qq);
949 }
950 bignum_rshift(u, u, d);
951 #ifdef DIV_DEBUG
952 Scm_Printf(SCM_CUROUT, "quot q="); Scm_DumpBignum(quotient, SCM_CUROUT);
953 Scm_Printf(SCM_CUROUT, "\nrem u="); Scm_DumpBignum(u, SCM_CUROUT);
954 Scm_Printf(SCM_CUROUT, "\n");
955 #endif
956 return u;
957 }
958
959 /* Fast path if divisor fits in a half word. Quotient remains in the
960 dividend's memory. Remainder returned. Quotient not normalized. */
bignum_sdiv(ScmBignum * dividend,u_long divisor)961 static u_long bignum_sdiv(ScmBignum *dividend, u_long divisor)
962 {
963 int n = dividend->size - 1;
964 u_long *pu = dividend->values;
965 u_long q0 = 0, r0 = 0, q1, r1;
966
967 for (; n > 0; n--) {
968 q1 = pu[n] / divisor + q0;
969 r1 = ((pu[n] % divisor) << HALF_BITS) + HI(pu[n-1]);
970 q0 = ((r1 / divisor) << HALF_BITS);
971 r0 = r1 % divisor;
972 pu[n] = q1;
973 pu[n-1] = (r0 << HALF_BITS) + LO(pu[n-1]);
974 }
975 q1 = pu[0] / divisor + q0;
976 r1 = pu[0] % divisor;
977 pu[0] = q1;
978 return r1;
979 }
980
981 /* assuming dividend is normalized. */
Scm_BignumDivSI(const ScmBignum * dividend,long divisor,long * remainder)982 ScmObj Scm_BignumDivSI(const ScmBignum *dividend, long divisor, long *remainder)
983 {
984 u_long dd = (divisor < 0)? -divisor : divisor;
985 u_long rr;
986 int d_sign = (divisor < 0)? -1 : 1;
987 ScmBignum *q;
988
989 if (dd < HALF_WORD) {
990 q = SCM_BIGNUM(Scm_BignumCopy(dividend));
991 rr = bignum_sdiv(q, dd);
992 } else {
993 ScmBignum *bv = SCM_BIGNUM(Scm_MakeBignumFromSI(dd));
994 ScmBignum *br;
995 q = make_bignum(dividend->size + 1);
996 br = bignum_gdiv(dividend, bv, q);
997 rr = br->values[0];
998 }
999 if (remainder) {
1000 *remainder = ((dividend->sign < 0)? -(signed long)rr : (signed long)rr);
1001 }
1002 q->sign = dividend->sign * d_sign;
1003 return Scm_NormalizeBignum(q);
1004 }
1005
1006 /* If we only need rem(bignum, fixnum), we don't need to copy bignum
1007 to keep quotient. */
Scm_BignumRemSI(const ScmBignum * dividend,long divisor)1008 long Scm_BignumRemSI(const ScmBignum *dividend, long divisor)
1009 {
1010 #if (SIZEOF_LONG == 4) && HAVE_UINT64_T
1011 /* ILP32 with 64bit integer - easy one. */
1012 uint64_t dd = (divisor < 0)? -divisor : divisor;
1013 int sign = SCM_BIGNUM_SIGN(dividend);
1014 int k = SCM_BIGNUM_SIZE(dividend) - 1;
1015 uint64_t m = 0;
1016 for (;k >= 0; k--) {
1017 uint64_t x = (uint64_t)dividend->values[k] + (m << WORD_BITS);
1018 m = x % dd;
1019 }
1020 return (long)m * sign;
1021 #else /* SIZEOF_LONG >= 4 || !HAVE_UINT64_T*/
1022 /* Here we need double-machine-word `rem` machine-word. */
1023 u_long dd = (divisor < 0)? -divisor : divisor;
1024 int sign = SCM_BIGNUM_SIGN(dividend);
1025 int k = SCM_BIGNUM_SIZE(dividend) - 1;
1026 int shift = WORD_BITS - Scm__HighestBitNumber(dd) - 1;
1027 u_long m = 0;
1028 for (;k >= 0; k--) {
1029 u_long x = dividend->values[k];
1030 int total_shift = 0;
1031 /* we calculate [m;x] modulo dd. since always m < dd it's safe to
1032 shift m first. */
1033 while (total_shift < WORD_BITS) {
1034 int s = ((total_shift + shift < WORD_BITS)?
1035 shift : WORD_BITS-total_shift);
1036 m = (m << s) | (x >> (WORD_BITS - s));
1037 x <<= s;
1038 m %= dd;
1039 total_shift += s;
1040 }
1041 }
1042 return (long)m * sign;
1043 #endif
1044 }
1045
1046 /* assuming dividend and divisor is normalized. returns quotient and
1047 remainder */
Scm_BignumDivRem(const ScmBignum * dividend,const ScmBignum * divisor)1048 ScmObj Scm_BignumDivRem(const ScmBignum *dividend, const ScmBignum *divisor)
1049 {
1050 /* special case */
1051 if (Scm_BignumAbsCmp(dividend, divisor) < 0) {
1052 return Scm_Cons(SCM_MAKE_INT(0), SCM_OBJ(dividend));
1053 }
1054
1055 ScmBignum *q = make_bignum(dividend->size - divisor->size + 1);
1056 ScmBignum *r = bignum_gdiv(dividend, divisor, q);
1057 q->sign = dividend->sign * divisor->sign;
1058 r->sign = dividend->sign;
1059
1060 return Scm_Cons(Scm_NormalizeBignum(q), Scm_NormalizeBignum(r));
1061 }
1062
1063 /*-----------------------------------------------------------------------
1064 * Logical (bitwise) operations
1065 */
1066
Scm_BignumAsh(const ScmBignum * x,int cnt)1067 ScmObj Scm_BignumAsh(const ScmBignum *x, int cnt)
1068 {
1069 if (cnt == 0) return SCM_OBJ(x);
1070 if (cnt > 0) {
1071 int rsize = SCM_BIGNUM_SIZE(x) + (cnt+WORD_BITS-1)/WORD_BITS;
1072 ScmBignum *r = make_bignum(rsize);
1073 return Scm_NormalizeBignum(bignum_lshift(r, x, cnt));
1074 } else {
1075 int rsize = SCM_BIGNUM_SIZE(x) + cnt/WORD_BITS;
1076 if (rsize < 1) {
1077 if (SCM_BIGNUM_SIGN(x) < 0) {
1078 return SCM_MAKE_INT(-1);
1079 } else {
1080 return SCM_MAKE_INT(0);
1081 }
1082 } else {
1083 if (SCM_BIGNUM_SIGN(x) < 0) {
1084 /* painful way */
1085 ScmObj r = Scm_Quotient(Scm_Add(SCM_OBJ(x), SCM_MAKE_INT(1)),
1086 Scm_Ash(SCM_MAKE_INT(1), -cnt),
1087 NULL);
1088 return Scm_Add(r, SCM_MAKE_INT(-1));
1089 } else {
1090 ScmBignum *r = make_bignum(rsize);
1091 return Scm_NormalizeBignum(bignum_rshift(r, x, -cnt));
1092 }
1093 }
1094 }
1095 }
1096
1097 /* internal routine for logand. z = x & y. assumes z has enough size.
1098 * assumes x and y are in 2's complement form (sign is ignored).
1099 */
bignum_and(ScmBignum * z,const ScmBignum * x,const ScmBignum * y,int commsize,int xsize,int ysize)1100 static ScmBignum *bignum_and(ScmBignum *z,
1101 const ScmBignum *x, const ScmBignum *y,
1102 int commsize, int xsize, int ysize)
1103 {
1104 int i;
1105 for (i = 0; i < commsize; i++) {
1106 z->values[i] = x->values[i] & y->values[i];
1107 }
1108 if (i < xsize) {
1109 for (; i < xsize; i++) z->values[i] = x->values[i];
1110 } else if (i < ysize) {
1111 for (; i < ysize; i++) z->values[i] = y->values[i];
1112 }
1113 return z;
1114 }
1115
Scm_BignumLogAnd(const ScmBignum * x,const ScmBignum * y)1116 ScmObj Scm_BignumLogAnd(const ScmBignum *x, const ScmBignum *y)
1117 {
1118 int xsize = SCM_BIGNUM_SIZE(x), xsign = SCM_BIGNUM_SIGN(x);
1119 int ysize = SCM_BIGNUM_SIZE(y), ysign = SCM_BIGNUM_SIGN(y);
1120 int zsize, minsize = min(xsize, ysize);
1121
1122 if (xsign > 0) {
1123 if (ysign > 0) {
1124 ScmBignum *z = bignum_and(make_bignum(minsize), x, y, minsize,
1125 0, 0);
1126 return Scm_NormalizeBignum(z);
1127 } else {
1128 ScmBignum *yy = SCM_BIGNUM(Scm_BignumComplement(y));
1129 ScmBignum *z = bignum_and(make_bignum(xsize), x, yy, minsize,
1130 xsize, 0);
1131 return Scm_NormalizeBignum(z);
1132 }
1133 } else {
1134 if (ysign > 0) {
1135 ScmBignum *xx = SCM_BIGNUM(Scm_BignumComplement(x));
1136 ScmBignum *z = bignum_and(make_bignum(ysize), xx, y, minsize,
1137 0, ysize);
1138 return Scm_NormalizeBignum(z);
1139 } else {
1140 ScmBignum *xx = SCM_BIGNUM(Scm_BignumComplement(x));
1141 ScmBignum *yy = SCM_BIGNUM(Scm_BignumComplement(y));
1142 zsize = max(xsize, ysize);
1143 ScmBignum *z = bignum_and(make_bignum(zsize), xx, yy, minsize,
1144 xsize, ysize);
1145 SCM_BIGNUM_SIGN(z) = -1;
1146 bignum_2scmpl(z);
1147 return Scm_NormalizeBignum(z);
1148 }
1149 }
1150 }
1151
1152 /* internal routine for logior. z = x | y. assumes z has enough size.
1153 * assumes x and y are in 2's complement form (sign is ignored).
1154 */
bignum_ior(ScmBignum * z,const ScmBignum * x,const ScmBignum * y,int commsize,int xsize,int ysize)1155 static ScmBignum *bignum_ior(ScmBignum *z,
1156 const ScmBignum *x, const ScmBignum *y,
1157 int commsize, int xsize, int ysize)
1158 {
1159 int i;
1160 for (i = 0; i < commsize; i++) {
1161 z->values[i] = x->values[i] | y->values[i];
1162 }
1163 if (i < xsize) {
1164 for (; i < xsize; i++) z->values[i] = x->values[i];
1165 } else if (i < ysize) {
1166 for (; i < ysize; i++) z->values[i] = y->values[i];
1167 }
1168 return z;
1169 }
1170
Scm_BignumLogIor(const ScmBignum * x,const ScmBignum * y)1171 ScmObj Scm_BignumLogIor(const ScmBignum *x, const ScmBignum *y)
1172 {
1173 int xsize = SCM_BIGNUM_SIZE(x), xsign = SCM_BIGNUM_SIGN(x);
1174 int ysize = SCM_BIGNUM_SIZE(y), ysign = SCM_BIGNUM_SIGN(y);
1175 int minsize = min(xsize, ysize);
1176
1177 if (xsign >= 0) {
1178 if (ysign >= 0) {
1179 int zsize = max(xsize, ysize);
1180 ScmBignum *z = bignum_ior(make_bignum(zsize), x, y, minsize,
1181 xsize, ysize);
1182 return Scm_NormalizeBignum(z);
1183 } else {
1184 ScmBignum *yy = SCM_BIGNUM(Scm_BignumComplement(y));
1185 ScmBignum *z = bignum_ior(make_bignum(ysize), x, yy, minsize,
1186 0, ysize);
1187 SCM_BIGNUM_SIGN(z) = -1;
1188 bignum_2scmpl(z);
1189 return Scm_NormalizeBignum(z);
1190 }
1191 } else {
1192 if (ysign >= 0) {
1193 ScmBignum *xx = SCM_BIGNUM(Scm_BignumComplement(x));
1194 ScmBignum *z = bignum_ior(make_bignum(xsize), xx, y, minsize,
1195 xsize, 0);
1196 SCM_BIGNUM_SIGN(z) = -1;
1197 bignum_2scmpl(z);
1198 return Scm_NormalizeBignum(z);
1199 } else {
1200 ScmBignum *xx = SCM_BIGNUM(Scm_BignumComplement(x));
1201 ScmBignum *yy = SCM_BIGNUM(Scm_BignumComplement(y));
1202 ScmBignum *z = bignum_ior(make_bignum(minsize), xx, yy, minsize,
1203 0, 0);
1204 SCM_BIGNUM_SIGN(z) = -1;
1205 bignum_2scmpl(z);
1206 return Scm_NormalizeBignum(z);
1207 }
1208 }
1209 }
1210
Scm_BignumLogXor(const ScmBignum * x,const ScmBignum * y)1211 ScmObj Scm_BignumLogXor(const ScmBignum *x, const ScmBignum *y)
1212 {
1213 /* TODO: more efficient implementation */
1214 ScmObj xandy = Scm_BignumLogAnd(x, y);
1215 ScmObj xory = Scm_BignumLogIor(x, y);
1216 return Scm_LogAnd(xory, Scm_LogNot(xandy));
1217 }
1218
Scm_BignumLogCount(const ScmBignum * b)1219 int Scm_BignumLogCount(const ScmBignum *b)
1220 {
1221 const ScmBignum *z = (SCM_BIGNUM_SIGN(b)>0)? b : SCM_BIGNUM(Scm_BignumComplement(b));
1222 int size = SCM_BIGNUM_SIZE(z) * SCM_WORD_BITS;
1223
1224 ScmBits *bits = (ScmBits*)z->values;
1225 if (b->sign > 0) {
1226 return Scm_BitsCount1(bits, 0, size);
1227 } else {
1228 return Scm_BitsCount0(bits, 0, size);
1229 }
1230 }
1231
1232
1233 /*-----------------------------------------------------------------------
1234 * Printing
1235 */
1236
Scm_BignumToString(const ScmBignum * b,int radix,int use_upper)1237 ScmObj Scm_BignumToString(const ScmBignum *b, int radix, int use_upper)
1238 {
1239 static const char ltab[] = "0123456789abcdefghijklmnopqrstuvwxyz";
1240 static const char utab[] = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
1241 const char *tab = use_upper? utab : ltab;
1242 ScmObj h = SCM_NIL, t = SCM_NIL;
1243 if (radix < SCM_RADIX_MIN || radix > SCM_RADIX_MAX)
1244 Scm_Error("radix out of range: %d", radix);
1245 ScmBignum *q = SCM_BIGNUM(Scm_BignumCopy(b));
1246 for (;q->size > 0;) {
1247 long rem = bignum_sdiv(q, radix);
1248 SCM_ASSERT(rem >= 0 && rem < radix);
1249 SCM_APPEND1(h, t, SCM_MAKE_CHAR(tab[rem]));
1250 for (; q->values[q->size-1] == 0 && q->size > 0; q->size--)
1251 ;
1252 }
1253 if (q->sign < 0) SCM_APPEND1(h, t, SCM_MAKE_CHAR('-'));
1254 return Scm_ListToString(Scm_ReverseX(h));
1255 }
1256
Scm_DumpBignum(const ScmBignum * b,ScmPort * out)1257 int Scm_DumpBignum(const ScmBignum *b, ScmPort *out)
1258 {
1259 Scm_Printf(out, "#<bignum ");
1260 if (b->sign < 0) SCM_PUTC('-', out);
1261 for (int i=(int)b->size-1; i>=0; i--) {
1262 Scm_Printf(out, "%08lx ", b->values[i]);
1263 }
1264 SCM_PUTC('>', out);
1265 return 0;
1266 }
1267
1268 /*-----------------------------------------------------------------------
1269 * Denormalized bignum API
1270 * These are provided for optimization of specific cases.
1271 */
1272
1273 /* Returns a bignum of specified size, initializing the least significant
1274 word by init. */
Scm_MakeBignumWithSize(int size,u_long init)1275 ScmBignum *Scm_MakeBignumWithSize(int size, u_long init)
1276 {
1277 ScmBignum *b = make_bignum(size);
1278 b->values[0] = init;
1279 return b;
1280 }
1281
1282 /* Calculate acc * coef + c and store the result to acc, if the result fits
1283 in acc. If acc's size is not enough, allocate new bignum, which is at
1284 least sizeincr words bigger than acc.
1285 Returns the bignum that has the result, without normalizing.
1286 Acc need not be normalized. */
Scm_BignumAccMultAddUI(ScmBignum * acc,u_long coef,u_long c)1287 ScmBignum *Scm_BignumAccMultAddUI(ScmBignum *acc, u_long coef, u_long c)
1288 {
1289 ScmBignum *r;
1290 u_int rsize = acc->size + 1;
1291 ALLOC_TEMP_BIGNUM(r, rsize);
1292 r->values[0] = c;
1293 bignum_mul_word(r, acc, coef, 0);
1294 if (r->values[rsize-1] == 0) {
1295 for (u_int i=0; i<acc->size; i++) {
1296 acc->values[i] = r->values[i];
1297 }
1298 return acc;
1299 } else {
1300 ScmBignum *rr;
1301 rr = make_bignum(rsize + 3); /* 3 is arbitrary size increment */
1302 rr->sign = acc->sign;
1303 for (u_int i=0; i<rsize; i++) {
1304 rr->values[i] = r->values[i];
1305 }
1306 return rr;
1307 }
1308 }
1309