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  return(newv);
368  }
369 
371 pointer x;
372 { pointer y;
373  register int size, i;
374  register eusinteger_t *xv, *yv;
375 
376  size=bigsize(x);
377  y=makebig(size);
378  xv=bigvec(x);
379  yv=bigvec(y);
380  for (i=0; i<size; i++) yv[i]=xv[i];
381 #ifdef SAFETY
382  take_care(y);
383 #endif
384  return(y);
385  }
386 
387 /* extend_big(b,newsize) allocates a new bignum which has the same value
388  as 'b' in the 'newsize' word. Sign is copied in the upper extra words.
389  Newsize should be greater than the original size of b.
390 */
391 pointer extend_big(pointer b, int newsize)
392 { pointer nb;
393  eusinteger_t *bv, *nbv;
394  int i,bsize;
395  nb=makebig(newsize);
396  bv=bigvec(b); nbv=bigvec(nb); bsize=bigsize(b);
397  for (i=0; i<bsize; i++) nbv[i]=bv[i];
398  if (big_sign(b)<0) {
399  for (i=bsize; i<newsize; i++) nbv[i]= MASK;
400  nbv[newsize-1] = -1;}
401  return(nb); }
402 
403 /*
404  Big_zerop(x) answers if bignum x is zero or not.
405  X may be any bignum.
406 */
408 pointer x;
409 { register eusinteger_t *xv;
410  register int i, size;
411  xv=bigvec(x); size=bigsize(x);
412  for (i=0; i<size; i++) if (xv[i] != 0) return(0);
413  return(1);
414 }
415 
416 /*
417  Big_sign(x) returns
418  something < 0 if x < -1
419  something >= 0 if x > 0.
420  X may be any bignum.
421 */
423 pointer x;
424 { pointer y;
425  y=x->c.bgnm.bv;
426  return(y->c.ivec.iv[vecsize(y)-1]);
427  }
428 
429 /* Big_compare(x, y) returns
430  -1 if x < y
431  0 if x = y
432  1 if x > y.
433  X and y may be any bignum.
434 */
435 
436 int big_compare(x, y)
437 register pointer x, y;
438 { register eusinteger_t *xv, *yv;
439  register int xsize,ysize,i,j=0;
440  eusinteger_t xsign, ysign;
441  xv=bigvec(x); yv=bigvec(y);
442  xsize=bigsize(x); ysize=bigsize(y);
443 
444  xsign=xv[xsize-1]; ysign=yv[ysize-1];
445  if (xsign<0 && ysign>=0) return(-1);
446  else if (xsign>=0 && ysign<0) return(1);
447  if (xsign<0) { /* ysign is also negative*/
448  /*assume both x and y are normalized*/
449  if (xsize>ysize) return(-1);
450  else if (xsize<ysize) return(1);
451  else { /* xsize=ysize*/
452  if (xv[xsize-1] < yv[xsize-1]) return(-1);
453  else if (xv[xsize-1]>yv[xsize-1]) return(1);
454  for (i=xsize-2; i>=0; i--) {
455  if (xv[i] < yv[i]) return(1);
456  else if (xv[i] > yv[i]) return(-1); }
457  return(0);} }
458  else { /* both x and y are positive */
459  if (xsize>ysize) return(1);
460  else if (xsize<ysize) return(-1);
461  else { /* xsize=ysize*/
462  for (i=xsize-1; i>=0; i--) {
463  if (xv[i] > yv[i]) return(1);
464  else if (xv[i] < yv[i]) return(-1); }
465  return(0);} }
466 }
467 
468 
469 /*
470  Complement_big(x) destructively takes
471  the 2's complement of bignum x.
472  X may be any bignum.
473 */
475 pointer x;
476 {
477  int size, i=0;
478  eusinteger_t *xv;
479  size=bigsize(x); xv=bigvec(x);
480 
481  while (i != size-1) {
482  if (xv[i] != 0) {
483  xv[i] = (-xv[i]) & MASK;
484  goto ONE;}
485  i++;}
486  if (xv[i] == ~MASK) {
487  xv[i] = 0;
488  stretch_big(x, 1);
489  xv=bigvec(x); /*reload*/ }
490  else xv[i] = -(xv[i]);
491  return;
492 
493 ONE:
494  for (;;) {
495  i++;
496  if (i==size-1) break;
497  xv[i] = (~xv[i]) & MASK;}
498  xv[i] = ~xv[i];
499  return;
500 }
501 
502 /*
503  big_minus(x) returns the 2's complement of bignum x.
504  X may be any bignum.
505 */
507 pointer x;
508 { int size, i;
509  eusinteger_t *xv, *yv;
510  pointer y;
511 
512  size=bigsize(x); xv=bigvec(x);
513  y = makebig(size);
514  yv=bigvec(y);
515  i=0;
516  while (i<size-1) {
517  if (xv[i] != 0) {
518  yv[i]= (-xv[i]) & MASK; goto ONE;}
519  i++;}
520  if (xv[i] == MSB) { /* 0x80000000 */
521  yv[i]= 0;
522  stretch_big(y, 1);
523  yv=bigvec(y);
524  }
525  else yv[i] = -(xv[i]);
526 #ifdef SAFETY
527  take_care(y);
528 #endif
529  return(y);
530 
531 ONE:
532  for (;;) {
533  i++;
534  if (i==size-1) break;
535  yv[i] = (~xv[i]) & MASK;
536  }
537  yv[i] = ~xv[i];
538 #ifdef SAFETY
539  take_care(y);
540 #endif
541  return(y);
542 }
543 
544 
545 /*
546  Add_int_big(c, x) destructively adds non-negative int c
547  to bignum x.
548  I should be non-negative.
549  X may be any bignum.
550 */
551 void add_int_big(c, x)
552 eusinteger_t c;
553 pointer x;
554 { register int size, i=0;
555  register eusinteger_t *xv;
556 
557  size=bigsize(x); xv=bigvec(x);
558  while (i<size-1) {
559  xv[i] += c;
560  if (xv[i] < 0) { /* carry */
561  c = 1;
562  xv[i] &= MASK;
563  i++; }
564  else return;} /* no carry propagation to upper bits */
565 
566  if (xv[i] >= 0) { /* MSB */
567  xv[i] += c;
568  if (xv[i] < 0) { /* overflow */
569  xv[i] &= MASK;
570  stretch_big(x, 1);
571  xv=bigvec(x);} }
572  else xv[i] += c;
573  return; }
574 
575 /*
576  Sub_int_big(c, x) destructively subtracts non-negative int c
577  from bignum x.
578  c should be non-negative.
579  X may be any bignum.
580 */
581 void sub_int_big(c, x)
582 eusinteger_t c;
583 pointer x;
584 { register int size, i=0;
585  register eusinteger_t *xv;
586  size=bigsize(x); xv=bigvec(x);
587  while (i<size-1) {
588  xv[i] -= c;
589  if (xv[i] < 0) { /* borrow */
590  c = 1; xv[i] &= MASK; i++;}
591  else return;}
592  if (xv[i] < 0) { /* MSB */
593  xv[i] -= c;
594  if (xv[i] >= 0) { /* if sign changes, underflow */
595  xv[i] &= MASK;
596  stretch_big(x, -2); } }
597  else xv[i] -= c;
598  return;}
599 
600 /*
601  Mul_int_big(i, x) destructively multiplies non-negative bignum x
602  by non-negative int i.
603  I should be non-negative.
604  X should be non-negative.
605 */
606 void mul_int_big(c, x)
607 eusinteger_t c;
608 pointer x;
609 { int size, i=0;
610  eusinteger_t *xv, h=0;
611 
612  size=bigsize(x); xv=bigvec(x);
613  while (i<size) {
614  extended_mul(c, xv[i], h, &h, &xv[i]);
615  i++;}
616  if (h > 0) stretch_big(x, h);
617  return;}
618 
619 /*
620  Div_int_big(c, x) destructively divides non-negative bignum x
621  by positive int c.
622  X will hold the remainder of the division.
623  Div_int_big(c, x) returns the remainder of the division.
624  c should be positive.
625  X should be non-negative.
626 */
628 eusinteger_t c;
629 pointer x;
630 { int i, size;
631  eusinteger_t *xv, r;
632  if (c == 0) error(E_USER,(pointer)"divide by zero in bignum div");
633  size=bigsize(x); xv=bigvec(x);
634  /* divide from MSB */
635  r = xv[size-1] % c;
636  xv[size-1] /= c;
637  i=size-2;
638  while (i>=0) {
639  extended_div(c, r, xv[i], &(xv[i]), &r);
640  i--;}
641  return(r); }
642 
643 /*
644  Big_plus(x, y) returns the sum of bignum x and bignum y.
645  X and y may be any bignum.
646 */
648 pointer x, y;
649 { pointer z;
650  int i, xsize, ysize;
651  eusinteger_t c, zlast, *xv, *yv, *zv;
652 
653  xsize=bigsize(x); ysize=bigsize(y);
654  if (xsize<ysize) { /*exchange x and y so that x is bigger than y */
655  z=x; x=y; y=z;
656  i=xsize; xsize=ysize; ysize=i;}
657 
658  /*fprintf(stderr, "big_plus xsize=%d ysize=%d\n", xsize,ysize);*/
659 
660  xv=bigvec(x); yv=bigvec(y);
661  z=copy_big(x); zv=bigvec(z);
662  c=0;
663  for (i=0; i<ysize; i++) {
664  zv[i]=zlast=zv[i]+yv[i]+c;
665  if (zlast<0) { c=1; zv[i] &= MASK;}
666  else c=0;}
667  i= ysize - 1;
668  if (ysize==xsize) {
669  if (xv[i]>=0 && yv[i]>=0 && zlast<0) stretch_big(z, 1);
670  else if (xv[i]<0 && yv[i]<0 && zlast>=0) {
671  stretch_big(z, -2);}
672  else if (zlast<0) zv[i] |= MSB;
673  return(z);}
674  else { /* xsize>ysize */
675  if (yv[ysize-1]>=0) c=1; else c= -1;
676  while (i<xsize-1) {
677  if (zlast<0) zv[i] &= MASK;
678  else return(z);
679  i++;
680  zv[i] += c;
681  zlast=zv[i];}
682  if (c>=0 && xv[i]>=0 && zv[i]<0) {
683  zv[i] &= MASK; stretch_big(z, 1);}
684  else if (c<0 && xv[i]<0 && zv[i]>=0) {
685  stretch_big(z, -2);}
686  return(z); }
687  }
688 
689 /*
690  Big_times(x0, y0) returns the product
691  of non-negative bignum x0 and non-negative bignum y0.
692  X0 and y0 should be non-negative.
693 */
695 pointer x, y;
696 { int xsize, ysize, zsize;
697  eusinteger_t *xv, *yv, *zv;
698  int i, j, k, m;
699  pointer z;
700  eusinteger_t hi, lo;
701 
702  xsize=bigsize(x); ysize=bigsize(y);
703  zsize=xsize+ysize;
704  z=makebig(zsize);
705  xv=bigvec(x); yv=bigvec(y); zv=bigvec(z);
706 
707  zv[0]=0;
708 
709  for (i=0; i<xsize; i++) {
710  m=xv[i];
711  hi=0;
712  for (j=0; j<ysize; j++) {
713  k=i+j;
714  extended_mul(m, yv[j], hi, &hi, &lo);
715  zv[k]+=lo;
716  if (zv[k] & MSB) { zv[k] &= MASK; hi++; }
717  if (j==ysize-1) {
718  while (hi>0) {
719  k++;
720  if (k>=xsize+ysize) error(E_USER,(pointer)"bignum mult overflow");
721  zv[k] += hi;
722  if (zv[k] & MSB) { zv[k] &= MASK; hi=1; }
723  else hi=0;}
724  }
725  }
726  }
727 #ifdef SAFETY
728  take_care(z);
729 #endif
730  return(z);}
731 
732 /*
733  Sub_int_big_big(c, x, y) destructively subtracts c*x from y.
734  C should be a non-negative int.
735  X should be a normalized non-negative bignum.
736  Y should be non-negative bignum and should be such that
737  y <= c*x.
738 */
739 
740 void sub_int_big_big(c, x, y)
741 eusinteger_t c;
742 pointer x, y;
743 { int i, j;
744  int xsize, ysize;
745  eusinteger_t *xv, *yv;
746  eusinteger_t hi, lo;
747 
748  xsize=bigsize(x); ysize=bigsize(y);
749  xv=bigvec(x); yv=bigvec(y);
750 
751  hi = 0;
752  for (i=0;i<xsize;i++) {
753  extended_mul(c, xv[i], hi, &hi, &lo);
754  yv[i] -= lo;
755  if (yv[i] & MSB) {
756  yv[i] &= MASK;
757  hi++;
758  }
759  if (i==xsize-1) {
760  while (hi > 0) {
761  i++;
762  yv[i] -= hi;
763  if (yv[i] & MSB) {
764  yv[i] &= MASK;
765  hi = 1;}
766  else break;}
767  break;}
768  }
769  y=normalize_big(y);
770 }
771 
772 /*
773  Get_standardizing_factor_and_normalize(x)
774  returns the standardizing factor of x.
775  As a side-effect, x will be normalized.
776  X should be a positive bignum.
777  If x is multiplied by the standardizing factor,
778  the most significant digit of x will be not less than 2**30,
779  i.e. the most significant bit of that digit will be set.
780 */
781 
783 pointer x;
784 { int size, s, newsize;
785  eusinteger_t *xv, i, j;
786 
787  size=bigsize(x); xv=bigvec(x);
788  s=size-1;
789  newsize=size;
790  while (xv[newsize-1]==0 && newsize>=0) newsize--;
791  x->c.bgnm.bv->c.ivec.length=makeint(newsize);
792  if (newsize==0) return(-1);
793  for (j = 1, i = xv[newsize-1]; (i <<= 1) >= 0; j <<= 1) ;
794  return(j);
795  }
796 
797 /*
798  Div_big_big(x0, y0) divides y0 by x0
799  and destructively places the remainder in y0.
800  X0 should be a normalized positive bignum,
801  whose most significant digit is not less than 2**30.
802  Y0 should be a non-negative bignum,
803  whose length must be equal to the length of x0
804  or one bigger than that.
805  Div_big_big(x0, y0) returns the quotient of the division,
806  which must be less than 2**31.
807 */
809 pointer x, y;
810 { int xsize, ysize, zsize;
811  eusinteger_t *xv, *yv, *zv;
812  eusinteger_t q,r; /*quotient, remainder*/
813 
814 
815  xsize=bigsize(x); ysize=bigsize(y);
816  xv=bigvec(x); yv=bigvec(y);
817 
818  if (xsize != ysize) {
819  if (yv[xsize] >= xv[xsize-1])
820  q = MASK-1;
821  /* This is the most critical point. */
822  /* The long division will fail, */
823  /* if the quotient can't lie in 31 bits. */
824  else {
825  extended_div(xv[xsize-1], yv[xsize], yv[xsize-1], &q, &r);
826  q -= 2;
827  /* You have to prove that 2 is enough, */
828  /* by doing elementary arithmetic, */
829  /* or refer to Knuth's dictionary. */
830  } }
831  else
832  q = yv[xsize-1] / xv[xsize-1] - 2;
833  if (q <= 0) q = 0;
834  else sub_int_big_big(q, x, y);
835  while (big_compare(x, y) <= 0) {
836  q++;
837  sub_int_big_big(1, x, y); }
838  return(q); }
839 
841 pointer x, y;
842 int i;
843 { pointer q, qq;
844  int xsize, ysize;
845  eusinteger_t *xv, *yv;
846 
847  xsize=bigsize(x); ysize=bigsize(y);
848 
849  if (i < 0) return(NULL);
850  if (i == 0) { /* x and y are of the same size */
851  i = div_big_big(y, x);
852  if (i == 0) return(NULL);
853  else {
854  qq = makebig(1);
855  bigvec(qq)[0] = i;
856 #ifdef SAFETY
857  take_care(qq);
858 #endif
859  return(qq);
860  }
861  }
862 
863  qq=makebig(i);
864 
865  /********************????????????????????????**************/
866 #ifdef SAFETY
867  take_care(qq);
868 #endif
869 
870  return(qq);}
871 
872 
873 
874 
875 
876 /*
877  Big_quotient_remainder(x0, y0, qp, rp)
878  sets the quotient and the remainder of the division of x0 by y0
879  to *qp and *rp respectively.
880  qp and rp should be place holders on the value stack to protect
881  from GC.
882  X0 should be a positive bignum.
883  Y0 should be a non-negative bignum.
884 */
885 void big_quotient_remainder(x0, y0, qp, rp)
886 pointer x0, y0, *qp, *rp;
887 {
888  context *ctx=euscontexts[thr_self()];
889  pointer x, y;
890  eusinteger_t i;
891  int l, m;
892 
893  x = copy_big(x0);
894  vpush(x);
895  y = y0;
897  mul_int_big(i, x);
898  mul_int_big(i, y);
899  l = bigsize(x);
900  m = bigsize(y);
901  *qp = big_quotient_remainder_auxiliary(x, y, l - m);
902  if (*qp == NULL) {
903  *qp = makebig1(0);
904  }
905  div_int_big(i, x);
906  div_int_big(i, y);
907  vpop();
908  *rp = x;
909 }
910 
911 
912 
914 pointer x;
915 {
916  return(bigsize(x));
917 }
918 
920 { int size, i, newsize, positive;
921  eusinteger_t *xv;
922 
923  size=bigsize(x); xv=bigvec(x);
924  newsize=size;
925  positive = gez(xv[size-1]);
926  for (i=size-1; i>=0; i--) {
927  if (positive && xv[i]==0) newsize--;
928  else if (!positive && xv[i]==-1 && size!=1) newsize--;
929  else break;}
930  if (newsize != size) {
931  x->c.bgnm.bv->c.ivec.length=makeint(newsize);
932  x->c.bgnm.size=makeint(newsize);
933  if (!positive) xv[newsize-1] |= MSB;}
934  return(x);
935  }
936 
938 pointer x;
939 { eusinteger_t y, msb3;
940  int bs;
941 
942  if (!isbignum(x)) return(x);
943  normalize_big(x);
944 
945  bs=bigsize(x);
946  if (bs==0) return(makeint(0));
947  else if (bs==1) {
948  y=x->c.bgnm.bv->c.ivec.iv[0];
949  msb3= (y >> (WORD_SIZE-3)) & 0x7;
950  if (msb3==0 || msb3==7) return(makeint(y));}
951  return(x);}
952 
954 { int size, i;
955  eusinteger_t *xv;
956  double d, e;
957  size=bigsize(x); xv=bigvec(x);
958 
959  for (i=0, d = 0.0, e = 1.0; i<size-1; i++) {
960  d += e * (double)(xv[i]);
961  e *= (double)((unsigned long)1 << (WORD_SIZE-1));
962  }
963  d += e * (double)(xv[size-1]);
964  return(d);
965 }
966 
968 { int i, j, exp, bsize, iz;
969  pointer b;
970  eusinteger_t *bv;
971  double sign, z;
972  extern double frexp(double, int *);
973  extern double ldexp(double, int);
974 
975  if (f<0.0) {
976  f= -f;
977  if (f<=MAXPOSFIXNUM+1) {
978  i=f;
979  return(makeint(-i)); }
980  else sign= -1;}
981  else if (f<=MAXPOSFIXNUM) { /* 0x1fffffff on 32bit machines */
982  i=f;
983  return(makeint(i));}
984  else sign= 1;
985 
986  /* here, f>MAXPOSFIXNUM */
987  z=frexp(f, &exp); /* f=z * 2^exp; 0.5<z<=1.0; exp>=30 */
988  /* z has at most 22bit or 52bit precision */
989  z= ldexp(z, 24);
990  bsize=(exp + WORD_SIZE - 2)/(WORD_SIZE - 1);
991  b=makebig(bsize);
992  exp -= 24;
993  iz=z;
994  bv=b->c.bgnm.bv->c.ivec.iv;
995  i=exp / (WORD_SIZE -1);
996  j=exp % (WORD_SIZE -1);
997  bv[i]= (iz << j) & MASK;
998  if (j+24 > 31) bv[i+1]= (iz >> (WORD_SIZE -1 - j)) & MASK;
999  if (sign<0) complement_big(b);
1000  return(b);
1001  }
1002 
1003 
context * euscontexts[MAXTHREAD]
Definition: eus.c:105
eusinteger_t iv[1]
Definition: eus.h:303
d
eusinteger_t div_int_big(eusinteger_t c, pointer x)
Definition: big.c:627
pointer makebig()
static char * rcsid
Definition: big.c:18
#define makeint(v)
Definition: sfttest.c:2
Definition: eus.h:522
eusinteger_t get_standardizing_factor_and_normalize(pointer x)
Definition: big.c:782
pointer length
Definition: eus.h:302
struct bignum bgnm
Definition: eus.h:422
pointer normalize_big(pointer x)
Definition: big.c:919
pointer bv
Definition: eus.h:376
pointer big_minus(pointer x)
Definition: big.c:506
Definition: eus.h:1002
pointer makevector(pointer, int)
Definition: makes.c:417
pointer copy_big(pointer x)
Definition: big.c:370
pointer makebig1()
struct intvector ivec
Definition: eus.h:414
union cell::cellunion c
void mul_int_big(eusinteger_t c, pointer x)
Definition: big.c:606
pointer extend_big(pointer b, int newsize)
Definition: big.c:391
long l
Definition: structsize.c:3
eusinteger_t div_big_big(pointer x, pointer y)
Definition: big.c:808
Definition: eus.h:379
#define ONE
Definition: snrm2.c:9
eusinteger_t big_zerop(pointer x)
Definition: big.c:407
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:375
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:740
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:551
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:647
pointer big_quotient_remainder_auxiliary(pointer x, pointer y, int i)
Definition: big.c:840
int big_compare(pointer x, pointer y)
Definition: big.c:436
pointer normalize_bignum(pointer x)
Definition: big.c:937
void sub_int_big(eusinteger_t c, pointer x)
Definition: big.c:581
#define NULL
Definition: transargv.c:8
pointer C_INTVECTOR
Definition: eus.c:146
void complement_big(pointer x)
Definition: big.c:474
unsigned int thr_self()
Definition: eus.c:25
eusinteger_t big_sign(pointer x)
Definition: big.c:422
pointer big_times(pointer x, pointer y)
Definition: big.c:694
double eusfloat_t
Definition: eus.h:20
int big_length(pointer x)
Definition: big.c:913
void big_quotient_remainder(pointer x0, pointer y0, pointer *qp, pointer *rp)
Definition: big.c:885
eusfloat_t big_to_float(pointer x)
Definition: big.c:953
pointer eusfloat_to_big(float f)
Definition: big.c:967


euslisp
Author(s): Toshihiro Matsui
autogenerated on Thu Jun 6 2019 20:00:43