Return to msqrt.c CVS log

Up to [MW Coherent from dump] / coherent / d / usr / lib / libmp

Annotation of coherent/d/usr/lib/libmp/msqrt.c, revision 1.1

1.1     ! root        1: #include "mprec.h"
        !             2: 
        !             3: 
        !             4: /*
        !             5:  *     Msqrt sets the mint pointed to by "b" to the greatest integer
        !             6:  *     which is less than or equal to the square root ot the mint
        !             7:  *     pointed to by "a".  It sets the mint pointed to by "r" to the
        !             8:  *     remainder (ie. r = a - b * b).  If "a" is negative, then mperr
        !             9:  *     is called with the appropriate error message.  Note that no 
        !            10:  *     assumption is made as to the distinction of "a", "b" and "r"
        !            11:  *     except that "b" must not equal "r".
        !            12:  */
        !            13: 
        !            14: void
        !            15: msqrt(a, b, r)
        !            16: register mint *a, *b, *r;
        !            17: {
        !            18:        mint    x, xs, tx, rem;
        !            19: 
        !            20:        if (!ispos(a))
        !            21:                mperr("square root of negative number");
        !            22: 
        !            23:        /* initialize temps */
        !            24:        minit(&x);
        !            25:        minit(&xs);
        !            26:        minit(&rem);
        !            27:        minit(&tx);
        !            28: 
        !            29: 
        !            30:        /* get initial guess */
        !            31:        mitom(BASE, &x);
        !            32:        spow(&x, a->len / 2, &x);
        !            33: 
        !            34:        /* loop till x * x <= a < (x + 1) * (x + 1) */
        !            35:        for (;;) {
        !            36:                mult(&x, &x, &xs);
        !            37:                msub(a, &xs, &rem);
        !            38:                madd(&x, &x, &tx);
        !            39:                if (ispos(&rem) && mcmp(&rem, &tx) <= 0)
        !            40:                        break;
        !            41:                madd(a, &xs, &xs);
        !            42:                mdiv(&xs, &tx, &x, &xs);
        !            43:        }
        !            44: 
        !            45:        /* throw away garbage and return results */
        !            46:        mpfree(b->val);
        !            47:        *b = x;
        !            48:        mpfree(r->val);
        !            49:        *r = rem;
        !            50:        mpfree(xs.val);
        !            51:        mpfree(tx.val);
        !            52: }

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.