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