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