![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to bundled CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10no / libmp|
Return to bundled CVS log | Up to [Research Unix] / researchv10no / libmp |
1.1 root 1: # To unbundle, sh this file
2: echo gcd.c 1>&2
3: sed 's/.//' >gcd.c <<'//GO.SYSIN DD gcd.c'
4: -#include "mp.h"
5: -mgcd(a,b,c) mint *a,*b,*c;
6: -{ mint *x,*y,*z,*w;
7: - x = itom(0), y = itom(0), z = itom(0), w = itom(0);
8: - move(a,x);
9: - move(b,y);
10: - while(y->len!=0){
11: - mdiv(x,y,w,z);
12: - move(y,x);
13: - move(z,y);
14: - }
15: - move(x,c);
16: - mfree(x);
17: - mfree(y);
18: - mfree(z);
19: - mfree(w);
20: -}
21: -invert(a, b, c) mint *a, *b, *c;
22: -{ mint x, y, z, w, Anew, Aold;
23: - int i = 0;
24: - x.len = y.len = z.len = w.len = Aold.len = 0;
25: - Anew.len = 1;
26: - Anew.val = xalloc(1, "invert");
27: - *Anew.val = 1;
28: - move(b, &x);
29: - move(a, &y);
30: - while(y.len != 0)
31: - { mdiv(&x, &y, &w, &z);
32: - move(&Anew, &x);
33: - mult(&w, &Anew, &Anew);
34: - madd(&Anew, &Aold, &Anew);
35: - move(&x, &Aold);
36: - move(&y, &x);
37: - move(&z, &y);
38: - i++;
39: - }
40: - move(&Aold, c);
41: - if( (i&01) == 0) msub(b, c, c);
42: - xfree(&x);
43: - xfree(&y);
44: - xfree(&z);
45: - xfree(&w);
46: - xfree(&Aold);
47: - xfree(&Anew);
48: -}
49: -
50: -lineq(a, b, x, y, u) /* ax + by = u */
51: -mint *a, *b, *x, *u, *y;
52: -{ mint *at, *bt, *xo, *yo, *q, *r, *z;
53: - static mint *zero, *one;
54: - int i;
55: - at = itom(0), bt = itom(0), xo = itom(0);
56: - yo = itom(0), q = itom(0), r = itom(0), z = itom(0);
57: - if(zero == 0) zero = itom(0);
58: - if(one == 0) one = itom(1);
59: - move(a, at);
60: - move(b, bt);
61: - move(zero, x);
62: - move(one, xo);
63: - move(zero, yo);
64: - move(one, y);
65: - if(bt->len == 0) {
66: - move(one, x);
67: - move(zero, y);
68: - move(a, u);
69: - goto out;
70: - }
71: - for(i = 0; ; i++) {
72: - mdiv(at, bt, q, r);
73: - if(r->len == 0)
74: - break;
75: - move(xo, z);
76: - move(x, xo);
77: - mult(q, xo, x);
78: - madd(z, x, x);
79: - move(yo, z);
80: - move(y, yo);
81: - mult(q, yo, y);
82: - madd(z, y, y);
83: - move(bt, at);
84: - move(r, bt);
85: - }
86: - move(bt, u);
87: - if(i & 1)
88: - y->len = -y->len;
89: - else
90: - x->len = -x->len;
91: -out:
92: - mfree(z), mfree(r), mfree(q), mfree(yo);
93: - mfree(xo), mfree(bt), mfree(at);
94: -}
95: //GO.SYSIN DD gcd.c
96: echo halloc.c 1>&2
97: sed 's/.//' >halloc.c <<'//GO.SYSIN DD halloc.c'
98: -/* list allocator */
99: -
100: -typedef struct hdr {
101: - struct hdr *next;
102: - struct hdr **list;
103: - short d[1];
104: -} hdr;
105: -extern char *malloc();
106: -
107: -hdr *h1, *h4, *h6;
108: -
109: -hfree(d)
110: -short *d;
111: -{ hdr *p;
112: - if(d == (short *)1)
113: - return;
114: - p = (hdr *)((char *)d - 2 * sizeof(hdr *));
115: - if(p->list == 0) {
116: - free((char *)p);
117: - return;
118: - }
119: - p->next = *p->list;
120: - *p->list = p;
121: -}
122: -
123: -short *
124: -halloc(n)
125: -{ short *x;
126: - hdr *p;
127: - if(n <= 0)
128: - return((short *)1); /* and unusable */
129: - if(n <= 2) {
130: - if(!h1)
131: - makelist(4, &h1);
132: - x = h1->d;
133: - h1 = h1->next;
134: - return(x);
135: - }
136: - if(n <= 8) {
137: - if(!h4)
138: - makelist(16, &h4);
139: - x = h4->d;
140: - h4 = h4->next;
141: - return(x);
142: - }
143: - if(n <= 64) {
144: - if(!h6)
145: - makelist(128, &h6);
146: - x = h6->d;
147: - h6 = h6->next;
148: - return(x);
149: - }
150: - p = (hdr *) malloc(sizeof(hdr) + (n - 1) * sizeof(short));
151: - p->list = 0;
152: - return(p->d);
153: -}
154: -
155: -makelist(n, p)
156: -hdr **p;
157: -{ int i, d;
158: - char *s;
159: - hdr *h;
160: - d = sizeof(hdr) + sizeof(short) * (n-1);
161: - d = (d + 3) & (~3);
162: - s = malloc(128 * d);
163: - h = *p = (hdr *)s;
164: - for(i = 0; i < 127; i++) {
165: - h->next = (hdr *)(s + d);
166: - h->list = p;
167: - s += d;
168: - h = (hdr *)s;
169: - }
170: - h->next = 0;
171: - h->list = p;
172: -}
173: //GO.SYSIN DD halloc.c
174: echo madd.c 1>&2
175: sed 's/.//' >madd.c <<'//GO.SYSIN DD madd.c'
176: -#include "mp.h"
177: -m_add(a,b,c) mint *a,*b,*c;
178: -{ register int carry,i;
179: - register int x;
180: - register short *cval;
181: - cval=xalloc(a->len+1,"m_add");
182: - carry=0;
183: - for(i=0;i<b->len;i++)
184: - { x=carry+a->val[i]+b->val[i];
185: - if(x&0100000)
186: - { carry=1;
187: - cval[i]=x&077777;
188: - }
189: - else
190: - { carry=0;
191: - cval[i]=x;
192: - }
193: - }
194: - for(;i<a->len;i++)
195: - { x=carry+a->val[i];
196: - if(x&0100000) cval[i]=x&077777;
197: - else
198: - { carry=0;
199: - cval[i]=x;
200: - }
201: -
202: - }
203: - if(carry==1)
204: - { cval[i]=1;
205: - c->len=i+1;
206: - }
207: - else c->len=a->len;
208: - c->val=cval;
209: - if(c->len==0) shfree(cval);
210: - return;
211: -}
212: -madd(a,b,c) mint *a,*b,*c;
213: -{ mint x,y,z;
214: - int sign;
215: - x.len=a->len;
216: - x.val=a->val;
217: - y.len=b->len;
218: - y.val=b->val;
219: - z.len=0;
220: - sign=1;
221: - if(x.len>=0)
222: - if(y.len>=0)
223: - if(x.len>=y.len) m_add(&x,&y,&z);
224: - else m_add(&y,&x,&z);
225: - else
226: - { y.len= -y.len;
227: - msub(&x,&y,&z);
228: - }
229: - else if(y.len<=0)
230: - { x.len = -x.len;
231: - y.len= -y.len;
232: - sign= -1;
233: - madd(&x,&y,&z);
234: - }
235: - else
236: - { x.len= -x.len;
237: - msub(&y,&x,&z);
238: - }
239: - xfree(c);
240: - c->val=z.val;
241: - c->len=sign*z.len;
242: - return;
243: -}
244: -m_sub(a,b,c) mint *a,*b,*c;
245: -{ register int x,i;
246: - register int borrow;
247: - short one;
248: - mint mone;
249: - one=1; mone.len= 1; mone.val= &one;
250: - c->val=xalloc(a->len,"m_sub");
251: - borrow=0;
252: - for(i=0;i<b->len;i++)
253: - { x=borrow+a->val[i]-b->val[i];
254: - if(x&0100000)
255: - { borrow= -1;
256: - c->val[i]=x&077777;
257: - }
258: - else
259: - { borrow=0;
260: - c->val[i]=x;
261: - }
262: - }
263: - for(;i<a->len;i++)
264: - { x=borrow+a->val[i];
265: - if(x&0100000) c->val[i]=x&077777;
266: - else
267: - { borrow=0;
268: - c->val[i]=x;
269: - }
270: - }
271: - if(borrow<0)
272: - { for(i=0;i<a->len;i++) c->val[i] ^= 077777;
273: - c->len=a->len;
274: - madd(c,&mone,c);
275: - }
276: - for(i=a->len-1;i>=0;--i)
277: - if(c->val[i]>0) {
278: - if(borrow==0) c->len=i+1;
279: - else c->len= -i-1;
280: - return;
281: - }
282: - shfree(c->val);
283: - return;
284: -}
285: -msub(a,b,c) mint *a,*b,*c;
286: -{ mint x,y,z;
287: - int sign;
288: - x.len=a->len;
289: - y.len=b->len;
290: - x.val=a->val;
291: - y.val=b->val;
292: - z.len=0;
293: - sign=1;
294: - if(x.len>=0)
295: - if(y.len>=0)
296: - if(x.len>=y.len) m_sub(&x,&y,&z);
297: - else
298: - { sign= -1;
299: - msub(&y,&x,&z);
300: - }
301: - else
302: - { y.len= -y.len;
303: - madd(&x,&y,&z);
304: - }
305: - else if(y.len<=0)
306: - { sign= -1;
307: - x.len= -x.len;
308: - y.len= -y.len;
309: - msub(&x, &y, &z);
310: - }
311: - else
312: - { x.len= -x.len;
313: - madd(&x,&y,&z);
314: - sign= -1;
315: - }
316: - if(a==c && x.len!=0) xfree(a);
317: - else if(b==c && y.len!=0) xfree(b);
318: - else xfree(c);
319: - c->val=z.val;
320: - c->len=sign*z.len;
321: - return;
322: -}
323: //GO.SYSIN DD madd.c
324: echo mdiv.c 1>&2
325: sed 's/.//' >mdiv.c <<'//GO.SYSIN DD mdiv.c'
326: -#include "mp.h"
327: -int mpdivdebug;
328: -mdiv(a,b,q,r) mint *a,*b,*q,*r;
329: -{ mint *x,*y, *z;
330: - int alen, blen;
331: - x = itom(0), y = itom(0);
332: - z = itom(0);
333: - move(a, x);
334: - move(b, y);
335: - alen = x->len;
336: - if(x->len<0) {x->len= -x->len;}
337: - blen = y->len;
338: - if(y->len<0) {y->len= -y->len;}
339: - xfree(q);
340: - xfree(r);
341: - m_div(x,y,q,r);
342: - if(mpdivdebug) {
343: - mult(y, q, z);
344: - madd(z, r, z);
345: - if(mcmp(z, x) != 0) {
346: - mout(a);
347: - mout(b);
348: - fatal("mdiv err");
349: - }
350: - }
351: - if(alen < 0) {
352: - mint o;
353: - short i = 1;
354: - o.len = 1;
355: - o.val = &i;
356: - if(r->len == 0) {
357: - if(blen > 0)
358: - q->len = -q->len;
359: - goto out;
360: - }
361: - msub(r, y, r);
362: - r->len = -r->len;
363: - madd(q, &o, q);
364: - if(blen > 0)
365: - q->len = -q->len;
366: - goto out;
367: - }
368: - if(blen < 0)
369: - q->len = -q->len;
370: -out:
371: - mfree(z);
372: - mfree(y);
373: - mfree(x);
374: -}
375: -#ifdef vax
376: -union zz {
377: - long a;
378: - struct {
379: - unsigned int lo:15, hi:15;
380: - } b;
381: -};
382: -#endif
383: -m_dsb(q,n,a,b) short *a,*b;
384: -{ long int qx, u;
385: - union zz x;
386: - int borrow,j;
387: - qx=q;
388: - x.a = 0;
389: - for(j = 0; j < n; j++) {
390: - x.a = qx * a[j] + x.b.hi;
391: - if((b[j] -= x.b.lo) < 0) {
392: - b[j] += (1 << 15);
393: - x.b.hi += 1;
394: - }
395: - }
396: - if((b[j] -= x.b.hi) >= 0)
397: - return(0);
398: - borrow=0;
399: - for(j=0;j<n;j++)
400: - { u=a[j]+b[j]+borrow;
401: - if(u & 0100000)
402: - borrow = 1;
403: - else borrow=0;
404: - b[j]=u&077777;
405: - }
406: - { return(1);}
407: -}
408: -m_trq(v1,v2,u1,u2,u3)
409: -{ long int d;
410: - long int x1;
411: - if(u1==v1) d=077777;
412: - else d=(u1*0100000L+u2)/v1;
413: - while(1)
414: - { x1=u1*0100000L+u2-v1*d;
415: - x1=x1*0100000L+u3-v2*d;
416: - if(x1<0) d=d-1;
417: - else {return(d);}
418: - }
419: -}
420: -m_div(a,b,q,r) mint *a,*b,*q,*r;
421: -{ mint u,v,x,w;
422: - short d,*qval;
423: - int qq,n,v1,v2,j;
424: - u.len=v.len=x.len=w.len=0;
425: - if(b->len==0) { fatal("mdiv divide by zero"); return;}
426: - if(b->len==1)
427: - { r->val=xalloc(1,"m_div1");
428: - sdiv(a,b->val[0],q,r->val);
429: - if(r->val[0]==0)
430: - { shfree(r->val);
431: - r->len=0;
432: - }
433: - else r->len=1;
434: - return;
435: - }
436: - if(a->len < b->len)
437: - { q->len=0;
438: - r->len=a->len;
439: - r->val=xalloc(r->len,"m_div2");
440: - for(qq=0;qq<r->len;qq++) r->val[qq]=a->val[qq];
441: - return;
442: - }
443: - x.len=1;
444: - x.val = &d;
445: - n=b->len;
446: - d=0100000L/(b->val[n-1]+1L);
447: - mult(a,&x,&u); /*subtle: relies on fact that mult allocates extra space */
448: - mult(b,&x,&v);
449: - v1=v.val[n-1];
450: - v2=v.val[n-2];
451: - qval=xalloc(a->len-n+1,"m_div3");
452: - for(j=a->len-n;j>=0;j--)
453: - { qq=m_trq(v1,v2,u.val[j+n],u.val[j+n-1],u.val[j+n-2]);
454: - if(m_dsb(qq,n,v.val,&(u.val[j]))) qq -= 1;
455: - qval[j]=qq;
456: - }
457: - x.len=n;
458: - x.val=u.val;
459: - mcan(&x);
460: - sdiv(&x,d,&w,(short *)&qq);
461: - r->len=w.len;
462: - r->val=w.val;
463: - q->val=qval;
464: - qq=a->len-n+1;
465: - if(qq>0 && qval[qq-1]==0) qq -= 1;
466: - q->len=qq;
467: - if(qq==0) shfree(qval);
468: - if(x.len!=0) xfree(&u);
469: - xfree(&v);
470: - return;
471: -}
472: //GO.SYSIN DD mdiv.c
473: echo mout.c 1>&2
474: sed 's/.//' >mout.c <<'//GO.SYSIN DD mout.c'
475: -#include <stdio.h>
476: -#include "mp.h"
477: -m_in(a,b,f) mint *a; FILE *f;
478: -{ mint x,y,ten;
479: - int sign,c;
480: - short qten,qy;
481: - xfree(a);
482: - sign=1;
483: - ten.len=1;
484: - ten.val= &qten;
485: - qten=b;
486: - x.len=0;
487: - y.len=1;
488: - y.val= &qy;
489: - while((c=getc(f))!=EOF)
490: - switch(c)
491: - {
492: - case '\\': getc(f);
493: - continue;
494: - case '\t':
495: - case '\n':
496: - goto done;
497: - case ' ':
498: - continue;
499: - case '-': sign = -sign;
500: - continue;
501: - default: if(c>='0' && c<= '9')
502: - { qy=c-'0';
503: - mult(&x,&ten,a);
504: - madd(a,&y,a);
505: - move(a,&x);
506: - continue;
507: - }
508: - else
509: - { (void) ungetc(c,stdin);
510: -done:
511: - a->len *= sign;
512: - mcan(a);
513: - return(0);
514: - }
515: - }
516: - return(EOF);
517: -}
518: -m_out(a,b,f) mint *a; FILE *f;
519: -{ int sign,xlen,i;
520: - short r;
521: - mint x;
522: - char *obuf;
523: - register char *bp;
524: - sign=1;
525: - xlen=a->len;
526: - if(xlen<0)
527: - { xlen= -xlen;
528: - sign= -1;
529: - }
530: - if(xlen==0)
531: - { fprintf(f,"0\n");
532: - return;
533: - }
534: - x.len=xlen;
535: - x.val=xalloc(xlen,"m_out");
536: - for(i=0;i<xlen;i++) x.val[i]=a->val[i];
537: - obuf=(char *)malloc(7*xlen);
538: - bp=obuf+7*xlen-1;
539: - *bp--=0;
540: - while(x.len>0)
541: - { for(i=0;i<10&&x.len>0;i++)
542: - { sdiv(&x,b,&x,&r);
543: - *bp--=r+'0';
544: - }
545: - if(x.len>0) *bp--=' ';
546: - }
547: - if(sign==-1) *bp--='-';
548: - fprintf(f,"%s\n",bp+1);
549: - free(obuf);
550: - xfree(&x);
551: - return;
552: -}
553: -sdiv(a,n,q,r) mint *a,*q; short *r;
554: -{ mint x,y;
555: - int sign;
556: - sign=1;
557: - x.len=a->len;
558: - x.val=a->val;
559: - if(n<0) {
560: - sign= -sign;
561: - n= -n;
562: - }
563: - if(x.len<0) {
564: - sign = -sign;
565: - x.len= -x.len;
566: - }
567: - else if(x.len == 0) {
568: - xfree(q);
569: - *r = 0;
570: - return;
571: - }
572: - s_div(&x,n,&y,r);
573: - xfree(q);
574: - q->val=y.val;
575: - q->len=sign*y.len;
576: - *r = *r*sign;
577: - return;
578: -}
579: -s_div(a,n,q,r) mint *a,*q; short *r;
580: -{ int qlen,i;
581: - long int x;
582: - short *qval;
583: - x=0;
584: - qlen=a->len;
585: - qval=xalloc(qlen,"s_div");
586: - for(i=qlen-1;i>=0;i--)
587: - {
588: - x=x*0100000L+a->val[i];
589: - qval[i]=x/n;
590: - x=x%n;
591: - }
592: - *r=x;
593: - if(qval[qlen-1]==0) qlen--;
594: - q->len=qlen;
595: - q->val=qval;
596: - if(qlen==0) shfree(qval);
597: - return;
598: -}
599: -min(a) mint *a;
600: -{
601: - return(m_in(a,10,stdin));
602: -}
603: -omin(a) mint *a;
604: -{
605: - return(m_in(a,8,stdin));
606: -}
607: -mout(a) mint *a;
608: -{
609: - m_out(a,10,stdout);
610: -}
611: -omout(a) mint *a;
612: -{
613: - m_out(a,8,stdout);
614: -}
615: -fmout(a,f) mint *a; FILE *f;
616: -{ m_out(a,10,f);
617: -}
618: -fmin(a,f) mint *a; FILE *f;
619: -{
620: - return(m_in(a,10,f));
621: -}
622: //GO.SYSIN DD mout.c
623: echo msqrt.c 1>&2
624: sed 's/.//' >msqrt.c <<'//GO.SYSIN DD msqrt.c'
625: -#include "mp.h"
626: -msqrt(a,b,r) mint *a,*b,*r;
627: -{ mint x,junk,y;
628: - int j;
629: - x.len=junk.len=y.len=0;
630: - if(a->len<0) fatal("msqrt: neg arg");
631: - if(a->len==0)
632: - { b->len=0;
633: - r->len=0;
634: - return(0);
635: - }
636: - if(a->len%2==1) x.len=(1+a->len)/2;
637: - else x.len=1+a->len/2;
638: - x.val=xalloc(x.len,"msqrt");
639: - for(j=0;j<x.len;x.val[j++]=0);
640: - if(a->len%2==1) x.val[x.len-1]=0400;
641: - else x.val[x.len-1]=1;
642: - xfree(b);
643: - xfree(r);
644: -loop:
645: - mdiv(a,&x,&y,&junk);
646: - xfree(&junk);
647: - madd(&x,&y,&y);
648: - sdiv(&y,2,&y,(short *)&j);
649: - if(mcmp(&x,&y)>0)
650: - { xfree(&x);
651: - move(&y,&x);
652: - xfree(&y);
653: - goto loop;
654: - }
655: - xfree(&y);
656: - move(&x,b);
657: - mult(&x,&x,&x);
658: - msub(a,&x,r);
659: - xfree(&x);
660: - return(r->len);
661: -}
662: -
663: -/* pathology: n<= 0 => r=rem=0, num <= 0, r=rem=0 */
664: -/* this is a lazy implementation, assuming n>=3 => root is a legal double */
665: -mroot(n, num, r, rem)
666: -mint *num, *r, *rem;
667: -{ extern double log(), exp();
668: - double z;
669: - int i;
670: - static mint *xn, *xn1, *top, *bot;
671: - static init;
672: - if(!init++) {
673: - xn = itom(0), xn1 = itom(0), top = itom(0), bot = itom(0);
674: - }
675: - if(n < 0 || num->len <= 0) {
676: - msub(r, r, r);
677: - move(r, rem);
678: - return;
679: - }
680: - if(n == 1) {
681: - move(num, r);
682: - msub(rem, rem, rem);
683: - return;
684: - }
685: - if(n == 2) {
686: - msqrt(num, r, rem);
687: - return;
688: - }
689: - /* compute approx bigger than root */
690: - for(z = 0, i = 0; i < num->len; i++)
691: - z = z/32768 + num->val[i];
692: - /* num = z * 2^15*(len-1) */
693: - z = exp((log(z) + 15*(num->len-1)*log(2.))/n);
694: - z = 1.01 * z + 1; /* make sure it is bigger than root */
695: - dtom(z, r);
696: - /* ((n-1)*x^n+num)/(n*x^(n-1))*/
697: - for(;;) {
698: - rpow(r, n - 1, xn1);
699: - mult(r, xn1, xn);
700: - imult(n-1, xn, top);
701: - madd(top, num, top);
702: - imult(n, xn1, bot);
703: - mdiv(top, bot, xn1, rem);
704: - switch(mcmp(xn1, r)) {
705: - case -1:
706: - move(xn1, r);
707: - continue;
708: - case 0:
709: - move(xn1, r);
710: - msub(num, xn, rem);
711: - return;
712: - case 1: /* since r was too large to start with */
713: - msub(num, xn, rem);
714: - return;
715: - }
716: - }
717: -}
718: -
719: -static
720: -imult(n, a, b)
721: -mint *a, *b;
722: -{ mint *x;
723: - x = itom(n);
724: - mult(x, a, b);
725: - msub(x, x,x);
726: -}
727: //GO.SYSIN DD msqrt.c
728: echo mult.c 1>&2
729: sed 's/.//' >mult.c <<'//GO.SYSIN DD mult.c'
730: -#include <mp.h>
731: -#undef unfast
732: -mult(a,b,c) mint *a,*b,*c;
733: -{ mint x,y,z;
734: - int sign;
735: - sign = 1;
736: - x.val=a->val;
737: - y.val=b->val;
738: - z.len=0;
739: - if(a->len<0)
740: - { x.len= -a->len;
741: - sign= -sign;
742: - }
743: - else x.len=a->len;
744: - if(b->len<0)
745: - { y.len= -b->len;
746: - sign= -sign;
747: - }
748: - else y.len=b->len;
749: - if(x.len == 0 || y.len == 0)
750: - z.len = 0;
751: - else if(x.len<y.len) m_mult(&y,&x,&z);
752: - else m_mult(&x,&y,&z);
753: - xfree(c);
754: - if(sign<0) c->len= -z.len;
755: - else c->len=z.len;
756: - c->val=z.val;
757: -}
758: -#define S2 x=a->val[j];
759: -#define S3 x=x*b->val[i-j];
760: -#ifdef unfast
761: -#define S4 tradd(&carry,&sum,x);
762: -#else
763: -#define S4 sum.xx += x; if(sum.yy.high & 0100000) {sum.yy.high &= 077777; carry++;}
764: -#endif
765: -#define S5 c->val[i]=sum.yy.low&077777;
766: -#define S6 sum.xx=sum.xx>>15;
767: -#define S7 sum.yy.high=carry;
768: -m_mult(a,b,c)
769: -register mint *a,*b,*c;
770: -{ register long x;
771: - union {long xx; struct half yy;} sum;
772: - int carry;
773: - int i,j;
774: - c->val=xalloc(a->len+b->len,"m_mult");
775: - sum.xx=0;
776: - for(i=0;i<b->len;i++)
777: - { carry=0;
778: - for(j=0;j<i+1;j++)
779: - { S2
780: - S3
781: - S4
782: - }
783: - S5
784: - S6
785: - S7
786: - }
787: - for(;i<a->len;i++)
788: - { carry=0;
789: - for(j=i-b->len+1;j<i+1;j++)
790: - { S2
791: - S3
792: - S4
793: - }
794: - S5
795: - S6
796: - S7
797: - }
798: - for(;i<a->len+b->len;i++)
799: - { carry=0;
800: - for(j=i-b->len+1;j<a->len;j++)
801: - { S2
802: - S3
803: - S4
804: - }
805: - S5
806: - S6
807: - S7
808: - }
809: - if(i == 0 || c->val[i-1]!=0)
810: - c->len=a->len+b->len;
811: - else c->len=a->len+b->len-1;
812: - return;
813: -}
814: -tradd(a,b,c)
815: -long c;
816: -int *a;
817: -register union {long xx; struct half yy;} *b;
818: -{
819: - b->xx += c;
820: - if(b->yy.high&0100000)
821: - { b->yy.high= b->yy.high&077777;
822: - *a += 1;
823: - }
824: - return;
825: -}
826: //GO.SYSIN DD mult.c
827: echo pow.c 1>&2
828: sed 's/.//' >pow.c <<'//GO.SYSIN DD pow.c'
829: -#include "mp.h"
830: -static
831: -mpow(a,b,c,d) mint *a,*b,*c,*d;
832: -{ int i,j,n;
833: - mint x,y;
834: - x.len=y.len=0;
835: - xfree(d);
836: - d->len=1;
837: - d->val=xalloc(1,"mpow");
838: - *d->val=1;
839: - for(j=0;j<b->len;j++)
840: - { n=b->val[b->len-j-1];
841: - for(i=0;i<15;i++)
842: - { mult(d,d,&x);
843: - mdiv(&x,c,&y,d);
844: - if((n=n<<1)&0100000)
845: - { mult(a,d,&x);
846: - mdiv(&x,c,&y,d);
847: - }
848: - }
849: - }
850: - xfree(&x);
851: - xfree(&y);
852: - return;
853: -}
854: -rpow(a,n,b) mint *a,*b;
855: -{ mint x,y;
856: - int i;
857: - x.len=1;
858: - x.val=xalloc(1,"rpow");
859: - *x.val=n;
860: - y.len=n*a->len+4;
861: - y.val=xalloc(y.len,"rpow2");
862: - for(i=0;i<y.len;i++) y.val[i]=0;
863: - y.val[y.len-1]=010000;
864: - xfree(b);
865: - mpow(a,&x,&y,b);
866: - xfree(&x);
867: - xfree(&y);
868: - return;
869: -}
870: //GO.SYSIN DD pow.c
871: echo util.c 1>&2
872: sed 's/.//' >util.c <<'//GO.SYSIN DD util.c'
873: -extern char *malloc();
874: -#include "stdio.h"
875: -#include "mp.h"
876: -move(a,b) mint *a,*b;
877: -{ int i,j;
878: - if(a == b)
879: - return;
880: - xfree(b);
881: - b->len=a->len;
882: - if((i=a->len)<0)
883: - i = -i;
884: - if(i==0)
885: - return;
886: - b->val=xalloc(i,"move");
887: - for(j=0;j<i;j++)
888: - b->val[j]=a->val[j];
889: - return;
890: -}
891: -dummy(){ }
892: -/*ARGSUSED*/
893: -short *xalloc(nint,s) char *s;
894: -{ short *i;
895: - extern short *halloc();
896: - i = halloc(nint);
897: -#ifdef DBG
898: - i=(short *)malloc(2*(unsigned)nint+4);
899: - if(dbg) fprintf(stderr, "%s: %o\n",s,i);
900: -#endif
901: - if(i!=NULL) return(i);
902: - fatal("mp: no free space");
903: - return(0);
904: -}
905: -fatal(s) char *s;
906: -{
907: - fprintf(stderr,"%s\n",s);
908: - (void) fflush(stdout);
909: - sleep(2);
910: - abort();
911: -}
912: -xfree(c) mint *c;
913: -{
914: -#ifdef DBG
915: - if(dbg) fprintf(stderr, "xfree ");
916: -#endif
917: - if(c->len==0) return;
918: - /*shfree(c->val);*/
919: - hfree(c->val);
920: - c->len=0;
921: - return;
922: -}
923: -
924: -mfree(a)
925: -mint *a;
926: -{
927: - xfree(a);
928: - hfree(a);
929: -}
930: -mcan(a) mint *a;
931: -{ int i,j;
932: - if((i=a->len)==0) return;
933: - else if(i<0) i= -i;
934: - for(j=i;j>0 && a->val[j-1]==0;j--);
935: - if(j==i) return;
936: - if(j==0) {
937: - xfree(a);
938: - return;
939: - }
940: - if(a->len > 0) a->len=j;
941: - else a->len = -j;
942: -}
943: -mint *itom(n)
944: -{ mint *a;
945: - a=(mint *)xalloc(2,"itom");
946: - if(n>0) {
947: - a->len=1;
948: - a->val=xalloc(1,"itom1");
949: - *a->val=n;
950: - return(a);
951: - }
952: - else if(n<0) {
953: - a->len = -1;
954: - a->val=xalloc(1,"itom2");
955: - *a->val= -n;
956: - return(a);
957: - }
958: - else {
959: - a->len=0;
960: - return(a);
961: - }
962: -}
963: -mcmp(a,b) mint *a,*b;
964: -{ mint c;
965: - int res;
966: - if(a->len < b->len)
967: - return(-1);
968: - if(a->len > b->len)
969: - return(1);
970: - c.len=0;
971: - msub(a,b,&c);
972: - res=c.len;
973: - xfree(&c);
974: - if(res < 0)
975: - return(-1);
976: - else if(res == 0)
977: - return(0);
978: - else
979: - return(1);
980: -}
981: -
982: -dtom(z, r)
983: -double z;
984: -mint *r;
985: -{ int i, sgn;
986: - static mint *c;
987: - if(!c) {
988: - c = itom(16384);
989: - madd(c, c, c);
990: - }
991: - if(z < 0) {
992: - sgn = 1;
993: - z = -z;
994: - }
995: - else
996: - sgn = 0;
997: - for(i = 0; z >= 32768; i++)
998: - z /= 32768;
999: - move(c, r);
1000: - r->len = 1;
1001: - r->val[0] = z;
1002: - while(--i >= 0) {
1003: - z -= r->val[0];
1004: - z *= 32768;
1005: - mult(r, c, r);
1006: - r->val[0] = z;
1007: - }
1008: - if(sgn)
1009: - r->len = -r->len;
1010: -}
1011: //GO.SYSIN DD util.c
1012: echo Makefile 1>&2
1013: sed 's/.//' >Makefile <<'//GO.SYSIN DD Makefile'
1014: -DESTDIR=/usr/lib
1015: -CFLAGS=
1016: -OBJS= pow.o gcd.o msqrt.o mdiv.o mout.o mult.o madd.o util.o halloc.o
1017: -
1018: -libmp.a: $(OBJS)
1019: - ar cr libmp.a $(OBJS)
1020: -
1021: -install: libmp.a
1022: - cp libmp.a ${DESTDIR}
1023: - ranlib ${DESTDIR}/libmp.a
1024: -
1025: -clean:
1026: - rm -f *.o libmp.a
1027: //GO.SYSIN DD Makefile
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.