OpenSolaris_b135/lib/libmp/common/msqrt.c

/*	Copyright (c) 1984, 1986, 1987, 1988, 1989 AT&T	*/
/*	  All Rights Reserved  	*/


/*
 * Copyright (c) 1980 Regents of the University of California.
 * All rights reserved.  The Berkeley software License Agreement
 * specifies the terms and conditions for redistribution.
 */
/* 	Portions Copyright(c) 1988, Sun Microsystems Inc.	*/
/*	All Rights Reserved					*/

/*
 * Copyright (c) 1997, by Sun Microsystems, Inc.
 * All rights reserved.
 */

#ident	"%Z%%M%	%I%	%E% SMI"	/* SVr4.0 1.1	*/

/* LINTLIBRARY */

#include <mp.h>
#include <sys/types.h>
#include "libmp.h"

int
mp_msqrt(MINT *a, MINT *b, MINT *r)
{
	MINT a0, x, junk, y;
	int j;

	a0.len = junk.len = y.len = 0;
	if (a->len < 0)
		_mp_fatal("mp_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 = _mp_xalloc(x.len, "mp_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;
	_mp_move(a, &a0);
	_mp_xfree(b);
	_mp_xfree(r);
loop:
	mp_mdiv(&a0, &x, &y, &junk);
	_mp_xfree(&junk);
	mp_madd(&x, &y, &y);
	mp_sdiv(&y, 2, &y, (short *) &j);
	if (mp_mcmp(&x, &y) > 0) {
		_mp_xfree(&x);
		_mp_move(&y, &x);
		_mp_xfree(&y);
		goto loop;
	}
	_mp_xfree(&y);
	_mp_move(&x, b);
	mp_mult(&x, &x, &x);
	mp_msub(&a0, &x, r);
	_mp_xfree(&x);
	_mp_xfree(&a0);
	return (r->len);
}