hcrypto: import libtommath v1.2.0

This commit is contained in:
Luke Howard
2020-04-12 18:37:13 +10:00
parent 7181c109d0
commit c403b66082
287 changed files with 28273 additions and 38374 deletions

View File

@@ -2,12 +2,13 @@
#include <tommath.h>
#include <time.h>
int sizes[] = {256, 512, 768, 1024, 1536, 2048, 3072, 4096};
static int sizes[] = {256, 512, 768, 1024, 1536, 2048, 3072, 4096};
int main(void)
{
char buf[2000];
int x, y;
size_t x;
mp_bool y;
mp_int q, p;
FILE *out;
clock_t t1;
@@ -16,69 +17,65 @@ int main(void)
mp_init_multi(&q, &p, NULL);
out = fopen("2kprime.1", "w");
for (x = 0; x < (int)(sizeof(sizes) / sizeof(sizes[0])); x++) {
top:
mp_2expt(&q, sizes[x]);
mp_add_d(&q, 3, &q);
z = -3;
if (out != NULL) {
for (x = 0; x < (sizeof(sizes) / sizeof(sizes[0])); x++) {
top:
mp_2expt(&q, sizes[x]);
mp_add_d(&q, 3uL, &q);
z = -3;
t1 = clock();
for(;;) {
mp_sub_d(&q, 4, &q);
z += 4;
t1 = clock();
for (;;) {
mp_sub_d(&q, 4uL, &q);
z += 4uL;
if (z > MP_MASK) {
printf("No primes of size %d found\n", sizes[x]);
break;
}
if ((clock() - t1) > CLOCKS_PER_SEC) {
printf(".");
fflush(stdout);
/* sleep((clock() - t1 + CLOCKS_PER_SEC/2)/CLOCKS_PER_SEC); */
t1 = clock();
}
/* quick test on q */
mp_prime_is_prime(&q, 1, &y);
if (y == MP_NO) {
continue;
}
/* find (q-1)/2 */
mp_sub_d(&q, 1uL, &p);
mp_div_2(&p, &p);
mp_prime_is_prime(&p, 3, &y);
if (y == MP_NO) {
continue;
}
/* test on q */
mp_prime_is_prime(&q, 3, &y);
if (y == MP_NO) {
continue;
}
if (z > MP_MASK) {
printf("No primes of size %d found\n", sizes[x]);
break;
}
if (clock() - t1 > CLOCKS_PER_SEC) {
printf("."); fflush(stdout);
// sleep((clock() - t1 + CLOCKS_PER_SEC/2)/CLOCKS_PER_SEC);
t1 = clock();
if (y == MP_NO) {
++sizes[x];
goto top;
}
/* quick test on q */
mp_prime_is_prime(&q, 1, &y);
if (y == 0) {
continue;
}
/* find (q-1)/2 */
mp_sub_d(&q, 1, &p);
mp_div_2(&p, &p);
mp_prime_is_prime(&p, 3, &y);
if (y == 0) {
continue;
}
/* test on q */
mp_prime_is_prime(&q, 3, &y);
if (y == 0) {
continue;
}
break;
}
if (y == 0) {
++sizes[x];
goto top;
}
mp_toradix(&q, buf, 10);
printf("\n\n%d-bits (k = %lu) = %s\n", sizes[x], z, buf);
fprintf(out, "%d-bits (k = %lu) = %s\n", sizes[x], z, buf); fflush(out);
mp_to_decimal(&q, buf, sizeof(buf));
printf("\n\n%d-bits (k = %lu) = %s\n", sizes[x], z, buf);
fprintf(out, "%d-bits (k = %lu) = %s\n", sizes[x], z, buf);
fflush(out);
}
fclose(out);
}
return 0;
}
/* $Source: /cvs/libtom/libtommath/etc/2kprime.c,v $ */
/* $Revision: 1.2 $ */
/* $Date: 2005/05/05 14:38:47 $ */

View File

@@ -1,10 +1,12 @@
/* Makes safe primes of a DR nature */
#include <tommath.h>
int sizes[] = { 1+256/DIGIT_BIT, 1+512/DIGIT_BIT, 1+768/DIGIT_BIT, 1+1024/DIGIT_BIT, 1+2048/DIGIT_BIT, 1+4096/DIGIT_BIT };
static int sizes[] = { 1+256/MP_DIGIT_BIT, 1+512/MP_DIGIT_BIT, 1+768/MP_DIGIT_BIT, 1+1024/MP_DIGIT_BIT, 1+2048/MP_DIGIT_BIT, 1+4096/MP_DIGIT_BIT };
int main(void)
{
int res, x, y;
mp_bool res;
int x, y;
char buf[4096];
FILE *out;
mp_int a, b;
@@ -13,52 +15,53 @@ int main(void)
mp_init(&b);
out = fopen("drprimes.txt", "w");
for (x = 0; x < (int)(sizeof(sizes)/sizeof(sizes[0])); x++) {
top:
printf("Seeking a %d-bit safe prime\n", sizes[x] * DIGIT_BIT);
mp_grow(&a, sizes[x]);
mp_zero(&a);
for (y = 1; y < sizes[x]; y++) {
a.dp[y] = MP_MASK;
}
if (out != NULL) {
for (x = 0; x < (int)(sizeof(sizes)/sizeof(sizes[0])); x++) {
top:
printf("Seeking a %d-bit safe prime\n", sizes[x] * MP_DIGIT_BIT);
mp_grow(&a, sizes[x]);
mp_zero(&a);
for (y = 1; y < sizes[x]; y++) {
a.dp[y] = MP_MASK;
}
/* make a DR modulus */
a.dp[0] = -1;
a.used = sizes[x];
/* make a DR modulus */
a.dp[0] = -1;
a.used = sizes[x];
/* now loop */
res = 0;
for (;;) {
a.dp[0] += 4;
if (a.dp[0] >= MP_MASK) break;
mp_prime_is_prime(&a, 1, &res);
if (res == 0) continue;
printf("."); fflush(stdout);
mp_sub_d(&a, 1, &b);
mp_div_2(&b, &b);
mp_prime_is_prime(&b, 3, &res);
if (res == 0) continue;
mp_prime_is_prime(&a, 3, &res);
if (res == 1) break;
}
/* now loop */
res = MP_NO;
for (;;) {
a.dp[0] += 4uL;
if (a.dp[0] >= MP_MASK) break;
mp_prime_is_prime(&a, 1, &res);
if (res == MP_NO) continue;
printf(".");
fflush(stdout);
mp_sub_d(&a, 1uL, &b);
mp_div_2(&b, &b);
mp_prime_is_prime(&b, 3, &res);
if (res == MP_NO) continue;
mp_prime_is_prime(&a, 3, &res);
if (res == MP_YES) break;
}
if (res != 1) {
printf("Error not DR modulus\n"); sizes[x] += 1; goto top;
} else {
mp_toradix(&a, buf, 10);
printf("\n\np == %s\n\n", buf);
fprintf(out, "%d-bit prime:\np == %s\n\n", mp_count_bits(&a), buf); fflush(out);
}
if (res != MP_YES) {
printf("Error not DR modulus\n");
sizes[x] += 1;
goto top;
} else {
mp_to_decimal(&a, buf, sizeof(buf));
printf("\n\np == %s\n\n", buf);
fprintf(out, "%d-bit prime:\np == %s\n\n", mp_count_bits(&a), buf);
fflush(out);
}
}
fclose(out);
}
fclose(out);
mp_clear(&a);
mp_clear(&b);
return 0;
}
/* $Source: /cvs/libtom/libtommath/etc/drprime.c,v $ */
/* $Revision: 1.2 $ */
/* $Date: 2005/05/05 14:38:47 $ */

View File

@@ -1,4 +1,5 @@
CFLAGS += -Wall -W -Wshadow -O3 -fomit-frame-pointer -funroll-loops -I../
LTM_CFLAGS += -Wall -W -Wextra -Wshadow -O3 -I../
LTM_CFLAGS += $(CFLAGS)
# default lib name (requires install with root)
# LIBNAME=-ltommath
@@ -8,43 +9,36 @@ LIBNAME=../libtommath.a
#provable primes
pprime: pprime.o
$(CC) pprime.o $(LIBNAME) -o pprime
$(CC) $(LTM_CFLAGS) pprime.o $(LIBNAME) -o pprime
# portable [well requires clock()] tuning app
tune: tune.o
$(CC) tune.o $(LIBNAME) -o tune
# same app but using RDTSC for higher precision [requires 80586+], coff based gcc installs [e.g. ming, cygwin, djgpp]
tune86: tune.c
nasm -f coff timer.asm
$(CC) -DX86_TIMER $(CFLAGS) tune.c timer.o $(LIBNAME) -o tune86
# for cygwin
tune86c: tune.c
nasm -f gnuwin32 timer.asm
$(CC) -DX86_TIMER $(CFLAGS) tune.c timer.o $(LIBNAME) -o tune86
$(CC) $(LTM_CFLAGS) tune.o $(LIBNAME) -o tune
./tune_it.sh
test_standalone: tune.o
# The benchmark program works as a testtool, too
$(CC) $(LTM_CFLAGS) tune.o $(LIBNAME) -o test
#make tune86 for linux or any ELF format
tune86l: tune.c
nasm -f elf -DUSE_ELF timer.asm
$(CC) -DX86_TIMER $(CFLAGS) tune.c timer.o $(LIBNAME) -o tune86l
# spits out mersenne primes
mersenne: mersenne.o
$(CC) mersenne.o $(LIBNAME) -o mersenne
$(CC) $(LTM_CFLAGS) mersenne.o $(LIBNAME) -o mersenne
# fines DR safe primes for the given config
# finds DR safe primes for the given config
drprime: drprime.o
$(CC) drprime.o $(LIBNAME) -o drprime
# fines 2k safe primes for the given config
$(CC) $(LTM_CFLAGS) drprime.o $(LIBNAME) -o drprime
# finds 2k safe primes for the given config
2kprime: 2kprime.o
$(CC) 2kprime.o $(LIBNAME) -o 2kprime
$(CC) $(LTM_CFLAGS) 2kprime.o $(LIBNAME) -o 2kprime
mont: mont.o
$(CC) mont.o $(LIBNAME) -o mont
$(CC) $(LTM_CFLAGS) mont.o $(LIBNAME) -o mont
clean:
rm -f *.log *.o *.obj *.exe pprime tune mersenne drprime tune86 tune86l mont 2kprime pprime.dat \
*.da *.dyn *.dpi *~
rm -f *.log *.o *.obj *.exe pprime tune mersenne drprime mont 2kprime pprime.dat \
tuning_list multiplying squaring test *.da *.dyn *.dpi *~
rm -rf .libs
.PHONY: tune

View File

@@ -8,7 +8,7 @@ CFLAGS += -I../
# -ax? specifies make code specifically for ? but compatible with IA-32
# -x? specifies compile solely for ? [not specifically IA-32 compatible]
#
# where ? is
# where ? is
# K - PIII
# W - first P4 [Williamette]
# N - P4 Northwood
@@ -28,15 +28,15 @@ LIBNAME=../libtommath.a
pprime: pprime.o
$(CC) pprime.o $(LIBNAME) -o pprime
# portable [well requires clock()] tuning app
tune: tune.o
$(CC) tune.o $(LIBNAME) -o tune
$(CC) $(CFLAGS) tune.o $(LIBNAME) -o tune
./tune_it.sh
# same app but using RDTSC for higher precision [requires 80586+], coff based gcc installs [e.g. ming, cygwin, djgpp]
tune86: tune.c
nasm -f coff timer.asm
$(CC) -DX86_TIMER $(CFLAGS) tune.c timer.o $(LIBNAME) -o tune86
# for cygwin
tune86c: tune.c
nasm -f gnuwin32 timer.asm
@@ -46,7 +46,7 @@ tune86c: tune.c
tune86l: tune.c
nasm -f elf -DUSE_ELF timer.asm
$(CC) -DX86_TIMER $(CFLAGS) tune.c timer.o $(LIBNAME) -o tune86l
# spits out mersenne primes
mersenne: mersenne.o
$(CC) mersenne.o $(LIBNAME) -o mersenne
@@ -54,7 +54,7 @@ mersenne: mersenne.o
# fines DR safe primes for the given config
drprime: drprime.o
$(CC) drprime.o $(LIBNAME) -o drprime
# fines 2k safe primes for the given config
2kprime: 2kprime.o
$(CC) 2kprime.o $(LIBNAME) -o 2kprime
@@ -62,6 +62,6 @@ drprime: drprime.o
mont: mont.o
$(CC) mont.o $(LIBNAME) -o mont
clean:
rm -f *.log *.o *.obj *.exe pprime tune mersenne drprime tune86 tune86l mont 2kprime pprime.dat *.il
rm -f *.log *.o *.obj *.exe pprime tune mersenne drprime tune86 tune86l mont 2kprime pprime.dat *.il tuning_list

View File

@@ -9,10 +9,11 @@ pprime: pprime.obj
mersenne: mersenne.obj
cl mersenne.obj ../tommath.lib
tune: tune.obj
cl tune.obj ../tommath.lib
mont: mont.obj
cl mont.obj ../tommath.lib

View File

@@ -5,140 +5,134 @@
#include <time.h>
#include <tommath.h>
int
is_mersenne (long s, int *pp)
static mp_err is_mersenne(long s, mp_bool *pp)
{
mp_int n, u;
int res, k;
mp_int n, u;
mp_err res;
int k;
*pp = 0;
*pp = MP_NO;
if ((res = mp_init (&n)) != MP_OKAY) {
return res;
}
if ((res = mp_init(&n)) != MP_OKAY) {
return res;
}
if ((res = mp_init (&u)) != MP_OKAY) {
goto LBL_N;
}
if ((res = mp_init(&u)) != MP_OKAY) {
goto LBL_N;
}
/* n = 2^s - 1 */
if ((res = mp_2expt(&n, s)) != MP_OKAY) {
goto LBL_MU;
}
if ((res = mp_sub_d (&n, 1, &n)) != MP_OKAY) {
goto LBL_MU;
}
/* set u=4 */
mp_set (&u, 4);
/* for k=1 to s-2 do */
for (k = 1; k <= s - 2; k++) {
/* u = u^2 - 2 mod n */
if ((res = mp_sqr (&u, &u)) != MP_OKAY) {
/* n = 2^s - 1 */
if ((res = mp_2expt(&n, (int)s)) != MP_OKAY) {
goto LBL_MU;
}
if ((res = mp_sub_d (&u, 2, &u)) != MP_OKAY) {
}
if ((res = mp_sub_d(&n, 1uL, &n)) != MP_OKAY) {
goto LBL_MU;
}
}
/* make sure u is positive */
while (u.sign == MP_NEG) {
if ((res = mp_add (&u, &n, &u)) != MP_OKAY) {
/* set u=4 */
mp_set(&u, 4uL);
/* for k=1 to s-2 do */
for (k = 1; k <= (s - 2); k++) {
/* u = u^2 - 2 mod n */
if ((res = mp_sqr(&u, &u)) != MP_OKAY) {
goto LBL_MU;
}
if ((res = mp_sub_d(&u, 2uL, &u)) != MP_OKAY) {
goto LBL_MU;
}
}
/* reduce */
if ((res = mp_reduce_2k (&u, &n, 1)) != MP_OKAY) {
goto LBL_MU;
}
}
/* make sure u is positive */
while (u.sign == MP_NEG) {
if ((res = mp_add(&u, &n, &u)) != MP_OKAY) {
goto LBL_MU;
}
}
/* if u == 0 then its prime */
if (mp_iszero (&u) == 1) {
mp_prime_is_prime(&n, 8, pp);
if (*pp != 1) printf("FAILURE\n");
}
/* reduce */
if ((res = mp_reduce_2k(&u, &n, 1uL)) != MP_OKAY) {
goto LBL_MU;
}
}
res = MP_OKAY;
LBL_MU:mp_clear (&u);
LBL_N:mp_clear (&n);
return res;
/* if u == 0 then its prime */
if (mp_iszero(&u) == MP_YES) {
mp_prime_is_prime(&n, 8, pp);
if (*pp != MP_YES) printf("FAILURE\n");
}
res = MP_OKAY;
LBL_MU:
mp_clear(&u);
LBL_N:
mp_clear(&n);
return res;
}
/* square root of a long < 65536 */
long
i_sqrt (long x)
static long i_sqrt(long x)
{
long x1, x2;
long x1, x2;
x2 = 16;
do {
x1 = x2;
x2 = x1 - ((x1 * x1) - x) / (2 * x1);
} while (x1 != x2);
x2 = 16;
do {
x1 = x2;
x2 = x1 - ((x1 * x1) - x) / (2 * x1);
} while (x1 != x2);
if (x1 * x1 > x) {
--x1;
}
if ((x1 * x1) > x) {
--x1;
}
return x1;
return x1;
}
/* is the long prime by brute force */
int
isprime (long k)
static int isprime(long k)
{
long y, z;
long y, z;
y = i_sqrt (k);
for (z = 2; z <= y; z++) {
if ((k % z) == 0)
return 0;
}
return 1;
y = i_sqrt(k);
for (z = 2; z <= y; z++) {
if ((k % z) == 0)
return 0;
}
return 1;
}
int
main (void)
int main(void)
{
int pp;
long k;
clock_t tt;
mp_bool pp;
long k;
clock_t tt;
k = 3;
k = 3;
for (;;) {
/* start time */
tt = clock ();
for (;;) {
/* start time */
tt = clock();
/* test if 2^k - 1 is prime */
if (is_mersenne (k, &pp) != MP_OKAY) {
printf ("Whoa error\n");
return -1;
}
/* test if 2^k - 1 is prime */
if (is_mersenne(k, &pp) != MP_OKAY) {
printf("Whoa error\n");
return -1;
}
if (pp == 1) {
/* count time */
tt = clock () - tt;
if (pp == MP_YES) {
/* count time */
tt = clock() - tt;
/* display if prime */
printf ("2^%-5ld - 1 is prime, test took %ld ticks\n", k, tt);
}
/* display if prime */
printf("2^%-5ld - 1 is prime, test took %ld ticks\n", k, (long)tt);
}
/* goto next odd exponent */
k += 2;
/* but make sure its prime */
while (isprime (k) == 0) {
/* goto next odd exponent */
k += 2;
}
}
return 0;
}
/* $Source: /cvs/libtom/libtommath/etc/mersenne.c,v $ */
/* $Revision: 1.3 $ */
/* $Date: 2006/03/31 14:18:47 $ */
/* but make sure its prime */
while (isprime(k) == 0) {
k += 2;
}
}
}

View File

@@ -1,50 +1,44 @@
/* tests the montgomery routines */
#include <tommath.h>
#include <stdlib.h>
#include <time.h>
int main(void)
{
mp_int modulus, R, p, pp;
mp_digit mp;
long x, y;
int x, y;
srand(time(NULL));
mp_init_multi(&modulus, &R, &p, &pp, NULL);
/* loop through various sizes */
for (x = 4; x < 256; x++) {
printf("DIGITS == %3ld...", x); fflush(stdout);
printf("DIGITS == %3d...", x);
fflush(stdout);
/* make up the odd modulus */
mp_rand(&modulus, x);
modulus.dp[0] |= 1;
/* make up the odd modulus */
mp_rand(&modulus, x);
modulus.dp[0] |= 1uL;
/* now find the R value */
mp_montgomery_calc_normalization(&R, &modulus);
mp_montgomery_setup(&modulus, &mp);
/* now find the R value */
mp_montgomery_calc_normalization(&R, &modulus);
mp_montgomery_setup(&modulus, &mp);
/* now run through a bunch tests */
for (y = 0; y < 1000; y++) {
mp_rand(&p, x/2); /* p = random */
mp_mul(&p, &R, &pp); /* pp = R * p */
mp_montgomery_reduce(&pp, &modulus, mp);
/* now run through a bunch tests */
for (y = 0; y < 1000; y++) {
mp_rand(&p, x/2); /* p = random */
mp_mul(&p, &R, &pp); /* pp = R * p */
mp_montgomery_reduce(&pp, &modulus, mp);
/* should be equal to p */
if (mp_cmp(&pp, &p) != MP_EQ) {
printf("FAILURE!\n");
exit(-1);
}
}
printf("PASSED\n");
}
/* should be equal to p */
if (mp_cmp(&pp, &p) != MP_EQ) {
printf("FAILURE!\n");
exit(-1);
}
}
printf("PASSED\n");
}
return 0;
return 0;
}
/* $Source: /cvs/libtom/libtommath/etc/mont.c,v $ */
/* $Revision: 1.2 $ */
/* $Date: 2005/05/05 14:38:47 $ */

View File

@@ -4,397 +4,408 @@
*
* Tom St Denis, tomstdenis@gmail.com, http://tom.gmail.com
*/
#include <stdlib.h>
#include <time.h>
#include "tommath.h"
int n_prime;
FILE *primes;
static int n_prime;
static FILE *primes;
/* fast square root */
static mp_digit
i_sqrt (mp_word x)
static mp_digit i_sqrt(mp_word x)
{
mp_word x1, x2;
mp_word x1, x2;
x2 = x;
do {
x1 = x2;
x2 = x1 - ((x1 * x1) - x) / (2 * x1);
} while (x1 != x2);
x2 = x;
do {
x1 = x2;
x2 = x1 - ((x1 * x1) - x) / (2u * x1);
} while (x1 != x2);
if (x1 * x1 > x) {
--x1;
}
if ((x1 * x1) > x) {
--x1;
}
return x1;
return x1;
}
/* generates a prime digit */
static void gen_prime (void)
static void gen_prime(void)
{
mp_digit r, x, y, next;
FILE *out;
mp_digit r, x, y, next;
FILE *out;
out = fopen("pprime.dat", "wb");
out = fopen("pprime.dat", "wb");
if (out != NULL) {
/* write first set of primes */
r = 3; fwrite(&r, 1, sizeof(mp_digit), out);
r = 5; fwrite(&r, 1, sizeof(mp_digit), out);
r = 7; fwrite(&r, 1, sizeof(mp_digit), out);
r = 11; fwrite(&r, 1, sizeof(mp_digit), out);
r = 13; fwrite(&r, 1, sizeof(mp_digit), out);
r = 17; fwrite(&r, 1, sizeof(mp_digit), out);
r = 19; fwrite(&r, 1, sizeof(mp_digit), out);
r = 23; fwrite(&r, 1, sizeof(mp_digit), out);
r = 29; fwrite(&r, 1, sizeof(mp_digit), out);
r = 31; fwrite(&r, 1, sizeof(mp_digit), out);
/* write first set of primes */
/* *INDENT-OFF* */
r = 3uL; fwrite(&r, 1uL, sizeof(mp_digit), out);
r = 5uL; fwrite(&r, 1uL, sizeof(mp_digit), out);
r = 7uL; fwrite(&r, 1uL, sizeof(mp_digit), out);
r = 11uL; fwrite(&r, 1uL, sizeof(mp_digit), out);
r = 13uL; fwrite(&r, 1uL, sizeof(mp_digit), out);
r = 17uL; fwrite(&r, 1uL, sizeof(mp_digit), out);
r = 19uL; fwrite(&r, 1uL, sizeof(mp_digit), out);
r = 23uL; fwrite(&r, 1uL, sizeof(mp_digit), out);
r = 29uL; fwrite(&r, 1uL, sizeof(mp_digit), out);
r = 31uL; fwrite(&r, 1uL, sizeof(mp_digit), out);
/* *INDENT-ON* */
/* get square root, since if 'r' is composite its factors must be < than this */
y = i_sqrt (r);
next = (y + 1) * (y + 1);
/* get square root, since if 'r' is composite its factors must be < than this */
y = i_sqrt(r);
next = (y + 1uL) * (y + 1uL);
for (;;) {
do {
r += 2; /* next candidate */
r &= MP_MASK;
if (r < 31) break;
for (;;) {
do {
r += 2uL; /* next candidate */
r &= MP_MASK;
if (r < 31uL) break;
/* update sqrt ? */
if (next <= r) {
++y;
next = (y + 1) * (y + 1);
}
/* update sqrt ? */
if (next <= r) {
++y;
next = (y + 1uL) * (y + 1uL);
}
/* loop if divisible by 3,5,7,11,13,17,19,23,29 */
if ((r % 3) == 0) {
x = 0;
continue;
}
if ((r % 5) == 0) {
x = 0;
continue;
}
if ((r % 7) == 0) {
x = 0;
continue;
}
if ((r % 11) == 0) {
x = 0;
continue;
}
if ((r % 13) == 0) {
x = 0;
continue;
}
if ((r % 17) == 0) {
x = 0;
continue;
}
if ((r % 19) == 0) {
x = 0;
continue;
}
if ((r % 23) == 0) {
x = 0;
continue;
}
if ((r % 29) == 0) {
x = 0;
continue;
}
/* loop if divisible by 3,5,7,11,13,17,19,23,29 */
if ((r % 3uL) == 0uL) {
x = 0uL;
continue;
}
if ((r % 5uL) == 0uL) {
x = 0uL;
continue;
}
if ((r % 7uL) == 0uL) {
x = 0uL;
continue;
}
if ((r % 11uL) == 0uL) {
x = 0uL;
continue;
}
if ((r % 13uL) == 0uL) {
x = 0uL;
continue;
}
if ((r % 17uL) == 0uL) {
x = 0uL;
continue;
}
if ((r % 19uL) == 0uL) {
x = 0uL;
continue;
}
if ((r % 23uL) == 0uL) {
x = 0uL;
continue;
}
if ((r % 29uL) == 0uL) {
x = 0uL;
continue;
}
/* now check if r is divisible by x + k={1,7,11,13,17,19,23,29} */
for (x = 30; x <= y; x += 30) {
if ((r % (x + 1)) == 0) {
x = 0;
break;
/* now check if r is divisible by x + k={1,7,11,13,17,19,23,29} */
for (x = 30uL; x <= y; x += 30uL) {
if ((r % (x + 1uL)) == 0uL) {
x = 0uL;
break;
}
if ((r % (x + 7uL)) == 0uL) {
x = 0uL;
break;
}
if ((r % (x + 11uL)) == 0uL) {
x = 0uL;
break;
}
if ((r % (x + 13uL)) == 0uL) {
x = 0uL;
break;
}
if ((r % (x + 17uL)) == 0uL) {
x = 0uL;
break;
}
if ((r % (x + 19uL)) == 0uL) {
x = 0uL;
break;
}
if ((r % (x + 23uL)) == 0uL) {
x = 0uL;
break;
}
if ((r % (x + 29uL)) == 0uL) {
x = 0uL;
break;
}
}
} while (x == 0uL);
if (r > 31uL) {
fwrite(&r, 1uL, sizeof(mp_digit), out);
printf("%9lu\r", r);
fflush(stdout);
}
if (r < 31uL) break;
}
if ((r % (x + 7)) == 0) {
x = 0;
break;
}
if ((r % (x + 11)) == 0) {
x = 0;
break;
}
if ((r % (x + 13)) == 0) {
x = 0;
break;
}
if ((r % (x + 17)) == 0) {
x = 0;
break;
}
if ((r % (x + 19)) == 0) {
x = 0;
break;
}
if ((r % (x + 23)) == 0) {
x = 0;
break;
}
if ((r % (x + 29)) == 0) {
x = 0;
break;
}
}
} while (x == 0);
if (r > 31) { fwrite(&r, 1, sizeof(mp_digit), out); printf("%9d\r", r); fflush(stdout); }
if (r < 31) break;
}
fclose(out);
fclose(out);
}
}
void load_tab(void)
static void load_tab(void)
{
primes = fopen("pprime.dat", "rb");
if (primes == NULL) {
gen_prime();
primes = fopen("pprime.dat", "rb");
}
fseek(primes, 0, SEEK_END);
fseek(primes, 0L, SEEK_END);
n_prime = ftell(primes) / sizeof(mp_digit);
}
mp_digit prime_digit(void)
static mp_digit prime_digit(void)
{
int n;
mp_digit d;
n = abs(rand()) % n_prime;
fseek(primes, n * sizeof(mp_digit), SEEK_SET);
fread(&d, 1, sizeof(mp_digit), primes);
fread(&d, 1uL, sizeof(mp_digit), primes);
return d;
}
/* makes a prime of at least k bits */
int
pprime (int k, int li, mp_int * p, mp_int * q)
static mp_err pprime(int k, int li, mp_int *p, mp_int *q)
{
mp_int a, b, c, n, x, y, z, v;
int res, ii;
static const mp_digit bases[] = { 2, 3, 5, 7, 11, 13, 17, 19 };
mp_int a, b, c, n, x, y, z, v;
mp_err res;
int ii;
static const mp_digit bases[] = { 2, 3, 5, 7, 11, 13, 17, 19 };
/* single digit ? */
if (k <= (int) DIGIT_BIT) {
mp_set (p, prime_digit ());
return MP_OKAY;
}
/* single digit ? */
if (k <= (int) MP_DIGIT_BIT) {
mp_set(p, prime_digit());
return MP_OKAY;
}
if ((res = mp_init (&c)) != MP_OKAY) {
return res;
}
if ((res = mp_init(&c)) != MP_OKAY) {
return res;
}
if ((res = mp_init (&v)) != MP_OKAY) {
goto LBL_C;
}
if ((res = mp_init(&v)) != MP_OKAY) {
goto LBL_C;
}
/* product of first 50 primes */
if ((res =
mp_read_radix (&v,
"19078266889580195013601891820992757757219839668357012055907516904309700014933909014729740190",
10)) != MP_OKAY) {
goto LBL_V;
}
/* product of first 50 primes */
if ((res =
mp_read_radix(&v,
"19078266889580195013601891820992757757219839668357012055907516904309700014933909014729740190",
10)) != MP_OKAY) {
goto LBL_V;
}
if ((res = mp_init (&a)) != MP_OKAY) {
goto LBL_V;
}
if ((res = mp_init(&a)) != MP_OKAY) {
goto LBL_V;
}
/* set the prime */
mp_set (&a, prime_digit ());
/* set the prime */
mp_set(&a, prime_digit());
if ((res = mp_init (&b)) != MP_OKAY) {
goto LBL_A;
}
if ((res = mp_init(&b)) != MP_OKAY) {
goto LBL_A;
}
if ((res = mp_init (&n)) != MP_OKAY) {
goto LBL_B;
}
if ((res = mp_init(&n)) != MP_OKAY) {
goto LBL_B;
}
if ((res = mp_init (&x)) != MP_OKAY) {
goto LBL_N;
}
if ((res = mp_init(&x)) != MP_OKAY) {
goto LBL_N;
}
if ((res = mp_init (&y)) != MP_OKAY) {
goto LBL_X;
}
if ((res = mp_init(&y)) != MP_OKAY) {
goto LBL_X;
}
if ((res = mp_init (&z)) != MP_OKAY) {
goto LBL_Y;
}
if ((res = mp_init(&z)) != MP_OKAY) {
goto LBL_Y;
}
/* now loop making the single digit */
while (mp_count_bits (&a) < k) {
fprintf (stderr, "prime has %4d bits left\r", k - mp_count_bits (&a));
fflush (stderr);
top:
mp_set (&b, prime_digit ());
/* now loop making the single digit */
while (mp_count_bits(&a) < k) {
fprintf(stderr, "prime has %4d bits left\r", k - mp_count_bits(&a));
fflush(stderr);
top:
mp_set(&b, prime_digit());
/* now compute z = a * b * 2 */
if ((res = mp_mul (&a, &b, &z)) != MP_OKAY) { /* z = a * b */
goto LBL_Z;
}
if ((res = mp_copy (&z, &c)) != MP_OKAY) { /* c = a * b */
goto LBL_Z;
}
if ((res = mp_mul_2 (&z, &z)) != MP_OKAY) { /* z = 2 * a * b */
goto LBL_Z;
}
/* n = z + 1 */
if ((res = mp_add_d (&z, 1, &n)) != MP_OKAY) { /* n = z + 1 */
goto LBL_Z;
}
/* check (n, v) == 1 */
if ((res = mp_gcd (&n, &v, &y)) != MP_OKAY) { /* y = (n, v) */
goto LBL_Z;
}
if (mp_cmp_d (&y, 1) != MP_EQ)
goto top;
/* now try base x=bases[ii] */
for (ii = 0; ii < li; ii++) {
mp_set (&x, bases[ii]);
/* compute x^a mod n */
if ((res = mp_exptmod (&x, &a, &n, &y)) != MP_OKAY) { /* y = x^a mod n */
goto LBL_Z;
/* now compute z = a * b * 2 */
if ((res = mp_mul(&a, &b, &z)) != MP_OKAY) { /* z = a * b */
goto LBL_Z;
}
/* if y == 1 loop */
if (mp_cmp_d (&y, 1) == MP_EQ)
continue;
/* now x^2a mod n */
if ((res = mp_sqrmod (&y, &n, &y)) != MP_OKAY) { /* y = x^2a mod n */
goto LBL_Z;
if ((res = mp_copy(&z, &c)) != MP_OKAY) { /* c = a * b */
goto LBL_Z;
}
if (mp_cmp_d (&y, 1) == MP_EQ)
continue;
/* compute x^b mod n */
if ((res = mp_exptmod (&x, &b, &n, &y)) != MP_OKAY) { /* y = x^b mod n */
goto LBL_Z;
if ((res = mp_mul_2(&z, &z)) != MP_OKAY) { /* z = 2 * a * b */
goto LBL_Z;
}
/* if y == 1 loop */
if (mp_cmp_d (&y, 1) == MP_EQ)
continue;
/* now x^2b mod n */
if ((res = mp_sqrmod (&y, &n, &y)) != MP_OKAY) { /* y = x^2b mod n */
goto LBL_Z;
/* n = z + 1 */
if ((res = mp_add_d(&z, 1uL, &n)) != MP_OKAY) { /* n = z + 1 */
goto LBL_Z;
}
if (mp_cmp_d (&y, 1) == MP_EQ)
continue;
/* compute x^c mod n == x^ab mod n */
if ((res = mp_exptmod (&x, &c, &n, &y)) != MP_OKAY) { /* y = x^ab mod n */
goto LBL_Z;
/* check (n, v) == 1 */
if ((res = mp_gcd(&n, &v, &y)) != MP_OKAY) { /* y = (n, v) */
goto LBL_Z;
}
/* if y == 1 loop */
if (mp_cmp_d (&y, 1) == MP_EQ)
continue;
if (mp_cmp_d(&y, 1uL) != MP_EQ)
goto top;
/* now compute (x^c mod n)^2 */
if ((res = mp_sqrmod (&y, &n, &y)) != MP_OKAY) { /* y = x^2ab mod n */
goto LBL_Z;
/* now try base x=bases[ii] */
for (ii = 0; ii < li; ii++) {
mp_set(&x, bases[ii]);
/* compute x^a mod n */
if ((res = mp_exptmod(&x, &a, &n, &y)) != MP_OKAY) { /* y = x^a mod n */
goto LBL_Z;
}
/* if y == 1 loop */
if (mp_cmp_d(&y, 1uL) == MP_EQ)
continue;
/* now x^2a mod n */
if ((res = mp_sqrmod(&y, &n, &y)) != MP_OKAY) { /* y = x^2a mod n */
goto LBL_Z;
}
if (mp_cmp_d(&y, 1uL) == MP_EQ)
continue;
/* compute x^b mod n */
if ((res = mp_exptmod(&x, &b, &n, &y)) != MP_OKAY) { /* y = x^b mod n */
goto LBL_Z;
}
/* if y == 1 loop */
if (mp_cmp_d(&y, 1uL) == MP_EQ)
continue;
/* now x^2b mod n */
if ((res = mp_sqrmod(&y, &n, &y)) != MP_OKAY) { /* y = x^2b mod n */
goto LBL_Z;
}
if (mp_cmp_d(&y, 1uL) == MP_EQ)
continue;
/* compute x^c mod n == x^ab mod n */
if ((res = mp_exptmod(&x, &c, &n, &y)) != MP_OKAY) { /* y = x^ab mod n */
goto LBL_Z;
}
/* if y == 1 loop */
if (mp_cmp_d(&y, 1uL) == MP_EQ)
continue;
/* now compute (x^c mod n)^2 */
if ((res = mp_sqrmod(&y, &n, &y)) != MP_OKAY) { /* y = x^2ab mod n */
goto LBL_Z;
}
/* y should be 1 */
if (mp_cmp_d(&y, 1uL) != MP_EQ)
continue;
break;
}
/* y should be 1 */
if (mp_cmp_d (&y, 1) != MP_EQ)
continue;
break;
}
/* no bases worked? */
if (ii == li)
goto top;
/* no bases worked? */
if (ii == li)
goto top;
{
char buf[4096];
{
char buf[4096];
mp_to_decimal(&n, buf, sizeof(buf));
printf("Certificate of primality for:\n%s\n\n", buf);
mp_to_decimal(&a, buf, sizeof(buf));
printf("A == \n%s\n\n", buf);
mp_to_decimal(&b, buf, sizeof(buf));
printf("B == \n%s\n\nG == %lu\n", buf, bases[ii]);
printf("----------------------------------------------------------------\n");
}
mp_toradix(&n, buf, 10);
printf("Certificate of primality for:\n%s\n\n", buf);
mp_toradix(&a, buf, 10);
printf("A == \n%s\n\n", buf);
mp_toradix(&b, buf, 10);
printf("B == \n%s\n\nG == %d\n", buf, bases[ii]);
printf("----------------------------------------------------------------\n");
}
/* a = n */
mp_copy(&n, &a);
}
/* a = n */
mp_copy (&n, &a);
}
/* get q to be the order of the large prime subgroup */
mp_sub_d(&n, 1uL, q);
mp_div_2(q, q);
mp_div(q, &b, q, NULL);
/* get q to be the order of the large prime subgroup */
mp_sub_d (&n, 1, q);
mp_div_2 (q, q);
mp_div (q, &b, q, NULL);
mp_exch(&n, p);
mp_exch (&n, p);
res = MP_OKAY;
LBL_Z:mp_clear (&z);
LBL_Y:mp_clear (&y);
LBL_X:mp_clear (&x);
LBL_N:mp_clear (&n);
LBL_B:mp_clear (&b);
LBL_A:mp_clear (&a);
LBL_V:mp_clear (&v);
LBL_C:mp_clear (&c);
return res;
res = MP_OKAY;
LBL_Z:
mp_clear(&z);
LBL_Y:
mp_clear(&y);
LBL_X:
mp_clear(&x);
LBL_N:
mp_clear(&n);
LBL_B:
mp_clear(&b);
LBL_A:
mp_clear(&a);
LBL_V:
mp_clear(&v);
LBL_C:
mp_clear(&c);
return res;
}
int
main (void)
int main(void)
{
mp_int p, q;
char buf[4096];
int k, li;
clock_t t1;
mp_int p, q;
char buf[4096];
int k, li;
clock_t t1;
srand (time (NULL));
load_tab();
srand(time(NULL));
load_tab();
printf ("Enter # of bits: \n");
fgets (buf, sizeof (buf), stdin);
sscanf (buf, "%d", &k);
printf("Enter # of bits: \n");
fgets(buf, sizeof(buf), stdin);
sscanf(buf, "%d", &k);
printf ("Enter number of bases to try (1 to 8):\n");
fgets (buf, sizeof (buf), stdin);
sscanf (buf, "%d", &li);
printf("Enter number of bases to try (1 to 8):\n");
fgets(buf, sizeof(buf), stdin);
sscanf(buf, "%d", &li);
mp_init (&p);
mp_init (&q);
mp_init(&p);
mp_init(&q);
t1 = clock ();
pprime (k, li, &p, &q);
t1 = clock () - t1;
t1 = clock();
pprime(k, li, &p, &q);
t1 = clock() - t1;
printf ("\n\nTook %ld ticks, %d bits\n", t1, mp_count_bits (&p));
printf("\n\nTook %d ticks, %d bits\n", t1, mp_count_bits(&p));
mp_toradix (&p, buf, 10);
printf ("P == %s\n", buf);
mp_toradix (&q, buf, 10);
printf ("Q == %s\n", buf);
mp_to_decimal(&p, buf, sizeof(buf));
printf("P == %s\n", buf);
mp_to_decimal(&q, buf, sizeof(buf));
printf("Q == %s\n", buf);
return 0;
return 0;
}
/* $Source: /cvs/libtom/libtommath/etc/pprime.c,v $ */
/* $Revision: 1.3 $ */
/* $Date: 2006/03/31 14:18:47 $ */

View File

@@ -34,4 +34,4 @@ _t_read:
sub eax,[time+4]
sbb edx,[time+0]
ret

View File

@@ -1,142 +1,542 @@
/* Tune the Karatsuba parameters
*
* Tom St Denis, tomstdenis@gmail.com
* Tom St Denis, tstdenis82@gmail.com
*/
#include <tommath.h>
#include "../tommath.h"
#include "../tommath_private.h"
#include <time.h>
#include <inttypes.h>
#include <errno.h>
/* how many times todo each size mult. Depends on your computer. For slow computers
* this can be low like 5 or 10. For fast [re: Athlon] should be 25 - 50 or so
*/
#define TIMES (1UL<<14UL)
/*
Please take in mind that both multiplicands are of the same size. The balancing
mechanism in mp_balance works well but has some overhead itself. You can test
the behaviour of it with the option "-o" followed by a (small) positive number 'x'
to generate ratios of the form 1:x.
*/
/* RDTSC from Scott Duplichan */
static ulong64 TIMFUNC (void)
{
#if defined __GNUC__
#if defined(__i386__) || defined(__x86_64__)
unsigned long long a;
__asm__ __volatile__ ("rdtsc\nmovl %%eax,%0\nmovl %%edx,4+%0\n"::"m"(a):"%eax","%edx");
return a;
#else /* gcc-IA64 version */
unsigned long result;
__asm__ __volatile__("mov %0=ar.itc" : "=r"(result) :: "memory");
while (__builtin_expect ((int) result == -1, 0))
__asm__ __volatile__("mov %0=ar.itc" : "=r"(result) :: "memory");
return result;
#endif
static uint64_t s_timer_function(void);
static void s_timer_start(void);
static uint64_t s_timer_stop(void);
static uint64_t s_time_mul(int size);
static uint64_t s_time_sqr(int size);
static void s_usage(char *s);
// Microsoft and Intel Windows compilers
#elif defined _M_IX86
__asm rdtsc
#elif defined _M_AMD64
return __rdtsc ();
#elif defined _M_IA64
#if defined __INTEL_COMPILER
#include <ia64intrin.h>
#endif
return __getReg (3116);
#else
#error need rdtsc function for this build
#endif
static uint64_t s_timer_function(void)
{
#if _POSIX_C_SOURCE >= 199309L
#define LTM_BILLION 1000000000
struct timespec ts;
/* TODO: Sets errno in case of error. Use? */
clock_gettime(CLOCK_MONOTONIC, &ts);
return (((uint64_t)ts.tv_sec) * LTM_BILLION + (uint64_t)ts.tv_nsec);
#else
clock_t t;
t = clock();
if (t < (clock_t)(0)) {
return (uint64_t)(0);
}
#ifndef X86_TIMER
return (uint64_t)(t);
#endif
}
/* generic ISO C timer */
ulong64 LBL_T;
void t_start(void) { LBL_T = TIMFUNC(); }
ulong64 t_read(void) { return TIMFUNC() - LBL_T; }
#else
extern void t_start(void);
extern ulong64 t_read(void);
#endif
ulong64 time_mult(int size, int s)
static uint64_t s_timer_tmp;
static void s_timer_start(void)
{
unsigned long x;
mp_int a, b, c;
ulong64 t1;
mp_init (&a);
mp_init (&b);
mp_init (&c);
mp_rand (&a, size);
mp_rand (&b, size);
if (s == 1) {
KARATSUBA_MUL_CUTOFF = size;
} else {
KARATSUBA_MUL_CUTOFF = 100000;
}
t_start();
for (x = 0; x < TIMES; x++) {
mp_mul(&a,&b,&c);
}
t1 = t_read();
mp_clear (&a);
mp_clear (&b);
mp_clear (&c);
return t1;
s_timer_tmp = s_timer_function();
}
static uint64_t s_timer_stop(void)
{
return s_timer_function() - s_timer_tmp;
}
ulong64 time_sqr(int size, int s)
static int s_check_result;
static int s_number_of_test_loops;
static int s_stabilization_extra;
static int s_offset = 1;
#define s_mp_mul(a, b, c) s_mp_mul_digs(a, b, c, (a)->used + (b)->used + 1)
static uint64_t s_time_mul(int size)
{
unsigned long x;
mp_int a, b;
ulong64 t1;
int x;
mp_err e;
mp_int a, b, c, d;
uint64_t t1;
mp_init (&a);
mp_init (&b);
if ((e = mp_init_multi(&a, &b, &c, &d, NULL)) != MP_OKAY) {
t1 = UINT64_MAX;
goto LTM_ERR;
}
mp_rand (&a, size);
if ((e = mp_rand(&a, size * s_offset)) != MP_OKAY) {
t1 = UINT64_MAX;
goto LTM_ERR;
}
if ((e = mp_rand(&b, size)) != MP_OKAY) {
t1 = UINT64_MAX;
goto LTM_ERR;
}
if (s == 1) {
KARATSUBA_SQR_CUTOFF = size;
} else {
KARATSUBA_SQR_CUTOFF = 100000;
}
s_timer_start();
for (x = 0; x < s_number_of_test_loops; x++) {
if ((e = mp_mul(&a,&b,&c)) != MP_OKAY) {
t1 = UINT64_MAX;
goto LTM_ERR;
}
if (s_check_result == 1) {
if ((e = s_mp_mul(&a,&b,&d)) != MP_OKAY) {
t1 = UINT64_MAX;
goto LTM_ERR;
}
if (mp_cmp(&c, &d) != MP_EQ) {
/* Time of 0 cannot happen (famous last words?) */
t1 = 0uLL;
goto LTM_ERR;
}
}
}
t_start();
for (x = 0; x < TIMES; x++) {
mp_sqr(&a,&b);
}
t1 = t_read();
mp_clear (&a);
mp_clear (&b);
return t1;
t1 = s_timer_stop();
LTM_ERR:
mp_clear_multi(&a, &b, &c, &d, NULL);
return t1;
}
int
main (void)
static uint64_t s_time_sqr(int size)
{
ulong64 t1, t2;
int x, y;
int x;
mp_err e;
mp_int a, b, c;
uint64_t t1;
for (x = 8; ; x += 2) {
t1 = time_mult(x, 0);
t2 = time_mult(x, 1);
printf("%d: %9llu %9llu, %9llu\n", x, t1, t2, t2 - t1);
if (t2 < t1) break;
}
y = x;
if ((e = mp_init_multi(&a, &b, &c, NULL)) != MP_OKAY) {
t1 = UINT64_MAX;
goto LTM_ERR;
}
for (x = 8; ; x += 2) {
t1 = time_sqr(x, 0);
t2 = time_sqr(x, 1);
printf("%d: %9llu %9llu, %9llu\n", x, t1, t2, t2 - t1);
if (t2 < t1) break;
}
printf("KARATSUBA_MUL_CUTOFF = %d\n", y);
printf("KARATSUBA_SQR_CUTOFF = %d\n", x);
if ((e = mp_rand(&a, size)) != MP_OKAY) {
t1 = UINT64_MAX;
goto LTM_ERR;
}
return 0;
s_timer_start();
for (x = 0; x < s_number_of_test_loops; x++) {
if ((e = mp_sqr(&a,&b)) != MP_OKAY) {
t1 = UINT64_MAX;
goto LTM_ERR;
}
if (s_check_result == 1) {
if ((e = s_mp_sqr(&a,&c)) != MP_OKAY) {
t1 = UINT64_MAX;
goto LTM_ERR;
}
if (mp_cmp(&c, &b) != MP_EQ) {
t1 = 0uLL;
goto LTM_ERR;
}
}
}
t1 = s_timer_stop();
LTM_ERR:
mp_clear_multi(&a, &b, &c, NULL);
return t1;
}
/* $Source: /cvs/libtom/libtommath/etc/tune.c,v $ */
/* $Revision: 1.3 $ */
/* $Date: 2006/03/31 14:18:47 $ */
struct tune_args {
int testmode;
int verbose;
int print;
int bncore;
int terse;
int upper_limit_print;
int increment_print;
} args;
static void s_run(const char *name, uint64_t (*op)(int), int *cutoff)
{
int x, count = 0;
uint64_t t1, t2;
if ((args.verbose == 1) || (args.testmode == 1)) {
printf("# %s.\n", name);
}
for (x = 8; x < args.upper_limit_print; x += args.increment_print) {
*cutoff = INT_MAX;
t1 = op(x);
if ((t1 == 0uLL) || (t1 == UINT64_MAX)) {
fprintf(stderr,"%s failed at x = INT_MAX (%s)\n", name,
(t1 == 0uLL)?"wrong result":"internal error");
exit(EXIT_FAILURE);
}
*cutoff = x;
t2 = op(x);
if ((t2 == 0uLL) || (t2 == UINT64_MAX)) {
fprintf(stderr,"%s failed (%s)\n", name,
(t2 == 0uLL)?"wrong result":"internal error");
exit(EXIT_FAILURE);
}
if (args.verbose == 1) {
printf("%d: %9"PRIu64" %9"PRIu64", %9"PRIi64"\n", x, t1, t2, (int64_t)t2 - (int64_t)t1);
}
if (t2 < t1) {
if (count == s_stabilization_extra) {
count = 0;
break;
} else if (count < s_stabilization_extra) {
count++;
}
} else if (count > 0) {
count--;
}
}
*cutoff = x - s_stabilization_extra * args.increment_print;
}
static long s_strtol(const char *str, char **endptr, const char *err)
{
const int base = 10;
char *_endptr;
long val;
errno = 0;
val = strtol(str, &_endptr, base);
if ((val > INT_MAX || val < 0) || (errno != 0)) {
fprintf(stderr, "Value %s not usable\n", str);
exit(EXIT_FAILURE);
}
if (_endptr == str) {
fprintf(stderr, "%s\n", err);
exit(EXIT_FAILURE);
}
if (endptr) *endptr = _endptr;
return val;
}
static int s_exit_code = EXIT_FAILURE;
static void s_usage(char *s)
{
fprintf(stderr,"Usage: %s [TvcpGbtrSLFfMmosh]\n",s);
fprintf(stderr," -T testmode, for use with testme.sh\n");
fprintf(stderr," -v verbose, print all timings\n");
fprintf(stderr," -c check results\n");
fprintf(stderr," -p print benchmark of final cutoffs in files \"multiplying\"\n");
fprintf(stderr," and \"squaring\"\n");
fprintf(stderr," -G [string] suffix for the filenames listed above\n");
fprintf(stderr," Implies '-p'\n");
fprintf(stderr," -b print benchmark of bncore.c\n");
fprintf(stderr," -t prints space (0x20) separated results\n");
fprintf(stderr," -r [64] number of rounds\n");
fprintf(stderr," -S [0xdeadbeef] seed for PRNG\n");
fprintf(stderr," -L [3] number of negative values accumulated until the result is accepted\n");
fprintf(stderr," -M [3000] upper limit of T-C tests/prints\n");
fprintf(stderr," -m [1] increment of T-C tests/prints\n");
fprintf(stderr," -o [1] multiplier for the second multiplicand\n");
fprintf(stderr," (Not for computing the cut-offs!)\n");
fprintf(stderr," -s 'preset' use values in 'preset' for printing.\n");
fprintf(stderr," 'preset' is a comma separated string with cut-offs for\n");
fprintf(stderr," ksm, kss, tc3m, tc3s in that order\n");
fprintf(stderr," ksm = karatsuba multiplication\n");
fprintf(stderr," kss = karatsuba squaring\n");
fprintf(stderr," tc3m = Toom-Cook 3-way multiplication\n");
fprintf(stderr," tc3s = Toom-Cook 3-way squaring\n");
fprintf(stderr," Implies '-p'\n");
fprintf(stderr," -h this message\n");
exit(s_exit_code);
}
struct cutoffs {
int KARATSUBA_MUL, KARATSUBA_SQR;
int TOOM_MUL, TOOM_SQR;
};
const struct cutoffs max_cutoffs =
{ INT_MAX, INT_MAX, INT_MAX, INT_MAX };
static void set_cutoffs(const struct cutoffs *c)
{
KARATSUBA_MUL_CUTOFF = c->KARATSUBA_MUL;
KARATSUBA_SQR_CUTOFF = c->KARATSUBA_SQR;
TOOM_MUL_CUTOFF = c->TOOM_MUL;
TOOM_SQR_CUTOFF = c->TOOM_SQR;
}
static void get_cutoffs(struct cutoffs *c)
{
c->KARATSUBA_MUL = KARATSUBA_MUL_CUTOFF;
c->KARATSUBA_SQR = KARATSUBA_SQR_CUTOFF;
c->TOOM_MUL = TOOM_MUL_CUTOFF;
c->TOOM_SQR = TOOM_SQR_CUTOFF;
}
int main(int argc, char **argv)
{
uint64_t t1, t2;
int x, i, j;
size_t n;
int printpreset = 0;
/*int preset[8];*/
char *endptr, *str;
uint64_t seed = 0xdeadbeef;
int opt;
struct cutoffs orig, updated;
FILE *squaring, *multiplying;
char mullog[256] = "multiplying";
char sqrlog[256] = "squaring";
s_number_of_test_loops = 64;
s_stabilization_extra = 3;
MP_ZERO_BUFFER(&args, sizeof(args));
args.testmode = 0;
args.verbose = 0;
args.print = 0;
args.bncore = 0;
args.terse = 0;
args.upper_limit_print = 3000;
args.increment_print = 1;
/* Very simple option parser, please treat it nicely. */
if (argc != 1) {
for (opt = 1; (opt < argc) && (argv[opt][0] == '-'); opt++) {
switch (argv[opt][1]) {
case 'T':
args.testmode = 1;
s_check_result = 1;
args.upper_limit_print = 1000;
args.increment_print = 11;
s_number_of_test_loops = 1;
s_stabilization_extra = 1;
s_offset = 1;
break;
case 'v':
args.verbose = 1;
break;
case 'c':
s_check_result = 1;
break;
case 'p':
args.print = 1;
break;
case 'G':
args.print = 1;
opt++;
if (opt >= argc) {
s_usage(argv[0]);
}
/* manual strcat() */
for (i = 0; i < 255; i++) {
if (mullog[i] == '\0') {
break;
}
}
for (j = 0; i < 255; j++, i++) {
mullog[i] = argv[opt][j];
if (argv[opt][j] == '\0') {
break;
}
}
for (i = 0; i < 255; i++) {
if (sqrlog[i] == '\0') {
break;
}
}
for (j = 0; i < 255; j++, i++) {
sqrlog[i] = argv[opt][j];
if (argv[opt][j] == '\0') {
break;
}
}
break;
case 'b':
args.bncore = 1;
break;
case 't':
args.terse = 1;
break;
case 'S':
opt++;
if (opt >= argc) {
s_usage(argv[0]);
}
str = argv[opt];
errno = 0;
seed = (uint64_t)s_strtol(argv[opt], NULL, "No seed given?\n");
break;
case 'L':
opt++;
if (opt >= argc) {
s_usage(argv[0]);
}
s_stabilization_extra = (int)s_strtol(argv[opt], NULL, "No value for option \"-L\"given");
break;
case 'o':
opt++;
if (opt >= argc) {
s_usage(argv[0]);
}
s_offset = (int)s_strtol(argv[opt], NULL, "No value for the offset given");
break;
case 'r':
opt++;
if (opt >= argc) {
s_usage(argv[0]);
}
s_number_of_test_loops = (int)s_strtol(argv[opt], NULL, "No value for the number of rounds given");
break;
case 'M':
opt++;
if (opt >= argc) {
s_usage(argv[0]);
}
args.upper_limit_print = (int)s_strtol(argv[opt], NULL, "No value for the upper limit of T-C tests given");
break;
case 'm':
opt++;
if (opt >= argc) {
s_usage(argv[0]);
}
args.increment_print = (int)s_strtol(argv[opt], NULL, "No value for the increment for the T-C tests given");
break;
case 's':
printpreset = 1;
args.print = 1;
opt++;
if (opt >= argc) {
s_usage(argv[0]);
}
str = argv[opt];
KARATSUBA_MUL_CUTOFF = (int)s_strtol(str, &endptr, "[1/4] No value for KARATSUBA_MUL_CUTOFF given");
str = endptr + 1;
KARATSUBA_SQR_CUTOFF = (int)s_strtol(str, &endptr, "[2/4] No value for KARATSUBA_SQR_CUTOFF given");
str = endptr + 1;
TOOM_MUL_CUTOFF = (int)s_strtol(str, &endptr, "[3/4] No value for TOOM_MUL_CUTOFF given");
str = endptr + 1;
TOOM_SQR_CUTOFF = (int)s_strtol(str, &endptr, "[4/4] No value for TOOM_SQR_CUTOFF given");
break;
case 'h':
s_exit_code = EXIT_SUCCESS;
/* FALLTHROUGH */
default:
s_usage(argv[0]);
}
}
}
/*
mp_rand uses the cryptographically secure
source of the OS by default. That is too expensive, too slow and
most important for a benchmark: it is not repeatable.
*/
s_mp_rand_jenkins_init(seed);
mp_rand_source(s_mp_rand_jenkins);
get_cutoffs(&orig);
updated = max_cutoffs;
if ((args.bncore == 0) && (printpreset == 0)) {
struct {
const char *name;
int *cutoff, *update;
uint64_t (*fn)(int);
} test[] = {
#define T_MUL_SQR(n, o, f) { #n, &o##_CUTOFF, &(updated.o), MP_HAS(S_MP_##o) ? f : NULL }
/*
The influence of the Comba multiplication cannot be
eradicated programmatically. It depends on the size
of the macro MP_WPARRAY in tommath.h which needs to
be changed manually (to 0 (zero)).
*/
T_MUL_SQR("Karatsuba multiplication", KARATSUBA_MUL, s_time_mul),
T_MUL_SQR("Karatsuba squaring", KARATSUBA_SQR, s_time_sqr),
T_MUL_SQR("Toom-Cook 3-way multiplying", TOOM_MUL, s_time_mul),
T_MUL_SQR("Toom-Cook 3-way squaring", TOOM_SQR, s_time_sqr),
#undef T_MUL_SQR
};
/* Turn all limits from bncore.c to the max */
set_cutoffs(&max_cutoffs);
for (n = 0; n < sizeof(test)/sizeof(test[0]); ++n) {
if (test[n].fn) {
s_run(test[n].name, test[n].fn, test[n].cutoff);
*test[n].update = *test[n].cutoff;
*test[n].cutoff = INT_MAX;
}
}
}
if (args.terse == 1) {
printf("%d %d %d %d\n",
updated.KARATSUBA_MUL,
updated.KARATSUBA_SQR,
updated.TOOM_MUL,
updated.TOOM_SQR);
} else {
printf("KARATSUBA_MUL_CUTOFF = %d\n", updated.KARATSUBA_MUL);
printf("KARATSUBA_SQR_CUTOFF = %d\n", updated.KARATSUBA_SQR);
printf("TOOM_MUL_CUTOFF = %d\n", updated.TOOM_MUL);
printf("TOOM_SQR_CUTOFF = %d\n", updated.TOOM_SQR);
}
if (args.print == 1) {
printf("Printing data for graphing to \"%s\" and \"%s\"\n",mullog, sqrlog);
multiplying = fopen(mullog, "w+");
if (multiplying == NULL) {
fprintf(stderr, "Opening file \"%s\" failed\n", mullog);
exit(EXIT_FAILURE);
}
squaring = fopen(sqrlog, "w+");
if (squaring == NULL) {
fprintf(stderr, "Opening file \"%s\" failed\n",sqrlog);
exit(EXIT_FAILURE);
}
for (x = 8; x < args.upper_limit_print; x += args.increment_print) {
set_cutoffs(&max_cutoffs);
t1 = s_time_mul(x);
set_cutoffs(&orig);
t2 = s_time_mul(x);
fprintf(multiplying, "%d: %9"PRIu64" %9"PRIu64", %9"PRIi64"\n", x, t1, t2, (int64_t)t2 - (int64_t)t1);
fflush(multiplying);
if (args.verbose == 1) {
printf("MUL %d: %9"PRIu64" %9"PRIu64", %9"PRIi64"\n", x, t1, t2, (int64_t)t2 - (int64_t)t1);
fflush(stdout);
}
set_cutoffs(&max_cutoffs);
t1 = s_time_sqr(x);
set_cutoffs(&orig);
t2 = s_time_sqr(x);
fprintf(squaring,"%d: %9"PRIu64" %9"PRIu64", %9"PRIi64"\n", x, t1, t2, (int64_t)t2 - (int64_t)t1);
fflush(squaring);
if (args.verbose == 1) {
printf("SQR %d: %9"PRIu64" %9"PRIu64", %9"PRIi64"\n", x, t1, t2, (int64_t)t2 - (int64_t)t1);
fflush(stdout);
}
}
printf("Finished. Data for graphing in \"%s\" and \"%s\"\n",mullog, sqrlog);
if (args.verbose == 1) {
set_cutoffs(&orig);
if (args.terse == 1) {
printf("%d %d %d %d\n",
KARATSUBA_MUL_CUTOFF,
KARATSUBA_SQR_CUTOFF,
TOOM_MUL_CUTOFF,
TOOM_SQR_CUTOFF);
} else {
printf("KARATSUBA_MUL_CUTOFF = %d\n", KARATSUBA_MUL_CUTOFF);
printf("KARATSUBA_SQR_CUTOFF = %d\n", KARATSUBA_SQR_CUTOFF);
printf("TOOM_MUL_CUTOFF = %d\n", TOOM_MUL_CUTOFF);
printf("TOOM_SQR_CUTOFF = %d\n", TOOM_SQR_CUTOFF);
}
}
}
exit(EXIT_SUCCESS);
}

View File

@@ -0,0 +1,107 @@
#!/bin/sh
die() {
echo "$1 failed"
echo "Exiting"
exit $2
}
# A linear congruential generator is sufficient for the purpose.
SEED=3735928559
LCG() {
SEED=$(((1103515245 * $SEED + 12345) % 2147483648))
echo $SEED
}
median() {
# read everything besides the header from file $1
# | cut-out the required column $2
# | sort all the entries numerically
# | show only the first $3 entries
# | show only the last entry
tail -n +2 $1 | cut -d' ' -f$2 | sort -n | head -n $3 | tail -n 1
}
MPWD=$(dirname $(readlink -f "$0"))
FILE_NAME="tuning_list"
TOMMATH_CUTOFFS_H="$MPWD/../tommath_cutoffs.h"
BACKUP_SUFFIX=".orig"
RNUM=0
#############################################################################
# It would be a good idea to isolate these processes (with e.g.: cpuset) #
# #
# It is not a good idea to e.g: watch high resolution videos while this #
# test are running if you do not have enough memory to avoid page faults. #
#############################################################################
# Number of rounds overall.
LIMIT=100
# Number of loops for each input.
RLOOPS=10
# Offset ( > 0 ) . Runs tests with asymmetric input of the form 1:OFFSET
# Please use another destination for TOMMATH_CUTOFFS_H if you change OFFSET, because the numbers
# with an offset different from 1 (one) are not usable as the general cut-off values
# in "tommath_cutoffs.h".
OFFSET=1
# Number ( >= 3 ) of positive results (TC-is-faster) accumulated until it is accepted.
# Due to the algorithm used to compute the median in this Posix compliant shell script
# the value needs to be 3 (three), not less, to keep the variation small.
LAG=3
# Keep the temporary file $FILE_NAME. Set to 0 (zero) to remove it at the end.
# The file is in a format fit to feed into R directly. If you do it and find the median
# of this program to be off by more than a couple: please contact the authors and report
# the numbers from this program and R and the standard deviation. This program is known
# to get larger errors if the standard deviation is larger than ~50.
KEEP_TEMP=1
echo "You might like to watch the numbers go up to $LIMIT but it will take a long time!"
# Might not have sufficient rights or disc full.
echo "km ks tc3m tc3s" > $FILE_NAME || die "Writing header to $FILE_NAME" $?
i=1
while [ $i -le $LIMIT ]; do
RNUM=$(LCG)
printf "\r%d" $i
"$MPWD"/tune -t -r $RLOOPS -L $LAG -S "$RNUM" -o $OFFSET >> $FILE_NAME || die "tune" $?
i=$((i + 1))
done
if [ $KEEP_TEMP -eq 0 ]; then
rm -v $FILE_NAME || die "Removing $KEEP_TEMP" $?
fi
echo "Writing cut-off values to \"$TOMMATH_CUTOFFS_H\"."
echo "In case of failure: a copy of \"$TOMMATH_CUTOFFS_H\" is in \"$TOMMATH_CUTOFFS_H$BACKUP_SUFFIX\""
cp -v $TOMMATH_CUTOFFS_H $TOMMATH_CUTOFFS_H$BACKUP_SUFFIX || die "Making backup copy of $TOMMATH_CUTOFFS_H" $?
cat << END_OF_INPUT > $TOMMATH_CUTOFFS_H || die "Writing header to $TOMMATH_CUTOFFS_H" $?
/* LibTomMath, multiple-precision integer library -- Tom St Denis */
/* SPDX-License-Identifier: Unlicense */
/*
Current values evaluated on an AMD A8-6600K (64-bit).
Type "make tune" to optimize them for your machine but
be aware that it may take a long time. It took 2:30 minutes
on the aforementioned machine for example.
*/
END_OF_INPUT
# The Posix shell does not offer an array data type so we create
# the median with 'standard tools'^TM
# read the file (without the first line) and count the lines
i=$(tail -n +2 $FILE_NAME | wc -l)
# our median point will be at $i entries
i=$(( (i / 2) + 1 ))
TMP=$(median $FILE_NAME 1 $i)
echo "#define MP_DEFAULT_KARATSUBA_MUL_CUTOFF $TMP"
echo "#define MP_DEFAULT_KARATSUBA_MUL_CUTOFF $TMP" >> $TOMMATH_CUTOFFS_H || die "(km) Appending to $TOMMATH_CUTOFFS_H" $?
TMP=$(median $FILE_NAME 2 $i)
echo "#define MP_DEFAULT_KARATSUBA_SQR_CUTOFF $TMP"
echo "#define MP_DEFAULT_KARATSUBA_SQR_CUTOFF $TMP" >> $TOMMATH_CUTOFFS_H || die "(ks) Appending to $TOMMATH_CUTOFFS_H" $?
TMP=$(median $FILE_NAME 3 $i)
echo "#define MP_DEFAULT_TOOM_MUL_CUTOFF $TMP"
echo "#define MP_DEFAULT_TOOM_MUL_CUTOFF $TMP" >> $TOMMATH_CUTOFFS_H || die "(tc3m) Appending to $TOMMATH_CUTOFFS_H" $?
TMP=$(median $FILE_NAME 4 $i)
echo "#define MP_DEFAULT_TOOM_SQR_CUTOFF $TMP"
echo "#define MP_DEFAULT_TOOM_SQR_CUTOFF $TMP" >> $TOMMATH_CUTOFFS_H || die "(tc3s) Appending to $TOMMATH_CUTOFFS_H" $?