factor.c

miracl-大数运算库,大家使用有什么问题请多多提意见

                    iv=p/MULT;
                    if (p%MULT>MULT/2) iv++;
                    interval=(long)iv*MULT;
                    p=interval+1;
                    marks(interval);
                    convert(iv,t);
                    lucas(fd,t,n,fp,fvw);
                    negify(fp,fp);
                    subtract(fvw,fu[p%MULT],q);
                    btch*=100;
                    i++;
                    continue;
                }
                pa=p;
                while ((lim1/p) > pa) pa*=p;
                convert((int)pa,t);
                lucas(b,t,n,fp,q);
                copy(q,b);
                decr(q,2,q);
            }
            else
            { /* phase 2 - looking for last large PRIME FACTOR of (p+1) */
                p+=2;
                pos=p%MULT;
                if (pos>MULT/2)
                { /* increment giant step */
                    iv++;
                    interval=(long)iv*MULT;
                    p=interval+1;
                    marks(interval);
                    pos=1;
                    copy(fvw,t);
                    mad(fvw,fd,fp,n,n,fvw);
                    negify(t,fp);
                }
                if (!cp[pos]) continue;

        /* if neither interval+/-pos is prime, don't bother */
                if (!plus[pos] && !minus[pos]) continue;

                subtract(fvw,fu[pos],t);
                mad(q,t,t,n,n,q);  /* batching gcds */
            }
            if (i++%btch==0)
            { /* try for a solution */
                if (!suppress)
                {
                     printf("\b\b\b\b\b\b\b\b%8ld",p);
                     fflush(stdout);
                }
                egcd(q,n,t);
                if (size(t)==1)
                {
                    if (p>lim2) 
                    {
                        if (!suppress)
                        {
                             printf("\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b");
                             fflush(stdout);
                        }
                        break;
                    }
                    else continue;
                }
                if (compare(t,n)==0)
                {
                    if (!suppress) 
                    {
                        printf("\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b");
                        printf("degenerate case\n");
                    }
                    break;
                }
                if (!suppress)
                {
                    printf("\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b");
                    if (isprime(t)) printf("PRIME FACTOR     ");
                    else            printf("COMPOSITE FACTOR ");
                }
                else if (!isprime(t)) printf("& ");
                cotnum(t,output);
                divide(n,t,n);
                if (isprime(n)) 
                {
                    if (!suppress) printf("PRIME FACTOR     ");
                    cotnum(n,output);
                    exit(0);
                }
                nt=ntrys;
                break;
            }
        } 
        if (nt>=ntrys) break;
    }
    gprime(0);
    mr_free(b);
    mr_free(q);
    mr_free(t);
    mr_free(fp);
    mr_free(fvw);
    mr_free(fd);
    mr_free(fn);
    for (m=1;m<=MULT/2;m+=2)
        if (igcd(MULT,m)==1) mr_free(fu[m]);
}


static big ak,t,ww,s1,d1,s2,d2;

void duplication(big sum,big diff,big x,big z)
{ /* double a point on the curve P(x,z)=2.P(x1,z1) */
    nres_modmult(sum,sum,t);
    nres_modmult(diff,diff,z);
    nres_modmult(t,z,x);          /* x = sum^2.diff^2 */
    nres_modsub(t,z,t);           /* t = sum^2-diff^2 */
    nres_modmult(ak,t,ww);
    nres_modadd(z,ww,z);           /* z = ak*t +diff^2 */
    nres_modmult(z,t,z);          /* z = z.t          */
}

void addition(big xd,big zd,big sm1,big df1,big sm2,big df2,big x,big z)
{ /* add two points on the curve P(x,z)=P(x1,z1)+P(x2,z2) *
   * given their difference P(xd,zd)                      */
    nres_modmult(df2,sm1,x);
    nres_modmult(df1,sm2,z);
    nres_modadd(z,x,t);
    nres_modsub(z,x,z);
    nres_modmult(t,t,x);
    nres_modmult(x,zd,x);     /* x = zd.[df1.sm2+sm1.df2]^2 */
    nres_modmult(z,z,z);
    nres_modmult(z,xd,z);     /* z = xd.[df1.sm2-sm1.df2]^2 */
}

void ellipse(big x,big z,int r,big x1,big z1,big x2,big z2)
{ /* calculate point r.P(x,z) on curve */
    int k,rr;
    k=1;
    rr=r;
    copy(x,x1);            
    copy(z,z1);
    nres_modadd(x1,z1,s1);
    nres_modsub(x1,z1,d1);
    duplication(s1,d1,x2,z2);  /* generate 2.P */
    while ((rr/=2)>1) k*=2;
    while (k>0)
    { /* use binary method */
        nres_modadd(x1,z1,s1);         /* form sums and differences */
        nres_modsub(x1,z1,d1);    /* x+z and x-z for P1 and P2 */
        nres_modadd(x2,z2,s2);
        nres_modsub(x2,z2,d2);
        if ((r&k)==0)
        { /* double P(x1,z1) mP to 2mP */
            addition(x,z,s1,d1,s2,d2,x2,z2);
            duplication(s1,d1,x1,z1);
        }
        else
        { /* double P(x2,z2) (m+1)P to (2m+2)P */
            addition(x,z,s1,d1,s2,d2,x1,z1);
            duplication(s2,d2,x2,z2);
        }    
        k/=2;
    }
}


int kurve=5;

void lenstra(int lim1,long lim2,int ncurves)
{  /*  factoring program using Lenstras Elliptic Curve method */
    int phase,m,nc,iv,pos,btch,u,v;
    long i,p,pa,interval;
    big q,x,z,a,x1,z1,x2,z2,xt,zt,fvw;
    q=mirvar(0);
    x=mirvar(0);
    z=mirvar(0);
    a=mirvar(0);
    x1=mirvar(0);
    z1=mirvar(0);
    x2=mirvar(0);
    z2=mirvar(0);
    s1=mirvar(0);
    d1=mirvar(0);
    s2=mirvar(0);
    d2=mirvar(0);
    ak=mirvar(0);
    xt=mirvar(0);
    zt=mirvar(0);
    fvw=mirvar(0);
    t=mirvar(0);
    ww=mirvar(0);
    gprime(lim1);

    for (m=1;m<=MULT/2;m+=2)
        if (igcd(MULT,m)==1)
        {
            fu[m]=mirvar(0);
            cp[m]=TRUE;
        }
        else cp[m]=FALSE;


    prepare_monty(n);
    for (nc=1;;)
    { /* try a new curve */
        kurve++;                     /* generating an elliptic curve */
        u=kurve*kurve-5;
        v=4*kurve;

        convert(u,x);  nres(x,x);
        convert(v,z);  nres(z,z);
        nres_modsub(z,x,a);   /* a=v-u */

        copy(x,t);
        nres_modmult(x,x,x);
        nres_modmult(x,t,x);  /* x=u^3 */

        copy(z,t);
        nres_modmult(z,z,z);
        nres_modmult(z,t,z);  /* z=v^3 */

        copy(a,t);
        nres_modmult(t,t,t);
        nres_modmult(t,a,t);  /* t=(v-u)^3 */

        convert(3*u,a); nres(a,a);
        convert(v,ak);  nres(ak,ak);
        nres_modadd(a,ak,a);
        nres_modmult(t,a,t);  /* t=(v-u)^3.(3u+v) */

        convert(u,a);  nres(a,a);
        copy(a,ak);
        nres_modmult(a,a,a);
        nres_modmult(a,ak,a);   /* a=u^3 */
        convert(v,ak); nres(ak,ak);
        nres_modmult(a,ak,a);   /* a=u^3.v */
        nres_premult(a,16,a);
        nres_moddiv(t,a,ak);     /* ak=(v-u)^3.(3u+v)/16u^3v */

        phase=1;
        p=0;
        i=0;
        btch=50;
        if (!suppress) printf("curve %3d phase 1 - trying all primes less than %d\n",nc,lim1);
        if (!suppress) printf("prime= %8ld",p);
        nc++;
        forever
        { /* main loop */
            if (phase==1)
            {
                p=mip->PRIMES[i];
                if (mip->PRIMES[i+1]==0)
                { /* now change gear */
                    phase=2;
                    if (!suppress) 
                    {
                        printf("\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b");
                        printf("          phase 2 - trying last prime less than %ld\n",
                            lim2);
                        printf("prime= %8ld",p);
                    }
                    copy(x,xt);
                    copy(z,zt);
                    nres_modadd(x,z,s2);
                    nres_modsub(x,z,d2);                    /*   P = (s2,d2) */
                    duplication(s2,d2,x,z);
                    nres_modadd(x,z,s1);
                    nres_modsub(x,z,d1);                    /* 2.P = (s1,d1) */

                    nres_moddiv(x1,z1,fu[1]);               /* fu[1] = x1/z1 */
                    
                    addition(x1,z1,s1,d1,s2,d2,x2,z2); /* 3.P = (x2,z2) */
                    for (m=5;m<=MULT/2;m+=2)
                    { /* calculate m.P = (x,z) and store fu[m] = x/z */
                        nres_modadd(x2,z2,s2);
                        nres_modsub(x2,z2,d2);
                        addition(x1,z1,s2,d2,s1,d1,x,z);
                        copy(x2,x1);
                        copy(z2,z1);
                        copy(x,x2);
                        copy(z,z2);
                        if (!cp[m]) continue;
                        copy(z2,fu[m]);
                        nres_moddiv(x2,fu[m],fu[m]);
                    }
                    ellipse(xt,zt,MULT,x,z,x2,z2);
                    nres_modadd(x,z,xt);
                    nres_modsub(x,z,zt);              /* MULT.P = (xt,zt) */
                    iv=(int)(p/MULT);
                    if (p%MULT>MULT/2) iv++;
                    interval=(long)iv*MULT; 
                    p=interval+1;
                    marks(interval);
                    ellipse(x,z,iv,x1,z1,x2,z2); /* (x1,z1) = iv.MULT.P */
                    nres_moddiv(x1,z1,fvw);                /* fvw = x1/z1 */
                    nres_modsub(fvw,fu[p%MULT],q);
                    btch*=100;
                    i++;
                    continue;
                }
                pa=p;
                while ((lim1/p) > pa) pa*=p;
                ellipse(x,z,(int)pa,x1,z1,x2,z2);
                copy(x1,x);
                copy(z1,z);
                copy(z,q);
            }
            else
            { /* phase 2 - looking for last large PRIME FACTOR of (p+1+d) */
                p+=2;
                pos=(int)(p%MULT);
                if (pos>MULT/2)
                { /* increment giant step */
                    iv++;
                    interval=(long)iv*MULT;
                    p=interval+1;
                    marks(interval);
                    pos=1;
                    nres_moddiv(x2,z2,fvw);
                    nres_modadd(x2,z2,s2);
                    nres_modsub(x2,z2,d2);
                    addition(x1,z1,s2,d2,xt,zt,x,z);
                    copy(x2,x1);
                    copy(z2,z1);
                    copy(x,x2);
                    copy(z,z2);
                }
                if (!cp[pos]) continue;

        /* if neither interval +/- pos is prime, don't bother */
                if (!plus[pos] && !minus[pos]) continue;

                nres_modsub(fvw,fu[pos],t);
                nres_modmult(q,t,q);
            }
            if (i++%btch==0)
            { /* try for a solution */
                if (!suppress) 
                {
                    printf("\b\b\b\b\b\b\b\b%8ld",p);
                    fflush(stdout);
                }
                egcd(q,n,t);
                if (size(t)==1)
                {
                    if (p>lim2) 
                    {
                        if (!suppress) 
                        {
                            printf("\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b");
                            fflush(stdout);
                        }
                        break;
                    }
                    else continue;
                }
                if (compare(t,n)==0)
                {
                    if (!suppress) 
                    {
                        printf("\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b");
                        printf("degenerate case\n");
                    }
                    break;
                }

                if (!suppress) 
                {
                    printf("\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b");
                    if (isprime(t)) printf("PRIME FACTOR     ");
                    else            printf("COMPOSITE FACTOR ");
                }
                else if (!isprime(t)) printf("& ");
                cotnum(t,output);
                divide(n,t,n);
                if (isprime(n))
                {
                    if (!suppress) printf("PRIME FACTOR     ");
                    cotnum(n,output);
                    exit(0);
                }
                nc=ncurves+1;
                break;
            }
        }
        if (nc>ncurves) break;
    } 


    gprime(0);
    mr_free(ww);
    mr_free(t);
    mr_free(q);
    mr_free(x);
    mr_free(z);
    mr_free(a);
    mr_free(x1);
    mr_free(z1);
    mr_free(x2);
    mr_free(z2);
    mr_free(s1);
    mr_free(d1);
    mr_free(s2);
    mr_free(d2);
    mr_free(ak);
    mr_free(xt);
    mr_free(zt);
    mr_free(fvw);
    for (m=1;m<=MULT/2;m+=2)
        if (igcd(MULT,m)==1) mr_free(fu[m]);
}


#define SSIZE 100000    /* Maximum sieve size            */

static big NN,TT,DD,RR,VV,PP,XX,YY,DG,IG,AA,BB;
static big *x,*y,*z,*w;
static unsigned int **EE,**G;
static int *epr,*r1,*r2,*rp,*b,*pr,*e,*hash;
static unsigned char *logp,*sieve;
static int mm,mlf,jj,nbts,nlp,lp,hmod,hmod2;
static BOOL partial;
static miracl *mip;