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