1 /*- 2 * %sccs.include.proprietary.c% 3 */ 4 5 #ifndef lint 6 static char sccsid[] = "@(#)msqrt.c 8.1 (Berkeley) 06/04/93"; 7 #endif /* not lint */ 8 9 #include <mp.h> 10 msqrt(a,b,r) MINT *a,*b,*r; 11 { MINT x,junk,y; 12 int j; 13 x.len=junk.len=y.len=0; 14 if(a->len<0) fatal("msqrt: neg arg"); 15 if(a->len==0) 16 { b->len=0; 17 r->len=0; 18 return(0); 19 } 20 if(a->len%2==1) x.len=(1+a->len)/2; 21 else x.len=1+a->len/2; 22 x.val=xalloc(x.len,"msqrt"); 23 for(j=0;j<x.len;x.val[j++]=0); 24 if(a->len%2==1) x.val[x.len-1]=0400; 25 else x.val[x.len-1]=1; 26 xfree(b); 27 xfree(r); 28 loop: 29 mdiv(a,&x,&y,&junk); 30 xfree(&junk); 31 madd(&x,&y,&y); 32 sdiv(&y,2,&y,(short *)&j); 33 if(mcmp(&x,&y)>0) 34 { xfree(&x); 35 move(&y,&x); 36 xfree(&y); 37 goto loop; 38 } 39 xfree(&y); 40 move(&x,b); 41 mult(&x,&x,&x); 42 msub(a,&x,r); 43 xfree(&x); 44 return(r->len); 45 } 46