vfactor source/factor.c

/**************************************************************
 *
 *  vfactor.c
 *
 *  AltiVec-based factoring program.
 *  Created as extension of original factor.c project at 
 *  Next Software, Inc.
 *
 *  This package is part of ongoing research in the
 *  Advanced Computation Group, Apple Computer.
 *  
 *  c. 1999 Apple Computer, Inc.
 *  All Rights Reserved.
 *
 *
 *************************************************************/
 
/*
    Disclaimer: IMPORTANT:  This Apple software is supplied to you by Apple Computer, Inc.
                ("Apple") in consideration of your agreement to the following terms, and your
                use, installation, modification or redistribution of this Apple software
                constitutes acceptance of these terms.  If you do not agree with these terms,
                please do not use, install, modify or redistribute this Apple software.
 
                In consideration of your agreement to abide by the following terms, and subject
                to these terms, Apple grants you a personal, non-exclusive license, under AppleÕs
                copyrights in this original Apple software (the "Apple Software"), to use,
                reproduce, modify and redistribute the Apple Software, with or without
                modifications, in source and/or binary forms; provided that if you redistribute
                the Apple Software in its entirety and without modifications, you must retain
                this notice and the following text and disclaimers in all such redistributions of
                the Apple Software.  Neither the name, trademarks, service marks or logos of
                Apple Computer, Inc. may be used to endorse or promote products derived from the
                Apple Software without specific prior written permission from Apple.  Except as
                expressly stated in this notice, no other rights or licenses, express or implied,
                are granted by Apple herein, including but not limited to any patent rights that
                may be infringed by your derivative works or by other works in which the Apple
                Software may be incorporated.
 
                The Apple Software is provided by Apple on an "AS IS" basis.  APPLE MAKES NO
                WARRANTIES, EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION THE IMPLIED
                WARRANTIES OF NON-INFRINGEMENT, MERCHANTABILITY AND FITNESS FOR A PARTICULAR
                PURPOSE, REGARDING THE APPLE SOFTWARE OR ITS USE AND OPERATION ALONE OR IN
                COMBINATION WITH YOUR PRODUCTS.
 
                IN NO EVENT SHALL APPLE BE LIABLE FOR ANY SPECIAL, INDIRECT, INCIDENTAL OR
                CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
                GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
                ARISING IN ANY WAY OUT OF THE USE, REPRODUCTION, MODIFICATION AND/OR DISTRIBUTION
                OF THE APPLE SOFTWARE, HOWEVER CAUSED AND WHETHER UNDER THEORY OF CONTRACT, TORT
                (INCLUDING NEGLIGENCE), STRICT LIABILITY OR OTHERWISE, EVEN IF APPLE HAS BEEN
                ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
 
/* include files */
 
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <events.h>
 
#define DO_RANDOM_CURVESEED 0 // use random curves (vs. iterate from start seed)
#define VERBOSE 1             // print curve info and show progress
 
#ifdef _WIN32 
 
#include <process.h>
 
#endif
 
#include <string.h>
#if __VEC__
#include "vgiants.h"
#else
#include "origgiants.h"
#endif
 
 
/* definitions */
 
#define D 100
#define NUM_PRIMES 6542 /* PrimePi[2^16]. */
 
 
/* compiler options */
 
#ifdef _WIN32
#pragma warning( disable : 4127 4706 ) /* disable conditional is constant warning */
#endif
 
 
/* global variables */
 
extern giant    scratch2;
int             pr[NUM_PRIMES];
giant           xr = NULL, xs = NULL, zs = NULL, zr = NULL, xorg = NULL,
                zorg = NULL, t1 = NULL, t2 = NULL, t3 = NULL, N = NULL,
                gg = NULL, An = NULL, Ad = NULL;
giant           xb[D+2], zb[D+2], xzb[D+2];
int             modmode = 0, Q, modcount = 0;
 
 
/* function prototypes */
 
void        ell_even(giant x1, giant z1, giant x2, giant z2, giant An,
                        giant Ad, giant N);
void        ell_odd(giant x1, giant z1, giant x2, giant z2, giant exor,
                        giant zor, giant N);
void        ell_mul(giant xx, giant zz, int n, giant An, giant Ad, giant N);
int         least_2(int n);
 
#if VERBOSE
void        dot(void);
#else
#define dot()
#endif
 
int         psi_rand();
 
 
/**************************************************************
 *
 *  Functions
 *
 **************************************************************/
 
 
int         
psi_rand(
    void
)
{
    unsigned short  hi;
    unsigned short  low;
    time_t          tp;
    int             result;
 
    time(&tp);
    low = (unsigned short)rand();
    hi = (unsigned short)rand();
    result = ((hi << 16) | low) ^ ((int)tp);
 
    return (result & 0x7fffffff);
}
 
 
static void
set_random_seed(
    void
)
{
    /* Start the random number generator at a new position. */
    time_t      tp;
 
    time(&tp);
    srand((int)tp + (int)TickCount());
}
 
 
static int
isprime(
    int     odd
)
{
    int     j;
    int     p;
 
    for (j=1; ; j++)
    {
        p = pr[j];
        if (p*p > odd)
            return(1);
        if (odd % p == 0)
            return(0);
    }
}
 
 
static int
primeq(
    int             p
)
{
    register int    j=3;
 
    if ((p&1)==0)
        return ((p==2)?1:0);
    if (j>=p)
        return (1);
    while ((p%j)!=0)
    {
        j+=2;
        if (j*j>p)
            return(1);
    }
    return(0);
}
 
 
static void
s_modg(
    giant       N,
    giant       t
)
{
    ++modcount;
    switch (modmode)
    {
        case 0:
            modg(N, t);
            break;
        case -1:
            mersennemod(Q, t);
            break;
        case 1:
            fermatmod(Q, t);
            break;
    }
}
 
 
static void
reset_mod(
    giant   x,
    giant   N
)
/* Perform a divide (by the discovered factor) and switch back
   to non-Fermat-non-Mersenne (i.e. normal) mod mode. */
{
    divg(x, N);
    modmode = 0;
}
 
void
ell_even(
    giant   x1,
    giant   z1,
    giant   x2,
    giant   z2,
    giant   An,
    giant   Ad,
    giant   N
)
{
    gtog(x1, t1);
    addg(z1, t1);
    squareg(t1);
    s_modg(N, t1);
    gtog(x1, t2);
    subg(z1, t2);
    squareg(t2);
    s_modg(N, t2);
    gtog(t1, t3);
    subg(t2, t3);
    gtog(t2, x2);
    mulg(t1, x2);
    gshiftleft(2, x2);
    s_modg(N, x2);
    mulg(Ad, x2);
    s_modg(N, x2);
    mulg(Ad, t2);
    gshiftleft(2, t2);
    s_modg(N, t2);
    gtog(t3, t1);
    mulg(An, t1);
    s_modg(N, t1);
    addg(t1, t2);
    mulg(t3, t2);
    s_modg(N, t2);
    gtog(t2,z2);
}
 
 
void
ell_odd(
    giant   x1,
    giant   z1,
    giant   x2,
    giant   z2,
    giant   exor,
    giant   zor,
    giant   N
)
{
    gtog(x1, t1);
    subg(z1, t1);
    gtog(x2, t2);
    addg(z2, t2);
    mulg(t1, t2);
    s_modg(N, t2);
    gtog(x1, t1);
    addg(z1, t1);
    gtog(x2, t3);
    subg(z2, t3);
    mulg(t3, t1);
    s_modg(N, t1);
    gtog(t2, x2);
    addg(t1, x2);
    squareg(x2);
    s_modg(N, x2);
    gtog(t2, z2);
    subg(t1, z2);
    squareg(z2);
    s_modg(N, z2);
    mulg(zor, x2);
    s_modg(N, x2);
    mulg(exor, z2);
    s_modg(N, z2);
}
 
 
void
ell_mul(
    giant           xx,
    giant           zz,
    int             n,
    giant           An,
    giant           Ad,
    giant           N
)
{
    unsigned int    c = (unsigned int)0x80000000;
 
    if (n==1)
        return;
    if (n==2)
    {
        ell_even(xx, zz, xx, zz, An, Ad, N);
        return;
    }
    gtog(xx, xorg);
    gtog(zz, zorg);
    ell_even(xx, zz, xs, zs, An, Ad, N);
 
    while((c&n) == 0)
    {
        c >>= 1;
    }
 
    c>>=1;
    do
    {
        if (c&n)
        {
            ell_odd(xs, zs, xx, zz, xorg, zorg, N);
            ell_even(xs, zs, xs, zs, An, Ad, N);
        }
        else
        {
            ell_odd(xx, zz, xs, zs, xorg, zorg, N);
            ell_even(xx, zz, xx, zz, An, Ad, N);
        }
        c >>= 1;
    } while(c);
}
 
 
 
/* From R. P. Brent, priv. comm. 1996:
Let s > 5 be a pseudo-random seed (called $\sigma$ in the Tech. Report),
 
    u/v = (s^2 - 5)/(4s)
 
Then starting point is (x_1, y_1) where
 
    x_1 = (u/v)^3
and
    a = (v-u)^3(3u+v)/(4u^3 v) - 2
*/
 
static void
choose12(
    giant   x,
    giant   z,
    int     k,
    giant   An,
    giant   Ad,
    giant   N
)
{
    itog(k, zs);
    gtog(zs, xs);
    squareg(xs);
    itog(5, t2);
    subg(t2, xs);
    s_modg(N, xs);
    addg(zs, zs);
    addg(zs, zs);
    s_modg(N, zs);
    gtog(xs, x);
    squareg(x);
    s_modg(N, x);
    mulg(xs, x);
    s_modg(N, x);
    gtog(zs, z);
    squareg(z);
    s_modg(N, z);
    mulg(zs, z);
    s_modg(N, z);
 
    /* Now for A. */
    gtog(zs, t2);
    subg(xs, t2);
    gtog(t2, t3);
    squareg(t2);
    s_modg(N, t2);
    mulg(t3, t2);
    s_modg(N, t2);  /* (v-u)^3. */
    gtog(xs, t3);
    addg(t3, t3);
    addg(xs, t3);
    addg(zs, t3);
    s_modg(N, t3);
    mulg(t3, t2);
    s_modg(N, t2);  /* (v-u)^3 (3u+v). */
    gtog(zs, t3);
    mulg(xs, t3);
    s_modg(N, t3);
    squareg(xs);
    s_modg(N, xs);
    mulg(xs, t3);
    s_modg(N, t3);
    addg(t3, t3);
    addg(t3, t3);
    s_modg(N, t3);
    gtog(t3, Ad);
    gtog(t2, An);  /* An/Ad is now A + 2. */
}
 
 
static void
ensure(
    int q
)
{
    int     nsh, j;
 
    N = newgiant(0);
    if(!q)
    {
        gread(N,stdin);
        q = bitlen(N) + 1;
    }
    nsh = q/4; /* Allowing (easily) enough space per giant,
                    since N is about 2^q, which is q bits, or
                   q/16 shorts.  But squaring, etc. is allowed,
                    so we need at least q/8, and we choose q/4
                    to be conservative. */
    if (!xr)
        xr = newgiant(nsh);
    if (!zr)
        zr = newgiant(nsh);
    if (!xs)
        xs = newgiant(nsh);
    if (!zs)
        zs = newgiant(nsh);
    if (!xorg)
        xorg = newgiant(nsh);
    if (!zorg)
        zorg = newgiant(nsh);
    if (!t1)
        t1 = newgiant(nsh);
    if (!t2)
        t2 = newgiant(nsh);
    if (!t3)
        t3 = newgiant(nsh);
    if (!gg)
        gg = newgiant(nsh);
    if (!An)
        An = newgiant(nsh);
    if (!Ad)
        Ad = newgiant(nsh);
    for (j=0;j<D+2;j++)
    {
        xb[j] = newgiant(nsh);
        zb[j] = newgiant(nsh);
        xzb[j] = newgiant(nsh);
    }
}
 
static int
bigprimeq(
    giant   x
)
{
    itog(1, t1);
    gtog(x, t2);
    subg(t1, t2);
    itog(5, t1);
    
    powermodg(t1, t2, x);
    if (isone(t1))
        return(1);
    return(0);
}
 
#if VERBOSE
void
dot(
    void
)
{
    printf(".");
    fflush(stdout);
}
#endif
/**************************************************************
 *
 *  Main Function
 *
 **************************************************************/
 
int main(
    int argc,
    char    *argv[]
)
{
    int     j, k, C, nshorts, cnt, count,
            limitbits = 0, pass, npr, rem;
    long    B;
    int     randmode = 0;
    long    cntseed;
 
 
    if (!strcmp(argv[argc-1], "-r"))
    {
        randmode = 1;
        if (argc > 4)
          /* This segment only takes effect in random mode. */
            limitbits = atoi(argv[argc-2]);
    }
    else
    {
        randmode = 0;
    }
 
    modmode = 0;
    if (argc > 2)
    {
        modmode = atoi(argv[1]);
        Q = atoi(argv[2]);
    }
    if (modmode==0)
        Q = 0;
    ensure(Q);
    if (modmode)
    {
        itog(1, N);
        gshiftleft(Q, N);
        itog(modmode, t1);
        addg(t1, N);
    }
    pr[0] = 2;
    for (k=0, npr=1;; k++)
    {
        if (primeq(3+2*k))
        {
            pr[npr++] = 3+2*k;
            if (npr >= NUM_PRIMES)
                break;
        }
    }
 
 
    if (randmode == 0)
    {
        printf("Sieving...\n");
        fflush(stdout);
        for (j=0; j < NUM_PRIMES; j++)
        {
            gtog(N, t1);
            rem = idivg(pr[j], t1);
            if (rem == 0)
            {
                printf("%d ", pr[j]);
                gtog(t1, N);
                if (isone(N))
                {
                    printf("\n");
                    goto exit;
                }
                else
                {
                    printf("* ");
                    fflush(stdout);
                }
                --j;
            }
        }
 
 
        if (bigprimeq(N))
        {
            gout(N);
            goto exit;
        }
 
        printf("\n");
        printf("Commencing Pollard rho...\n");
        fflush(stdout);
        itog(1, gg);
        itog(3, t1); itog(3, t2);
 
        for (j=0; j < 15000; j++)
        {
            if((j%100) == 0)
            {
                dot();
                gcdg(N, gg);
                if (!isone(gg))
                    break;
            }
            squareg(t1);
            iaddg(2, t1);
            s_modg(N, t1);
            squareg(t2);
            iaddg(2, t2);
            s_modg(N, t2);
            squareg(t2);
            iaddg(2, t2);
            s_modg(N, t2);
            gtog(t2, t3);
            subg(t1, t3);
            absg(t3);
            mulg(t3, gg);
            s_modg(N, gg);
        }
        gcdg(N, gg);
 
        if ((gcompg(N,gg) != 0) && (!isone(gg)))
        {
            fprintf(stdout,"\n");
            gout(gg);
            reset_mod(gg, N);
            if (isone(N))
            {
                printf("\n");
                goto exit;
            }
            else
            {
                printf("* ");
                fflush(stdout);
            }
            if (bigprimeq(N))
            {
                gout(N);
                goto exit;
            }
        }
 
        printf("\n");
        printf("Commencing Pollard (p-1)...\n");
        fflush(stdout);
        itog(1, gg);
        itog(3, t1);
        for (j=0; j< NUM_PRIMES; j++)
        {
            cnt = (int)(8*log(2.0)/log(1.0*pr[j]));
            if (cnt < 2)
                cnt = 1;
            for (k=0; k< cnt; k++)
            {
                powermod(t1, pr[j], N);
            }
            itog(1, t2);
            subg(t1, t2);
            mulg(t2, gg);
            s_modg(N, gg);
 
            if (j % 100 == 0)
            {
                dot();
                gcdg(N, gg);
                if  (!isone(gg))
                    break;
            }
        }
        gcdg(N, gg);
        if ((gcompg(N,gg) != 0) && (!isone(gg)))
        {
            fprintf(stdout,"\n");
            gout(gg);
            reset_mod(gg, N);
            if (isone(N))
            {
                printf("\n");
                goto exit;
            }
            else
            {
                printf("* ");
                fflush(stdout);
            }
            if (bigprimeq(N))
            {
                gout(N);
                goto exit;
            }
        }
    } /* This is the end of (randmode == 0) */
 
    printf("\n");
    printf("Commencing ECM...\n");
    fflush(stdout);
 
    if (randmode)
        set_random_seed();
    pass = 0;
    
    
#if !DO_RANDOM_CURVESEED
    cntseed = 1438858426;       // initial curve seed point
#endif
 
    
    while (++pass)
    {
        if (randmode == 0)
        {
            if (pass <= 3)
            {
                B = 1000;
            }
            else if (pass <= 10)
            {
                B = 10000;
            }
            else if (pass <= 100)
            {
                B = 100000L;
            } else
            {
                B = 1000000L;
            }
        }
        else
        {
            B = 2000000L;
        }
        C = 50*((int)B);
 
        /* Next, choose curve with order divisible by 16 and choose
         *  a point (xr/zr) on said curve.
         */
 
        /* Order-div-12 case. 
         * cnt = 8020345;   Brent's parameter for stage one discovery
         * of 27-digit factor of F_13.
         */
#if DO_RANDOM_CURVESEED
        cnt = psi_rand(); /* cnt = 8020345; */
#else
        cnt = cntseed++;
#endif
        choose12(xr, zr, cnt, An, Ad, N);
#if VERBOSE      
        printf("Choosing curve %d, with s = %d, B = %d, C = %d:\n", pass,cnt, B, C);   fflush(stdout);
#endif      
        cnt = 0;
        nshorts = 1;
        count = 0;
        for (j=0;j<nshorts;j++)
        {
            ell_mul(xr, zr, 1<<16, An, Ad, N);
            ell_mul(xr, zr, 3*3*3*3*3*3*3*3*3*3*3, An, Ad, N);
            ell_mul(xr, zr, 5*5*5*5*5*5*5, An, Ad, N);
            ell_mul(xr, zr, 7*7*7*7*7*7, An, Ad, N);
            ell_mul(xr, zr, 11*11*11*11, An, Ad, N);
            ell_mul(xr, zr, 13*13*13*13, An, Ad, N);
            ell_mul(xr, zr, 17*17*17, An, Ad, N);
        }
        k = 19;
        while (k<B)
        {
            if (isprime(k))
            {
                ell_mul(xr, zr, k, An, Ad, N);
                if (k<100)
                    ell_mul(xr, zr, k, An, Ad, N);
                if (cnt++ %100==0)
                    dot();
            }
            k += 2;
        }
        count = 0;
 
        gtog(zr, gg);
        gcdg(N, gg);
        if ((!isone(gg))&&(bitlen(gg)>limitbits))
        {
            fprintf(stdout,"\n");
            gwriteln(gg, stdout);
            fflush(stdout);
            reset_mod(gg, N);
            if (isone(N))
            {
                printf("\n");
                goto exit;
            }
            else
            {
                printf("* ");
                fflush(stdout);
            }
            if (bigprimeq(N))
            {
                 gout(N);
                 goto exit;
            }
            continue;
        }
        else
        {
#if VERBOSE
            printf("\n");           
            fflush(stdout);
#endif          
        }
 
        /* Continue;  Invoke, to test Stage 1 only. */
        k = ((int)B)/D;
        gtog(xr, xb[0]);
        gtog(zr, zb[0]);
        ell_mul(xb[0], zb[0], k*D+1 , An, Ad, N);
        gtog(xr, xb[D+1]);
        gtog(zr, zb[D+1]);
        ell_mul(xb[D+1], zb[D+1], (k+2)*D+1 , An, Ad, N);
 
        for (j=1; j <= D; j++)
        {
            gtog(xr, xb[j]);
            gtog(zr, zb[j]);
            ell_mul(xb[j], zb[j], 2*j , An, Ad, N);
            gtog(zb[j], xzb[j]);
            mulg(xb[j], xzb[j]);
            s_modg(N, xzb[j]);
        }
        modcount = 0;
#if VERBOSE
        printf("\nCommencing second stage, curve %d...\n",pass); fflush(stdout);
#endif      
        count = 0;
        itog(1, gg);
 
        while (1)
        {
            gtog(zb[0], xzb[0]);
            mulg(xb[0], xzb[0]);
            s_modg(N, xzb[0]);
            mulg(zb[0], gg);
            s_modg(N,gg); /* Accumulate. */
            for (j = 1; j < D; j++)
            {
                if (!isprime(k*D+1+ 2*j))
                    continue;
 
                /* Next, accumulate (xa - xb)(za + zb) - xa za + xb zb. */
                gtog(xb[0], t1);
                subg(xb[j], t1);
                gtog(zb[0], t2);
                addg(zb[j], t2);
                mulg(t1, t2);
                s_modg(N, t1);
                subg(xzb[0], t2);
                addg(xzb[j], t2);
                s_modg(N, t2);
                --modcount;
                mulg(t2, gg);
                s_modg(N, gg);
                if((++count)%1000==0)
                    dot();
            }
 
            k += 2;
            if(k*D > C)
                break;
            gtog(xb[D+1], xs);
            gtog(zb[D+1], zs);
            ell_odd(xb[D], zb[D], xb[D+1], zb[D+1], xb[0], zb[0], N);
            gtog(xs, xb[0]);
            gtog(zs, zb[0]);
        }
 
        gcdg(N, gg);
        if((!isone(gg))&&(bitlen(gg)>limitbits))
        {
            fprintf(stdout,"\n");
            gwriteln(gg, stdout);
            fflush(stdout);
            reset_mod(gg, N);
            if (isone(N))
            {
                printf("\n");
                goto exit;
            }
            else
            {
                printf("* ");
                fflush(stdout);
            }
            if (bigprimeq(N))
            {
                gout(N);
                goto exit;
            }
            continue;
        }
 
#if VERBOSE
        printf("\n");
        fflush(stdout);
#endif      
    }
 
exit:
 
    return 0;
}


Copyright © 2003 Apple Computer, Inc. All Rights Reserved. Terms of Use | Privacy Policy | Updated: 2003-01-14