big.c
Go to the documentation of this file.
1 /*
2  big.c
3 
4  bignum routines for EusLisp
5  1995 (c) Toshihiro Matsui, Electrotechnical Laboraotory
6  Skeltal codes are borrowed from KCL's big.c and earith.c.
7 
8  EusLisp's bignum is represented by an integer-vector.
9  A value is placed in a vector from LSB to MSB.
10  Sign bit is not used in the first to (MSB -1)th words.
11  Sign is held in the sign bit of the MSB word.
12  example: #x1122334455667788
13  bv[0]= #x55667788
14  bv[1]= #x22446688
15 
16 */
17 
18 static char *rcsid="@(#)$Id$";
19 #include "eus.h"
20 
21 #define ltz(x) ((x)<0)
22 #define gez(x) ((x)>=0)
23 
24 extern pointer makebig();
25 
26 #if (WORD_SIZE == 64)
27 
28 void extended_mul(d, q, r, hp, lp)
29 eusinteger_t d, q, r;
30 eusinteger_t *hp, *lp;
31 {
32  eusinteger_t ld, hd, lq, hq, m, lm, hm, lr;
33  eusinteger_t hz, lz;
34 
35 
36  ld = (d & 0x00000000ffffffff);
37  hd = (d >> 32);
38  hd = (hd & 0x000000007fffffff);
39  lq = (q & 0x00000000ffffffff);
40  hq = (q >> 32);
41  hq = (hq & 0x000000007fffffff);
42  lr = r;
43 
44  m = ld * lq;
45  hm = (m & 0x8000000000000000);
46  hm = (hm >> 32);
47  hm = (hm & 0x0000000080000000);
48  m = (m & 0x7fffffffffffffff);
49 
50  m = m + lr;
51 
52  lz = (m & 0x00000000ffffffff);
53 
54  m = (m >> 32);
55  m = (m & 0x00000000ffffffff);
56  m = m + hm;
57  m = m + (lq * hd);
58  hm = (m & 0x8000000000000000);
59  hm = (hm >> 31);
60  hm = (hm & 0x00000001ffffffff);
61  m = (m & 0x7fffffffffffffff);
62 
63  m = m + (hq * ld);
64 
65  lm = (m & 0x000000007fffffff);
66 
67  lz = lz + (lm << 32);
68 
69  m = (m >> 31);
70  m = (m & 0x00000001ffffffff);
71  m = m + hm;
72 
73  hz = m + (hq * hd * 2);
74 
75 
76  *hp = hz;
77  *lp = lz;
78 
79 }
80 
81 void extended_div(d, h, l, qp, rp)
82 eusinteger_t d, h, l;
83 eusinteger_t *qp, *rp;
84 {
85  eusinteger_t lh, ld, ll, ans;
86  int i;
87 
88  ld = d;
89  lh = h;
90  ll = l;
91 
92 
93  ans = 0;
94  for(i=63;i>0;i--) {
95  ans = (ans << 1);
96  lh = (lh << 1);
97  ll = (ll << 1);
98  if (ll < 0) lh += 1;
99  if (lh >= d) {
100  lh -= ld;
101  ans += 1;
102  }
103  if (lh < 0) {
104  lh = (0x7fffffffffffffff - ld) +
105  (lh & 0x7fffffffffffffff) + 1;
106  ans += 1;
107  }
108  }
109 
110 
111  *qp = ans;
112  *rp = lh;
113 }
114 
115 #else /* WORD_SIZE==32 */
116 #if (!Solaris2 || !sun4)
117 
118 void extended_mul(d, q, r, hp, lp)
119 eusinteger_t d, q, r;
120 eusinteger_t *hp, *lp;/* ???? */
121 {
122  long int ld, hd, lq, hq, m, lm, hm, lr;
123  int hz, lz;
124 
125 
126  ld = (d & 0x0000ffff);
127  hd = (d >> 16);
128  hd = (hd & 0x00007fff);
129  lq = (q & 0x0000ffff);
130  hq = (q >> 16);
131  hq = (hq & 0x00007fff);
132  lr = r;
133 
134  m = ld * lq;
135  hm = (m & 0x80000000);
136  hm = (hm >> 16);
137  hm = (hm & 0x00008000);
138  m = (m & 0x7fffffff);
139 
140  m = m + lr;
141 
142  lz = (m & 0x0000ffff);
143 
144  m = (m >> 16);
145  m = (m & 0x0000ffff);
146  m = m + hm;
147  m = m + (lq * hd);
148  hm = (m & 0x80000000);
149  hm = (hm >> 15);
150  hm = (hm & 0x0001ffff);
151  m = (m & 0x7fffffff);
152 
153  m = m + (hq * ld);
154 
155  lm = (m & 0x00007fff);
156 
157  lz = lz + (lm << 16);
158 
159  m = (m >> 15);
160  m = (m & 0x0001ffff);
161  m = m + hm;
162 
163  hz = m + (hq * hd * 2);
164 
165 
166  *hp = hz;
167  *lp = lz;
168 
169 }
170 
171 void extended_div(d, h, l, qp, rp)
172 eusinteger_t d, h, l;
173 eusinteger_t *qp, *rp;
174 {
175  long int lh, ld, ll;
176  int i,ans;
177 
178  ld = d;
179  lh = h;
180  ll = l;
181 
182 
183  ans = 0;
184  for(i=31;i>0;i--) {
185  ans = (ans << 1);
186  lh = (lh << 1);
187  ll = (ll << 1);
188  if (ll < 0) lh += 1;
189  if (lh >= d) {
190  lh -= ld;
191  ans += 1;
192  }
193  if (lh < 0) {
194  lh = (0x7fffffff - ld) + (lh & 0x7fffffff) + 1;
195  ans += 1;
196  }
197  }
198 
199 
200  *qp = ans;
201  *rp = lh;
202 }
203 
204 #else /* !(sun4 && Solaris2) */
205 
206 /* Compile with -O option */
207 #define u_int unsigned int
208 
209 #define DEBUG
210 #undef DEBUG
211 #ifdef DEBUG
212 extended_mul(d,q,r,hp,lp)
213  u_int d,q,r,*hp,*lp;
214  {
215  printf("ext.mul (%x %x %x)",d,q,r);fflush(stdout);
216  low_extended_mul(d,q,r,hp,lp);
217  printf(" =%x:%08x\n",*hp,*lp);fflush(stdout);
218  }
219 static low_extended_mul(d, q, r, hp, lp)
220 #else
221 extended_mul(d,q,r,hp,lp)
222 #endif
223 register u_int d, q, r;
224 register u_int *hp, *lp;
225 {
226  asm(" save %sp, -0x60, %sp");
227  asm(" mov %i0,%o0"); /* d*q+r */
228  asm(" call .umul,2"); /* delayed inst. */
229  asm(" mov %i1,%o1");
230  asm(" addcc %o0,%i2,%o0");
231  asm(" addx %o1, 0,%o1");
232 
233  asm(" sll %o1, 2,%o1"); /* ((H<<1) & 0x3fffffff) */
234  asm(" srl %o1, 1,%o1");
235  asm(" srl %o0, 31,%i2"); /* | (L>>31) -> %o1 */
236  asm(" or %o1,%i2,%o1");
237 
238  asm(" sll %o0, 1,%o0"); /* L & 0x7ffffffff */
239  asm(" srl %o0, 1,%o0");
240 
241  asm(" st %o0,[%i4]");
242  asm(" st %o1,[%i3]");
243  asm(" ret");
244  asm(" restore");
245 }
246 
247 extended_div(d, h, l, qp, rp)
248 register u_int d, h, l;
249 u_int *qp, *rp;
250 {
251  register u_int q,one;
252 #ifdef DEBUG
253 printf("ext.div (%x %x:%x)",d,h,l);fflush(stdout);
254 #endif
255  l=l<<1; one=1<<31; q=0;
256  if (h>=d) { q+= one;h-=d; }
257  h<<=1; h|=l>>31; l<<=1; one>>=1;
258  if (h>=d) { q+= one;h-=d; }
259  h<<=1; h|=l>>31; l<<=1; one>>=1;
260  if (h>=d) { q+= one;h-=d; }
261  h<<=1; h|=l>>31; l<<=1; one>>=1;
262  if (h>=d) { q+= one;h-=d; }
263  h<<=1; h|=l>>31; l<<=1; one>>=1;
264  if (h>=d) { q+= one;h-=d; }
265  h<<=1; h|=l>>31; l<<=1; one>>=1;
266  if (h>=d) { q+= one;h-=d; }
267  h<<=1; h|=l>>31; l<<=1; one>>=1;
268  if (h>=d) { q+= one;h-=d; }
269  h<<=1; h|=l>>31; l<<=1; one>>=1;
270  if (h>=d) { q+= one;h-=d; }
271  h<<=1; h|=l>>31; l<<=1; one>>=1;
272 
273  if (h>=d) { q+= one;h-=d; }
274  h<<=1; h|=l>>31; l<<=1; one>>=1;
275  if (h>=d) { q+= one;h-=d; }
276  h<<=1; h|=l>>31; l<<=1; one>>=1;
277  if (h>=d) { q+= one;h-=d; }
278  h<<=1; h|=l>>31; l<<=1; one>>=1;
279  if (h>=d) { q+= one;h-=d; }
280  h<<=1; h|=l>>31; l<<=1; one>>=1;
281  if (h>=d) { q+= one;h-=d; }
282  h<<=1; h|=l>>31; l<<=1; one>>=1;
283  if (h>=d) { q+= one;h-=d; }
284  h<<=1; h|=l>>31; l<<=1; one>>=1;
285  if (h>=d) { q+= one;h-=d; }
286  h<<=1; h|=l>>31; l<<=1; one>>=1;
287  if (h>=d) { q+= one;h-=d; }
288  h<<=1; h|=l>>31; l<<=1; one>>=1;
289 
290  if (h>=d) { q+= one;h-=d; }
291  h<<=1; h|=l>>31; l<<=1; one>>=1;
292  if (h>=d) { q+= one;h-=d; }
293  h<<=1; h|=l>>31; l<<=1; one>>=1;
294  if (h>=d) { q+= one;h-=d; }
295  h<<=1; h|=l>>31; l<<=1; one>>=1;
296  if (h>=d) { q+= one;h-=d; }
297  h<<=1; h|=l>>31; l<<=1; one>>=1;
298  if (h>=d) { q+= one;h-=d; }
299  h<<=1; h|=l>>31; l<<=1; one>>=1;
300  if (h>=d) { q+= one;h-=d; }
301  h<<=1; h|=l>>31; l<<=1; one>>=1;
302  if (h>=d) { q+= one;h-=d; }
303  h<<=1; h|=l>>31; l<<=1; one>>=1;
304  if (h>=d) { q+= one;h-=d; }
305  h<<=1; h|=l>>31; l<<=1; one>>=1;
306 
307  if (h>=d) { q+= one;h-=d; }
308  h<<=1; h|=l>>31; l<<=1; one>>=1;
309  if (h>=d) { q+= one;h-=d; }
310  h<<=1; h|=l>>31; l<<=1; one>>=1;
311  if (h>=d) { q+= one;h-=d; }
312  h<<=1; h|=l>>31; l<<=1; one>>=1;
313  if (h>=d) { q+= one;h-=d; }
314  h<<=1; h|=l>>31; l<<=1; one>>=1;
315  if (h>=d) { q+= one;h-=d; }
316  h<<=1; h|=l>>31; l<<=1; one>>=1;
317  if (h>=d) { q+= one;h-=d; }
318  h<<=1; h|=l>>31; l<<=1; one>>=1;
319  if (h>=d) { q+= one;h-=d; }
320  h<<=1; h|=l>>31; l<<=1; one>>=1;
321  if (h>=d) { q+= one;h-=d; }
322  h<<=1; h|=l>>31; l<<=1; one>>=1;
323 
324  *qp=q;
325  *rp=h>>1;
326 
327 /* l=l<<1; one=1<<31;
328  for(q=0;one;one>>=1)
329  {
330  if (h>=d) { q+= one;h-=d; }
331  h<<=1;
332  h|=l>>31;
333  l<<=1;
334  }
335  *qp=q;
336  *rp=h>>1;
337 */
338 #ifdef DEBUG
339 printf(" =%x:%x\n",*qp,*rp);fflush(stdout);
340 #endif
341 }
342 
343 #endif /* !Solaris2 */
344 #endif
345 
346 
347 
348 /****************************************************************/
349 /* stretch_big(x, i)
350 /* 'x' is a bignum and 'i' is integer.
351 /* Allocate bigvector whose size is one word bigger than x,
352 /* and put 'i' in the MSB word.
353 */
354 
355 inline pointer stretch_big(x, i)
356 pointer x;
357 eusinteger_t i;
358 { pointer bn=x->c.bgnm.bv;
359  int vlen=vecsize(bn);
360  pointer newv, oldv;
361  register int j;
362 
363  newv=makevector(C_INTVECTOR, vlen+1);
364  for (j=0; j<vlen; j++) newv->c.ivec.iv[j]=bn->c.ivec.iv[j];
365  newv->c.ivec.iv[vlen]=i;
366  pointer_update(x->c.bgnm.bv, newv);
367  x->c.bgnm.size=makeint(vlen+1);
368  return(newv);
369  }
370 
372 pointer x;
373 { pointer y;
374  register int size, i;
375  register eusinteger_t *xv, *yv;
376 
377  size=bigsize(x);
378  y=makebig(size);
379  xv=bigvec(x);
380  yv=bigvec(y);
381  for (i=0; i<size; i++) yv[i]=xv[i];
382 #ifdef SAFETY
383  take_care(y);
384 #endif
385  return(y);
386  }
387 
388 /* extend_big(b,newsize) allocates a new bignum which has the same value
389  as 'b' in the 'newsize' word. Sign is copied in the upper extra words.
390  Newsize should be greater than the original size of b.
391 */
392 pointer extend_big(pointer b, int newsize)
393 { pointer nb;
394  eusinteger_t *bv, *nbv;
395  int i,bsize;
396  nb=makebig(newsize);
397  bv=bigvec(b); nbv=bigvec(nb); bsize=bigsize(b);
398  for (i=0; i<bsize; i++) nbv[i]=bv[i];
399  if (big_sign(b)<0) {
400  for (i=bsize; i<newsize; i++) nbv[i]= MASK;
401  nbv[newsize-1] = -1;}
402  return(nb); }
403 
404 /*
405  Big_zerop(x) answers if bignum x is zero or not.
406  X may be any bignum.
407 */
409 pointer x;
410 { register eusinteger_t *xv;
411  register int i, size;
412  xv=bigvec(x); size=bigsize(x);
413  for (i=0; i<size; i++) if (xv[i] != 0) return(0);
414  return(1);
415 }
416 
417 /*
418  Big_sign(x) returns
419  something < 0 if x < -1
420  something >= 0 if x > 0.
421  X may be any bignum.
422 */
424 pointer x;
425 { pointer y;
426  y=x->c.bgnm.bv;
427  return(y->c.ivec.iv[vecsize(y)-1]);
428  }
429 
430 /* Big_compare(x, y) returns
431  -1 if x < y
432  0 if x = y
433  1 if x > y.
434  X and y may be any bignum.
435 */
436 
437 int big_compare(x, y)
438 register pointer x, y;
439 { register eusinteger_t *xv, *yv;
440  register int xsize,ysize,i,j=0;
441  eusinteger_t xsign, ysign;
442  xv=bigvec(x); yv=bigvec(y);
443  xsize=bigsize(x); ysize=bigsize(y);
444 
445  xsign=xv[xsize-1]; ysign=yv[ysize-1];
446  if (xsign<0 && ysign>=0) return(-1);
447  else if (xsign>=0 && ysign<0) return(1);
448  if (xsign<0) { /* ysign is also negative*/
449  /*assume both x and y are normalized*/
450  if (xsize>ysize) return(-1);
451  else if (xsize<ysize) return(1);
452  else { /* xsize=ysize*/
453  if (xv[xsize-1] < yv[xsize-1]) return(-1);
454  else if (xv[xsize-1]>yv[xsize-1]) return(1);
455  for (i=xsize-2; i>=0; i--) {
456  if (xv[i] < yv[i]) return(1);
457  else if (xv[i] > yv[i]) return(-1); }
458  return(0);} }
459  else { /* both x and y are positive */
460  if (xsize>ysize) return(1);
461  else if (xsize<ysize) return(-1);
462  else { /* xsize=ysize*/
463  for (i=xsize-1; i>=0; i--) {
464  if (xv[i] > yv[i]) return(1);
465  else if (xv[i] < yv[i]) return(-1); }
466  return(0);} }
467 }
468 
469 
470 /*
471  Complement_big(x) destructively takes
472  the 2's complement of bignum x.
473  X may be any bignum.
474 */
476 pointer x;
477 {
478  int size, i=0;
479  eusinteger_t *xv;
480  size=bigsize(x); xv=bigvec(x);
481 
482  while (i != size-1) {
483  if (xv[i] != 0) {
484  xv[i] = (-xv[i]) & MASK;
485  goto ONE;}
486  i++;}
487  if (xv[i] == ~MASK) {
488  xv[i] = 0;
489  stretch_big(x, 1);
490  xv=bigvec(x); /*reload*/ }
491  else xv[i] = -(xv[i]);
492  return;
493 
494 ONE:
495  for (;;) {
496  i++;
497  if (i==size-1) break;
498  xv[i] = (~xv[i]) & MASK;}
499  xv[i] = ~xv[i];
500  return;
501 }
502 
503 /*
504  big_minus(x) returns the 2's complement of bignum x.
505  X may be any bignum.
506 */
508 pointer x;
509 { int size, i;
510  eusinteger_t *xv, *yv;
511  pointer y;
512 
513  size=bigsize(x); xv=bigvec(x);
514  y = makebig(size);
515  yv=bigvec(y);
516  i=0;
517  while (i<size-1) {
518  if (xv[i] != 0) {
519  yv[i]= (-xv[i]) & MASK; goto ONE;}
520  i++;}
521  if ((unsigned long)xv[i] == MSB) { /* 0x80000000 */
522  yv[i]= 0;
523  stretch_big(y, 1);
524  yv=bigvec(y);
525  }
526  else yv[i] = -(xv[i]);
527 #ifdef SAFETY
528  take_care(y);
529 #endif
530  return(y);
531 
532 ONE:
533  for (;;) {
534  i++;
535  if (i==size-1) break;
536  yv[i] = (~xv[i]) & MASK;
537  }
538  yv[i] = ~xv[i];
539 #ifdef SAFETY
540  take_care(y);
541 #endif
542  return(y);
543 }
544 
545 
546 /*
547  Add_int_big(c, x) destructively adds non-negative int c
548  to bignum x.
549  I should be non-negative.
550  X may be any bignum.
551 */
552 void add_int_big(c, x)
553 eusinteger_t c;
554 pointer x;
555 { register int size, i=0;
556  register eusinteger_t *xv;
557 
558  size=bigsize(x); xv=bigvec(x);
559  while (i<size-1) {
560  xv[i] += c;
561  if (xv[i] < 0) { /* carry */
562  c = 1;
563  xv[i] &= MASK;
564  i++; }
565  else return;} /* no carry propagation to upper bits */
566 
567  if (xv[i] >= 0) { /* MSB */
568  xv[i] += c;
569  if (xv[i] < 0) { /* overflow */
570  xv[i] &= MASK;
571  stretch_big(x, 1);
572  xv=bigvec(x);} }
573  else xv[i] += c;
574  return; }
575 
576 /*
577  Sub_int_big(c, x) destructively subtracts non-negative int c
578  from bignum x.
579  c should be non-negative.
580  X may be any bignum.
581 */
582 void sub_int_big(c, x)
583 eusinteger_t c;
584 pointer x;
585 { register int size, i=0;
586  register eusinteger_t *xv;
587  size=bigsize(x); xv=bigvec(x);
588  while (i<size-1) {
589  xv[i] -= c;
590  if (xv[i] < 0) { /* borrow */
591  c = 1; xv[i] &= MASK; i++;}
592  else return;}
593  if (xv[i] < 0) { /* MSB */
594  xv[i] -= c;
595  if (xv[i] >= 0) { /* if sign changes, underflow */
596  xv[i] &= MASK;
597  stretch_big(x, -2); } }
598  else xv[i] -= c;
599  return;}
600 
601 /*
602  Mul_int_big(i, x) destructively multiplies non-negative bignum x
603  by non-negative int i.
604  I should be non-negative.
605  X should be non-negative.
606 */
607 void mul_int_big(c, x)
608 eusinteger_t c;
609 pointer x;
610 { int size, i=0;
611  eusinteger_t *xv, h=0;
612 
613  size=bigsize(x); xv=bigvec(x);
614  while (i<size) {
615  extended_mul(c, xv[i], h, &h, &xv[i]);
616  i++;}
617  if (h > 0) stretch_big(x, h);
618  return;}
619 
620 /*
621  Div_int_big(c, x) destructively divides non-negative bignum x
622  by positive int c.
623  X will hold the remainder of the division.
624  Div_int_big(c, x) returns the remainder of the division.
625  c should be positive.
626  X should be non-negative.
627 */
629 eusinteger_t c;
630 pointer x;
631 { int i, size;
632  eusinteger_t *xv, r;
633  if (c == 0) error(E_USER,(pointer)"divide by zero in bignum div");
634  size=bigsize(x); xv=bigvec(x);
635  /* divide from MSB */
636  r = xv[size-1] % c;
637  xv[size-1] /= c;
638  i=size-2;
639  while (i>=0) {
640  extended_div(c, r, xv[i], &(xv[i]), &r);
641  i--;}
642  return(r); }
643 
644 /*
645  Big_plus(x, y) returns the sum of bignum x and bignum y.
646  X and y may be any bignum.
647 */
649 pointer x, y;
650 { pointer z;
651  int i, xsize, ysize;
652  eusinteger_t c, zlast, *xv, *yv, *zv;
653 
654  xsize=bigsize(x); ysize=bigsize(y);
655  if (xsize<ysize) { /*exchange x and y so that x is bigger than y */
656  z=x; x=y; y=z;
657  i=xsize; xsize=ysize; ysize=i;}
658 
659  /*fprintf(stderr, "big_plus xsize=%d ysize=%d\n", xsize,ysize);*/
660 
661  xv=bigvec(x); yv=bigvec(y);
662  z=copy_big(x); zv=bigvec(z);
663  c=0;
664  for (i=0; i<ysize; i++) {
665  zv[i]=zlast=zv[i]+yv[i]+c;
666  if (zlast<0) { c=1; zv[i] &= MASK;}
667  else c=0;}
668  i= ysize - 1;
669  if (ysize==xsize) {
670  if (xv[i]>=0 && yv[i]>=0 && zlast<0) stretch_big(z, 1);
671  else if (xv[i]<0 && yv[i]<0 && zlast>=0) {
672  stretch_big(z, -2);}
673  else if (zlast<0) zv[i] |= MSB;
674  return(z);}
675  else { /* xsize>ysize */
676  if (yv[ysize-1]>=0) c=1; else c= -1;
677  while (i<xsize-1) {
678  if (zlast<0) zv[i] &= MASK;
679  else return(z);
680  i++;
681  zv[i] += c;
682  zlast=zv[i];}
683  if (c>=0 && xv[i]>=0 && zv[i]<0) {
684  zv[i] &= MASK; stretch_big(z, 1);}
685  else if (c<0 && xv[i]<0 && zv[i]>=0) {
686  stretch_big(z, -2);}
687  return(z); }
688  }
689 
690 /*
691  Big_times(x0, y0) returns the product
692  of non-negative bignum x0 and non-negative bignum y0.
693  X0 and y0 should be non-negative.
694 */
696 pointer x, y;
697 { int xsize, ysize, zsize;
698  eusinteger_t *xv, *yv, *zv;
699  int i, j, k, m;
700  pointer z;
701  eusinteger_t hi, lo;
702 
703  xsize=bigsize(x); ysize=bigsize(y);
704  zsize=xsize+ysize;
705  z=makebig(zsize);
706  xv=bigvec(x); yv=bigvec(y); zv=bigvec(z);
707 
708  zv[0]=0;
709 
710  for (i=0; i<xsize; i++) {
711  m=xv[i];
712  hi=0;
713  for (j=0; j<ysize; j++) {
714  k=i+j;
715  extended_mul(m, yv[j], hi, &hi, &lo);
716  zv[k]+=lo;
717  if (zv[k] & MSB) { zv[k] &= MASK; hi++; }
718  if (j==ysize-1) {
719  while (hi>0) {
720  k++;
721  if (k>=xsize+ysize) error(E_USER,(pointer)"bignum mult overflow");
722  zv[k] += hi;
723  if (zv[k] & MSB) { zv[k] &= MASK; hi=1; }
724  else hi=0;}
725  }
726  }
727  }
728 #ifdef SAFETY
729  take_care(z);
730 #endif
731  return(z);}
732 
733 /*
734  Sub_int_big_big(c, x, y) destructively subtracts c*x from y.
735  C should be a non-negative int.
736  X should be a normalized non-negative bignum.
737  Y should be non-negative bignum and should be such that
738  y <= c*x.
739 */
740 
741 void sub_int_big_big(c, x, y)
742 eusinteger_t c;
743 pointer x, y;
744 { int i, j;
745  int xsize, ysize;
746  eusinteger_t *xv, *yv;
747  eusinteger_t hi, lo;
748 
749  xsize=bigsize(x); ysize=bigsize(y);
750  xv=bigvec(x); yv=bigvec(y);
751 
752  hi = 0;
753  for (i=0;i<xsize;i++) {
754  extended_mul(c, xv[i], hi, &hi, &lo);
755  yv[i] -= lo;
756  if (yv[i] & MSB) {
757  yv[i] &= MASK;
758  hi++;
759  }
760  if (i==xsize-1) {
761  while (hi > 0) {
762  i++;
763  yv[i] -= hi;
764  if (yv[i] & MSB) {
765  yv[i] &= MASK;
766  hi = 1;}
767  else break;}
768  break;}
769  }
770  y=normalize_big(y);
771 }
772 
773 /*
774  Get_standardizing_factor_and_normalize(x)
775  returns the standardizing factor of x.
776  As a side-effect, x will be normalized.
777  X should be a positive bignum.
778  If x is multiplied by the standardizing factor,
779  the most significant digit of x will be not less than 2**30,
780  i.e. the most significant bit of that digit will be set.
781 */
782 
784 pointer x;
785 { int size, s, newsize;
786  eusinteger_t *xv, i, j;
787 
788  size=bigsize(x); xv=bigvec(x);
789  s=size-1;
790  newsize=size;
791  while (xv[newsize-1]==0 && newsize>=0) newsize--;
792  x->c.bgnm.bv->c.ivec.length=makeint(newsize);
793  if (newsize==0) return(-1);
794  for (j = 1, i = xv[newsize-1]; (i <<= 1) >= 0; j <<= 1) ;
795  return(j);
796  }
797 
798 /*
799  Div_big_big(x0, y0) divides y0 by x0
800  and destructively places the remainder in y0.
801  X0 should be a normalized positive bignum,
802  whose most significant digit is not less than 2**30.
803  Y0 should be a non-negative bignum,
804  whose length must be equal to the length of x0
805  or one bigger than that.
806  Div_big_big(x0, y0) returns the quotient of the division,
807  which must be less than 2**31.
808 */
810 pointer x, y;
811 { int xsize, ysize, zsize;
812  eusinteger_t *xv, *yv, *zv;
813  eusinteger_t q,r; /*quotient, remainder*/
814 
815 
816  xsize=bigsize(x); ysize=bigsize(y);
817  xv=bigvec(x); yv=bigvec(y);
818 
819  if (xsize != ysize) {
820  if (yv[xsize] >= xv[xsize-1])
821  q = MASK-1;
822  /* This is the most critical point. */
823  /* The long division will fail, */
824  /* if the quotient can't lie in 31 bits. */
825  else {
826  extended_div(xv[xsize-1], yv[xsize], yv[xsize-1], &q, &r);
827  q -= 2;
828  /* You have to prove that 2 is enough, */
829  /* by doing elementary arithmetic, */
830  /* or refer to Knuth's dictionary. */
831  } }
832  else
833  q = yv[xsize-1] / xv[xsize-1] - 2;
834  if (q <= 0) q = 0;
835  else sub_int_big_big(q, x, y);
836  while (big_compare(x, y) <= 0) {
837  q++;
838  sub_int_big_big(1, x, y); }
839  return(q); }
840 
842 pointer x, y;
843 int i;
844 { pointer q, qq;
845  int xsize, ysize;
846  eusinteger_t *xv, *yv;
847 
848  xsize=bigsize(x); ysize=bigsize(y);
849 
850  if (i < 0) return(NULL);
851  if (i == 0) { /* x and y are of the same size */
852  i = div_big_big(y, x);
853  if (i == 0) return(NULL);
854  else {
855  qq = makebig(1);
856  bigvec(qq)[0] = i;
857 #ifdef SAFETY
858  take_care(qq);
859 #endif
860  return(qq);
861  }
862  }
863 
864  qq=makebig(i);
865 
866  /********************????????????????????????**************/
867 #ifdef SAFETY
868  take_care(qq);
869 #endif
870 
871  return(qq);}
872 
873 
874 
875 
876 
877 /*
878  Big_quotient_remainder(x0, y0, qp, rp)
879  sets the quotient and the remainder of the division of x0 by y0
880  to *qp and *rp respectively.
881  qp and rp should be place holders on the value stack to protect
882  from GC.
883  X0 should be a positive bignum.
884  Y0 should be a non-negative bignum.
885 */
886 void big_quotient_remainder(x0, y0, qp, rp)
887 pointer x0, y0, *qp, *rp;
888 {
889  context *ctx=euscontexts[thr_self()];
890  pointer x, y;
891  eusinteger_t i;
892  int l, m;
893 
894  x = copy_big(x0);
895  vpush(x);
896  y = y0;
898  mul_int_big(i, x);
899  mul_int_big(i, y);
900  l = bigsize(x);
901  m = bigsize(y);
902  *qp = big_quotient_remainder_auxiliary(x, y, l - m);
903  if (*qp == NULL) {
904  *qp = makebig1(0);
905  }
906  div_int_big(i, x);
907  div_int_big(i, y);
908  vpop();
909  *rp = x;
910 }
911 
912 
913 
915 pointer x;
916 {
917  return(bigsize(x));
918 }
919 
921 { int size, i, newsize, positive;
922  eusinteger_t *xv;
923 
924  size=bigsize(x); xv=bigvec(x);
925  newsize=size;
926  positive = gez(xv[size-1]);
927  for (i=size-1; i>=0; i--) {
928  if (positive && xv[i]==0) newsize--;
929  else if (!positive && xv[i]==-1 && size!=1) newsize--;
930  else break;}
931  if (newsize != size) {
932  x->c.bgnm.bv->c.ivec.length=makeint(newsize);
933  x->c.bgnm.size=makeint(newsize);
934  if (!positive) xv[newsize-1] |= MSB;}
935  return(x);
936  }
937 
939 pointer x;
940 { eusinteger_t y, msb3;
941  int bs;
942 
943  if (!isbignum(x)) return(x);
944  normalize_big(x);
945 
946  bs=bigsize(x);
947  if (bs==0) return(makeint(0));
948  else if (bs==1) {
949  y=x->c.bgnm.bv->c.ivec.iv[0];
950  msb3= (y >> (WORD_SIZE-3)) & 0x7;
951  if (msb3==0 || msb3==7) return(makeint(y));}
952  return(x);}
953 
955 { int size, i;
956  eusinteger_t *xv;
957  double d, e;
958  size=bigsize(x); xv=bigvec(x);
959 
960  for (i=0, d = 0.0, e = 1.0; i<size-1; i++) {
961  d += e * (double)(xv[i]);
962  e *= (double)((unsigned long)1 << (WORD_SIZE-1));
963  }
964  d += e * (double)(xv[size-1]);
965  return(d);
966 }
967 
969 { int i, j, exp, bsize, iz;
970  pointer b;
971  eusinteger_t *bv;
972  double sign, z;
973  extern double frexp(double, int *);
974  extern double ldexp(double, int);
975 
976  if (f<0.0) {
977  f= -f;
978  if (f<=MAXPOSFIXNUM+1) {
979  i=f;
980  return(makeint(-i)); }
981  else sign= -1;}
982  else if (f<=MAXPOSFIXNUM) { /* 0x1fffffff on 32bit machines */
983  i=f;
984  return(makeint(i));}
985  else sign= 1;
986 
987  /* here, f>MAXPOSFIXNUM */
988  z=frexp(f, &exp); /* f=z * 2^exp; 0.5<z<=1.0; exp>=30 */
989  /* z has at most 22bit or 52bit precision */
990  z= ldexp(z, 24);
991  bsize=(exp + WORD_SIZE - 2)/(WORD_SIZE - 1);
992  b=makebig(bsize);
993  exp -= 24;
994  iz=z;
995  bv=b->c.bgnm.bv->c.ivec.iv;
996  i=exp / (WORD_SIZE -1);
997  j=exp % (WORD_SIZE -1);
998  bv[i]= (iz << j) & MASK;
999  if (j+24 > 31) bv[i+1]= (iz >> (WORD_SIZE -1 - j)) & MASK;
1000  if (sign<0) complement_big(b);
1001  return(b);
1002  }
1003 
1004 
context * euscontexts[MAXTHREAD]
Definition: eus.c:105
eusinteger_t iv[1]
Definition: eus.h:305
d
eusinteger_t div_int_big(eusinteger_t c, pointer x)
Definition: big.c:628
pointer makebig()
static char * rcsid
Definition: big.c:18
#define makeint(v)
Definition: sfttest.c:2
Definition: eus.h:524
eusinteger_t get_standardizing_factor_and_normalize(pointer x)
Definition: big.c:783
pointer length
Definition: eus.h:304
struct bignum bgnm
Definition: eus.h:424
pointer normalize_big(pointer x)
Definition: big.c:920
pointer bv
Definition: eus.h:378
pointer big_minus(pointer x)
Definition: big.c:507
Definition: eus.h:1006
pointer makevector(pointer, int)
Definition: makes.c:417
pointer copy_big(pointer x)
Definition: big.c:371
pointer makebig1()
struct intvector ivec
Definition: eus.h:416
union cell::cellunion c
void mul_int_big(eusinteger_t c, pointer x)
Definition: big.c:607
pointer extend_big(pointer b, int newsize)
Definition: big.c:392
long l
Definition: structsize.c:3
eusinteger_t div_big_big(pointer x, pointer y)
Definition: big.c:809
Definition: eus.h:381
#define ONE
Definition: snrm2.c:9
eusinteger_t big_zerop(pointer x)
Definition: big.c:408
void extended_mul(eusinteger_t d, eusinteger_t q, eusinteger_t r, eusinteger_t *hp, eusinteger_t *lp)
Definition: big.c:28
pointer size
Definition: eus.h:377
short s
Definition: structsize.c:2
pointer stretch_big(pointer x, eusinteger_t i)
Definition: big.c:355
void sub_int_big_big(eusinteger_t c, pointer x, pointer y)
Definition: big.c:741
static low_extended_mul(extended_mul(d, extended_mul(q, extended_mul(r, extended_mul(hp, u_int *lp)
Definition: big.c:219
void add_int_big(eusinteger_t c, pointer x)
Definition: big.c:552
pointer error(enum errorcode ec,...) pointer error(va_alist) va_dcl
Definition: eus.c:297
long eusinteger_t
Definition: eus.h:19
void extended_div(eusinteger_t d, eusinteger_t h, eusinteger_t l, eusinteger_t *qp, eusinteger_t *rp)
Definition: big.c:81
pointer big_plus(pointer x, pointer y)
Definition: big.c:648
pointer big_quotient_remainder_auxiliary(pointer x, pointer y, int i)
Definition: big.c:841
int big_compare(pointer x, pointer y)
Definition: big.c:437
pointer normalize_bignum(pointer x)
Definition: big.c:938
void sub_int_big(eusinteger_t c, pointer x)
Definition: big.c:582
#define NULL
Definition: transargv.c:8
pointer C_INTVECTOR
Definition: eus.c:146
void complement_big(pointer x)
Definition: big.c:475
unsigned int thr_self()
Definition: eus.c:25
eusinteger_t big_sign(pointer x)
Definition: big.c:423
pointer big_times(pointer x, pointer y)
Definition: big.c:695
double eusfloat_t
Definition: eus.h:21
int big_length(pointer x)
Definition: big.c:914
void big_quotient_remainder(pointer x0, pointer y0, pointer *qp, pointer *rp)
Definition: big.c:886
eusfloat_t big_to_float(pointer x)
Definition: big.c:954
pointer eusfloat_to_big(float f)
Definition: big.c:968


euslisp
Author(s): Toshihiro Matsui
autogenerated on Mon Feb 28 2022 22:18:27