Return to msqrt.c CVS log

Up to [Research Unix] / researchv10no / libmp

Annotation of researchv10no/libmp/msqrt.c, revision 1.1.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.