2N/A/* Copyright (c) 1984, 1986, 1987, 1988, 1989 AT&T */
2N/A/* All Rights Reserved */
2N/A
2N/A
2N/A/*
2N/A * Copyright (c) 1980 Regents of the University of California.
2N/A * All rights reserved. The Berkeley software License Agreement
2N/A * specifies the terms and conditions for redistribution.
2N/A */
2N/A/* Portions Copyright(c) 1988, Sun Microsystems Inc. */
2N/A/* All Rights Reserved */
2N/A
2N/A/*
2N/A * Copyright (c) 1997, by Sun Microsystems, Inc.
2N/A * All rights reserved.
2N/A */
2N/A
2N/A#ident "%Z%%M% %I% %E% SMI" /* SVr4.0 1.1 */
2N/A
2N/A/* LINTLIBRARY */
2N/A
2N/A#include <mp.h>
2N/A#include <sys/types.h>
2N/A#include "libmp.h"
2N/A
2N/Aint
2N/Amp_msqrt(MINT *a, MINT *b, MINT *r)
2N/A{
2N/A MINT a0, x, junk, y;
2N/A int j;
2N/A
2N/A a0.len = junk.len = y.len = 0;
2N/A if (a->len < 0)
2N/A _mp_fatal("mp_msqrt: neg arg");
2N/A if (a->len == 0) {
2N/A b->len = 0;
2N/A r->len = 0;
2N/A return (0);
2N/A }
2N/A if (a->len % 2 == 1)
2N/A x.len = (1 + a->len) / 2;
2N/A else
2N/A x.len = 1 + a->len / 2;
2N/A x.val = _mp_xalloc(x.len, "mp_msqrt");
2N/A for (j = 0; j < x.len; x.val[j++] = 0)
2N/A ;
2N/A if (a->len % 2 == 1)
2N/A x.val[x.len - 1] = 0400;
2N/A else
2N/A x.val[x.len - 1] = 1;
2N/A _mp_move(a, &a0);
2N/A _mp_xfree(b);
2N/A _mp_xfree(r);
2N/Aloop:
2N/A mp_mdiv(&a0, &x, &y, &junk);
2N/A _mp_xfree(&junk);
2N/A mp_madd(&x, &y, &y);
2N/A mp_sdiv(&y, 2, &y, (short *) &j);
2N/A if (mp_mcmp(&x, &y) > 0) {
2N/A _mp_xfree(&x);
2N/A _mp_move(&y, &x);
2N/A _mp_xfree(&y);
2N/A goto loop;
2N/A }
2N/A _mp_xfree(&y);
2N/A _mp_move(&x, b);
2N/A mp_mult(&x, &x, &x);
2N/A mp_msub(&a0, &x, r);
2N/A _mp_xfree(&x);
2N/A _mp_xfree(&a0);
2N/A return (r->len);
2N/A}