xref: /illumos-gate/usr/src/lib/libmp/common/msqrt.c (revision 7c478bd9)
1*7c478bd9Sstevel@tonic-gate /*	Copyright (c) 1984, 1986, 1987, 1988, 1989 AT&T	*/
2*7c478bd9Sstevel@tonic-gate /*	  All Rights Reserved  	*/
3*7c478bd9Sstevel@tonic-gate 
4*7c478bd9Sstevel@tonic-gate 
5*7c478bd9Sstevel@tonic-gate /*
6*7c478bd9Sstevel@tonic-gate  * Copyright (c) 1980 Regents of the University of California.
7*7c478bd9Sstevel@tonic-gate  * All rights reserved.  The Berkeley software License Agreement
8*7c478bd9Sstevel@tonic-gate  * specifies the terms and conditions for redistribution.
9*7c478bd9Sstevel@tonic-gate  */
10*7c478bd9Sstevel@tonic-gate /* 	Portions Copyright(c) 1988, Sun Microsystems Inc.	*/
11*7c478bd9Sstevel@tonic-gate /*	All Rights Reserved					*/
12*7c478bd9Sstevel@tonic-gate 
13*7c478bd9Sstevel@tonic-gate /*
14*7c478bd9Sstevel@tonic-gate  * Copyright (c) 1997, by Sun Microsystems, Inc.
15*7c478bd9Sstevel@tonic-gate  * All rights reserved.
16*7c478bd9Sstevel@tonic-gate  */
17*7c478bd9Sstevel@tonic-gate 
18*7c478bd9Sstevel@tonic-gate #ident	"%Z%%M%	%I%	%E% SMI"	/* SVr4.0 1.1	*/
19*7c478bd9Sstevel@tonic-gate 
20*7c478bd9Sstevel@tonic-gate /* LINTLIBRARY */
21*7c478bd9Sstevel@tonic-gate 
22*7c478bd9Sstevel@tonic-gate #include <mp.h>
23*7c478bd9Sstevel@tonic-gate #include <sys/types.h>
24*7c478bd9Sstevel@tonic-gate #include "libmp.h"
25*7c478bd9Sstevel@tonic-gate 
26*7c478bd9Sstevel@tonic-gate int
mp_msqrt(MINT * a,MINT * b,MINT * r)27*7c478bd9Sstevel@tonic-gate mp_msqrt(MINT *a, MINT *b, MINT *r)
28*7c478bd9Sstevel@tonic-gate {
29*7c478bd9Sstevel@tonic-gate 	MINT a0, x, junk, y;
30*7c478bd9Sstevel@tonic-gate 	int j;
31*7c478bd9Sstevel@tonic-gate 
32*7c478bd9Sstevel@tonic-gate 	a0.len = junk.len = y.len = 0;
33*7c478bd9Sstevel@tonic-gate 	if (a->len < 0)
34*7c478bd9Sstevel@tonic-gate 		_mp_fatal("mp_msqrt: neg arg");
35*7c478bd9Sstevel@tonic-gate 	if (a->len == 0) {
36*7c478bd9Sstevel@tonic-gate 		b->len = 0;
37*7c478bd9Sstevel@tonic-gate 		r->len = 0;
38*7c478bd9Sstevel@tonic-gate 		return (0);
39*7c478bd9Sstevel@tonic-gate 	}
40*7c478bd9Sstevel@tonic-gate 	if (a->len % 2 == 1)
41*7c478bd9Sstevel@tonic-gate 		x.len = (1 + a->len) / 2;
42*7c478bd9Sstevel@tonic-gate 	else
43*7c478bd9Sstevel@tonic-gate 		x.len = 1 + a->len / 2;
44*7c478bd9Sstevel@tonic-gate 	x.val = _mp_xalloc(x.len, "mp_msqrt");
45*7c478bd9Sstevel@tonic-gate 	for (j = 0; j < x.len; x.val[j++] = 0)
46*7c478bd9Sstevel@tonic-gate 		;
47*7c478bd9Sstevel@tonic-gate 	if (a->len % 2 == 1)
48*7c478bd9Sstevel@tonic-gate 		x.val[x.len - 1] = 0400;
49*7c478bd9Sstevel@tonic-gate 	else
50*7c478bd9Sstevel@tonic-gate 		x.val[x.len - 1] = 1;
51*7c478bd9Sstevel@tonic-gate 	_mp_move(a, &a0);
52*7c478bd9Sstevel@tonic-gate 	_mp_xfree(b);
53*7c478bd9Sstevel@tonic-gate 	_mp_xfree(r);
54*7c478bd9Sstevel@tonic-gate loop:
55*7c478bd9Sstevel@tonic-gate 	mp_mdiv(&a0, &x, &y, &junk);
56*7c478bd9Sstevel@tonic-gate 	_mp_xfree(&junk);
57*7c478bd9Sstevel@tonic-gate 	mp_madd(&x, &y, &y);
58*7c478bd9Sstevel@tonic-gate 	mp_sdiv(&y, 2, &y, (short *) &j);
59*7c478bd9Sstevel@tonic-gate 	if (mp_mcmp(&x, &y) > 0) {
60*7c478bd9Sstevel@tonic-gate 		_mp_xfree(&x);
61*7c478bd9Sstevel@tonic-gate 		_mp_move(&y, &x);
62*7c478bd9Sstevel@tonic-gate 		_mp_xfree(&y);
63*7c478bd9Sstevel@tonic-gate 		goto loop;
64*7c478bd9Sstevel@tonic-gate 	}
65*7c478bd9Sstevel@tonic-gate 	_mp_xfree(&y);
66*7c478bd9Sstevel@tonic-gate 	_mp_move(&x, b);
67*7c478bd9Sstevel@tonic-gate 	mp_mult(&x, &x, &x);
68*7c478bd9Sstevel@tonic-gate 	mp_msub(&a0, &x, r);
69*7c478bd9Sstevel@tonic-gate 	_mp_xfree(&x);
70*7c478bd9Sstevel@tonic-gate 	_mp_xfree(&a0);
71*7c478bd9Sstevel@tonic-gate 	return (r->len);
72*7c478bd9Sstevel@tonic-gate }
73