|
|
1.1 root 1: #ifndef lint
2: static char sccsid[] = "@(#)msqrt.c 5.2 (Berkeley) 3/2/87";
3: #endif not lint
4:
5: #include <mp.h>
6: msqrt(a,b,r) MINT *a,*b,*r;
7: { MINT x,junk,y;
8: int j;
9: x.len=junk.len=y.len=0;
10: if(a->len<0) fatal("msqrt: neg arg");
11: if(a->len==0)
12: { b->len=0;
13: r->len=0;
14: return(0);
15: }
16: if(a->len%2==1) x.len=(1+a->len)/2;
17: else x.len=1+a->len/2;
18: x.val=xalloc(x.len,"msqrt");
19: for(j=0;j<x.len;x.val[j++]=0);
20: if(a->len%2==1) x.val[x.len-1]=0400;
21: else x.val[x.len-1]=1;
22: xfree(b);
23: xfree(r);
24: loop:
25: mdiv(a,&x,&y,&junk);
26: xfree(&junk);
27: madd(&x,&y,&y);
28: sdiv(&y,2,&y,(short *)&j);
29: if(mcmp(&x,&y)>0)
30: { xfree(&x);
31: move(&y,&x);
32: xfree(&y);
33: goto loop;
34: }
35: xfree(&y);
36: move(&x,b);
37: mult(&x,&x,&x);
38: msub(a,&x,r);
39: xfree(&x);
40: return(r->len);
41: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.