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
26 int
mp_msqrt(MINT * a,MINT * b,MINT * r)27 mp_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);
54 loop:
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