factor.c

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

/*
 *   Program to factor Integers
 *
 *   Current parameters assume a large 32-bit address space is available.
 *
 *   This program is cobbled together from the implementations of various
 *   factoring algorithms better described in:-
 *
 *   brute.c/cpp     - brute force division by small primes
 *   brent.c/cpp     - Pollard Rho method as improved by Brent
 *   pollard.c/cpp   - Pollard (p-1) method
 *   williams.c/cpp  - Williams (p+1) method
 *   lenstra.c/cpp   - Lenstra's Elliptic Curve method
 *   qsieve.c/cpp    - The Multiple polynomial quadratic sieve
 *
 *   Note that the .cpp C++ implementations are easier to follow
 *
 *   NOTE: The quadratic sieve program requires a lot of memory for 
 *   bigger numbers. It may fail if you system cannot provide the memory
 *   requested.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "miracl.h"

#define LIMIT 15000
#define BTRIES 1000
#define MULT 2310      /* 2*3*5*7*11    */
#define NEXT 13        /* .. next prime */
#define mr_min(a,b) ((a) < (b)? (a) : (b))

static big *fu;
static BOOL *cp,*plus,*minus;
big n;
FILE *output;
static BOOL suppress=FALSE;
static int PADDING;
static miracl *mip;

void brute(void)
{ /* find factors by brute force division */
    big x,y;
    int m,p;
    gprime(LIMIT);
    x=mirvar(0);
    y=mirvar(0);
    m=0;
    p=mip->PRIMES[0];
    forever
    { /* try division by each prime in turn */
        if (subdiv(n,p,y)==0)
        { /* factor found */
            copy(y,n);
            if (!suppress) printf("PRIME FACTOR     ");
            fprintf(output,"%d\n",p);
            if (size(n)==1) exit(0);
            continue;
        }
        if (size(y)<=p) 
        { /* must be prime */
            if (!suppress) printf("PRIME FACTOR     ");
            cotnum(n,output);
            exit(0);
        }
        p=mip->PRIMES[++m];
        if (p==0) break;
    }
    if (isprime(n)) 
    {
        if (!suppress) printf("PRIME FACTOR     ");
        cotnum(n,output);
        exit(0);
    }
    mr_free(x);
    mr_free(y);
    gprime(0);
}

void brent(void)
{  /*  factoring program using Brents method */
    long k,r,i,m,iter;
    big x,y,ys,z,q,c3,t;
    x=mirvar(0);
    y=mirvar(0);
    ys=mirvar(0);
    z=mirvar(0);
    q=mirvar(0);
    c3=mirvar(3);
    t=mirvar(0);
    m=10;
    r=1;
    iter=0;
    do
    {
        if (!suppress) printf("iterations=%5ld",iter);
        convert(1,q);
        do
        {
            copy(y,x);
            for (i=1;i<=r;i++)
                mad(y,y,c3,n,n,y);
            k=0;
            do
            {
                iter++;
                if (iter>BTRIES)
                {
                    if (!suppress) 
                    {
                         printf("\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b");
                         fflush(stdout);
                    }
                    mr_free(t);
                    mr_free(x);
                    mr_free(y);
                    mr_free(ys);
                    mr_free(z);
                    mr_free(q);
                    mr_free(c3);             
                    return;
                }
                if (iter%10==0) if (!suppress) 
                {
                    printf("\b\b\b\b\b%5ld",iter);
                    fflush(stdout);
                }
                copy(y,ys);
                for (i=1;i<=mr_min(m,r-k);i++)
                {
                    mad(y,y,c3,n,n,y);
                    subtract(y,x,z);
                    mad(z,q,q,n,n,q);
                }
                egcd(q,n,z);
                k+=m;
            } while (k<r && size(z)==1);
            r*=2;
        } while (size(z)==1);
        if (compare(z,n)==0) do 
        { /* back-track */
            mad(ys,ys,c3,n,n,ys);
            subtract(ys,x,z);
        } while (egcd(z,n,z)==1);
        if (!suppress) 
        {
            printf("\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b");
            fflush(stdout);
        }
        forever
        {
            if (isprime(z))
            {
                if (!suppress) printf("PRIME FACTOR     ");
                cotnum(z,output);
            }
            else
            { /* if end of factorisation - pass it on... */
                if (compare(z,n)==0) return;
              /* you will have to come back to it */
                if (!suppress) printf("COMPOSITE FACTOR ");
                else printf("& ");
                cotnum(z,output);
            }
       
            divide(n,z,n);
            divide(y,n,n);

            copy(n,t);
            divide(t,z,z);
            if (size(t)==0) continue;
            break;

        }
        if (size(n)==1) exit(0);
    } while (!isprime(n));      
    if (!suppress) printf("PRIME FACTOR     ");
    cotnum(n,output);
    exit(0);
}

void marks(long start)
{ /* mark non-primes in this interval. Note    *
   * that those < NEXT are dealt with already  */
    int i,pr,j,k;
    for (j=1;j<=MULT/2;j+=2) plus[j]=minus[j]=TRUE;
    for (i=0;;i++)
    { /* mark in both directions */
        pr=mip->PRIMES[i];
        if (pr<NEXT) continue;
        if ((long)pr*pr>start) break;
        k=pr-start%pr;
        for (j=k;j<=MULT/2;j+=pr)
            plus[j]=FALSE;
        k=start%pr;
        for (j=k;j<=MULT/2;j+=pr)
            minus[j]=FALSE;
    }        
}

void pollard(int lim1,long lim2)
{ /* factoring program using Pollards (p-1) method */
    long i,p,pa,interval;
    int phase,m,pos,btch,iv;
    big t,b,bw,bvw,bd,bp,q;
    t=mirvar(0);
    b=mirvar(0);
    q=mirvar(0);
    bw=mirvar(0);
    bvw=mirvar(0);
    bd=mirvar(0);
    bp=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;
    phase=1;
    p=0;
    btch=50;
    i=0;
    convert(3,b);
    if (!suppress) printf("phase 1 - trying all primes less than %d\n",lim1);
    if (!suppress) printf("prime= %8ld",p);
    forever
    {
        if (phase==1)
        { /* looking for all factors of (p-1) < lim1 */
            p=mip->PRIMES[i];
            if (mip->PRIMES[i+1]==0)
            {
                 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);
                 }
                 power(b,8,n,bw);
                 convert(1,t);
                 copy(b,bp);
                 copy(b,fu[1]);
                 for (m=3;m<=MULT/2;m+=2)
                 { /* store fu[m] = b^(m*m) */
                     mad(t,bw,bw,n,n,t);
                     mad(bp,t,t,n,n,bp);
                     if (cp[m]) copy(bp,fu[m]);
                 }
                 power(b,MULT,n,t);
                 power(t,MULT,n,t);
                 mad(t,t,t,n,n,bd);        /* bd=b^(2*MULT*MULT) */
                 iv=p/MULT;
                 if (p%MULT>MULT/2) iv++;
                 interval=(long)iv*MULT;
                 p=interval+1;
                 marks(interval);
                 power(t,2*iv-1,n,bw);
                 power(t,iv,n,bvw);
                 power(bvw,iv,n,bvw);      /* bvw = b^(MULT*MULT*iv*iv) */
                 subtract(bvw,fu[p%MULT],q);
                 btch*=100;
                 i++;
                 continue;
            }
            pa=p;
            while ((lim1/p) > pa) pa*=p;
            power(b,(int)pa,n,b);
            decr(b,1,q);
        }
        else
        { /* looking for last PRIME FACTOR of (p-1) < lim2 */
            p+=2;
            pos=p%MULT;
            if (pos>MULT/2) 
            { /* increment giant step */
                iv++;
                interval=(long)iv*MULT;
                p=interval+1;
                marks(interval);
                pos=1;
                mad(bw,bd,bd,n,n,bw);
                mad(bvw,bw,bw,n,n,bvw);
            }
            if (!cp[pos]) continue;

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

            subtract(bvw,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);
            }
            break;
        }
    }
    gprime(0);
    mr_free(t);
    mr_free(b);
    mr_free(q);
    mr_free(bw);
    mr_free(bvw);
    mr_free(bd);
    mr_free(bp);
    for (m=1;m<=MULT/2;m+=2)
        if (igcd(MULT,m)==1) mr_free(fu[m]);
}

void williams(int lim1,long lim2,int ntrys)
{  /*  factoring program using Williams (p+1) method */
    int k,phase,m,nt,iv,pos,btch;
    long i,p,pa,interval;
    big b,q,fp,fvw,fd,fn,t;
    b=mirvar(0);
    q=mirvar(0);
    t=mirvar(0);
    fp=mirvar(0);
    fvw=mirvar(0);
    fd=mirvar(0);
    fn=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;
    for (nt=0,k=3;k<10;k++)
    { /* try more than once for p+1 condition (may be p-1) */
        convert(k,b);              /* try b=3,4,5..        */
        convert((k*k-4),t);
        if (egcd(t,n,t)!=1) continue; /* check (b*b-4,n)!=0 */
        nt++;
        phase=1;
        p=0;
        btch=50;
        i=0;
        if (!suppress) printf("phase 1 - trying all primes less than %d\n",lim1);
        if (!suppress) printf("prime= %8ld",p);
        forever
        { /* main loop */
            if (phase==1)
            { /* looking for all factors of p+1 < lim1 */
                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(b,fu[1]);
                    copy(b,fp);
                    mad(b,b,b,n,n,fd);
                    decr(fd,2,fd);
                    negify(b,t);
                    mad(fd,b,t,n,n,fn);
                    for (m=5;m<=MULT/2;m+=2)
                    { /* store fu[m] = Vm(b) */
                        negify(fp,t);
                        mad(fn,fd,t,n,n,t);
                        copy(fn,fp);
                        copy(t,fn);
                        if (!cp[m]) continue;
                        copy(t,fu[m]);
                    }
                    convert(MULT,t);
                    lucas(b,t,n,fp,fd);