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 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