Next Previous

Retired Document

Important: This sample code may not represent best practices for current development. The project may use deprecated symbols and illustrate technologies and techniques that are no longer recommended.

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;

Next Previous