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

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.