![[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.