[unix-history] / usr / src / libmp / msqrt.c

#include <mp.h>
msqrt(a,b,r) MINT *a,*b,*r;
{	MINT x,junk,y;
	int j;
	x.len=junk.len=y.len=0;
	if(a->len<0) fatal("msqrt: neg arg");
	if(a->len==0)
	{	b->len=0;
		r->len=0;
		return(0);
	}
	if(a->len%2==1) x.len=(1+a->len)/2;
	else x.len=1+a->len/2;
	x.val=xalloc(x.len,"msqrt");
	for(j=0;j<x.len;x.val[j++]=0);
	if(a->len%2==1) x.val[x.len-1]=0400;
	else x.val[x.len-1]=1;
	xfree(b);
	xfree(r);
loop:
	mdiv(a,&x,&y,&junk);
	xfree(&junk);
	madd(&x,&y,&y);
	sdiv(&y,2,&y,(short *)&j);
	if(mcmp(&x,&y)>0)
	{	xfree(&x);
		move(&y,&x);
		xfree(&y);
		goto loop;
	}
	xfree(&y);
	move(&x,b);
	mult(&x,&x,&x);
	msub(a,&x,r);
	xfree(&x);
	return(r->len);
}
Commit	Line	Data
168f834f PW	1	#include <mp.h>
	2	msqrt(a,b,r) MINT a,b,*r;
	3	{ MINT x,junk,y;
	4	int j;
	5	x.len=junk.len=y.len=0;
	6	if(a->len<0) fatal("msqrt: neg arg");
	7	if(a->len==0)
	8	{ b->len=0;
	9	r->len=0;
	10	return(0);
	11	}
	12	if(a->len%2==1) x.len=(1+a->len)/2;
	13	else x.len=1+a->len/2;
	14	x.val=xalloc(x.len,"msqrt");
	15	for(j=0;j<x.len;x.val[j++]=0);
	16	if(a->len%2==1) x.val[x.len-1]=0400;
	17	else x.val[x.len-1]=1;
	18	xfree(b);
	19	xfree(r);
	20	loop:
	21	mdiv(a,&x,&y,&junk);
	22	xfree(&junk);
	23	madd(&x,&y,&y);
	24	sdiv(&y,2,&y,(short *)&j);
	25	if(mcmp(&x,&y)>0)
	26	{ xfree(&x);
	27	move(&y,&x);
	28	xfree(&y);
	29	goto loop;
	30	}
	31	xfree(&y);
	32	move(&x,b);
	33	mult(&x,&x,&x);
	34	msub(a,&x,r);
	35	xfree(&x);
	36	return(r->len);
	37	}