factor.c

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

        incr(DG,1,TT);
        subdiv(TT,4,TT);
        powmod(DD,TT,DG,BB);
        negify(DD,TT);
        mad(BB,BB,TT,DG,TT,TT);
        negify(TT,TT);
        premult(BB,2,AA);
        xgcd(AA,DG,AA,AA,AA);
        mad(AA,TT,TT,DG,DG,AA);
        multiply(AA,DG,TT);
        add(BB,TT,BB);        /* BB^2 = DD mod DG^2 */
        multiply(DG,DG,AA);   /* AA = DG*DG         */
        xgcd(DG,DD,IG,IG,IG); /* IG = 1/DG mod DD  */

        r1[0]=r2[0]=0;
        for (k=1;k<=mm;k++) 
        { /* find roots of quadratic mod each prime */
            s=subdiv(BB,epr[k],TT);
            r=subdiv(AA,epr[k],TT);
            r=invers(r,epr[k]);     /* r = 1/AA mod p */
            s1=(epr[k]-s+rp[k]);
            s2=(epr[k]-s+epr[k]-rp[k]);
            r1[k]=smul(s1,r,epr[k]);
            r2[k]=smul(s2,r,epr[k]);
        }

        for (ptr=(-NS);ptr<NS;ptr++)
        { /* sieve over next period */
            la=(long)ptr*SSIZE;
            SV=(unsigned int *)sieve;
            for (i=0;i<SSIZE/sizeof(int);i++) *SV++=0;
            for (k=1;k<=mm;k++)
            { /* sieving with each prime */
                epri=epr[k];
                logpi=logp[k];
                r=(int)(la%epr[k]);
                s1=(r1[k]-r)%epri;
                if (s1<0) s1+=epri;
                s2=(r2[k]-r)%epri;
                if (s2<0) s2+=epri;
            /* these loops are time-critical */
   
                for (j=s1;j<SSIZE;j+=epri) sieve[j]+=logpi;
                if (s1==s2) continue;
                for (j=s2;j<SSIZE;j+=epri) sieve[j]+=logpi;
            }

            for (a=0;a<SSIZE;a++)
            { /* main loop - look for fully factored residues */
                if (sieve[a]<threshold) continue;
                lptr=la+a;
                lgconv(lptr,TT);
                S=0;
                multiply(AA,TT,TT);           /* TT = AAx + BB      */
                add(TT,BB,TT);
                mad(TT,IG,TT,DD,DD,PP);       /* PP = (AAx + BB)/G  */
                if (size(PP)<0) add(PP,DD,PP);
                mad(PP,PP,PP,DD,DD,VV);       /* VV = PP^2 mod kN  */
                absol(TT,TT);
                if (compare(TT,RR)<0) S=1;    /* check for -ve VV */
                if (S==1) subtract(DD,VV,VV);
                copy(VV,TT);
                e[0]=S;
                for (k=1;k<=mm;k++) e[k]=0;
                if (!factored(lptr,TT)) continue;
                if (gotcha())
                { /* factors found! */
                    egcd(TT,NN,PP);
                    if (!suppress)
                    {
                        if (isprime(PP)) printf("PRIME FACTOR     ");
                        else             printf("COMPOSITE FACTOR ");
                    }
                    else if (!isprime(PP)) printf("& ");

                    cotnum(PP,output);
                    divide(NN,PP,NN);
                    if (!suppress)
                    {
                        if (isprime(NN)) printf("PRIME FACTOR     ");
                        else             printf("COMPOSITE FACTOR ");
                    }
                    else if (!isprime(NN)) printf("& ");

                    cotnum(NN,output);
                    exit(0);
                }
            }
        }
    }
}


/* Code to parse formula in command line
   This code isn't mine, but its public domain
   Shamefully I forget the source
 
  NOTE: It may be necessary on some platforms to change the operators * and #
*/

#define TIMES '*'
#define RAISE '#'


int digits(void)
{ /* size of n */
    int d;
    t=mirvar(0);
    d=0;
    copy(n,t);
    while (size(t)!=0)
    {
        subdiv(t,10,t);
        d++;
    }
    mr_free(t);
    return d;
}

static char *s;

void eval_power (big oldn,big n,char op)
{
        if (op) power(oldn,size(n),n,n);
}

void eval_product (big oldn,big n,char op)
{
        switch (op)
        {
        case TIMES:
                multiply(n,oldn,n);
                break;
        case '/':
                copy(oldn,t);
                divide(t,n,t);
                copy(t,n);
                break;
        case '%':
                copy(oldn,t);
                divide(t,n,n);
                copy(t,n);
        }
}

void eval_sum (big oldn,big n,char op)
{
        switch (op)
        {
        case '+':
                add(n,oldn,n);
                break;
        case '-':
                subtract(oldn,n,n);
        }
}

void eval (void)
{
        big oldn[3];
        big n;
        int i;
        char oldop[3];
        char op;
        char minus;
        n=mirvar(0);
        for (i=0;i<3;i++)
        {
            oldop[i]=0;
            oldn[i]=mirvar(0);
        }
LOOP:
        while (*s==' ')
        s++;
        if (*s=='-')    /* Unary minus */
        {
        s++;
        minus=1;
        }
        else
        minus=0;
        while (*s==' ')
        s++;
        if (*s=='(' || *s=='[' || *s=='{')    /* Number is subexpression */
        {
        s++;
        eval ();
        copy(t,n);     
        }
        else            /* Number is decimal value */
        {
        for (i=0;s[i]>='0' && s[i]<='9';i++)
                ;
        if (!i)         /* No digits found */
        {
                printf ("Error - invalid number\n");
                exit (20);
        }
        op=s[i];
        s[i]=0;
        lgconv(atol(s),n);
        s+=i;
        *s=op;
        }
        if (minus) negify(n,n);
        do
        op=*s++;
        while (op==' ');
        if (op==0 || op==')' || op==']' || op=='}')
        {
        eval_power (oldn[2],n,oldop[2]);
        eval_product (oldn[1],n,oldop[1]);
        eval_sum (oldn[0],n,oldop[0]);
        copy(n,t);
        mr_free(n);
        for (i=0;i<2;i++) mr_free(oldn[i]);
        return;
        }
        else
        {
        if (op==RAISE)
        {
                eval_power (oldn[2],n,oldop[2]);
                copy(n,oldn[2]);
                oldop[2]=RAISE;
        }
        else
        {
                if (op==TIMES || op=='/' || op=='%')
                {
                eval_power (oldn[2],n,oldop[2]);
                oldop[2]=0;
                eval_product (oldn[1],n,oldop[1]);
                copy(n,oldn[1]);
                oldop[1]=op;
                }
                else
                {
                if (op=='+' || op=='-')
                {
                        eval_power (oldn[2],n,oldop[2]);
                        oldop[2]=0;
                        eval_product (oldn[1],n,oldop[1]);
                        oldop[1]=0;
                        eval_sum (oldn[0],n,oldop[0]);
                        copy(n,oldn[0]);
                        oldop[0]=op;
                }
                else    /* Error - invalid operator */
                {
                        printf ("Error - invalid operator\n");
                        exit (20);
                }
                }
        }
        }
        goto LOOP;
}

int main(argc,argv)
int argc;
char **argv;
{
    FILE *ifile;
    int ip,b,d=250;
    argv++;argc--;
    if (argc<1)
    {
      printf("Incorrect Usage\n");
      printf("factor <number>\n");
      printf("OR\n");
      printf("factor -f <formula>\n");
      printf("e.g. factor 999999999999999999999999999999999999999999999999999997\n");
      printf("or   factor -f 10#100-19\n\n");
      printf("To suppress the commentary, use flag -s\n");
      printf("To input from a file, use flag -i <filename>\n");
      printf("To output to a file, use flag -o <filename>\n");
      printf("To set max. number size, set -dn, where n is number of decimal digits\n");
      printf("(Default is -d150). Must be first flag and before number.\n");
      printf("e.g. factor -d200 -f 10#200-1 -s -o factors.dat\n\n");
      printf("Freeware from Shamus Software, Dublin, Ireland\n");
      printf("Full C source code and MIRACL multiprecision library available\n");
      printf("Email to [e-mail] for details\n");
      printf("Web page http://indigo.ie/~mscott\n");
      printf("Source code from ftp://ftp.computing.dcu.ie/pub/crypto/miracl.zip\n");
      return 0;
    }

    b=(d*45)/100;
#ifndef MR_NOFULLWIDTH
    mip=mirsys(-b,0);
#else
    mip=mirsys(-b,MAXBASE);
#endif

    mip->NTRY=100;
    n=mirvar(0);

    ip=0;
    output=stdout;

    while (ip<argc)
    {
        if (strncmp(argv[ip],"-d",2)==0)
        {
            d=atoi(argv[ip++]+2);
            mr_free(n);
            mirexit();
            if (d<20) d=20;
            b=(d*45)/100;
#ifndef MR_NOFULLWIDTH
            mip=mirsys(-b,0);
#else
            mip=mirsys(-b,MAXBASE);
#endif
            mip->NTRY=100;
            n=mirvar(0);
            continue;

        }
        if (strcmp(argv[ip],"-f")==0)
        {
            ip++;
            s=argv[ip++];
            t=mirvar(0);
            eval();
            copy(t,n);
            mr_free(t);
            cotnum(n,stdout);
            continue;
        }
        if (strcmp(argv[ip],"-s")==0)
        {
            ip++;
            suppress=TRUE;
            continue;
        }
        if (strcmp(argv[ip],"-i")==0)
        {
            ip++;
            ifile=fopen(argv[ip++],"rt");
            if (ifile==NULL) break;
            cinnum(n,ifile);
            cotnum(n,stdout);
            continue;
        }
        if (strcmp(argv[ip],"-o")==0)
        {
            ip++;
            output=fopen(argv[ip++],"wt");
            continue;
        }
        cinstr(n,argv[ip++]);
    }

    if (size(n)==0) 
    {
        printf("No number to factor!\n");
        return 0;
    }
    if (size(n)<0)
    {
        printf("Positive numbers only!\n");
        return 0;
    }
    if (isprime(n))
    {
        cotnum(n,output);
        printf("this number is prime!\n");
        return 0;
    }

    if (!suppress) printf("first trying brute force division by small primes\n");
    brute();
    if (!suppress) printf("now trying %d iterations of brent's method\n",BTRIES);
    brent();
    fu= (big *)mr_alloc((1+MULT/2),sizeof(big));
    cp=(BOOL *)mr_alloc((1+MULT/2),sizeof(BOOL));
    plus=(BOOL *)mr_alloc((1+MULT/2),sizeof(BOOL));
    minus=(BOOL *)mr_alloc((1+MULT/2),sizeof(BOOL));

    if (digits()>25)
    {
        if (!suppress) printf("now trying william's (p+1) method\n");
        williams(10000,1000000L,1);
        if (!suppress) printf("now trying pollard's (p-1) method\n");
        pollard(100000,5000000L);
    }
    if (digits()>35)
    {
        if (!suppress) printf("now trying lenstra's method using 10 curves\n");
        lenstra(20000,2000000L,10);
        if (digits()>64) 
        {
             if (!suppress) printf("now trying 80 more curves\n");
             lenstra(20000,2000000L,80);
        }
        if (digits()>72)
        {
             if (!suppress) printf("trying 300 last curves\n");
             lenstra(50000,5000000,300);
        } 
    }  

    mr_free(minus);
    mr_free(plus);
    mr_free(cp);
    mr_free(fu);
    if (digits()<100)
    {
        if (!suppress) printf("finally - the multiple polynomial quadratic sieve - with large prime (*)\n");
        qsieve(digits());
    }
    if (!suppress) printf("I give up          \nCOMPOSITE FACTOR ");
    else printf("& ");
    cotnum(n,output);
    return 0;
}