2 /* --------------------------------------------------------------------------
3 * Yet another implementation of Integer
5 * Copyright (c) Glasgow University, 1999.
6 * All rights reserved. See NOTICE for details and conditions of use etc...
7 * ------------------------------------------------------------------------*/
14 #include "sainteger.h"
17 /* --------------------------------------------------------------------------
19 * ------------------------------------------------------------------------*/
21 typedef unsigned char uchar;
22 typedef unsigned short ush;
25 static int maxused_add ( B*, B* );
26 static int maxused_sub ( B*, B* );
27 static int maxused_mul ( B*, B* );
28 static int maxused_qrm ( B*, B* );
29 static int maxused_neg ( B* );
31 static int ucmp ( B*, B* );
32 static void uadd ( B*, B*, B* );
33 static void usub ( B*, B*, B* );
34 static void umul ( B*, B*, B* );
35 static void uqrm ( B*, B*, B*, B* );
37 /*#define DEBUG_SAINTEGER*/
38 /*#define DEBUG_SAINTEGER_UQRM*/
41 #ifdef DEBUG_SAINTEGER
42 #define myassert(zzzz) assert(zzzz)
44 #define myassert(zzzz) /* */
48 /* --------------------------------------------------------------------------
50 * ------------------------------------------------------------------------*/
55 printf ( "sign=%2d used=%d size=%d ", x->sign, x->used, x->size );
56 for (i = x->used-1; i >= 0; i--)
57 printf ( "%2x ", (int)(x->stuff[i]) );
62 static int sane ( B* x )
66 if (x->sign == 0 && x->used != 0) return 0;
67 if (x->sign != -1 && x->sign != 0 && x->sign != 1) return 0;
69 if (x->used < 0) return 0;
70 if (x->size < 0) return 0;
71 if (x->used > x->size) return 0;
72 if (x->used == 0) return 1;
73 if (x->stuff[x->used-1] == 0) return 0;
74 for (i = 0; i < x->used; i++)
75 if (x->stuff[i] >= B_BASE) return 0;
86 static void u_renormalise ( B* b )
88 while (b->used > 0 && b->stuff[b->used-1] == 0) b->used--;
89 if (b->used == 0) b->sign = 0; else b->sign = 1;
93 void do_renormalise ( B* b )
95 while (b->used > 0 && b->stuff[b->used-1] == 0) b->used--;
96 if (b->used == 0) b->sign = 0;
99 /* --------------------------------------------------------------------------
101 * ------------------------------------------------------------------------*/
103 static int maxused_add ( B* x, B* y )
107 return 1 + (x->used > y->used ? x->used : y->used);
110 static int maxused_sub ( B* x, B* y )
114 return 1 + (x->used > y->used ? x->used : y->used);
117 static int maxused_mul ( B* x, B* y )
121 return x->used + y->used;
124 static int maxused_qrm ( B* x, B* y )
128 return (x->used > y->used ? x->used : y->used);
131 static int maxused_neg ( B* x )
138 /* quick, safe approx */
139 static int maxused_fromInt ( int sizeof_int )
141 if (B_BASE == 256) return sizeof_int;
142 if (B_BASE >= 16) return 2 * sizeof_int;
143 if (B_BASE >= 4) return 4 * sizeof_int;
144 /* (B_BASE >= 2) */ return 8 * sizeof_int;
148 static int maxused_fromStr ( char* str )
151 if (*str == '-') str++;
152 while (isdigit((int)(*str))) { str++; nd++; };
154 if (B_BASE >= 100) return ((nd+1) / 2);
155 if (B_BASE >= 10) return nd;
156 /* (B_BASE >= 2)*/ return 4 * nd;
160 int size_add ( B* x, B* y )
162 return sizeof(B) + maxused_add(x,y);
165 int size_sub ( B* x, B* y )
167 return sizeof(B) + maxused_sub(x,y);
170 int size_mul ( B* x, B* y )
172 return sizeof(B) + maxused_mul(x,y);
175 int size_qrm ( B* x, B* y )
177 return sizeof(B) + maxused_qrm(x,y);
180 int size_neg ( B* x )
182 return sizeof(B) + maxused_neg(x);
185 int size_fromInt ( void )
187 int sizeof_int = sizeof(int);
188 return sizeof(B) + maxused_fromInt ( sizeof_int );
191 int size_fromWord ( void )
193 int sizeof_word = sizeof(unsigned int);
194 return sizeof(B) + maxused_fromInt ( sizeof_word );
197 int size_fromStr ( char* str )
199 return sizeof(B) + maxused_fromStr ( str );
202 int size_fltmantissa ( void )
204 return sizeof(B) + sizeof(float);
207 int size_dblmantissa ( void )
209 return sizeof(B) + sizeof(double);
213 /* --------------------------------------------------------------------------
215 * ------------------------------------------------------------------------*/
217 void do_fromInt ( int n, int sizeRes, B* res )
220 res->size = sizeRes - sizeof(B);
221 res->sign = res->used = 0;
222 if (n == 0) { myassert(sane(res)); return; };
223 if (n < 0) res->sign = -1; else res->sign = 1;
227 res->stuff[res->used] = (uchar)(n % B_BASE);
234 void do_fromWord ( unsigned int n, int sizeRes, B* res )
237 res->size = sizeRes - sizeof(B);
238 res->sign = res->used = 0;
239 if (n == 0) { myassert(sane(res)); return; };
243 res->stuff[res->used] = (uchar)(n % B_BASE);
250 /* NOTE: This only works currectly if B_BASE >= 10 */
251 void do_fromStr ( char* str, int sizeRes, B* res )
253 int sign, d, t, j, carry;
255 res->size = sizeRes - sizeof(B);
256 res->sign = res->used = 0;
258 if (*str == '-') { str++; sign = -1; };
260 while (isdigit((int)(*str))) {
262 /* multiply res by 10 */
264 for (j = 0; j < res->used; j++) {
265 t = 10 * res->stuff[j] + carry;
266 res->stuff[j] = t % B_BASE;
269 myassert(carry < B_BASE);
271 res->stuff[res->used++] = carry;
278 for (j = 0; j < res->used; j++) {
279 carry += res->stuff[j];
280 res->stuff[j] = carry % B_BASE;
282 if (carry == 0) break;
285 res->stuff[res->used++] = carry;
292 int do_toInt ( B* x )
295 if (x->sign == 0) return 0;
297 for (i = x->used-1; i >= 0; i--) {
299 res = res * B_BASE + d;
301 if (x->sign < 0) res = -res;
305 unsigned int do_toWord ( B* x )
309 if (x->sign == 0) return 0;
311 for (i = x->used-1; i >= 0; i--) {
313 res = res * B_BASE + d;
318 float do_toFloat ( B* x )
322 if (x->sign == 0) return 0.0;
324 for (i = x->used-1; i >= 0; i--) {
326 res = res * B_BASE_FLT + d;
328 if (x->sign < 0) res = -res;
332 double do_toDouble ( B* x )
336 if (x->sign == 0) return 0.0;
338 for (i = x->used-1; i >= 0; i--) {
340 res = res * B_BASE_FLT + d;
342 if (x->sign < 0) res = -res;
347 /* --------------------------------------------------------------------------
349 * ------------------------------------------------------------------------*/
351 /* A helper for signed + and -. sdiff(x,y) ignores the signs of x and y
352 sets p to the signed value abs(x)-abs(y).
354 static void sdiff ( B* x, B* y, B* res )
359 myassert(res->size == maxused_sub(x,y));
361 if (t == 0) { res->sign = res->used = 0; return; }
374 int do_getsign ( B* x )
380 void do_neg ( B* x, int sizeRes, B* res )
384 res->size = sizeRes - sizeof(B);
386 for (i = 0; i < x->used; i++)
387 res->stuff[i] = x->stuff[i];
388 res->sign = - (x->sign);
391 void do_add ( B* x, B* y, int sizeRes, B* res )
395 res->size = sizeRes - sizeof(B);
396 res->used = res->sign = 0;
398 if ( (x->sign >= 0 && y->sign >= 0) ||
399 (x->sign < 0 && y->sign < 0)) {
400 /* same sign; add magnitude and clone sign */
402 if (x->sign < 0 && res->sign != 0) res->sign = -1;
405 /* signs differ; employ sdiff */
406 if (x->sign >= 0 && y->sign < 0) {
409 myassert(x->sign < 0 && y->sign >= 0);
415 void do_sub ( B* x, B* y, int sizeRes, B* res )
419 res->size = sizeRes - sizeof(B);
420 res->used = res->sign = 0;
422 if ( (x->sign >= 0 && y->sign < 0) ||
423 (x->sign < 0 && y->sign >= 0)) {
424 /* opposite signs; add magnitudes and clone sign of x */
426 myassert(res->sign != 0);
427 if (x->sign < 0) res->sign = -1;
430 /* signs are the same; employ sdiff */
431 if (x->sign >= 0 && y->sign >= 0) {
434 myassert(x->sign < 0 && y->sign < 0);
441 void do_mul ( B* x, B* y, int sizeRes, B* res )
445 res->size = sizeRes - sizeof(B);
446 res->used = res->sign = 0;
448 if (x->sign == 0 || y->sign == 0) {
449 res->sign = res->used = 0;
454 if (x->sign != y->sign) res->sign = -1;
459 void do_qrm ( B* x, B* y, int sizeRes, B* q, B* r )
464 q->size = r->size = sizeRes - sizeof(B);
465 q->used = r->used = q->sign = r->sign = 0;
468 fprintf(stderr, "do_qrm: division by zero -- exiting now!\n");
474 q->used = r->used = q->sign = r->sign = 0;
475 myassert(sane(q)); myassert(sane(r));
480 if (x->sign != y->sign && q->sign != 0) q->sign = -1;
481 if (x->sign == -1 && r->sign != 0) r->sign = -1;
483 myassert(sane(q)); myassert(sane(r));
486 int do_cmp ( B* x, B* y )
492 if (x->sign < y->sign) return -1;
493 if (x->sign > y->sign) return 1;
494 myassert(x->sign == y->sign);
495 if (x->sign == 0) return 0;
496 if (x->sign == 1) return ucmp(x,y); else return ucmp(y,x);
500 /* --------------------------------------------------------------------------
502 * ------------------------------------------------------------------------*/
504 static int ucmp ( B* x, B* y )
509 if (x->used < y->used) return -1;
510 if (x->used > y->used) return 1;
511 for (i = x->used-1; i >= 0; i--) {
512 if (x->stuff[i] < y->stuff[i]) return -1;
513 if (x->stuff[i] > y->stuff[i]) return 1;
520 static void uadd ( B* x, B* y, B* res )
527 myassert (res->size == maxused_add(x,y));
528 res->used = res->size;
529 res->stuff[res->used-1] = 0;
531 if (x->used > y->used) {
540 for (i = 0; i < n; i++) {
541 t = x->stuff[i] + y->stuff[i] + c;
543 res->stuff[i] = t-B_BASE;
551 for (i = n; i < longer->used; i++) {
552 t = longer->stuff[i] + c;
554 res->stuff[i] = t-B_BASE;
561 myassert(res->used == longer->used+1);
562 res->stuff[longer->used] = c;
570 static void usub ( B* x, B* y, B* res )
575 myassert (x->used >= y->used);
576 myassert (res->size == maxused_sub(x,y));
579 for (i = 0; i < y->used; i++) {
580 t = x->stuff[i] - y->stuff[i] - b;
582 res->stuff[i] = t + B_BASE;
590 for (i = y->used; i < x->used; i++) {
593 res->stuff[i] = t + B_BASE;
607 void umul ( B* x, B* y, B* res )
613 myassert(res->size == maxused_mul(x,y));
615 for (j = 0; j < y->used; j++) res->stuff[j] = 0;
617 for (i = 0; i < x->used; i++) {
619 for (j = 0; j < y->used; j++) {
620 carry += res->stuff[i+j] + x->stuff[i]*y->stuff[j];
621 res->stuff[i+j] = carry % B_BASE;
623 myassert (carry < B_BASE);
625 res->stuff[i+y->used] = carry;
628 res->used = x->used+y->used;
634 static void uqrm ( B* dend, B* isor, B* dres, B* mres )
636 int i, j, t, vh, toolarge, delta, carry, scaleup;
637 uchar *dend_stuff, *isor_stuff, *tmp;
639 myassert(sane(isor));
640 myassert(sane(dend));
641 myassert(isor->used > 0); // against division by zero
643 myassert(dres->size == maxused_qrm(isor,dend));
644 myassert(mres->size == maxused_qrm(isor,dend));
646 if (dend->used < isor->used) {
647 // Result of division must be zero, since dividend has
648 // fewer digits than the divisor. Remainder is the
649 // original dividend.
651 mres->used = dend->used;
652 for (j = 0; j < mres->used; j++) mres->stuff[j] = dend->stuff[j];
653 u_renormalise(dres); u_renormalise(mres);
654 myassert(sane(dres));
655 myassert(sane(mres));
659 if (isor->used == 1) {
661 // Simple case; divisor is a single digit
663 for (j = dend->used-1; j >= 0; j--) {
664 carry += dend->stuff[j];
665 dres->stuff[j] = carry/isor->stuff[0];
666 carry = B_BASE*(carry%isor->stuff[0]);
669 dres->used = dend->used;
672 // Remainder is the final carry value
676 mres->stuff[0] = carry;
678 u_renormalise(dres); u_renormalise(mres);
679 myassert(sane(dres));
680 myassert(sane(mres));
685 // Complex case: both dividend and divisor have two or more digits.
686 myassert(isor->used >= 2);
687 myassert(dend->used >= 2);
689 // Allocate space for a copy of both dividend and divisor, since we
690 // need to mess with them. Also allocate tmp as a place to hold
691 // values of the form quotient_digit * divisor.
692 dend_stuff = malloc ( sizeof(uchar)*(dend->used+1) );
693 isor_stuff = malloc ( sizeof(uchar)*isor->used );
694 tmp = malloc ( sizeof(uchar)*(isor->used+1) );
695 myassert (dend_stuff && isor_stuff && tmp);
697 // Calculate a scaling-up factor, and multiply both divisor and
698 // dividend by it. Doing this reduces the number of corrections
699 // needed to the quotient-digit-estimates made in the loop below,
700 // and thus speeds up division, but is not actually needed to
701 // get the correct results. The scaleup factor should not increase
702 // the number of digits needed to represent either the divisor
703 // (since the factor is derived from it) or the dividend (since
704 // we already gave it a new leading zero).
705 scaleup = B_BASE / (1 + isor->stuff[isor->used-1]);
706 myassert (1 <= scaleup && scaleup <= B_BASE/2);
709 // Don't bother to multiply; just copy.
710 for (j = 0; j < dend->used; j++) dend_stuff[j] = dend->stuff[j];
711 for (j = 0; j < isor->used; j++) isor_stuff[j] = isor->stuff[j];
713 // Extend dividend with leading zero.
714 dend_stuff[dend->used] = tmp[isor->used] = 0;
718 for (j = 0; j < isor->used; j++) {
719 t = scaleup * isor->stuff[j] + carry;
720 isor_stuff[j] = t % B_BASE;
723 myassert (carry == 0);
726 for (j = 0; j < dend->used; j++) {
727 t = scaleup * dend->stuff[j] + carry;
728 dend_stuff[j] = t % B_BASE;
731 dend_stuff[dend->used] = carry;
735 // For each quotient digit ...
736 for (i = dend->used; i >= isor->used; i--) {
738 myassert (i <= dend->used);
739 myassert (isor->used >= 2);
741 #if DEBUG_SAINTEGER_UQRM
742 printf("\n---------\nqdigit %d\n", i );
743 printf("dend_stuff is ");
744 for (j = dend->used; j>= 0; j--) printf("%d ",dend_stuff[j]);
747 // Make a guess vh of the quotient digit
748 vh = (B_BASE*B_BASE*dend_stuff[i] + B_BASE*dend_stuff[i-1] + dend_stuff[i-2])
750 (B_BASE*isor_stuff[isor->used-1] + isor_stuff[isor->used-2]);
751 if (vh > B_BASE-1) vh = B_BASE-1;
752 #if DEBUG_SAINTEGER_UQRM
753 printf("guess formed from %d %d %d %d %d\n",
754 dend_stuff[i], dend_stuff[i-1] , dend_stuff[i-2],
755 isor_stuff[isor->used-1], isor_stuff[isor->used-2]);
756 printf("guess is %d\n", vh );
758 // Check if vh is too large (by 1). Calculate vh * isor into tmp
759 // and see if it exceeds the same length prefix of dend. If so,
760 // vh needs to be decremented.
762 for (j = 0; j < isor->used; j++) {
763 t = vh * isor_stuff[j] + carry;
767 tmp[isor->used] = carry;
768 delta = i - isor->used;
769 #if DEBUG_SAINTEGER_UQRM
770 printf("final carry is %d\n", carry);
771 printf("vh * isor is " );
772 for (j = isor->used; j >=0; j--) printf("%d ",tmp[j]);printf("\n");
773 printf("delta = %d\n", delta );
776 for (j = isor->used; j >= 0; j--) {
777 #if DEBUG_SAINTEGER_UQRM
778 printf ( "(%d,%d) ", (int)(tmp[j]), (int)(dend_stuff[j+delta]) );
780 if (tmp[j] > dend_stuff[j+delta]) {toolarge=1; break;};
781 if (tmp[j] < dend_stuff[j+delta]) break;
784 // If we did guess too large, decrement vh and subtract a copy of
785 // isor from tmp. This had better not go negative!
787 #if DEBUG_SAINTEGER_UQRM
788 printf ( "guess too large\n" );
792 for (j = 0; j < isor->used; j++) {
793 if (carry + isor_stuff[j] > tmp[j]) {
794 tmp[j] = (B_BASE + tmp[j]) - isor_stuff[j] - carry;
797 tmp[j] = tmp[j] - isor_stuff[j] - carry;
801 //if (carry > 0) {pp(isor);pp(dend);};
802 //myassert(carry == 0);
804 myassert(tmp[isor->used] > 0);
807 #if DEBUG_SAINTEGER_UQRM
808 printf("after adjustment of tmp ");
809 for (j = isor->used; j >=0; j--) printf("%d ",tmp[j]);
814 // Now vh really is the i'th quotient digit.
815 // Subtract (tmp << delta) from
818 for (j = 0; j <= isor->used; j++) {
819 if (carry + tmp[j] > dend_stuff[j+delta]) {
820 dend_stuff[j+delta] = (B_BASE+dend_stuff[j+delta]) - tmp[j] - carry;
823 dend_stuff[j+delta] = dend_stuff[j+delta] - tmp[j] - carry;
829 #if DEBUG_SAINTEGER_UQRM
830 printf("after final sub ");
831 for(j=dend->used; j>=0; j--) printf("%d ", dend_stuff[j]);
835 // park vh in the result array
836 #if DEBUG_SAINTEGER_UDIV
837 printf("[%d] <- %d\n", i-isor->used, vh );
839 dres->stuff[i-isor->used] = vh;
843 // Now we've got all the quotient digits. Zap leading zeroes.
844 dres->used = dend->used - isor->used + 1;
846 myassert(sane(dres));
848 // The remainder is in dend_stuff. Copy, divide by the original scaling
849 // factor, and zap leading zeroes.
850 mres->used = dend->used;
851 for (j = 0; j < dend->used; j++) mres->stuff[j] = dend_stuff[j];
853 myassert(sane(mres));
857 for (j = mres->used-1; j >= 0; j--) {
858 carry += mres->stuff[j];
859 mres->stuff[j] = carry/scaleup;
860 carry = B_BASE*(carry%scaleup);
862 myassert (carry == 0);
864 myassert(sane(mres));
873 /* --------------------------------------------------------------------------
875 * ------------------------------------------------------------------------*/
878 int main ( int argc, char** argv )
881 B *bi, *bj, *bk, *bm;
883 for (i = -10007; i <= 10007; i++) {
884 printf ( "i = %d\n", i );
886 t = size_fromInt(); bi = malloc(t); myassert(bi);
887 do_fromInt(i, t, bi);
889 t = do_toInt(bi); myassert(i == t);
891 for (j = -10007; j <= 10007; j++) {
893 t = size_fromInt(); bj = malloc(t); myassert(bj);
894 do_fromInt(j, t, bj);
896 t = do_toInt(bj); myassert(j == t);
899 t = size_add(bi,bj); bk = malloc(t); myassert(bk);
903 pp(bi); pp(bj); pp(bk);
910 t = size_sub(bi,bj); bk = malloc(t); myassert(bk);
914 pp(bi); pp(bj); pp(bk);
921 t = size_mul(bi,bj); bk = malloc(t); myassert(bk);
925 pp(bi); pp(bj); pp(bk);
933 bk = malloc(t); myassert(bk);
934 bm = malloc(t); myassert(bm);
935 do_qrm(bi,bj,t,bk,bm);
954 int main ( int argc, char** argv )
956 B *a, *b, *c, *d, *e;
957 a = fromInt(1); b=fromInt(9); pp(a); pp(b);
958 c = mkB( maxused_uqrm(a,b) );
959 d = mkB( maxused_uqrm(a,b) );
960 e = mkB( maxused_uadd(a,b) );
962 //uqrm(a,b,c,d); pp(c); pp(d);
968 /*-------------------------------------------------------------------------*/