/*      Copyright (c) 1984, 1986, 1987, 1988, 1989 AT&T */
/*        All Rights Reserved   */


/*
 * Copyright (c) 1980 Regents of the University of California.
 * All rights reserved.  The Berkeley software License Agreement
 * specifies the terms and conditions for redistribution.
 */
/*      Portions Copyright(c) 1988, Sun Microsystems Inc.       */
/*      All Rights Reserved                                     */

/*
 * Copyright (c) 1997, by Sun Microsystems, Inc.
 * All rights reserved.
 */

#ident  "%Z%%M% %I%     %E% SMI"        /* SVr4.0 1.1   */

/* LINTLIBRARY */

#include <mp.h>
#include <sys/types.h>
#include "libmp.h"

int
mp_msqrt(MINT *a, MINT *b, MINT *r)
{
        MINT a0, x, junk, y;
        int j;

        a0.len = junk.len = y.len = 0;
        if (a->len < 0)
                _mp_fatal("mp_msqrt: neg arg");
        if (a->len == 0) {
                b->len = 0;
                r->len = 0;
                return (0);
        }
        if (a->len % 2 == 1)
                x.len = (1 + a->len) / 2;
        else
                x.len = 1 + a->len / 2;
        x.val = _mp_xalloc(x.len, "mp_msqrt");
        for (j = 0; j < x.len; x.val[j++] = 0)
                ;
        if (a->len % 2 == 1)
                x.val[x.len - 1] = 0400;
        else
                x.val[x.len - 1] = 1;
        _mp_move(a, &a0);
        _mp_xfree(b);
        _mp_xfree(r);
loop:
        mp_mdiv(&a0, &x, &y, &junk);
        _mp_xfree(&junk);
        mp_madd(&x, &y, &y);
        mp_sdiv(&y, 2, &y, (short *) &j);
        if (mp_mcmp(&x, &y) > 0) {
                _mp_xfree(&x);
                _mp_move(&y, &x);
                _mp_xfree(&y);
                goto loop;
        }
        _mp_xfree(&y);
        _mp_move(&x, b);
        mp_mult(&x, &x, &x);
        mp_msub(&a0, &x, r);
        _mp_xfree(&x);
        _mp_xfree(&a0);
        return (r->len);
}