/* > ecdl2K-108.c * Purpose: Fast arithmetic for computing discrete logs on elliptic curves. * Copyright: Robert J. Harley, 1997-1999. * Contact: Robert.Harley@inria.fr * Legalese: This source code is subject to the GNU General Public Licence v2. */ /* This is Rob's 64-bit ECDL program for Unix, version 1.1.1. * * Consult the following URL for more information: * http://cristal.inria.fr/~harley/ecdl7/ */ /*** Quickstart ***/ /* Do this: cc -O2 ecdl2K-108.c -o ecdl.exe echo Machine name is `hostname` mkdir `hostname` cd `hostname` nice ../ecdl.exe mail by my@email.address on `hostname` > ecdl.log & * * Then go sign up, using the exact same email address, at this form: * * http://cristal.inria.fr/bin/ecdldb * */ /*** Detailed instructions ***/ /** 1. Signing up ** * * Go to this form to sign up: * * http://cristal.inria.fr/bin/ecdldb * * Give your email address (or at least something that identifies you * uniquely) and a passphrase you just invented. Keep a safe copy of * the passphrase! * * Note: If you prefer to stay anonymous and not sign up, that's fine. * **/ /** 2. Compiling ** * * Compile with something like: * gcc -O2 -freg-struct-return ecdl2K-108.c -o ecdl.exe * or: * cc -O5 ecdl2K-108.c -o ecdl * * You can test the binary like this: * ./ecdl.exe test by me@here on bla * * The test output should match this sample except for the reported * times and rates: > ./ecdl.exe test by me@here on bla This is Rob's ECDL program for Unix64, version 1.1.1 Mode ...... test Email ..... me@here Machine ... bla Starting at Sun Dec 12 15:34:13 1999 Generating 32 new starting points. Computing iterations... Test point found at: Fri Nov 26 17:49:14 1999 TEST-L0|Unix64|1.1.1| TEST-L1|EF7DC6C12ACAF01D48CE4D44EB6|5163373649C868AD3265FACCF99| TEST-L2|me@here| TEST-L3|000000004887| TEST-L4|bla| Iterations used = 18567. Total iterations = 594144. Rate = 204368 per second. Test point found at: Fri Nov 26 17:49:28 1999 TEST-L0|Unix64|1.1.1| TEST-L1|B0152E85E4E22B6EE3C7120FC46|026063744AFC9DB003098294BD6| TEST-L2|me@here| TEST-L3|00000001BACA| TEST-L4|bla| Iterations used = 113354. Total iterations = 3627328. Rate = 205715 per second. Test point found at: Fri Nov 26 17:49:33 1999 TEST-L0|Unix64|1.1.1| TEST-L1|607BB600DECAF410A7B62ED3DD3|F43D17AAF642BBF4EBFAA9C1984| TEST-L2|me@here| TEST-L3|000000023C69| TEST-L4|bla| Iterations used = 146537. Total iterations = 4689184. Rate = 205606 per second. [...and so on...] * If the program complains that it cannot find sendmail, you can * compile with -DSENDMAIL='"/usr/sbin/sendmail"' or something similar. * * Note: Beware if you compile on Digital Unix and run on Linux. The * program uses popen() to run sendmail; when compiled on Digital Unix * popen() looks for /sbin/sh whereas Linux has /bin/sh instead. * You will have to create an appropriate link or else use batch mode * (see section 4 below). * **/ /** 3. Optimising ** * * Good speeds on Alpha are roughly 70 k iterations per second for a * 300 MHz EV4, 200 k for 500 MHz EV5 and 290 k for 500 MHz EV6. Best * results are obtained with gcc 2.95 or Compaq cc >= 5.7. A 333 MHz * UltraSPARC-IIi running Solaris 7 gets 140 k with Sun's cc 5.0. * * You can compile with various optimisation flags to maximise speed * before starting on the real calculations. Try adding -DPROD=2 (or 3 * or 4) to the compilation flags to use a different implementation of * speed critical functions. It might be a little faster; then again, * maybe not! * * You can also try tweaking the value of PARAL. Just add -DPARAL=40, * for instance. The default is 32 and useful values are roughly in * the range 20 to 50. Note that the test output will no longer * exactly match the sample given above. * **/ /** 4. Starting up ** * * When you're ready to start up: on each machine you want to run on * make a directory for that machine, "cd" into it and go. * * echo Machine name is `hostname` * mkdir `hostname` * cd `hostname` * unlimit cputime || (echo Using ulimit instead; ulimit -t unlimited) * nice +15 ../ecdl.exe [...arguments, see below...] > ecdl.log & * * Note: Some shells need "nice -15" instead. * Note: You can keep an eye on the output using "tail -f ecdl.log". * Note: If you have a machine with several CPUs, make a directory for * each CPU and start up an ecdl process in each one. * * The arguments have the following syntax: * ecdl by on * * : Either "test", "batch", "mail", "alt" or "http". * : Something to identify you uniquely e.g., your email address. * : The machine name e.g., as given by "hostname". * * Note: For give the exact same string you used when signing up, * since otherwise the Web pages will not be able to show correct * statistics for your contribution. * * Examples: * ecdl mail by John.Smith@example.com on avalon * ecdl batch by anonymous on 'Mystery machine' * * Test mode quickly produces some (useless) test output to stdout along * with info on iterations per second. This allows you to check that * your binary is OK, to compare compiler optimisation flags etc. * * The other modes produce things called "distinguished points" every * few hours. Each distinguished point looks something like this: Distinguished point found at: Wed Nov 24 03:52:02 1999 ECC-L0|Unix64|1.1.1| ECC-L1|6B445FC02D031AEA01F28B7AF92|ABD2FF8960D196A4AB43719A64F| ECC-L2|me@here| ECC-L3|000015AE4C35| ECC-L4|bla| * These points are valuable! The ECDL project needs to collect * roughly 1.3 million of them to solve Certicom's ECC2K-108 problem. * The points are always appended to a file called "dist.points" and * also written to stdout along with some verbose info. * * When running in "mail" mode each point is automatically emailed * (using "sendmail") to ecdl2K-108@cristal.inria.fr as well. INRIA's * excessively aggressive spam filters might refuse the email; in that * case switch to "alt" mode which sends to ecdl2K-108@rupture.net instead. * * When running in "http" mode each point is automatically sent (using * the HTTP POST protocol) to the Web server on cristal.inria.fr. If * your machine is behind a firewall, you'll need to specify a Web * proxy to use for sending the points via HTTP. Specify it on the * command line as follows: * * ecdl http by on proxy : * * where is the host name of the proxy, e.g. * "myproxy.mydomain.com", and is the port number of the * HTTP proxy on this host, e.g. "3128". * * For machines unable to send email or make HTTP connections, * for instance not permanently connected to the Net, use "batch" mode. * Then every few days do: * * mv dist.points dist.points.nov.25 * * or similar using the correct date, and email the "dist.points.nov.25" * file manually. When the next distinguished point is found, a new * "dist.points" file will be created for it and further points will * be appended to this new file until the next time you "mv" it. * * Mode: test batch mail alt http * ---------------------------------------------------------------------- * * Writes info on stdout x x x x x * * Writes points to "dist.points" x x x x * * Sends points via sendmail x x * * Sends points via http x * **/ /** 5. Shutting down (and saved states) ** * * When the program is running it saves its internal state in a file * called "saved.state" from time to time. If you want to stop the * program, send it a SIGINT or SIGTERM signal so it can save its * state before exiting. If it is running in the foreground, you can * send SIGINT by typing Ctrl-C. If it is running in the background * send SIGTERM by "killall ecdl.exe". If you don't have "killall" * installed, find the PID (first column of the "ps" output) and use "kill": * * ps x | grep ecdl | grep -v grep * kill * * When you restart in the same directory: * * nice +15 ../ecdl.exe [...arguments...] >> ecdl.log & * * the program will look for the saved state file. If no saved state * is found it will generate a new state and work from there. If one * is found, the program reads it and continues where it left off. * Note that it would be wasteful to use a given saved state on more * than one CPU because work would just be duplicated. * * The reason each machine should have its own directory is to avoid * the dist.points and saved.state files getting mixed up! Note that * test mode does not mess with these files. * **/ /*== #includes =============================================================*/ /** Ansi includes. **/ #include /* For strerror(), strcmp() and strlen(). */ #include /* For signal(), SIGINT and SIGTERM. */ #include /* For errno. */ #include /* For time and date: time() and ctime(). */ #include /* For malloc(). */ #include /** Unix includes. **/ /* For getpid(), access() and timing. */ #include /* For timing: getrusage() and struct rusage. */ #include #include /* Socket stuff */ #include #include #include #include /*== Types =================================================================*/ /* Abbreviation used when bit-length does not matter. */ typedef unsigned int uint; /* Integers. */ typedef unsigned char u8; typedef unsigned int u32; typedef unsigned long u64; typedef struct { u64 hi, lo; } u128; /* Polynomials over the two-element field GF(2) = Z/2Z. */ typedef struct { u64 hi, lo; } poly128; /* Operating modes. */ typedef enum { Bad, Test, Batch, Mail, Alt, Http } modeType; /*== #defines ==============================================================*/ /*-= Stuff that you can change =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=*/ /* The place to find sendmail. Often it's in /usr/sbin/sendmail but * then /usr/lib/sendmail is a link to it, so this works fairly universally. */ #ifndef SENDMAIL #define SENDMAIL "/usr/lib/sendmail" #endif /* Name of device for reading entropy from (it isn't strictly necessary). */ #ifndef RANDOM_DEV_FILENAME #define RANDOM_DEV_FILENAME "/dev/urandom" #endif /* Name of file for writing distinguished points to (as well as stdout). */ #ifndef DIST_POINTS_FILENAME #define DIST_POINTS_FILENAME "dist.points" #endif /* Name of file for saving state after each distinguished point or when * stopping due to a SIGINT or SIGTERM. It is read when restarting. */ #ifndef SAVED_STATE_FILENAME #define SAVED_STATE_FILENAME "saved.state" #endif /* The state is checked and saved every 2^SAVED_STATE_SHIFT iterations * and also after a distinguished point is found. */ #ifndef SAVED_STATE_SHIFT #define SAVED_STATE_SHIFT (27) #endif #if SAVED_STATE_SHIFT < 20 || SAVED_STATE_SHIFT > 63 #error Please define SAVED_STATE_SHIFT in the range 20..63. #endif /* Number of iterations to run "in parallel". This is done so that * PARAL inversions can be replaced by one inversion and 3*PARAL-3 mults. * Try something roughly in the range 20 to 50. */ #ifndef PARAL #define PARAL (32) #endif #if PARAL < 1 #error Please define PARAL to be greater than 0. #endif /* Define PROD to choose a different implementation of the critical * GF2ProductNNxNN() and product() functions. The default is 1. * Try 2, 3 or 4. They might be faster. Then again maybe not! */ #ifndef PROD #define PROD 1 #endif #if PROD < 1 || PROD > 4 #error Please define PROD in the range 1..4. #endif /* Keyword for giving "inline" hint. * Adjust for your compiler or leave it out. */ #if defined(__GNUC__) || defined(__DECC) #define INLINE __inline #else #define INLINE #endif /*-= Stuff that must not be changed =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/ /* Version 0.9.0: 31 September 1999 * - Does 187 k iters/sec with gcc 2.95+Linux+500 MHz Alpha EV5 (PARAL=20). * Version 0.99.0: 24 November 1999 * - Does 204 k iters/sec with gcc 2.95+Linux+500 MHz Alpha EV5 (PARAL=32). * Version 1.0.0: 3 December 1999 * Version 1.1.0: 12 December 1999 * - Added HTTP mode (Xavier Leroy). * - Print time at start-up. * - Proper solution for http mode: use HTTP/1.0 (Roy T. Fielding). * Version 1.1.1: 24 March 1999 * - Some minor tweaks. */ #define VARIANT "Unix64" #define VERSION "1.1.1" /* These are used a few places for types poly128 and u128. */ #define ZERO { 0, 0 } #define ONE { 0, 1 } /* The number of cases used in pseudo-random function. */ #define CASES (7) /* Population count used to decide if a point is distinguished. * Expect a given popc to occur one time in 2^108/binomial(109, popc) * if it's odd (even ones never occur). */ /* 1 time in 1422267 or so (for testing). */ #define TEST_POP_COUNT (29) /* 1 time in 1407700775 or so (for real). */ #define DIST_POP_COUNT (23) /* For clarity. */ #define square(x) (squareNTimes((x), 1)) /*-= Data defining which ECDL problem to solve =-=-=-=-=-=-=-=-=-=-=-=-=-=-=*/ /* Field is GF(2^109) represented as (Z/2Z)[t] / (t^109+t^9+t^2+t+1) * (t^109+t^9+t^2+t+1 is irreducible in (Z/2Z)[t]). */ /* Curve from Certicom ECC2K-108 problem is y^2 + x*y = x^3 + x^2 + 1 * over GF(2^109). * * Group order is G = 2^109+1-w^109-conj(w)^109 where w = (1-sqrt(-7))/2. * Use subgroup of order g = G/2 = 324518553658426701487448656461467 (prime). * * Frobenius is multiplication by w i.e., by * 138423345589698157369693034392981 modulo g. */ #define A { 0, 1 } #define B { 0, 1 } #define ORDER { 0x0FFFFFFFFFFF, 0xFFA621B02C383E9B } /* Points from Certicom ECC2K-108 problem: */ #define XP { 0x0478C46CC963, 0x38CED91565E17257 } #define YP { 0x1E7965E4A3AF, 0xB73A48FC9AB790E9 } #define XQ { 0x1FF0CE5EC618, 0x93F2119C3077C59E } #define YQ { 0x1F20E9B010AC, 0x691C9B87B438241D } /*== Function declarations =================================================*/ /*-= Short little functions =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/ static INLINE int ge(u128 x, u128 y); static INLINE u128 diff(u128 x, u128 y); static INLINE u128 sumMod(u128 x, u128 y, u128 m); static INLINE u128 doubleMod(u128 x, u128 m); static INLINE int equal(poly128 x, poly128 y); static INLINE poly128 xor(poly128 x, poly128 y); static INLINE u64 topBit(u64 x); static INLINE u128 squareMod(u128 x, u128 m); /*-= Main and friends =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/ static modeType parse(int argc, char *argv[]); void catchSignal(int sigNum); int main(int argc, char *argv[]); /*-= Stuff for checking that a point is distinguished =-=-=-=-=-=-=-=-=-=-=-*/ static uint popc(poly128 n); static poly128 changeBasis(poly128 x); /*-= Arithmetic in the field GF(2^109) -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/ static poly128 squareNTimes(poly128 x, uint n); #if PROD == 1 || PROD == 2 static u64 GF2Product54x54(u64 *ph, u64 x, u64 y); #else static u64 GF2Product56x56(u64 *ph, u64 x, u64 y); #endif static poly128 product(poly128 x, poly128 y); static poly128 inverse(poly128 y); static INLINE poly128 quotient(poly128 x, poly128 y); /*-= Arithmetic mod the group order =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/ static u128 productMod(u128 x, u128 y, u128 m); static u128 powerMod(u128 x, u64 y, u128 m); static u128 findMultiplier(const u64 *pCases); /*-= Arithmetic on elliptic curve over the field =-=-=-=-=-=-=-=-=-=-=-=-=-=*/ static int ellipticDouble (poly128 *px2, poly128 *py2, poly128 x, poly128 y, int z); static int ellipticSum ( poly128 *px3, poly128 *py3 , poly128 x1, poly128 y1, int z1, poly128 x2, poly128 y2, int z2 ); static int ellipticProduct (poly128 *px2, poly128 *py2, poly128 x, poly128 y, int z, u128 fac); /*-= Various file manipulation thingies =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/ static int checkUVXY(u128 u, u128 v, poly128 x, poly128 y); static int syncAndCheck (u64 *pCases, u128 *pu, u128 *pv, poly128 x, poly128 y); static void reportDistinguished ( modeType mode, char *argv[] , u64 iters, u128 u, u128 v, u64 total, double dStart ); static void u64ToBytes(u64 data, u8 *p); static u64 bytesToU64(u8 *p); static void writeState (u64 *itersT, u128 *uT, u128 *vT, poly128 *xT, poly128 *yT); static uint readState (u64 *itersT, u128 *uT, u128 *vT, poly128 *xT, poly128 *yT); /*-= Stuff for http mode =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=*/ static int httpPost(char * payload); static int httpParseProxy(char * proxy); static int httpInit(void); /*== Global variables ======================================================*/ static int gotSignal = 0; /*== Function definitions ==================================================*/ /*-= Short little functions =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/ /*-- ge --------------------------------------------------------------------*/ /* Return whether x >= y, for 0 <= x, y < 2^128. */ static INLINE int ge(u128 x, u128 y) { return x.hi > y.hi || x.hi == y.hi && x.lo >= y.lo; } /* end function ge */ /*-- diff ------------------------------------------------------------------*/ /* Return x-y for 0 <= y <= x < 2^128. */ static INLINE u128 diff(u128 x, u128 y) { u64 c; c = x.lo < y.lo; x.hi -= y.hi; x.lo -= y.lo; x.hi -= c; return x; } /* end function diff */ /*-- sumMod ----------------------------------------------------------------*/ /* Return sum x+y modulo m, 0 <= x, y < m <= 2^127. */ static INLINE u128 sumMod(u128 x, u128 y, u128 m) { x.lo += y.lo; x.hi += y.hi; x.hi += x.lo < y.lo; if (ge(x, m)) x = diff(x, m); return x; } /* end function sumMod */ /*-- doubleMod -------------------------------------------------------------*/ /* Return double x<<1 modulo m, where 0 <= x < m <= 2^127. */ static INLINE u128 doubleMod(u128 x, u128 m) { x.hi <<= 1; x.hi |= x.lo>>63; x.lo <<= 1; if (ge(x, m)) x = diff(x, m); return x; } /* end function doubleMod */ /*-- equal -----------------------------------------------------------------*/ /* Return whether x == y, for polys over Z/2Z. */ static INLINE int equal(poly128 x, poly128 y) { return !(x.hi ^ y.hi | x.lo ^ y.lo); } /* end function equal */ /*-- xor -------------------------------------------------------------------*/ /* Return x XOR y, in other words add them as polys over Z/2Z. */ static INLINE poly128 xor(poly128 x, poly128 y) { x.hi ^= y.hi; x.lo ^= y.lo; return x; } /* end function xor */ /*-- topBit ----------------------------------------------------------------*/ /* Return top bit of x, or 0 if x = 0. */ static INLINE u64 topBit(u64 x) { u64 t, y, z; z = x>>1; t = x & 0xFFFFFFFF00000000; y = x | z; if (t) x = t; z = y>>2; y |= z; t = x & 0xFFFF0000FFFF0000; z = y>>4; if (t) x = t; y |= z; z = y>>1; t = x & 0xFF00FF00FF00FF00; if (t) x = t; return x & ~z; } /* end function topBit */ /*-- squareMod -------------------------------------------------------------*/ /* Return square x*x modulo m, for 0 <= x < m <= 2^127. */ static INLINE u128 squareMod(u128 x, u128 m) { return productMod(x, x, m); } /* end function squareMod */ /*-= Main and friends =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/ /*-- parse -----------------------------------------------------------------*/ /* Parse command line and return mode Bad, Batch, Test, Mail, Alt or Http. */ static modeType parse(int argc, char *argv[]) { char ch, *str; modeType mode; if (argc != 6 && argc != 8) return Bad; str = argv[1]; if (!strcmp(str, "test")) mode = Test; else if (!strcmp(str, "batch")) mode = Batch; else if (!strcmp(str, "mail")) mode = Mail; else if (!strcmp(str, "alt")) mode = Alt; else if (!strcmp(str, "http")) mode = Http; else return Bad; if (strcmp(argv[2], "by") ) return Bad; for (str = argv[3]; (ch = *str) != '\0'; ++str) if (ch == '|') return Bad; if (strcmp(argv[4], "on")) return Bad; for (str = argv[5]; (ch = *str) != '\0'; ++str) if (ch == '|') return Bad; if ( argc == 8 && ( mode != Http || strcmp(argv[6], "proxy") || httpParseProxy(argv[7]) == -1 ) ) { return Bad; } /* end if */ return mode; } /* end function parse */ /*-- catchSignal -----------------------------------------------------------*/ /* Signal catcher for SIGINT or SIGTERM. * Allows program to stop cleanly and save state. * User can send the signal by Ctrl-C or "kill ", for instance. */ void catchSignal(int sigNum) { gotSignal = 1; } /* end function catchSignal */ /*-- main ------------------------------------------------------------------*/ int main(int argc, char *argv[]) { /* Constant stuff. */ const int zP = 1, zQ = 1; const u128 order = ORDER; const poly128 xP = XP, yP = YP, xQ = XQ, yQ = YQ; /* These get initialised once and for all. */ uint distPopc; modeType mode; int zj; poly128 xj, yj; /* This changes a few times, then stays fixed in the main loop. */ u128 seed; /* These are variable. */ u64 total; double dStart; /* For timing. */ uint pops[PARAL]; u64 casesTT[PARAL][CASES], itersT[PARAL]; u128 uT[PARAL], vT[PARAL]; poly128 denT[PARAL], x1T[PARAL], x2T[PARAL], y1T[PARAL], y2T[PARAL]; /** Simple sanity check. **/ if (sizeof(u32) != 4 || sizeof(u64) != 8) { puts("Fatal error: *** the size of u32 or u64 is wrong. ***"); return 1; } /* end if */ /** Parse command line. **/ mode = parse(argc, argv); if (mode == Bad) { printf( "Syntax: %s by on " "[proxy :]\n" "See top of source code or README file for details. Stopping.\n" , argv[0] ); return 1; } /* end if */ /** See if we can find sendmail with execute permission. **/ if ((mode == Mail || mode == Alt) && access(SENDMAIL, X_OK)) { printf( "Fatal error: *** sendmail program " SENDMAIL " does not exist\n" " or is not executable ***\n" " (errno = %d: %s).\n" , errno, strerror(errno) ); return 1; } /* end if */ /** Initialize HTTP stuff if http mode selected */ if (mode == Http && httpInit() == -1) { puts("Fatal error: *** http initialisation failed. ***"); return 1; } /* end if */ /** Put up banner. **/ printf( "This is Rob's ECDL program for " VARIANT ", version " VERSION "\n" "Mode ...... %s\n" "Email ..... %s\n" "Machine ... %s\n" , argv[1], argv[3], argv[5] ); { time_t now; /* For time and date. */ now = time(NULL); printf("Starting at %s", now == -1 ? "(unknown time)\n" : ctime(&now)); } /* end block */ fflush(stdout); /** Get random seed. **/ if (mode == Test) { /* Fixed seed for testing. */ seed.hi = 0; seed.lo = 0x0123456789ABCDEF; } else { /* Get really random seed. */ FILE *handle; handle = fopen(RANDOM_DEV_FILENAME, "rb"); if (handle) { int status; fread((void *)&seed, sizeof(seed), (size_t)1, handle); status = fclose(handle); if (status == EOF) { printf( "Warning: fclose() of random device " RANDOM_DEV_FILENAME " failed\n" " (errno = %d: %s).\n" , errno, strerror(errno) ); } /* end if */ } else { printf( "Warning: random device " RANDOM_DEV_FILENAME " not found\n" " (errno = %d: %s). Making do without...\n" , errno, strerror(errno) ); } /* end if */ /* Randomise seed some more (if random device was not found, * seed is stack junk which is OK but may not be random enough!) */ seed.hi ^= getpid(); seed.lo ^= time(NULL); /* Ensure seed < order. */ seed.hi &= 0x00000FFFFFFFFFFF; if (ge(seed, order)) seed = diff(seed, order); } /* end if/else */ /** Initialisation. **/ /* Read starting points from the saved state or generate new random ones. */ { uint i; if (mode == Test) i = 0; else i = readState(itersT, uT, vT, x1T, y1T); if (i < PARAL) { /* Generate random points. */ printf("Generating %u new starting points.\n", (uint)PARAL-i); for (/* i was set above */ ; i < PARAL; ++i) { int z, zu, zv; poly128 xu, xv, yu, yv; /* Euler gamma times 2^108. */ const u128 gamma = { 0x0000093C467E37DB, 0xC7A4D1BE3F81015 }; itersT[i] = 0; do { /* Ensure points are non-zero. */ uT[i] = seed; zu = ellipticProduct(&xu, &yu, xP, yP, zP, seed); seed = sumMod(seed, gamma, order); /* Not very random but will do! */ vT[i] = seed; zv = ellipticProduct(&xv, &yv, xQ, yQ, zQ, seed); seed = sumMod(seed, gamma, order); /* Not very random but will do! */ z = ellipticSum(&x1T[i], &y1T[i], xu, yu, zu, xv, yv, zv); } while (!z); } /* end for (i) */ } /* end if */ /* Initialise counts of cases to 0. */ for (i = 0; i < PARAL; ++i) { uint j; for (j = 0; j < CASES; ++j) casesTT[i][j] = 0; } /* end for (i) */ } /* end block */ /* Amount to jump by after finding a distinguished point. * Note: This is random so that processes accidentally started with * identical saved state files will diverge after a while. * Seed is fixed from here on. It may be zero, that's OK. */ zj = ellipticProduct(&xj, &yj, xQ, yQ, zQ, seed); /* Timing (measured in seconds of user time for this process). */ { struct rusage ru; getrusage(RUSAGE_SELF, &ru); dStart = (double)ru.ru_utime.tv_sec+(double)ru.ru_utime.tv_usec*0.000001; } /* end block */ /** Main loop. **/ /* Sketch of main loop structure: * for ever { * for (all PARAL points) { * if (this one is distinguished) { output it; jump to new one; } * } * for (all PARAL points) { * get denominator for this one; * } * get inverses of all PARAL denominators at once; * for (all PARAL points) { * apply pseudo-random function to this one using inverse(denominator); * } * } */ puts("Computing iterations..."); fflush(stdout); distPopc = mode == Test ? TEST_POP_COUNT : DIST_POP_COUNT; signal(SIGINT, catchSignal); signal(SIGTERM, catchSignal); for (total = 0; !gotSignal; total += PARAL) { uint i; int someWereDistinguished = 0; /** Check for distinguished points. **/ for (i = 0; i < PARAL; ++i) { pops[i] = popc(changeBasis(x1T[i])); if (pops[i] == distPopc) { /* Point is distinguished. */ int z; /* Synchronize the point and check that it is OK. * I.e., update u and v using cases and check [u]P+[v]Q = (x:y:1)... */ if (!syncAndCheck(&casesTT[i][0], &uT[i], &vT[i], x1T[i], y1T[i])) { puts("Fatal error: *** bad distinguished point detected! ***"); return 2; } /* end if */ /* Report it. */ reportDistinguished( mode, argv , itersT[i], uT[i], vT[i], total, dStart ); /* Restart this point by jumping to a new point. */ itersT[i] = 0; z = 1; do { /* Ensure point is non-zero. */ vT[i] = sumMod(vT[i], seed, order); z = ellipticSum(&x1T[i], &y1T[i], x1T[i], y1T[i], z, xj, yj, zj); } while (!z); /* Note this for use immediately below. */ someWereDistinguished = 1; } /* end if (distinguished point) */ } /* end for (i) */ /** Save state after a bunch of iterations or a distinguished point. **/ if ( total && (total & ((1UL<>6); if (m >= 7) m -= 7; /* Reduce modulo 7. */ /* Don't update u and v every time. Just count cases that occur * so that they can be synced up with x and y when needed. */ ++casesTT[i][m]; m += 2; if (m == 3) m = 1; x2T[i] = squareNTimes(x1T[i], m); denT[i] = xor(x1T[i], x2T[i]); y2T[i] = squareNTimes(y1T[i], m); } /* end for (i) */ /** Invert PARAL denominators with 1 inversion and 3*PARAL-3 mults. **/ { poly128 ix, prod, q, prodT[PARAL]; prod = denT[0]; for (i = 1; i < PARAL; ++i) { prodT[i] = prod; prod = product(prod, denT[i]); } /* end for */ q = inverse(prod); for (i = PARAL; --i; ) { ix = product(q, prodT[i]); q = product(q, denT[i]); denT[i] = ix; } /* end for */ denT[0] = q; } /* end block */ /** Get new points. **/ for (i = 0; i < PARAL; ++i) { /* Always in general case. */ poly128 nx, y1 = y1T[i]; poly128 lam = product(xor(y1, y2T[i]), denT[i]); poly128 t = xor(xor(square(lam), lam), x2T[i]); const poly128 a = A; t.lo ^= a.lo; /* t = xor(t, a) with a = 1. */ x1T[i] = nx = xor(t, x1T[i]); y1T[i] = xor(xor(product(lam, t), nx), y1); ++itersT[i]; } /* end for (i) */ } /* end main loop */ puts("Received signal."); /** Save state after receiving signal (except when testing). **/ { uint i; /* Sync and check all points. */ for (i = 0; i < PARAL; ++i) { if (!syncAndCheck(&casesTT[i][0], &uT[i], &vT[i], x1T[i], y1T[i])) { puts("Fatal error: *** bad point detected after signal! ***"); return 4; } /* end if */ } /* end for */ /* Save state. */ if (mode != Test) writeState(itersT, uT, vT, x1T, y1T); } /* end block */ { time_t now; /* For time and date. */ now = time(NULL); printf("Stopping at %s", now == -1 ? "(unknown time)\n" : ctime(&now)); } /* end block */ puts("Bye!"); return 0; } /* end function main */ /*-= Stuff for checking that a point is distinguished =-=-=-=-=-=-=-=-=-=-=-*/ /*-- popc ------------------------------------------------------------------*/ /* Count bits set in a 109-bit poly over Z/2Z. */ static uint popc(poly128 n) { u64 nh, nl; const u64 mask1 = 0x5555555555555555 , mask2 = 0x3333333333333333 , mask3 = 0x0F0F0F0F0F0F0F0F ; /* Partial popc. */ nh = n.hi; nh -= nh>>1 & mask1; /* Count 2-by-2. */ nh = (nh & mask2)+(nh>>2 & mask2); /* Count 4-by-4. */ nh += nh>>4; /* Count 8-by-8. */ nh &= mask3; /* Partial popc. */ nl = n.lo; nl -= nl>>1 & mask1; nl = (nl & mask2)+(nl>>2 & mask2); nl += nl>>4; nl &= mask3; /* Accumulate pieces. */ nh += nl; nh += nh>> 8; nh += nh>>16; nh += nh>>32; return nh & 255; } /* end function popc */ /*-- changeBasis -----------------------------------------------------------*/ /* Go from basis { t^k | 0 <= k < 109 } to basis { (t+1)^2^k | 0 <= k < 109 } * (both over Z/2Z). */ static poly128 changeBasis(poly128 x) { static const poly128 matrix[218] = { { 0, 0 } , { 0x1FFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF } , { 0x1FFFFFFFFFFF, 0xFFFFFFFFFFFFFFFE } , { 0x000000000000, 0x0000000000000001 } , { 0, 0 } , { 0x1FFFFFFFFFFF, 0xFFFFFFFFFFFFFFFD } , { 0x09EA1A34BEEC, 0x0A908559BDEDA6AB } , { 0x1615E5CB4113, 0xF56F7AA642125956 } , { 0, 0 } , { 0x1FFFFFFFFFFF, 0xFFFFFFFFFFFFFFFB } , { 0x0536810828A6, 0x034361502FC46FD6 } , { 0x1AC97EF7D759, 0xFCBC9EAFD03B902D } , { 0, 0 } , { 0x13D434697DD8, 0x15210AB37BDB4D56 } , { 0x1748B0C11920, 0x5DC6E0E909A60AC2 } , { 0x049C84A864F8, 0x48E7EA5A727D4794 } , { 0, 0 } , { 0x1FFFFFFFFFFF, 0xFFFFFFFFFFFFFFF7 } , { 0x1E31E07B050A, 0x23A4555030FAD343 } , { 0x01CE1F84FAF5, 0xDC5BAAAFCF052CB4 } , { 0, 0 } , { 0x0A6D0210514C, 0x0686C2A05F88DFAC } , { 0x02C794F5E66A, 0xAA1EE7366F735E1B } , { 0x08AA96E5B726, 0xAC98259630FB81B7 } , { 0, 0 } , { 0x07A868D2FBB0, 0x2A421566F7B69AAD } , { 0x03469E40B570, 0xCE15C90071B8CCC4 } , { 0x04EEF6924EC0, 0xE457DC66860E5669 } , { 0, 0 } , { 0x0E9161823240, 0xBB8DC1D2134C1585 } , { 0x0A5E6CC8E146, 0x486897C20ACD96E5 } , { 0x04CF0D4AD306, 0xF3E5561019818360 } , { 0, 0 } , { 0x1FFFFFFFFFFF, 0xFFFFFFFFFFFFFFEF } , { 0x1DE5DF8084BB, 0x9970280CA3BF99C1 } , { 0x021A207F7B44, 0x668FD7F35C40662E } , { 0, 0 } , { 0x1C63C0F60A14, 0x4748AAA061F5A687 } , { 0x1F5501B79A31, 0xD77E3E4D9FBDDCC1 } , { 0x0336C1419025, 0x903694EDFE487A46 } , { 0, 0 } , { 0x14DA0420A298, 0x0D0D8540BF11BF58 } , { 0x0952396CA3D7, 0x5297C756F9CBB7EC } , { 0x1D883D4C014F, 0x5F9A421646DA08B4 } , { 0, 0 } , { 0x058F29EBCCD5, 0x543DCE6CDEE6BC36 } , { 0x1BFA18AA70FE, 0x45B2D97ACD057310 } , { 0x1E753141BC2B, 0x118F171613E3CF26 } , { 0, 0 } , { 0x0F50D1A5F760, 0x54842ACDEF6D355A } , { 0x192AAFD26253, 0x0FD989FC56604C4E } , { 0x167A7E779533, 0x5B5DA331B90D7914 } , { 0, 0 } , { 0x068D3C816AE1, 0x9C2B9200E3719988 } , { 0x0095EAC2D27F, 0xAD933BEFB6D89E6C } , { 0x0618D643B89E, 0x31B8A9EF55A907E4 } , { 0, 0 } , { 0x1D22C3046481, 0x771B83A426982B0A } , { 0x0A5E542E870A, 0x9F98AD6E4E573666 } , { 0x177C972AE38B, 0xE8832ECA68CF1D6C } , { 0, 0 } , { 0x14BCD991C28C, 0x90D12F84159B2DCA } , { 0x09C527449B58, 0x9E98DABF2E783164 } , { 0x1D79FED559D4, 0x0E49F53B3BE31CAE } , { 0, 0 } , { 0x1FFFFFFFFFFF, 0xFFFFFFFFFFFFFFDF } , { 0x1C9FE648868C, 0xB844216D22FDFD76 } , { 0x036019B77973, 0x47BBDE92DD0202A9 } , { 0, 0 } , { 0x1BCBBF010977, 0x32E05019477F3383 } , { 0x00C44DACA9C6, 0xF0CE03AB2356535F } , { 0x1B0FF2ADA0B1, 0xC22E53B2642960DC } , { 0, 0 } , { 0x18C781EC1428, 0x8E915540C3EB4D0F } , { 0x1003982774E6, 0xEC619F3152DD5203 } , { 0x08C419CB60CE, 0x62F0CA7191361F0C } , { 0, 0 } , { 0x1EAA036F3463, 0xAEFC7C9B3F7BB983 } , { 0x0BE2DA6A67E3, 0x9960FC5632ADBEF9 } , { 0x1548D9055380, 0x379C80CD0DD6077A } , { 0, 0 } , { 0x09B408414530, 0x1A1B0A817E237EB1 } , { 0x1CF911FD1F88, 0x9BF5DF89D921E2F2 } , { 0x154D19BC5AB8, 0x81EED508A7029C43 } , { 0, 0 } , { 0x12A472D947AE, 0xA52F8EADF3976FD8 } , { 0x1CE48176618A, 0x5063CF7B362CFC89 } , { 0x0E40F3AF2624, 0xF54C41D6C5BB9351 } , { 0, 0 } , { 0x0B1E53D799AA, 0xA87B9CD9BDCD786C } , { 0x0048A61FFC6D, 0x7BBE85CFB0EC51FA } , { 0x0B56F5C865C7, 0xD3C519160D212996 } , { 0, 0 } , { 0x17F43154E1FC, 0x8B65B2F59A0AE621 } , { 0x19E3B34684B4, 0xAE95923514E90A5D } , { 0x0E1782126548, 0x25F020C08EE3EC7C } , { 0, 0 } , { 0x1EA1A34BEEC0, 0xA908559BDEDA6AB4 } , { 0x128CFFAC3CF8, 0x8C2931BE4DB5CD19 } , { 0x0C2D5CE7D238, 0x25216425936FA7AD } , { 0, 0 } , { 0x12555FA4C4A6, 0x1FB313F8ACC0989D } , { 0x0EEFF28924AA, 0x2F8873511C93544F } , { 0x1CBAAD2DE00C, 0x303B60A9B053CCD2 } , { 0, 0 } , { 0x0D1A7902D5C3, 0x38572401C6E33310 } , { 0x0C5C19AD60B0, 0x5E11225D09CBC60F } , { 0x014660AFB573, 0x6646065CCF28F51F } , { 0, 0 } , { 0x012BD585A4FF, 0x5B2677DF6DB13CD8 } , { 0x01C38C1277D0, 0x060B23FCF132BC82 } , { 0x00E85997D32F, 0x5D2D54239C83805A } , { 0, 0 } , { 0x1A458608C902, 0xEE3707484D305615 } , { 0x179DB2113F0F, 0x7814A56C277A9EE3 } , { 0x0DD83419F60D, 0x9623A2246A4AC8F6 } , { 0, 0 } , { 0x14BCA85D0E15, 0x3F315ADC9CAE6CCC } , { 0x157C5B345C1C, 0x643A6D3330FA7A39 } , { 0x01C0F3695209, 0x5B0B37EFAC5416F5 } , { 0, 0 } , { 0x0979B3238519, 0x21A25F082B365B95 } , { 0x01035D70D551, 0xB4E6A8BE5B6FF9C9 } , { 0x087AEE535048, 0x9544F7B67059A25C } , { 0, 0 } , { 0x138A4E8936B1, 0x3D31B57E5CF062C8 } , { 0x1C2F11203F3F, 0xCAB69770DF433D07 } , { 0x0FA55FA9098E, 0xF787220E83B35FCF } , { 0, 0 } , { 0x1FFFFFFFFFFF, 0xFFFFFFFFFFFFFFBF } , { 0x01CDE8DE6EB8, 0x69E4FFE27FD9AA80 } , { 0x1E3217219147, 0x961B001D8026553F } , { 0, 0 } , { 0x193FCC910D19, 0x708842DA45FBFAED } , { 0x127CE34869DD, 0x0640B805225B3C14 } , { 0x0B432FD964C4, 0x76C8FADF67A0C6F9 } , { 0, 0 } , { 0x17977E0212EE, 0x65C0A0328EFE6707 } , { 0x14C6D4CA755E, 0x9F9692672567C9A5 } , { 0x0351AAC867B0, 0xFA563255AB99AEA2 } , { 0, 0 } , { 0x01889B59538D, 0xE19C075646ACA6BE } , { 0x0A9A331E30C9, 0xEFA2FE395AA1F939 } , { 0x0B12A8476344, 0x0E3EF96F1C0D5F87 } , { 0, 0 } , { 0x118F03D82851, 0x1D22AA8187D69A1F } , { 0x095FB83B634D, 0x284C5684E280D9AC } , { 0x18D0BBE34B1C, 0x356EFC05655643B3 } , { 0, 0 } , { 0x0007304EE9CD, 0xD8C33E62A5BAA407 } , { 0x00765EFB7E85, 0x388546D3D82D749F } , { 0x00716EB59748, 0xE04678B17D97D098 } , { 0, 0 } , { 0x1D5406DE68C7, 0x5DF8F9367EF77307 } , { 0x07BA78C41EED, 0x02B4696829DF40AF } , { 0x1AEE7E1A762A, 0x5F4C905E572833A8 } , { 0, 0 } , { 0x17C5B4D4CFC7, 0x32C1F8AC655B7DF2 } , { 0x1D457D6E6970, 0xC19185E5B0A43683 } , { 0x0A80C9BAA6B7, 0xF3507D49D5FF4B71 } , { 0, 0 } , { 0x136810828A60, 0x34361502FC46FD62 } , { 0x1F9F0759C297, 0x1F06E18064DC124E } , { 0x0CF717DB48F7, 0x2B30F482989AEF2C } , { 0, 0 } , { 0x19F223FA3F11, 0x37EBBF13B243C5E5 } , { 0x17B146F4B00F, 0xA9CBE82CE76AB97D } , { 0x0E43650E8F1E, 0x9E20573F55297C98 } , { 0, 0 } , { 0x0548E5B28F5D, 0x4A5F1D5BE72EDFB1 } , { 0x07974CC047A9, 0x517EC6F34F38005C } , { 0x02DFA972C8F4, 0x1B21DBA8A816DFED } , { 0, 0 } , { 0x19C902ECC314, 0xA0C79EF66C59F913 } , { 0x0544FBA4F1D8, 0xE3F79DAFDEE1E43F } , { 0x1C8DF94832CC, 0x43300359B2B81D2C } , { 0, 0 } , { 0x163CA7AF3355, 0x50F739B37B9AF0D8 } , { 0x0408E42CECEE, 0xD1B4C9D22E18F5E8 } , { 0x12344383DFBB, 0x8143F06155820530 } , { 0, 0 } , { 0x00914C3FF8DA, 0xF77D0B9F61D8A3F4 } , { 0x186E7ABA658A, 0x7FF0C89316AA6668 } , { 0x18FF36859D50, 0x888DC30C7772C59C } , { 0, 0 } , { 0x0FE862A9C3F9, 0x16CB65EB3415CC43 } , { 0x1A79D9C93DA7, 0x2891C56BC83CF966 } , { 0x1591BB60FE5E, 0x3E5AA080FC293525 } , { 0, 0 } , { 0x13C7668D0969, 0x5D2B246A29D214BB } , { 0x08A697345AA9, 0xBB2DDE902816E691 } , { 0x1B61F1B953C0, 0xE606FAFA01C4F22A } , { 0, 0 } , { 0x1D434697DD81, 0x5210AB37BDB4D569 } , { 0x146EE982BE06, 0x27B2F0E06A2904E5 } , { 0x092DAF156387, 0x75A25BD7D79DD18C } , { 0, 0 } , { 0x0519FF5879F1, 0x1852637C9B6B9A33 } , { 0x1CF716F88427, 0xA0B5B4D16100D523 } , { 0x19EEE9A0FDD6, 0xB8E7D7ADFA6B4F10 } , { 0, 0 } , { 0x04AABF49894C, 0x3F6627F15981313B } , { 0x0E63E6B11399, 0xC8C3EC99F86F51E9 } , { 0x0AC959F89AD5, 0xF7A5CB68A1EE60D2 } , { 0, 0 } , { 0x1DDFE5124954, 0x5F10E6A23926A89E } , { 0x1A9C19BC4158, 0x11BBB1FBC1EAE3F5 } , { 0x0743FCAE080C, 0x4EAB5759F8CC4B6B } , { 0, 0 } , { 0x1A34F205AB86, 0x70AE48038DC66620 } , { 0x0553811DFA1A, 0xC53F7E1CF16331C3 } , { 0x1F677318519C, 0xB591361F7CA557E3 } , { 0, 0 } , { 0x18B8335AC160, 0xBC2244BA13978C1E } , { 0x0E4FB2233CDC, 0x78EC2FD59BFF7FF6 } , { 0x16F78179FDBC, 0xC4CE6B6F8868F3E8 } , { 0, 0 } , { 0x0257AB0B49FE, 0xB64CEFBEDB6279B0 } } ; uint i; u64 t, yh,yl; const poly128 *p; yh = yl = 0; p = &matrix[0]; #define FOR_BLOCK \ { yh ^= p[t & 3].hi; yl ^= p[t & 3].lo; t >>= 2; p += 4; \ yh ^= p[t & 3].hi; yl ^= p[t & 3].lo; t >>= 2; p += 4; \ } t = x.lo; for (i = 16; i; --i) FOR_BLOCK t = x.hi; for (i = 11; i; --i) FOR_BLOCK /* Last bit. */ yh ^= p[t & 1].hi; yl ^= p[t & 1].lo; #undef FOR_BLOCK { poly128 res; res.hi = yh; res.lo = yl; return res; } /* end block */ } /* end function changeBasis */ /*-= Arithmetic in the field GF(2^109) -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/ /*-- squareNTimes ----------------------------------------------------------*/ /* Square x, n times, in the field GF(2^109) i.e., as a poly in (Z/2Z)[t] * reduced modulo t^109+t^9+t^2+t+1. Degree x < 109. */ static poly128 squareNTimes(poly128 x, uint n) { static const int tab2[256] = { 0, 1, 4, 5, 16, 17, 20, 21, 64, 65 , 68, 69, 80, 81, 84, 85, 256, 257, 260, 261 , 272, 273, 276, 277, 320, 321, 324, 325, 336, 337 , 340, 341, 1024, 1025, 1028, 1029, 1040, 1041, 1044, 1045 , 1088, 1089, 1092, 1093, 1104, 1105, 1108, 1109, 1280, 1281 , 1284, 1285, 1296, 1297, 1300, 1301, 1344, 1345, 1348, 1349 , 1360, 1361, 1364, 1365, 4096, 4097, 4100, 4101, 4112, 4113 , 4116, 4117, 4160, 4161, 4164, 4165, 4176, 4177, 4180, 4181 , 4352, 4353, 4356, 4357, 4368, 4369, 4372, 4373, 4416, 4417 , 4420, 4421, 4432, 4433, 4436, 4437, 5120, 5121, 5124, 5125 , 5136, 5137, 5140, 5141, 5184, 5185, 5188, 5189, 5200, 5201 , 5204, 5205, 5376, 5377, 5380, 5381, 5392, 5393, 5396, 5397 , 5440, 5441, 5444, 5445, 5456, 5457, 5460, 5461, 16384, 16385 , 16388, 16389, 16400, 16401, 16404, 16405, 16448, 16449, 16452, 16453 , 16464, 16465, 16468, 16469, 16640, 16641, 16644, 16645, 16656, 16657 , 16660, 16661, 16704, 16705, 16708, 16709, 16720, 16721, 16724, 16725 , 17408, 17409, 17412, 17413, 17424, 17425, 17428, 17429, 17472, 17473 , 17476, 17477, 17488, 17489, 17492, 17493, 17664, 17665, 17668, 17669 , 17680, 17681, 17684, 17685, 17728, 17729, 17732, 17733, 17744, 17745 , 17748, 17749, 20480, 20481, 20484, 20485, 20496, 20497, 20500, 20501 , 20544, 20545, 20548, 20549, 20560, 20561, 20564, 20565, 20736, 20737 , 20740, 20741, 20752, 20753, 20756, 20757, 20800, 20801, 20804, 20805 , 20816, 20817, 20820, 20821, 21504, 21505, 21508, 21509, 21520, 21521 , 21524, 21525, 21568, 21569, 21572, 21573, 21584, 21585, 21588, 21589 , 21760, 21761, 21764, 21765, 21776, 21777, 21780, 21781, 21824, 21825 , 21828, 21829, 21840, 21841, 21844, 21845 }; static const u32 tab4[256] = { 0x00000000, 0x00000001, 0x00000010, 0x00000011, 0x00000100, 0x00000101 , 0x00000110, 0x00000111, 0x00001000, 0x00001001, 0x00001010, 0x00001011 , 0x00001100, 0x00001101, 0x00001110, 0x00001111, 0x00010000, 0x00010001 , 0x00010010, 0x00010011, 0x00010100, 0x00010101, 0x00010110, 0x00010111 , 0x00011000, 0x00011001, 0x00011010, 0x00011011, 0x00011100, 0x00011101 , 0x00011110, 0x00011111, 0x00100000, 0x00100001, 0x00100010, 0x00100011 , 0x00100100, 0x00100101, 0x00100110, 0x00100111, 0x00101000, 0x00101001 , 0x00101010, 0x00101011, 0x00101100, 0x00101101, 0x00101110, 0x00101111 , 0x00110000, 0x00110001, 0x00110010, 0x00110011, 0x00110100, 0x00110101 , 0x00110110, 0x00110111, 0x00111000, 0x00111001, 0x00111010, 0x00111011 , 0x00111100, 0x00111101, 0x00111110, 0x00111111, 0x01000000, 0x01000001 , 0x01000010, 0x01000011, 0x01000100, 0x01000101, 0x01000110, 0x01000111 , 0x01001000, 0x01001001, 0x01001010, 0x01001011, 0x01001100, 0x01001101 , 0x01001110, 0x01001111, 0x01010000, 0x01010001, 0x01010010, 0x01010011 , 0x01010100, 0x01010101, 0x01010110, 0x01010111, 0x01011000, 0x01011001 , 0x01011010, 0x01011011, 0x01011100, 0x01011101, 0x01011110, 0x01011111 , 0x01100000, 0x01100001, 0x01100010, 0x01100011, 0x01100100, 0x01100101 , 0x01100110, 0x01100111, 0x01101000, 0x01101001, 0x01101010, 0x01101011 , 0x01101100, 0x01101101, 0x01101110, 0x01101111, 0x01110000, 0x01110001 , 0x01110010, 0x01110011, 0x01110100, 0x01110101, 0x01110110, 0x01110111 , 0x01111000, 0x01111001, 0x01111010, 0x01111011, 0x01111100, 0x01111101 , 0x01111110, 0x01111111, 0x10000000, 0x10000001, 0x10000010, 0x10000011 , 0x10000100, 0x10000101, 0x10000110, 0x10000111, 0x10001000, 0x10001001 , 0x10001010, 0x10001011, 0x10001100, 0x10001101, 0x10001110, 0x10001111 , 0x10010000, 0x10010001, 0x10010010, 0x10010011, 0x10010100, 0x10010101 , 0x10010110, 0x10010111, 0x10011000, 0x10011001, 0x10011010, 0x10011011 , 0x10011100, 0x10011101, 0x10011110, 0x10011111, 0x10100000, 0x10100001 , 0x10100010, 0x10100011, 0x10100100, 0x10100101, 0x10100110, 0x10100111 , 0x10101000, 0x10101001, 0x10101010, 0x10101011, 0x10101100, 0x10101101 , 0x10101110, 0x10101111, 0x10110000, 0x10110001, 0x10110010, 0x10110011 , 0x10110100, 0x10110101, 0x10110110, 0x10110111, 0x10111000, 0x10111001 , 0x10111010, 0x10111011, 0x10111100, 0x10111101, 0x10111110, 0x10111111 , 0x11000000, 0x11000001, 0x11000010, 0x11000011, 0x11000100, 0x11000101 , 0x11000110, 0x11000111, 0x11001000, 0x11001001, 0x11001010, 0x11001011 , 0x11001100, 0x11001101, 0x11001110, 0x11001111, 0x11010000, 0x11010001 , 0x11010010, 0x11010011, 0x11010100, 0x11010101, 0x11010110, 0x11010111 , 0x11011000, 0x11011001, 0x11011010, 0x11011011, 0x11011100, 0x11011101 , 0x11011110, 0x11011111, 0x11100000, 0x11100001, 0x11100010, 0x11100011 , 0x11100100, 0x11100101, 0x11100110, 0x11100111, 0x11101000, 0x11101001 , 0x11101010, 0x11101011, 0x11101100, 0x11101101, 0x11101110, 0x11101111 , 0x11110000, 0x11110001, 0x11110010, 0x11110011, 0x11110100, 0x11110101 , 0x11110110, 0x11110111, 0x11111000, 0x11111001, 0x11111010, 0x11111011 , 0x11111100, 0x11111101, 0x11111110, 0x11111111 }; u64 xh,xl; xh = x.hi; xl = x.lo; /* Two at a time. */ for ( ; n >= 2; n -= 2) { u64 t, yh, yl; t = (u64)tab4[xh>>40 /* & 31 */ ]<<32 | (u64)tab4[xh>>32 & 255]; yh = (u64)tab4[xh>>24 & 255]<<32 | (u64)tab4[xh>>16 & 255]; yl = (u64)tab4[xh>> 8 & 255]<<32 | (u64)tab4[xh & 255]; yh = yh ^ t>>45 ^ t>>44 ^ t>>43 ^ t>>36; yl = yl ^ t<<19 ^ t<<20 ^ t<<21 ^ t<<28; t = yh; yh = yl; yl = (u64)tab4[xl>>56]<<32 | (u64)tab4[xl>>48 & 255]; yh = yh ^ t>>45 ^ t>>44 ^ t>>43 ^ t>>36; yl = yl ^ t<<19 ^ t<<20 ^ t<<21 ^ t<<28; t = yh; yh = yl; yl = (u64)tab4[xl>>40 & 255]<<32 | (u64)tab4[xl>>32 & 255]; yh = yh ^ t>>45 ^ t>>44 ^ t>>43 ^ t>>36; yl = yl ^ t<<19 ^ t<<20 ^ t<<21 ^ t<<28; t = yh; yh = yl; yl = (u64)tab4[xl>>24 & 255]<<32 | (u64)tab4[xl>>16 & 255]; yh = yh ^ t>>45 ^ t>>44 ^ t>>43 ^ t>>36; yl = yl ^ t<<19 ^ t<<20 ^ t<<21 ^ t<<28; t = yh; yh = yl; yl = (u64)tab4[xl>> 8 & 255]<<32 | (u64)tab4[xl & 255]; yh = yh ^ t>>45 ^ t>>44 ^ t>>43 ^ t>>36; yl = yl ^ t<<19 ^ t<<20 ^ t<<21 ^ t<<28; /* Reduce modulo t^109+t^9+t^2+t+1. */ t = yh>>45; xh = yh ^ t<<45; xl = yl ^ t ^ t<<1 ^ t<<2 ^ t<<9; } /* end for */ /* Last one (if n was odd). */ if (n) { u64 t0,t1,t2,t3, tmp; t3 = (u64)tab2[xh>>40 /* & 31 */ ]<<16 | (u64)tab2[xh>>32 & 255]; t2 = (u64)tab2[xh>>24 & 255]<<48 | (u64)tab2[xh>>16 & 255]<<32 | (u64)tab2[xh>> 8 & 255]<<16 | (u64)tab2[xh & 255] ; t1 = (u64)tab2[xl>>56]<<48 | (u64)tab2[xl>>48 & 255]<<32 | (u64)tab2[xl>>40 & 255]<<16 | (u64)tab2[xl>>32 & 255] ; t0 = (u64)tab2[xl>>24 & 255]<<48 | (u64)tab2[xl>>16 & 255]<<32 | (u64)tab2[xl>> 8 & 255]<<16 | (u64)tab2[xl & 255] ; /* Reduce modulo t^109+t^9+t^2+t+1. */ t1 ^= t3<<19 | t2>>45; t0 ^= t2<<19; t1 ^= t3<<20 | t2>>44; t0 ^= t2<<20; t1 ^= t3<<21 | t2>>43; t0 ^= t2<<21; t1 ^= t3<<28 | t2>>36; t0 ^= t2<<28; tmp = t1>>45; t1 ^= tmp<<45; t0 ^= tmp ^ tmp<<1 ^ tmp<<2 ^ tmp<<9; xh = t1; xl = t0; } /* end if */ { poly128 r; r.hi = xh; r.lo = xl; return r; } /* end block */ } /* end function squareNTimes */ /*# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # */ #if PROD == 1 || PROD == 2 /*-- GF2Product54x54 -------------------------------------------------------*/ /* Multiply y by low 54 bits of x, as polys over Z/2Z, degree y < 54. * Speed-critical auxiliary function used for product(). * Returns low 64 bits of result, puts high 43 bits in *ph. */ static u64 GF2Product54x54(u64 *ph, u64 x, u64 y) { u64 tab[16]; { u64 e, y1,y2,y3; y1 = y<<1; y2 = y<<2; y3 = y<<3; /* Gray code walk through table. */ #if PROD == 1 e = 0; tab[0] = e; tab[8] = e ^ y3; e ^= y; tab[1] = e; tab[9] = e ^ y3; e ^= y1; tab[3] = e; tab[11] = e ^ y3; e ^= y; tab[2] = e; tab[10] = e ^ y3; e ^= y2; tab[6] = e; tab[14] = e ^ y3; e ^= y; tab[7] = e; tab[15] = e ^ y3; e ^= y1; tab[5] = e; tab[13] = e ^ y3; e ^= y; tab[4] = e; tab[12] = e ^ y3; #else /* PROD == 2 */ e = 0; tab[0] = 0; e ^= y; tab[1] = y; e ^= y1; tab[3] = e; e ^= y; tab[2] = y1; e ^= y2; tab[6] = e; e ^= y; tab[7] = e; e ^= y1; tab[5] = e; e ^= y; tab[4] = y2; e ^= y3; tab[12] = e; e ^= y; tab[13] = e; e ^= y1; tab[15] = e; e ^= y; tab[14] = e; e ^= y2; tab[10] = e; e ^= y; tab[11] = e; e ^= y1; tab[9] = e; e ^= y; tab[8] = y3; #endif } /* end block */ { u64 a,b,c,d,e,f,g; a = tab[x & 15] ^ tab[x>> 4 & 15]<<4; b = tab[x>> 8 & 15] ^ tab[x>>12 & 15]<<4; c = tab[x>>16 & 15] ^ tab[x>>20 & 15]<<4; d = tab[x>>24 & 15] ^ tab[x>>28 & 15]<<4; e = tab[x>>32 & 15] ^ tab[x>>36 & 15]<<4; f = tab[x>>40 & 15] ^ tab[x>>44 & 15]<<4; g = tab[x>>48 & 15] ^ tab[x>>52 & 3]<<4; *ph = b>>56 ^ c>>48 ^ d>>40 ^ e>>32 ^ f>>24 ^ g>>16; return a ^ b<<8 ^ c<<16 ^ d<<24 ^ e<<32 ^ f<<40 ^ g<<48; } /* end block */ } /* end function GF2Product54x54 */ /*# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # */ #else /* PROD == 3 || PROD == 4 */ /*-- GF2Product56x56 -------------------------------------------------------*/ /* Multiply y by low 56 bits of x, as polys over Z/2Z, degree y < 56. * Speed-critical auxiliary function used for product(). * Returns low 64 bits of result, puts high 47 bits in *ph. */ static u64 GF2Product56x56(u64 *ph, u64 x, u64 y) { u64 tab[16]; { u64 e, y1,y2,y3; y1 = y<<1; y2 = y<<2; y3 = y<<3; /* Gray code walk through table. */ #if PROD == 3 e = 0; tab[0] = e; tab[8] = e ^ y3; e ^= y; tab[1] = e; tab[9] = e ^ y3; e ^= y1; tab[3] = e; tab[11] = e ^ y3; e ^= y; tab[2] = e; tab[10] = e ^ y3; e ^= y2; tab[6] = e; tab[14] = e ^ y3; e ^= y; tab[7] = e; tab[15] = e ^ y3; e ^= y1; tab[5] = e; tab[13] = e ^ y3; e ^= y; tab[4] = e; tab[12] = e ^ y3; #else /* PROD == 4 */ e = 0; tab[0] = 0; e ^= y; tab[1] = y; e ^= y1; tab[3] = e; e ^= y; tab[2] = y1; e ^= y2; tab[6] = e; e ^= y; tab[7] = e; e ^= y1; tab[5] = e; e ^= y; tab[4] = y2; e ^= y3; tab[12] = e; e ^= y; tab[13] = e; e ^= y1; tab[15] = e; e ^= y; tab[14] = e; e ^= y2; tab[10] = e; e ^= y; tab[11] = e; e ^= y1; tab[9] = e; e ^= y; tab[8] = y3; #endif } /* end block */ { u64 a,b,c,d,e,f,g; a = tab[x & 15] ^ tab[x>> 4 & 15]<<4; b = tab[x>> 8 & 15] ^ tab[x>>12 & 15]<<4; c = tab[x>>16 & 15] ^ tab[x>>20 & 15]<<4; d = tab[x>>24 & 15] ^ tab[x>>28 & 15]<<4; e = tab[x>>32 & 15] ^ tab[x>>36 & 15]<<4; f = tab[x>>40 & 15] ^ tab[x>>44 & 15]<<4; g = tab[x>>48 & 15] ^ tab[x>>52 & 15]<<4; *ph = b>>56 ^ c>>48 ^ d>>40 ^ e>>32 ^ f>>24 ^ g>>16; return a ^ b<<8 ^ c<<16 ^ d<<24 ^ e<<32 ^ f<<40 ^ g<<48; } /* end block */ } /* end function GF2Product56x56 */ /*# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # */ #endif /*-- product ---------------------------------------------------------------*/ /* Return product of x and y in the field GF(2^109) i.e., as polys in * (Z/2Z)[t] reduced modulo t^109+t^9+t^2+t+1. Degree x, y < 109. */ static poly128 product(poly128 x, poly128 y) { u64 lh,ll, mh,ml, hl,hh, tmp, t0,t1,t2,t3, xh, yh; #if PROD == 1 || PROD == 2 xh = x.hi<<10 | x.lo>>54; yh = y.hi<<10 | y.lo>>54; hl = GF2Product54x54(&hh, xh, yh & 0x003FFFFFFFFFFFFF); ml = GF2Product54x54(&mh, xh ^ x.lo, (yh ^ y.lo) & 0x003FFFFFFFFFFFFF); ml ^= hl; mh ^= hh; ll = GF2Product54x54(&lh, x.lo, y.lo & 0x003FFFFFFFFFFFFF); ml ^= ll; mh ^= lh; t0 = ll; t1 = lh | hl<<44; t2 = hl>>20 | hh<<44; t3 = hh>>20; t0 ^= ml<<54; t1 ^= ml>>10 | mh<<54; t2 ^= mh>>10; /* Handle 109th bit. */ { u64 xb, yb; xb = x.hi>>44; yb = y.hi>>44; if (xb) { t1 ^= y.lo<<44; t2 ^= y.lo>>20; t2 ^= y.hi<<44; t3 ^= y.hi>>20; } if (yb) { t1 ^= x.lo<<44; t2 ^= x.lo>>20; t2 ^= x.hi<<44; t3 ^= x.hi>>20; } if (xb & yb) t3 ^= 1UL<<24; } /* end block */ #else xh = x.hi<<8 | x.lo>>56; yh = y.hi<<8 | y.lo>>56; hl = GF2Product56x56(&hh, xh, yh); ll = GF2Product56x56(&lh, x.lo, y.lo & 0x00FFFFFFFFFFFFFF); ml = GF2Product56x56(&mh, xh ^ x.lo, yh ^ y.lo & 0x00FFFFFFFFFFFFFF); ml ^= ll; mh ^= lh; ml ^= hl; mh ^= hh; t0 = ll; t1 = lh | hl<<48; t2 = hl>>16 | hh<<48; t3 = hh>>16; t0 ^= ml<<56; t1 ^= ml>>8 | mh<<56; t2 ^= mh>>8; #endif /* Reduce modulo t^109+t^9+t^2+t+1. */ t1 ^= t3<<19 | t2>>45; t0 ^= t2<<19; t1 ^= t3<<20 | t2>>44; t0 ^= t2<<20; t1 ^= t3<<21 | t2>>43; t0 ^= t2<<21; t1 ^= t3<<28 | t2>>36; t0 ^= t2<<28; tmp = t1>>45; t1 ^= tmp<<45; t0 ^= tmp ^ tmp<<1 ^ tmp<<2 ^ tmp<<9; { poly128 r; r.hi = t1; r.lo = t0; return r; } /* end block */ } /* end function product */ /*-- inverse ---------------------------------------------------------------*/ /* Invert y in the field GF(2^109) i.e., as a poly in (Z/2Z)[t] * reduced modulo t^109+t^9+t^2+t+1. Degree y < 109, y != 0. * Note: 1/(1+t+t^2+t^9+O(t^10)) = 1+t+t^3+t^4+t^6+t^7+O(t^10) */ static poly128 inverse(poly128 y) { u64 ah,al, bh,bl, sh,th, t, uh,ul, vh,vl; static const int shiftTab[64] = { 6,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0 , 5,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0 }; static const int xorTab[64] = { 0, 27, 54, 45, 44, 55, 26, 1, 24, 3, 46, 53, 52, 47, 2, 25, 48, 43 , 6, 29, 28, 7, 42, 49, 40, 51, 30, 5, 4, 31, 50, 41, 32, 59, 22, 13 , 12, 23, 58, 33, 56, 35, 14, 21, 20, 15, 34, 57, 16, 11, 38, 61, 60, 39 , 10, 17, 8, 19, 62, 37, 36, 63, 18, 9 }; static const int xorTab2[64] = { 0, 13889, 27778, 23235, 22724, 28293, 13382, 519, 12360, 1545 , 23754, 27275, 26764, 24269, 1038, 12879, 24720, 22225, 3090, 14931 , 14420, 3605, 21718, 25239, 20696, 26265, 15450, 2587, 2076, 15965 , 25758, 21215, 16608, 30369, 11362, 6691, 6180, 11877, 29862, 17127 , 28840, 18153, 7210, 10859, 10348, 7725, 17646, 29359, 8304, 5681 , 19698, 31411, 30900, 20213, 5174, 8823, 4152, 9849, 31930, 19195 , 18684, 32445, 9342, 4671 }; const u64 mh = 1UL<<45, ml = 519; /* Maintain: u.y = a and v.y = b modulo t^109+t^9+t^2+t+1. */ ah = y.hi; al = y.lo; uh = 0; ul = 1; while (!(al & 1)) { t = ul & 1; uh ^= t<<45; if (t) ul ^= ml; al = al>>1 | ah<<63; ah >>= 1; ul = ul>>1 | uh<<63; uh >>= 1; } /* end while */ bh = mh; bl = ml; vh = 0; vl = 0; do { do { do { bh ^= ah; bl ^= al; vh ^= uh; vl ^= ul; sh = shiftTab[bl & 63]; th = 64-sh; t = vl & 63; vh ^= (u64)xorTab[t]<<45; vl ^= xorTab2[t]; bl = bl>>sh | bh<>= sh; vl = vl>>sh | vh<>= sh; while (!(bl & 1)) { t = vl & 1; vh ^= t<<45; if (t) vl ^= ml; bl = bl>>1 | bh<<63; bh >>= 1; vl = vl>>1 | vh<<63; vh >>= 1; } /* end while */ } while (ah < bh || ah == bh && al < bl); if (al == bl && ah == bh) break; do { ah ^= bh; al ^= bl; uh ^= vh; ul ^= vl; sh = shiftTab[al & 63]; th = 64-sh; t = ul & 63; uh ^= (u64)xorTab[t]<<45; ul ^= xorTab2[t]; al = al>>sh | ah<>= sh; ul = ul>>sh | uh<>= sh; while (!(al & 1)) { t = ul & 1; uh ^= t<<45; if (t) ul ^= ml; al = al>>1 | ah<<63; ah >>= 1; ul = ul>>1 | uh<<63; uh >>= 1; } /* end while */ } while (ah > bh || ah == bh && al > bl); } while (al != bl); } while (ah != bh); /* Now a (and b) must equal 1. */ /* Reduce u modulo t^109+t^9+t^2+t+1. */ t = uh>>45; uh ^= t<<45; ul ^= t ^ t<<1 ^ t<<2 ^ t<<9; { poly128 r; r.hi = uh; r.lo = ul; return r; } /* end block */ } /* end function inverse */ /*-- quotient --------------------------------------------------------------*/ /* Divide x by y in the field GF(2^109) i.e., as polys in (Z/2Z)[t] * reduced modulo t^109+t^9+t^2+t+1. Degree x, y < 109, y != 0. */ static INLINE poly128 quotient(poly128 x, poly128 y) { return product(x, inverse(y)); } /* end function quotient */ /*-= Arithmetic mod the group order =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/ /*-- productMod ------------------------------------------------------------*/ /* Return product x*y modulo m, 0 <= x, y < m <= 2^127. * Could this possibly be any slower? =:-) */ static u128 productMod(u128 x, u128 y, u128 m) { u64 b, w; u128 z; const u128 uZero = ZERO; w = y.hi; if (!w) w = y.lo; if (!w) return uZero; b = topBit(w); z = x; while (b >>= 1) { z = doubleMod(z, m); if (w & b) z = sumMod(z, x, m); } /* end while */ if (y.hi) { b = 1UL<<63; w = y.lo; do { z = doubleMod(z, m); if (w & b) z = sumMod(z, x, m); b >>= 1; } while (b); } /* end if */ return z; } /* end function productMod */ /*-- powerMod --------------------------------------------------------------*/ /* Return power x^y modulo m, 0 <= x < m <= 2^127, 0 <= y < 2^64. * Note that 0^0 returns 1 (even when m = 1). */ static u128 powerMod(u128 x, u64 y, u128 m) { u64 b; u128 z; const u128 uOne = ONE; if (!y) return uOne; b = topBit(y); z = x; while (b >>= 1) { z = squareMod(z, m); if (y & b) z = productMod(z, x, m); } /* end while */ return z; } /* end function powerMod */ /*-- findMultiplier --------------------------------------------------------*/ /* Compute multiplier Prod_i pow(multipliers[i], pCases[i]) modulo order. */ static u128 findMultiplier(const u64 *pCases) { uint i; u128 f; const u128 order = ORDER; /* Here multiplier[i] = (Frobenius to power i) + 1 modulo group order. */ static const u128 multipliers[CASES] = { { 0x6D325F3CBE9, 0x20245E2659822D94 } , { 0x6D325F3CBE9, 0x20245E2659822D96 } , { 0xB868E249C44, 0x9EDF28ED4BE9F47A } , { 0x92CDA0C3416, 0xDF81C389D2B6110D } , { 0x21FBDC2FB8D, 0xA169935F671A66B6 } , { 0xFC609AA935F, 0xE1B24FAC1A1EC1D4 } , { 0xB868E249C44, 0x9EDF28ED4BE9F46A } }; f = powerMod(multipliers[0], pCases[0], order); for (i = 1; i < CASES; ++i) { f = productMod(f, powerMod(multipliers[i], pCases[i], order), order); } /* end for */ return f; } /* end function findMultiplier */ /*-= Arithmetic on elliptic curve over the field =-=-=-=-=-=-=-=-=-=-=-=-=-=*/ /*-- ellipticDouble --------------------------------------------------------*/ /* Given a point (x:y:z) on y^2 + x*y = x^3 + a*x^2 + b, compute * its double (x2:y2:z2) by the group law. * This puts x2 in *px2, y2 in *py2 and returns z2. * The point at infinity is represented by (0:1:0). * Finite points are represented by (x:y:1). */ static int ellipticDouble ( poly128 *px2, poly128 *py2, poly128 x, poly128 y, int z ) { poly128 lam, x2,y2; const poly128 a = A, zero = ZERO, one = ONE; if (!z || equal(x, zero)) { *px2 = zero; *py2 = one; return 0; } lam = xor(x, quotient(y, x)); x2 = xor(xor(square(lam), lam), a); y2 = xor(xor(product(lam, xor(x, x2)), x2), y); *px2 = x2; *py2 = y2; return 1; } /* end function ellipticDouble */ /*-- ellipticSum -----------------------------------------------------------*/ /* Given points (x1:y1:z1) and (x2:y2:z2) on y^2 + x*y = x^3 + a*x^2 + b, * compute their sum (x3:y3:z3) by the group law. * This puts x3 in *px3, y3 in *py3 and returns z3. * The point at infinity is represented by (0:1:0). * Finite points are represented by (x:y:1). */ static int ellipticSum ( poly128 *px3, poly128 *py3 , poly128 x1, poly128 y1, int z1, poly128 x2, poly128 y2, int z2 ) { poly128 lam, t, x3,y3; const poly128 a = A, zero = ZERO, one = ONE; if (!z1) { *px3 = x2; *py3 = y2; return z2; } if (!z2) { *px3 = x1; *py3 = y1; return z1; } if (equal(x1, x2)) { if (equal(y1, y2)) return ellipticDouble(px3, py3, x1, y1, z1); *px3 = zero; *py3 = one; return 0; } /* end if */ lam = quotient(xor(y1, y2), xor(x1, x2)); t = xor(xor(xor(square(lam), lam), x2), a); x3 = xor(t, x1); y3 = xor(xor(product(lam, t), x3), y1); *px3 = x3; *py3 = y3; return 1; } /* end function ellipticSum */ /*-- ellipticProduct -------------------------------------------------------*/ /* Given a point (x:y:z) on y^2 + x*y = x^3 + a*x^2 + b, compute * its fac-th multiple (x2:y2:z2) by the group law. * This puts x2 in *px2, y2 in *py2 and returns z2. * The point at infinity is represented by (0:1:0). * Finite points are represented by (x:y:1). */ static int ellipticProduct ( poly128 *px2, poly128 *py2, poly128 x, poly128 y, int z, u128 fac ) { int z2; u64 f, m; poly128 x2,y2; const poly128 zero = ZERO, one = ONE; /* Get most significant 64-bit half of fac. */ f = fac.hi; if (!f) f = fac.lo; /* Quick exit if fac was 0. */ if (!f) { *px2 = zero; *py2 = one; return 0; } /* Run through bits of most significant half. */ m = topBit(f); x2 = x; y2 = y; z2 = z; while (m >>= 1) { z2 = ellipticDouble(&x2, &y2, x2, y2, z2); if (f & m) z2 = ellipticSum(&x2, &y2, x, y, z, x2, y2, z2); } /* end while */ /* Run through bits of least significant half, if necessary. */ if (fac.hi) { m = 1UL<<63; f = fac.lo; do { z2 = ellipticDouble(&x2, &y2, x2, y2, z2); if (f & m) z2 = ellipticSum(&x2, &y2, x, y, z, x2, y2, z2); m >>= 1; } while (m); } /* end if */ *px2 = x2; *py2 = y2; return z2; } /* end function ellipticProduct */ /*-= Various file manipulation thingies =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/ /*-- checkUVXY -------------------------------------------------------------*/ /* Check that u,v,x,y are in range and that [u]P+[v]Q = (x:y:1). */ static int checkUVXY(u128 u, u128 v, poly128 x, poly128 y) { int zc, zu, zv; poly128 xc,yc, xu,yu, xv,yv; const int zP = 1, zQ = 1; const poly128 xP = XP, xQ = XQ, yP = YP, yQ = YQ; const u128 order = ORDER; if (ge(u, order)) return 0; if (ge(v, order)) return 0; if (x.hi > 0x00001FFFFFFFFFFF) return 0; if (y.hi > 0x00001FFFFFFFFFFF) return 0; zu = ellipticProduct(&xu, &yu, xP, yP, zP, u); zv = ellipticProduct(&xv, &yv, xQ, yQ, zQ, v); zc = ellipticSum(&xc, &yc, xu, yu, zu, xv, yv, zv); if (!equal(xc, x) || !equal(yc, y) || zc != 1) return 0; return 1; } /* end function checkUVXY */ /*-- syncAndCheck ----------------------------------------------------------*/ /* Synchronize u and v with x and y, then check that everything is OK. * Updates u and v using cases and resets pCases[]. * Then checks u,v,x,y are in range and that [u]P+[v]Q = (x:y:1). * Note: used before reporting a distinguished point or saving state. */ static int syncAndCheck ( u64 *pCases, u128 *pu, u128 *pv, poly128 x, poly128 y ) { uint i; u128 f; const u128 order = ORDER; f = findMultiplier(pCases); for (i = 0; i < CASES; ++i) pCases[i] = 0; *pu = productMod(*pu, f, order); *pv = productMod(*pv, f, order); return checkUVXY(*pu, *pv, x, y); } /* end function syncAndCheck */ /*-- reportDistinguished ---------------------------------------------------*/ /* Report distinguished points back to headquarters. */ static void reportDistinguished ( modeType mode, char *argv[] , u64 iters, u128 u, u128 v, u64 total, double dStart ) { const char format[] = "ECC-L0|" VARIANT "|" VERSION "|\n" "ECC-L1|%011lX%016lX|%011lX%016lX|\n" "ECC-L2|%s|\n" "ECC-L3|%012lX|\n" "ECC-L4|%s|\n" , test_format[] = "TEST-L0|" VARIANT "|" VERSION "|\n" "TEST-L1|%011lX%016lX|%011lX%016lX|\n" "TEST-L2|%s|\n" "TEST-L3|%012lX|\n" "TEST-L4|%s|\n" ; int status; FILE *handle; time_t now; /* For time and date. */ static char *buf = NULL; /** Allocate buffer first time around. **/ if (!buf) { buf = malloc( strlen(format) /* Upper bound on size of fixed text. */ + 27 + 27 /* Two 27-digit hex numbers on line 1. */ + strlen(argv[3]) /* The user ID on line 2. */ + 12 /* One 12-digit hex number on line 3. */ + strlen(argv[5]) /* The machine ID on line 4. */ + 1 /* The terminating '\0'. */ ); } /* end if */ /** Output the distinguished point to stdout and the points file. **/ if (mode == Test) handle = NULL; else { handle = fopen(DIST_POINTS_FILENAME, "a"); if (!handle) { printf( "Warning: couldn't fopen() file " DIST_POINTS_FILENAME " for appending\n" " (errno = %d: %s). Following point is not saved!\n" , errno, strerror(errno) ); } /* end if */ } /* end if */ putchar('\n'); if (handle) putc('\n', handle); /* Time and date. */ now = time(NULL); { char *str1, *str2; str1 = mode == Test ? "Test" : "Distinguished"; str2 = now == -1 ? "(unknown time)\n" : ctime(&now); printf("%s point found at: %s", str1, str2); if (handle) fprintf(handle, "%s point found at: %s", str1, str2); } /* end if/else */ /* The distinguished point. */ printf( mode == Test ? test_format : format , u.hi, u.lo, v.hi, v.lo, argv[3], iters, argv[5] ); if (handle) { /* Note: mode != Test. */ fprintf(handle, format, u.hi, u.lo, v.hi, v.lo, argv[3], iters, argv[5]); } /* end if */ if (handle) { status = fclose(handle); if (status == EOF) { printf( "Warning: fclose() of file " DIST_POINTS_FILENAME " failed\n" " (errno = %d: %s).\n" , errno, strerror(errno) ); } /* end if */ } /* end if */ /** Output extra verbosity to stdout. **/ /* Iteration counts. */ printf( "Iterations used = %lu.\n" "Total iterations = %lu.\n" , iters, total ); /* Iteration rate. */ { double dNow, dUserTime; /* For timing. */ struct rusage ru; getrusage(RUSAGE_SELF, &ru); dNow = (double)ru.ru_utime.tv_sec+(double)ru.ru_utime.tv_usec*0.000001; dUserTime = dNow-dStart; if (dUserTime) printf("Rate = %g per second.\n", (double)total/dUserTime); } /* end block */ /** Send the distinguished point by email. **/ if (mode == Mail || mode == Alt) { handle = popen(SENDMAIL " -t", "w"); if (handle == NULL) { puts("Warning: couldn't pipe to " SENDMAIL ", send by hand!"); } else { fprintf( handle , "To: %s\n\n" , mode == Mail ? "ecdl2K-108@cristal.inria.fr" : "ecdl2K-108@rupture.net" ); /* The distinguished point. Note: mode != Test. */ fprintf( handle, format , u.hi, u.lo, v.hi, v.lo, argv[3], iters, argv[5] ); status = fflush(handle); if (status == EOF) { printf( "Warning: fflush() of pipe failed\n" " (errno = %d: %s).\n" , errno, strerror(errno) ); } /* end if */ status = pclose(handle); if (status) { if (status == -1) { printf( "Warning: pclose() failed\n" " (errno = %d: %s).\n" , errno, strerror(errno) ); } else printf("Warning: " SENDMAIL " returned status %d.\n", status); } /* end if */ } /* end if/else */ } /* end if (mode == Mail || mode == Alt) */ /** Send the distinguished point by HTTP. **/ if (mode == Http) { if (!buf) { puts( "Warning: couldn't allocate temp buffer for HTTP sending, " "send by hand!" ); } else { sprintf(buf, format, u.hi, u.lo, v.hi, v.lo, argv[3], iters, argv[5]); if (httpPost(buf) == -1) { puts( "Warning: an error occurred while sending point via HTTP, " "send by hand!" ); } /* end if */ } /* end if/else */ } /* end if (mode == Http) */ fflush(stdout); } /* end function reportDistinguished */ /*-- u64ToBytes ------------------------------------------------------------*/ static void u64ToBytes(u64 data, u8 *p) { p[0] = data; p[1] = data>> 8; p[2] = data>>16; p[3] = data>>24; p[4] = data>>32; p[5] = data>>40; p[6] = data>>48; p[7] = data>>56; } /* end function u64ToBytes */ /*-- bytesToU64 ------------------------------------------------------------*/ static u64 bytesToU64(u8 *p) { return p[0] | (u64)p[1]<<8 | (u64)p[2]<<16 | (u64)p[3]<<24 | (u64)p[4]<<32 | (u64)p[5]<<40 | (u64)p[6]<<48 | (u64)p[7]<<56 ; } /* end function bytesToU64 */ /*-- writeState ------------------------------------------------------------*/ /* Write state out to file in a portable binary format. */ static void writeState ( u64 *itersT, u128 *uT, u128 *vT, poly128 *xT, poly128 *yT ) { FILE *handle; time_t now; now = time(NULL); printf( "Saving state to " SAVED_STATE_FILENAME " at %s" , now == -1 ? "(unknown time)\n" : ctime(&now) ); handle = fopen(SAVED_STATE_FILENAME, "wb"); if (handle) { uint i; int status; u8 buf[72]; for (i = 0; i < PARAL; ++i) { #define OUT128(offset, X) \ u64ToBytes(X.lo, buf+offset); u64ToBytes(X.hi, buf+offset+8) /* ; */ u64ToBytes(itersT[i], buf); OUT128( 8, uT[i]); OUT128(24, vT[i]); OUT128(40, xT[i]); OUT128(56, yT[i]); #undef OUT128 if (fwrite((void *)buf, sizeof(buf), 1, handle) != 1) { printf( "Warning: couldn't fwrite() to saved state file\n" " (errno = %d: %s).\n" , errno, strerror(errno) ); break; } /* end if */ } /* end for */ status = fclose(handle); if (status == EOF) { printf( "Warning: fclose() of saved state file failed\n" " (errno = %d: %s).\n" , errno, strerror(errno) ); } /* end if */ } else printf( "Warning: couldn't fopen() saved state file for writing\n" " (errno = %d: %s).\n" , errno, strerror(errno) ); fflush(stdout); } /* end function writeState */ /*-- readState -------------------------------------------------------------*/ /* Read state back from file (if it exists) in a portable binary format * and return number of valid items (i.e., u,v,x,y data) found. * Note: There is nothing to check for the itersT[] array. */ static uint readState ( u64 *itersT, u128 *uT, u128 *vT, poly128 *xT, poly128 *yT ) { uint i; FILE *handle; i = 0; handle = fopen(SAVED_STATE_FILENAME, "rb"); if (handle) { int status; u8 buf[72]; for (i = 0; i < PARAL; ++i) { if (fread((void *)buf, sizeof(buf), 1, handle) != 1) { if (feof(handle)) puts("Warning: short saved state file."); else { printf( "Warning: couldn't fread() from saved state file\n" " (errno = %d: %s).\n" , errno, strerror(errno) ); } /* end if/else */ break; } /* end if */ #define IN128(offset, X) \ X.lo = bytesToU64(buf+offset); X.hi = bytesToU64(buf+offset+8) /* ; */ itersT[i] = bytesToU64(buf); IN128( 8, uT[i]); IN128(24, vT[i]); IN128(40, xT[i]); IN128(56, yT[i]); #undef IN128 if (!checkUVXY(uT[i], vT[i], xT[i], yT[i])) { puts("Warning: invalid starting point found in saved state!"); break; } /* end if */ } /* end for */ printf( "Read %u of %u starting points from saved state.\n" , i, (uint)PARAL ); status = fclose(handle); if (status == EOF) { printf( "Warning: fclose() of saved state file failed\n" " (errno = %d: %s).\n" , errno, strerror(errno) ); } /* end if */ } else printf( "Warning: couldn't fopen() saved state file for reading\n" " (errno = %d: %s).\n" "Note: 'No such file' is normal for the first run.\n" , errno, strerror(errno) ); fflush(stdout); return i; } /* end function readState */ /*-= Stuff for http mode =-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=*/ /*-- httpPost --------------------------------------------------------------*/ /* Configuration */ #define HTTP_URL "http://cristal.inria.fr/bin/ecdl7" #define HTTP_SERVER_HOST "cristal.inria.fr" #define HTTP_SERVER_PORT 80 /* INET address of HTTP server or proxy */ static struct sockaddr_in http_addr; /* Name of HTTP server or proxy */ static char * http_host = HTTP_SERVER_HOST; /* Port number on HTTP server or proxy */ static int http_port = HTTP_SERVER_PORT; /* Send a string back to the base using the HTTP POST protocol */ static int httpPost(char * payload) { int s, retcode; FILE * f, * r; /* Connect to HTTP server or proxy */ s = socket(AF_INET, SOCK_STREAM, 0); if (s == -1) { printf("Warning: socket() failed (errno = %d: %s).\n", errno, strerror(errno)); return -1; } retcode = connect(s, (struct sockaddr *) &http_addr, sizeof(http_addr)); if (retcode == -1) { printf("Warning: couldn't connect to host %s port %d (errno = %d: %s).\n", http_host, http_port, errno, strerror(errno)); return -1; } /* Prepare to send data */ f = fdopen(s, "wb"); if (f == NULL) { close(s); printf("Warning: fdopen(socket, \"wb\") failed (errno = %d: %s).\n", errno, strerror(errno)); return -1; } /* Send POST request + data */ fprintf(f, "POST " HTTP_URL " HTTP/1.0\r\n"); fprintf(f, "Host: %s:%d\r\n", HTTP_SERVER_HOST, HTTP_SERVER_PORT); fprintf(f, "Content-length: %lu\r\n", (unsigned long)strlen(payload)); fprintf(f, "\r\n"); fputs(payload, f); fflush(f); if (ferror(f)) { fclose(f); printf("Warning: error while sending data via HTTP to %s.\n", http_host); return -1; } /* Say that we're done sending data */ shutdown(s, 1); /* Read response from server or proxy */ r = fdopen(s, "rb"); if (r == NULL) { fclose(f); printf("Warning: fdopen(socket, \"rb\") failed (errno = %d: %s).\n", errno, strerror(errno)); return -1; } /* First line of response is the status line, with format HTTP-Version SP Status-Code SP Reason-Phrase CRLF. */ retcode = -1; fscanf(r, "%*s %d ", &retcode); /* Discard remainder of response */ fclose(f); fclose(r); /* Check status code -- should be 2xx */ if (! (retcode >= 200 && retcode < 300)) { printf("Warning: Web server or proxy %s returned an error " "(status code = %d)\n", http_host, retcode); return -1; } return 0; } /*-- httpParseProxy --------------------------------------------------------*/ /* Parse proxy specification and store proxy host name and port number * in http_host and http_port */ static int httpParseProxy(char * proxy) { http_host = malloc(strlen(proxy)); if (http_host == NULL) { printf("Error: cannot allocate proxy host name\n"); return -1; } if (sscanf(proxy, "%[^:]:%d", http_host, &http_port) != 2) { printf("Error: bad proxy specification `%s'\n", proxy); return -1; } return 0; } /*-- httpInit --------------------------------------------------------------*/ /* Initialize the HTTP machinery. Determine INET address of * server or proxy and store it in http_addr */ static int httpInit(void) { struct hostent * entry; entry = gethostbyname(http_host); if (entry == NULL) { printf("Error: unknown host %s\n", http_host); return -1; } memset(&http_addr, 0, sizeof(struct sockaddr_in)); http_addr.sin_family = AF_INET; if ( entry->h_length < 0 || (size_t)entry->h_length > sizeof(http_addr.sin_addr.s_addr) ) { puts("Error: address buffer overflow!"); return -1; } /* end if */ memcpy(&http_addr.sin_addr.s_addr, entry->h_addr, entry->h_length); http_addr.sin_port = htons(http_port); return 0; } /*== end of file ecdl2K-108.c ==============================================*/