Return to msqrt.c CVS log

Up to [Research Unix] / researchv10no / libmp

Annotation of researchv10no/libmp/msqrt.c, revision 1.1

1.1     ! root        1: #include "mp.h"
        !             2: msqrt(a,b,r) mint *a,*b,*r;
        !             3: {      mint x,junk,y;
        !             4:        int j;
        !             5:        x.len=junk.len=y.len=0;
        !             6:        if(a->len<0) fatal("msqrt: neg arg");
        !             7:        if(a->len==0)
        !             8:        {       b->len=0;
        !             9:                r->len=0;
        !            10:                return(0);
        !            11:        }
        !            12:        if(a->len%2==1) x.len=(1+a->len)/2;
        !            13:        else x.len=1+a->len/2;
        !            14:        x.val=xalloc(x.len,"msqrt");
        !            15:        for(j=0;j<x.len;x.val[j++]=0);
        !            16:        if(a->len%2==1) x.val[x.len-1]=0400;
        !            17:        else x.val[x.len-1]=1;
        !            18:        xfree(b);
        !            19:        xfree(r);
        !            20: loop:
        !            21:        mdiv(a,&x,&y,&junk);
        !            22:        xfree(&junk);
        !            23:        madd(&x,&y,&y);
        !            24:        sdiv(&y,2,&y,(short *)&j);
        !            25:        if(mcmp(&x,&y)>0)
        !            26:        {       xfree(&x);
        !            27:                move(&y,&x);
        !            28:                xfree(&y);
        !            29:                goto loop;
        !            30:        }
        !            31:        xfree(&y);
        !            32:        move(&x,b);
        !            33:        mult(&x,&x,&x);
        !            34:        msub(a,&x,r);
        !            35:        xfree(&x);
        !            36:        return(r->len);
        !            37: }
        !            38: 
        !            39: /* pathology: n<= 0 => r=rem=0, num <= 0, r=rem=0 */
        !            40: /* this is a lazy implementation, assuming n>=3 => root is a legal double */
        !            41: mroot(n, num, r, rem)
        !            42: mint *num, *r, *rem;
        !            43: {      extern double log(), exp();
        !            44:        double z;
        !            45:        int i;
        !            46:        static mint *xn, *xn1, *top, *bot;
        !            47:        static init;
        !            48:        if(!init++) {
        !            49:                xn = itom(0), xn1 = itom(0), top = itom(0), bot = itom(0);
        !            50:        }
        !            51:        if(n < 0 || num->len <= 0) {
        !            52:                msub(r, r, r);
        !            53:                move(r, rem);
        !            54:                return;
        !            55:        }
        !            56:        if(n == 1) {
        !            57:                move(num, r);
        !            58:                msub(rem, rem, rem);
        !            59:                return;
        !            60:        }
        !            61:        if(n == 2) {
        !            62:                msqrt(num, r, rem);
        !            63:                return;
        !            64:        }
        !            65:        /* compute approx bigger than root */
        !            66:        for(z = 0, i = 0; i < num->len; i++)
        !            67:                z = z/32768 + num->val[i];
        !            68:        /* num = z * 2^15*(len-1) */
        !            69:        z = exp((log(z) + 15*(num->len-1)*log(2.))/n);
        !            70:        z = 1.01 * z + 1;       /* make sure it is bigger than root */
        !            71:        dtom(z, r);
        !            72:        /* ((n-1)*x^n+num)/(n*x^(n-1))*/
        !            73:        for(;;) {
        !            74:                rpow(r, n - 1, xn1);
        !            75:                mult(r, xn1, xn);
        !            76:                imult(n-1, xn, top);
        !            77:                madd(top, num, top);
        !            78:                imult(n, xn1, bot);
        !            79:                mdiv(top, bot, xn1, rem);
        !            80:                switch(mcmp(xn1, r)) {
        !            81:                case -1:
        !            82:                        move(xn1, r);
        !            83:                        continue;
        !            84:                case 0:
        !            85:                        move(xn1, r);
        !            86:                        msub(num, xn, rem);
        !            87:                        return;
        !            88:                case 1: /* since r was too large to start with */
        !            89:                        msub(num, xn, rem);
        !            90:                        return;
        !            91:                }
        !            92:        }
        !            93: }
        !            94: 
        !            95: static
        !            96: imult(n, a, b)
        !            97: mint *a, *b;
        !            98: {      mint *x;
        !            99:        x = itom(n);
        !           100:        mult(x, a, b);
        !           101:        msub(x, x,x);
        !           102: }

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