3
0
Fork 0

added tomsfastmath-0.01

This commit is contained in:
Tom St Denis 2004-08-25 02:43:43 +00:00 committed by Steffen Jaeckel
commit 5e92ed2a59
75 changed files with 11069 additions and 0 deletions

7
LICENSE Normal file
View File

@ -0,0 +1,7 @@
TomsFastMath is public domain.
Note some ideas were borrowed from LibTomMath and OpenSSL. All of the code is original or ported
from LibTomMath [no code was ported from OpenSSL]. As such the origins and status of this code
are both public domain.
-- Tom St Denis

5
SPONSORS Normal file
View File

@ -0,0 +1,5 @@
Development of TomsFastMath was sponsored by three groups. Two companies that use LTC and LTM commercially
and one individual who decided he wanted to help out by being generous.
Thanks goes to them [though they wished to remain anonymous] and people like them.

6
TODO Normal file
View File

@ -0,0 +1,6 @@
1. Write more documentation ;-)
2. Ports to PPC and MIPS
3. Fix any lingering bugs, add additional requested functionality.
NOTE: The library is still fairly new. I've tested it quite a bit but that doesn't mean surprises
can't happen. Please test the results you get for correctness.

2
changes.txt Normal file
View File

@ -0,0 +1,2 @@
August 25th, 2004
TFM 0.01 -- Initial Release

50
comba_mult_gen.c Normal file
View File

@ -0,0 +1,50 @@
/* program emits a NxN comba multiplier */
#include <stdio.h>
int main(int argc, char **argv)
{
int N, x, y, z;
N = atoi(argv[1]);
/* print out preamble */
printf(
"void fp_mul_comba%d(fp_int *A, fp_int *B, fp_int *C)\n"
"{\n"
" fp_word t;\n"
" fp_digit c0, c1, c2, at[%d];\n"
"\n"
" memcpy(at, A->dp, %d * sizeof(fp_digit));\n"
" memcpy(at+%d, B->dp, %d * sizeof(fp_digit));\n"
" COMBA_START;\n"
"\n"
" COMBA_CLEAR;\n", N, N+N, N, N, N, N);
/* now do the rows */
for (x = 0; x < (N+N-1); x++) {
printf(
" /* %d */\n", x);
if (x > 0) {
printf(
" COMBA_FORWARD;\n");
}
for (y = 0; y < N; y++) {
for (z = 0; z < N; z++) {
if ((y+z)==x) {
printf(" MULADD(at[%d], at[%d]); ", y, z+N);
}
}
}
printf(
"\n"
" COMBA_STORE(C->dp[%d]);\n", x);
}
printf(
" COMBA_STORE2(C->dp[%d]);\n"
" C->used = %d;\n"
" C->sign = A->sign ^ B->sign;\n"
" fp_clamp(C);\n"
" COMBA_FINI;\n"
"}\n\n\n", N+N-1, N+N, N+N);
return 0;
}

54
comba_sqr_gen.c Normal file
View File

@ -0,0 +1,54 @@
/* Generates squaring comba code... it learns it knows our secrets! */
#include <stdio.h>
int main(int argc, char **argv)
{
int x, y, z, N;
N = atoi(argv[1]);
printf(
"void fp_sqr_comba%d(fp_int *A, fp_int *B)\n"
"{\n"
" fp_word t;\n"
" fp_digit *a, b[%d], c0, c1, c2;\n"
"\n"
" a = A->dp;\n"
" COMBA_START; \n"
"\n"
" /* clear carries */\n"
" CLEAR_CARRY;\n"
"\n"
" /* output 0 */\n"
" SQRADD(a[0],a[0]);\n"
" COMBA_STORE(b[0]);\n", N, N+N);
for (x = 1; x < N+N-1; x++) {
printf(
"\n /* output %d */\n"
" CARRY_FORWARD;\n ", x);
for (y = 0; y < N; y++) {
for (z = 0; z < N; z++) {
if (y<=z && (y+z)==x) {
if (y == z) {
printf("SQRADD(a[%d], a[%d]); ", y, y);
} else {
printf("SQRADD2(a[%d], a[%d]); ", y, z);
}
}
}
}
printf("\n COMBA_STORE(b[%d]);\n", x);
}
printf(" COMBA_STORE2(b[%d]);\n", N+N-1);
printf(
" COMBA_FINI;\n"
"\n"
" B->used = %d;\n"
" B->sign = FP_ZPOS;\n"
" memcpy(B->dp, b, %d * sizeof(fp_digit));\n"
" fp_clamp(B);\n"
"}\n\n\n", N+N, N+N);
return 0;
}

144
demo/stest.c Normal file
View File

@ -0,0 +1,144 @@
/* A simple static test program. */
#include <tfm.h>
#ifdef GBA_MODE
#include <gba.h>
#define DISPLAY(x) modetxt_puts(vfb, x, 1)
#endif
#ifndef DISPLAY
#define DISPLAY(x) printf(x)
#endif
#ifdef GBA_MODE
int c_main(void)
#else
int main(void)
#endif
{
fp_int a,b,c,d,e,f;
fp_digit dp;
fp_init(&a);
fp_init(&b);
fp_init(&c);
fp_init(&d);
fp_init(&e);
fp_init(&f);
#ifdef GBA_MODE
install_common();
modetxt_init();
modetxt_gotoxy(0,0);
#endif
/* test multiplication */
fp_read_radix(&a, "3453534534535345345341230891273", 10);
fp_read_radix(&b, "2394873294871238934718923" , 10);
fp_read_radix(&c, "8270777629674273015508507050766235312931312159028658979", 10);
fp_mul(&a, &b, &d);
if (fp_cmp(&c, &d)) {
DISPLAY("mul failed\n");
return 0;
} else {
DISPLAY("mul passed\n");
}
/* test multiplication */
fp_read_radix(&a, "30481290320498235987349712308523652378643912563478232907782361237864278207235782364578264891274789264278634289739", 10);
fp_read_radix(&b, "48761478126387263782638276327836287632836278362837627838736278362923698724823749238732" , 10);
fp_read_radix(&c, "1486312771227034563307950634490737985563993459700941115664257275795366623795590136120579100118233580357115074068815507257715906295105536107921754177810976863679300283932188006885811950341132768970948", 10);
fp_mul(&a, &b, &d);
if (fp_cmp(&c, &d)) {
DISPLAY("mul failed\n");
return 0;
} else {
DISPLAY("mul passed\n");
}
/* test multiplication */
fp_read_radix(&a, "115792089237316195423570985008687907853269984665640564039457584007913129639935", 10);
fp_read_radix(&b, "174224571863520493293247799005065324265471" , 10);
fp_read_radix(&c, "20173827172553973356686868531273530268200710714389071377794102651988800859098544338487575161443744102709980552583184385", 10);
fp_mul(&a, &b, &d);
if (fp_cmp(&c, &d)) {
DISPLAY("mul failed\n");
return 0;
} else {
DISPLAY("mul passed\n");
}
/* test squaring */
fp_read_radix(&a, "298723982748923478923473927489237289347238947238947238947238972893", 10);
fp_read_radix(&b, "89236017869379132235512787068367546521309689412262624434964313994127411682542855190667724226920696163962644836740110835385588789449" , 10);
fp_sqr(&a, &c);
if (fp_cmp(&c, &b)) {
DISPLAY("sqr failed\n");
return 0;
} else {
DISPLAY("sqr passed\n");
}
fp_read_radix(&a, "397823894238973128942895123894327123941724927848927349274897238978927593487012378490184789429812734982738972389", 10);
fp_read_radix(&b, "158263850827461677491961439999264901067636282938352531932899298293270945997930087353471903166601507321298827087008336951419604640736464667188494668962822678461626245753696845719301945679092882499787869509090904187704367321" , 10);
fp_sqr(&a, &c);
if (fp_cmp(&c, &b)) {
DISPLAY("sqr failed\n");
return 0;
} else {
DISPLAY("sqr passed\n");
}
fp_read_radix(&a, "13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084095", 10);
fp_read_radix(&b, "179769313486231590772930519078902473361797697894230657273430081157732675805500963132708477322407536021120113879871393357658789768814416622492847430639474097562152033539671286128252223189553839160721441767298250321715263238814402734379959506792230903356495130620869925267845538430714092411695463462326211969025" , 10);
fp_sqr(&a, &c);
if (fp_cmp(&c, &b)) {
DISPLAY("sqr failed\n");
return 0;
} else {
DISPLAY("sqr passed\n");
}
/* montgomery reductions */
fp_read_radix(&a, "234892374892374893489123428937892781237863278637826327367637836278362783627836783678363", 10);
fp_read_radix(&b, "4447823492749823749234123489273987393983289319382762756425425425642727352327452374521", 10);
fp_read_radix(&c, "2396271882990732698083317035605836523697277786556053771759862552557086442129695099100", 10);
fp_montgomery_setup(&b, &dp);
fp_montgomery_reduce(&a, &b, dp);
if (fp_cmp(&a, &c)) {
DISPLAY("mont failed\n");
return 0;
} else {
DISPLAY("mont passed\n");
}
fp_read_radix(&a, "2348923748923748934891234456645654645645684576353428937892781237863278637826327367637836278362783627836783678363", 10);
fp_read_radix(&b, "444782349274982374923412348927398739398328931938276275642542542564272735232745237452123424324324444121111119", 10);
fp_read_radix(&c, "45642613844554582908652603086180267403823312390990082328515008314514368668691233331246183943400359349283420", 10);
fp_montgomery_setup(&b, &dp);
fp_montgomery_reduce(&a, &b, dp);
if (fp_cmp(&a, &c)) {
DISPLAY("mont failed\n");
return 0;
} else {
DISPLAY("mont passed\n");
}
fp_read_radix(&a, "234823424242342923748923748934891234456645654645645684576353424972378234762378623891236834132352375235378462378489378927812378632786378263273676378362783627555555555539568389052478124618461834763837685723645827529034853490580134568947341278498542893481762349723907847892983627836783678363", 10);
fp_read_radix(&b, "44478234927456563455982374923412348927398739398328931938276275642485623481638279025465891276312903262837562349056234783648712314678120389173890128905425242424239784256427", 10);
fp_read_radix(&c, "33160865265453361650564031464519042126185632333462754084489985719613480783282357410514898819797738034600484519472656152351777186694609218202276509271061460265488348645081", 10);
fp_montgomery_setup(&b, &dp);
fp_montgomery_reduce(&a, &b, dp);
if (fp_cmp(&a, &c)) {
DISPLAY("mont failed\n");
return 0;
} else {
DISPLAY("mont passed\n");
}
return 0;
}

537
demo/test.c Normal file
View File

@ -0,0 +1,537 @@
/* TFM demo program */
#include <tfm.h>
void draw(fp_int *a)
{
int x;
printf("%d, %d, ", a->used, a->sign);
for (x = a->used - 1; x >= 0; x--) {
printf("%08lx ", a->dp[x]);
}
printf("\n");
}
int myrng(unsigned char *dst, int len, void *dat)
{
int x;
for (x = 0; x < len; x++) dst[x] = rand() & 0xFF;
return len;
}
/* 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
// 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
}
char cmd[4096], buf[4096];
int main(void)
{
fp_int a,b,c,d,e,f;
fp_digit fp;
int n, err;
unsigned long expt_n, add_n, sub_n, mul_n, div_n, sqr_n, mul2d_n, div2d_n, gcd_n, lcm_n, inv_n,
div2_n, mul2_n, add_d_n, sub_d_n, mul_d_n, t, cnt, rr, ix;
ulong64 t1, t2;
fp_zero(&b); fp_zero(&c); fp_zero(&d); fp_zero(&e); fp_zero(&f);
fp_zero(&a); draw(&a);
/* test set and simple shifts */
printf("Testing mul/div 2\n");
fp_set(&a, 1); draw(&a);
for (n = 0; n <= DIGIT_BIT; n++) {
fp_mul_2(&a, &a); printf("(%d) ", fp_count_bits(&a));
draw(&a);
}
for (n = 0; n <= (DIGIT_BIT + 1); n++) {
fp_div_2(&a, &a);
draw(&a);
}
fp_set(&a, 1);
/* test lshd/rshd */
printf("testing lshd/rshd\n");
fp_lshd(&a, 3); draw(&a);
fp_rshd(&a, 3); draw(&a);
/* test more complicated shifts */
printf("Testing mul/div 2d\n");
fp_mul_2d(&a, DIGIT_BIT/2, &a); draw(&a);
fp_div_2d(&a, DIGIT_BIT/2, &a, NULL); draw(&a);
fp_mul_2d(&a, DIGIT_BIT + DIGIT_BIT/2, &a); draw(&a);
fp_div_2d(&a, DIGIT_BIT + DIGIT_BIT/2, &a, NULL); draw(&a);
/* test neg/abs */
printf("testing neg/abs\n");
fp_neg(&a, &a); draw(&a);
fp_neg(&a, &a); draw(&a);
fp_neg(&a, &a); draw(&a);
fp_abs(&a, &a); draw(&a);
/* test comparisons */
fp_set(&b, 3);
fp_set(&c, 4); fp_neg(&c, &c);
fp_set(&d, 1);
printf("Testing compares\n%d, %d, %d, %d\n", fp_cmp(&a, &b), fp_cmp(&a, &c), fp_cmp(&a, &d), fp_cmp(&b, &c));
/* test add/sub */
printf("Testing add/sub \n");
fp_set(&a, ((fp_digit)1)<<(DIGIT_BIT-1)); draw(&a);
fp_set(&b, ((fp_digit)1)<<(DIGIT_BIT-2));
fp_add(&a, &b, &a); draw(&a);
fp_add(&a, &b, &a); draw(&a);
fp_add(&a, &b, &a); draw(&a);
printf("sub...\n");
printf("cmp returns: %d, ", fp_cmp(&a, &b)); fp_sub(&a, &b, &a); draw(&a);
printf("cmp returns: %d, ", fp_cmp(&a, &b)); fp_sub(&a, &b, &a); draw(&a);
printf("cmp returns: %d, ", fp_cmp(&a, &b)); fp_sub(&a, &b, &a); draw(&a);
printf("cmp returns: %d, ", fp_cmp(&a, &b)); fp_sub(&a, &b, &a); draw(&a);
printf("cmp returns: %d, ", fp_cmp(&a, &b)); fp_sub(&a, &b, &a); draw(&a);
printf("cmp returns: %d, ", fp_cmp(&a, &b)); fp_sub(&a, &b, &a); draw(&a);
/* test mul_d */
printf("Testing mul_d and div_d\n");
fp_set(&a, 1);
fp_mul_d(&a, ((fp_digit)1)<<(DIGIT_BIT/2), &a); draw(&a);
fp_mul_d(&a, ((fp_digit)1)<<(DIGIT_BIT/2), &a); draw(&a);
fp_mul_d(&a, ((fp_digit)1)<<(DIGIT_BIT/2), &a); draw(&a);
printf("div_d\n");
fp_div_d(&a, ((fp_digit)1)<<(DIGIT_BIT/2), &a, NULL); draw(&a);
fp_div_d(&a, ((fp_digit)1)<<(DIGIT_BIT/2), &a, NULL); draw(&a);
fp_div_d(&a, ((fp_digit)1)<<(DIGIT_BIT/2), &a, NULL); draw(&a);
/* testing read radix */
printf("Testing read_radix\n");
fp_read_radix(&a, "123456789012345678901234567890", 16); draw(&a);
/* test mont */
printf("Montgomery test\n");
fp_set(&a, 1);
fp_lshd(&a, 4);
fp_add_d(&a, 1, &a);
fp_montgomery_setup(&a, &fp);
fp_montgomery_calc_normalization(&b, &a);
fp_read_radix(&d, "123456789123", 16);
for (n = 0; n < 100000; n++) {
fp_add_d(&d, 1, &d); fp_sqrmod(&d, &a, &d);
fp_mul(&d, &b, &c);
fp_montgomery_reduce(&c, &a, fp);
if (fp_cmp(&c, &d) != FP_EQ) {
printf("Failed mont %d\n", n);
draw(&a);
draw(&d);
draw(&c);
return EXIT_FAILURE;
}
}
printf("Passed.\n");
/* test for size */
for (ix = 8*DIGIT_BIT; ix < 10*DIGIT_BIT; ix++) {
printf("Testing (not safe-prime): %9d bits \r", ix); fflush(stdout);
err = fp_prime_random_ex(&a, 8, ix, (rand()&1)?TFM_PRIME_2MSB_OFF:TFM_PRIME_2MSB_ON, myrng, NULL);
if (err != FP_OKAY) {
printf("failed with err code %d\n", err);
return EXIT_FAILURE;
}
if (fp_count_bits(&a) != ix) {
printf("Prime is %d not %d bits!!!\n", fp_count_bits(&a), ix);
return EXIT_FAILURE;
}
}
printf("\n\n");
#if 0
/* do some timings... */
printf("Addition:\n");
for (t = 2; t <= FP_SIZE/2; t += 2) {
fp_zero(&a);
fp_zero(&b);
fp_zero(&c);
for (ix = 0; ix < t; ix++) {
a.dp[ix] = ix;
b.dp[ix] = ix;
}
a.used = t;
b.used = t;
t2 = -1;
for (ix = 0; ix < 2500; ++ix) {
t1 = TIMFUNC();
fp_add(&a, &b, &c); fp_add(&a, &b, &c);
fp_add(&a, &b, &c); fp_add(&a, &b, &c);
fp_add(&a, &b, &c); fp_add(&a, &b, &c);
fp_add(&a, &b, &c); fp_add(&a, &b, &c);
t2 = (TIMFUNC() - t1)>>3;
if (t1<t2) { --ix; t2 = t1; }
}
printf("%5lu-bit: %9llu\n", t * DIGIT_BIT, t2);
}
printf("Multiplication:\n");
for (t = 2; t <= FP_SIZE/2; t += 2) {
fp_zero(&a);
fp_zero(&b);
fp_zero(&c);
for (ix = 0; ix < t; ix++) {
a.dp[ix] = ix;
b.dp[ix] = ix;
}
a.used = t;
b.used = t;
t2 = -1;
for (ix = 0; ix < 10000; ++ix) {
t1 = TIMFUNC();
fp_mul(&a, &b, &c); fp_mul(&a, &b, &c);
fp_mul(&a, &b, &c); fp_mul(&a, &b, &c);
t2 = (TIMFUNC() - t1)>>2;
if (t1<t2) { --ix; t2 = t1; }
}
printf("%5lu-bit: %9llu\n", t * DIGIT_BIT, t2);
}
//#else
printf("Squaring:\n");
for (t = 2; t <= FP_SIZE/2; t += 2) {
fp_zero(&a);
fp_zero(&b);
for (ix = 0; ix < t; ix++) {
a.dp[ix] = ix;
}
a.used = t;
t2 = -1;
for (ix = 0; ix < 10000; ++ix) {
t1 = TIMFUNC();
fp_sqr(&a, &b); fp_sqr(&a, &b);
fp_sqr(&a, &b); fp_sqr(&a, &b);
t2 = (TIMFUNC() - t1)>>2;
if (t1<t2) { --ix; t2 = t1; }
}
printf("%5lu-bit: %9llu\n", t * DIGIT_BIT, t2);
}
//#else
printf("Montgomery:\n");
for (t = 2; t <= (FP_SIZE/2)-2; t += 2) {
fp_zero(&a);
for (ix = 0; ix < t; ix++) {
a.dp[ix] = ix | 1;
}
a.used = t;
fp_montgomery_setup(&a, &fp);
fp_sub_d(&a, 3, &b);
fp_sqr(&b, &b);
fp_copy(&b, &c);
fp_copy(&b, &d);
t2 = -1;
for (ix = 0; ix < 10000; ++ix) {
t1 = TIMFUNC();
fp_montgomery_reduce(&c, &a, &fp);
fp_montgomery_reduce(&d, &a, &fp);
t2 = (TIMFUNC() - t1)>>1;
fp_copy(&b, &c);
fp_copy(&b, &d);
if (t1<t2) { --ix; t2 = t1; }
}
printf("%5lu-bit: %9llu\n", t * DIGIT_BIT, t2);
}
//#else
printf("Exptmod:\n");
for (t = 512/DIGIT_BIT; t <= (FP_SIZE/2)-2; t += t) {
fp_zero(&a);
fp_zero(&b);
fp_zero(&c);
for (ix = 0; ix < t; ix++) {
a.dp[ix] = ix+1;
b.dp[ix] = (fp_digit)rand() * (fp_digit)rand();
c.dp[ix] = ix;
}
a.used = t;
b.used = t;
c.used = t;
t2 = -1;
for (ix = 0; ix < 50; ++ix) {
t1 = TIMFUNC();
fp_exptmod(&c, &b, &a, &d);
fp_exptmod(&c, &b, &a, &d);
t2 = (TIMFUNC() - t1)>>1;
fp_copy(&b, &c);
fp_copy(&b, &d);
if (t1<t2) { t2 = t1; --ix; }
}
printf("%5lu-bit: %9llu\n", t * DIGIT_BIT, t2);
}
#endif
div2_n = mul2_n = inv_n = expt_n = lcm_n = gcd_n = add_n =
sub_n = mul_n = div_n = sqr_n = mul2d_n = div2d_n = cnt = add_d_n = sub_d_n= mul_d_n = 0;
for (;;) {
printf("%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu/%4lu ", add_n, sub_n, mul_n, div_n, sqr_n, mul2d_n, div2d_n, gcd_n, lcm_n, expt_n, inv_n, div2_n, mul2_n, add_d_n, sub_d_n, mul_d_n);
fgets(cmd, 4095, stdin);
cmd[strlen(cmd)-1] = 0;
printf("%s ]\r",cmd); fflush(stdout);
if (!strcmp(cmd, "mul2d")) { ++mul2d_n;
fgets(buf, 4095, stdin); fp_read_radix(&a, buf, 64);
fgets(buf, 4095, stdin); sscanf(buf, "%d", &rr);
fgets(buf, 4095, stdin); fp_read_radix(&b, buf, 64);
fp_mul_2d(&a, rr, &a);
a.sign = b.sign;
if (fp_cmp(&a, &b) != FP_EQ) {
printf("mul2d failed, rr == %d\n",rr);
draw(&a);
draw(&b);
return 0;
}
} else if (!strcmp(cmd, "div2d")) { ++div2d_n;
fgets(buf, 4095, stdin); fp_read_radix(&a, buf, 64);
fgets(buf, 4095, stdin); sscanf(buf, "%d", &rr);
fgets(buf, 4095, stdin); fp_read_radix(&b, buf, 64);
fp_div_2d(&a, rr, &a, &e);
a.sign = b.sign;
if (a.used == b.used && a.used == 0) { a.sign = b.sign = FP_ZPOS; }
if (fp_cmp(&a, &b) != FP_EQ) {
printf("div2d failed, rr == %d\n",rr);
draw(&a);
draw(&b);
return 0;
}
} else if (!strcmp(cmd, "add")) { ++add_n;
fgets(buf, 4095, stdin); fp_read_radix(&a, buf, 64);
fgets(buf, 4095, stdin); fp_read_radix(&b, buf, 64);
fgets(buf, 4095, stdin); fp_read_radix(&c, buf, 64);
fp_copy(&a, &d);
fp_add(&d, &b, &d);
if (fp_cmp(&c, &d) != FP_EQ) {
printf("add %lu failure!\n", add_n);
draw(&a);draw(&b);draw(&c);draw(&d);
return 0;
}
/* test the sign/unsigned storage functions */
rr = fp_signed_bin_size(&c);
fp_to_signed_bin(&c, (unsigned char *)cmd);
memset(cmd+rr, rand()&255, sizeof(cmd)-rr);
fp_read_signed_bin(&d, (unsigned char *)cmd, rr);
if (fp_cmp(&c, &d) != FP_EQ) {
printf("fp_signed_bin failure!\n");
draw(&c);
draw(&d);
return 0;
}
rr = fp_unsigned_bin_size(&c);
fp_to_unsigned_bin(&c, (unsigned char *)cmd);
memset(cmd+rr, rand()&255, sizeof(cmd)-rr);
fp_read_unsigned_bin(&d, (unsigned char *)cmd, rr);
if (fp_cmp_mag(&c, &d) != FP_EQ) {
printf("fp_unsigned_bin failure!\n");
draw(&c);
draw(&d);
return 0;
}
} else if (!strcmp(cmd, "sub")) { ++sub_n;
fgets(buf, 4095, stdin); fp_read_radix(&a, buf, 64);
fgets(buf, 4095, stdin); fp_read_radix(&b, buf, 64);
fgets(buf, 4095, stdin); fp_read_radix(&c, buf, 64);
fp_copy(&a, &d);
fp_sub(&d, &b, &d);
if (fp_cmp(&c, &d) != FP_EQ) {
printf("sub %lu failure!\n", sub_n);
draw(&a);draw(&b);draw(&c);draw(&d);
return 0;
}
} else if (!strcmp(cmd, "mul")) {
fgets(buf, 4095, stdin); fp_read_radix(&a, buf, 64);
fgets(buf, 4095, stdin); fp_read_radix(&b, buf, 64);
fgets(buf, 4095, stdin); fp_read_radix(&c, buf, 64);
//continue;
fp_copy(&a, &d);
fp_mul(&d, &b, &d); ++mul_n;
if (fp_cmp(&c, &d) != FP_EQ) {
printf("mul %lu failure!\n", mul_n);
draw(&a);draw(&b);draw(&c);draw(&d);
return 0;
}
} else if (!strcmp(cmd, "div")) {
fgets(buf, 4095, stdin); fp_read_radix(&a, buf, 64);
fgets(buf, 4095, stdin); fp_read_radix(&b, buf, 64);
fgets(buf, 4095, stdin); fp_read_radix(&c, buf, 64);
fgets(buf, 4095, stdin); fp_read_radix(&d, buf, 64);
// continue;
fp_div(&a, &b, &e, &f); ++div_n;
if (fp_cmp(&c, &e) != FP_EQ || fp_cmp(&d, &f) != FP_EQ) {
printf("div %lu failure!\n", div_n);
draw(&a);draw(&b);draw(&c);draw(&d); draw(&e); draw(&f);
return 0;
}
} else if (!strcmp(cmd, "sqr")) {
fgets(buf, 4095, stdin); fp_read_radix(&a, buf, 64);
fgets(buf, 4095, stdin); fp_read_radix(&b, buf, 64);
// continue;
fp_copy(&a, &c);
fp_sqr(&c, &c); ++sqr_n;
if (fp_cmp(&b, &c) != FP_EQ) {
printf("sqr %lu failure!\n", sqr_n);
draw(&a);draw(&b);draw(&c);
return 0;
}
} else if (!strcmp(cmd, "gcd")) {
fgets(buf, 4095, stdin); fp_read_radix(&a, buf, 64);
fgets(buf, 4095, stdin); fp_read_radix(&b, buf, 64);
fgets(buf, 4095, stdin); fp_read_radix(&c, buf, 64);
// continue;
fp_copy(&a, &d);
fp_gcd(&d, &b, &d); ++gcd_n;
d.sign = c.sign;
if (fp_cmp(&c, &d) != FP_EQ) {
printf("gcd %lu failure!\n", gcd_n);
draw(&a);draw(&b);draw(&c);draw(&d);
return 0;
}
} else if (!strcmp(cmd, "lcm")) {
fgets(buf, 4095, stdin); fp_read_radix(&a, buf, 64);
fgets(buf, 4095, stdin); fp_read_radix(&b, buf, 64);
fgets(buf, 4095, stdin); fp_read_radix(&c, buf, 64);
//continue;
fp_copy(&a, &d);
fp_lcm(&d, &b, &d); ++lcm_n;
d.sign = c.sign;
if (fp_cmp(&c, &d) != FP_EQ) {
printf("lcm %lu failure!\n", lcm_n);
draw(&a);draw(&b);draw(&c);draw(&d);
return 0;
}
} else if (!strcmp(cmd, "expt")) {
fgets(buf, 4095, stdin); fp_read_radix(&a, buf, 64);
fgets(buf, 4095, stdin); fp_read_radix(&b, buf, 64);
fgets(buf, 4095, stdin); fp_read_radix(&c, buf, 64);
fgets(buf, 4095, stdin); fp_read_radix(&d, buf, 64);
// continue;
fp_copy(&a, &e);
fp_exptmod(&e, &b, &c, &e); ++expt_n;
if (fp_cmp(&d, &e) != FP_EQ) {
printf("expt %lu failure!\n", expt_n);
draw(&a);draw(&b);draw(&c);draw(&d); draw(&e);
return 0;
}
} else if (!strcmp(cmd, "invmod")) {
fgets(buf, 4095, stdin); fp_read_radix(&a, buf, 64);
fgets(buf, 4095, stdin); fp_read_radix(&b, buf, 64);
fgets(buf, 4095, stdin); fp_read_radix(&c, buf, 64);
//continue;
fp_invmod(&a, &b, &d);
#if 1
fp_mulmod(&d,&a,&b,&e); ++inv_n;
if (fp_cmp_d(&e, 1) != FP_EQ) {
#else
if (fp_cmp(&d, &c) != FP_EQ) {
#endif
printf("inv [wrong value from MPI?!] failure\n");
draw(&a);draw(&b);draw(&c);draw(&d);
return 0;
}
} else if (!strcmp(cmd, "div2")) { ++div2_n;
fgets(buf, 4095, stdin); fp_read_radix(&a, buf, 64);
fgets(buf, 4095, stdin); fp_read_radix(&b, buf, 64);
fp_div_2(&a, &c);
if (fp_cmp(&c, &b) != FP_EQ) {
printf("div_2 %lu failure\n", div2_n);
draw(&a);
draw(&b);
draw(&c);
return 0;
}
} else if (!strcmp(cmd, "mul2")) { ++mul2_n;
fgets(buf, 4095, stdin); fp_read_radix(&a, buf, 64);
fgets(buf, 4095, stdin); fp_read_radix(&b, buf, 64);
fp_mul_2(&a, &c);
if (fp_cmp(&c, &b) != FP_EQ) {
printf("mul_2 %lu failure\n", mul2_n);
draw(&a);
draw(&b);
draw(&c);
return 0;
}
} else if (!strcmp(cmd, "add_d")) { ++add_d_n;
fgets(buf, 4095, stdin); fp_read_radix(&a, buf, 64);
fgets(buf, 4095, stdin); sscanf(buf, "%d", &ix);
fgets(buf, 4095, stdin); fp_read_radix(&b, buf, 64);
fp_add_d(&a, ix, &c);
if (fp_cmp(&b, &c) != FP_EQ) {
printf("add_d %lu failure\n", add_d_n);
draw(&a);
draw(&b);
draw(&c);
printf("d == %d\n", ix);
return 0;
}
} else if (!strcmp(cmd, "sub_d")) { ++sub_d_n;
fgets(buf, 4095, stdin); fp_read_radix(&a, buf, 64);
fgets(buf, 4095, stdin); sscanf(buf, "%d", &ix);
fgets(buf, 4095, stdin); fp_read_radix(&b, buf, 64);
fp_sub_d(&a, ix, &c);
if (fp_cmp(&b, &c) != FP_EQ) {
printf("sub_d %lu failure\n", sub_d_n);
draw(&a);
draw(&b);
draw(&c);
printf("d == %d\n", ix);
return 0;
}
} else if (!strcmp(cmd, "mul_d")) { ++mul_d_n;
fgets(buf, 4095, stdin); fp_read_radix(&a, buf, 64);
fgets(buf, 4095, stdin); sscanf(buf, "%d", &ix);
fgets(buf, 4095, stdin); fp_read_radix(&b, buf, 64);
fp_mul_d(&a, ix, &c);
if (fp_cmp(&b, &c) != FP_EQ) {
printf("mul_d %lu failure\n", sub_d_n);
draw(&a);
draw(&b);
draw(&c);
printf("d == %d\n", ix);
return 0;
}
}
}
}

BIN
doc/tfm.pdf Normal file

Binary file not shown.

35
fp_2expt.c Normal file
View File

@ -0,0 +1,35 @@
/* TomsFastMath, a fast ISO C bignum library.
*
* This project is meant to fill in where LibTomMath
* falls short. That is speed ;-)
*
* This project is public domain and free for all purposes.
*
* Tom St Denis, tomstdenis@iahu.ca
*/
#include <tfm.h>
/* computes a = 2**b */
void fp_2expt(fp_int *a, int b)
{
int z;
/* zero a as per default */
fp_zero (a);
if (b < 0) {
return;
}
z = b / DIGIT_BIT;
if (z >= FP_SIZE) {
return;
}
/* set the used count of where the bit will go */
a->used = z + 1;
/* put the single bit in its place */
a->dp[z] = ((fp_digit)1) << (b % DIGIT_BIT);
}

39
fp_add.c Normal file
View File

@ -0,0 +1,39 @@
/* TomsFastMath, a fast ISO C bignum library.
*
* This project is meant to fill in where LibTomMath
* falls short. That is speed ;-)
*
* This project is public domain and free for all purposes.
*
* Tom St Denis, tomstdenis@iahu.ca
*/
#include <tfm.h>
void fp_add(fp_int *a, fp_int *b, fp_int *c)
{
int sa, sb;
/* get sign of both inputs */
sa = a->sign;
sb = b->sign;
/* handle two cases, not four */
if (sa == sb) {
/* both positive or both negative */
/* add their magnitudes, copy the sign */
c->sign = sa;
s_fp_add (a, b, c);
} else {
/* one positive, the other negative */
/* subtract the one with the greater magnitude from */
/* the one of the lesser magnitude. The result gets */
/* the sign of the one with the greater magnitude. */
if (fp_cmp_mag (a, b) == FP_LT) {
c->sign = sb;
s_fp_sub (b, a, c);
} else {
c->sign = sa;
s_fp_sub (a, b, c);
}
}
}

18
fp_add_d.c Normal file
View File

@ -0,0 +1,18 @@
/* TomsFastMath, a fast ISO C bignum library.
*
* This project is meant to fill in where LibTomMath
* falls short. That is speed ;-)
*
* This project is public domain and free for all purposes.
*
* Tom St Denis, tomstdenis@iahu.ca
*/
#include <tfm.h>
/* c = a + b */
void fp_add_d(fp_int *a, fp_digit b, fp_int *c)
{
fp_int tmp;
fp_set(&tmp, b);
fp_add(a,&tmp,c);
}

19
fp_addmod.c Normal file
View File

@ -0,0 +1,19 @@
/* TomsFastMath, a fast ISO C bignum library.
*
* This project is meant to fill in where LibTomMath
* falls short. That is speed ;-)
*
* This project is public domain and free for all purposes.
*
* Tom St Denis, tomstdenis@iahu.ca
*/
#include <tfm.h>
/* d = a + b (mod c) */
int fp_addmod(fp_int *a, fp_int *b, fp_int *c, fp_int *d)
{
fp_int tmp;
fp_zero(&tmp);
fp_add(a, b, &tmp);
return fp_mod(&tmp, c, d);
}

27
fp_cmp.c Normal file
View File

@ -0,0 +1,27 @@
/* TomsFastMath, a fast ISO C bignum library.
*
* This project is meant to fill in where LibTomMath
* falls short. That is speed ;-)
*
* This project is public domain and free for all purposes.
*
* Tom St Denis, tomstdenis@iahu.ca
*/
#include <tfm.h>
int fp_cmp(fp_int *a, fp_int *b)
{
if (a->sign == FP_NEG && b->sign == FP_ZPOS) {
return FP_LT;
} else if (a->sign == FP_ZPOS && b->sign == FP_NEG) {
return FP_GT;
} else {
/* compare digits */
if (a->sign == FP_NEG) {
/* if negative compare opposite direction */
return fp_cmp_mag(b, a);
} else {
return fp_cmp_mag(a, b);
}
}
}

34
fp_cmp_d.c Normal file
View File

@ -0,0 +1,34 @@
/* TomsFastMath, a fast ISO C bignum library.
*
* This project is meant to fill in where LibTomMath
* falls short. That is speed ;-)
*
* This project is public domain and free for all purposes.
*
* Tom St Denis, tomstdenis@iahu.ca
*/
#include <tfm.h>
/* compare against a single digit */
int fp_cmp_d(fp_int *a, fp_digit b)
{
/* compare based on sign */
if ((b && a->used == 0) || a->sign == FP_NEG) {
return FP_LT;
}
/* compare based on magnitude */
if (a->used > 1) {
return FP_GT;
}
/* compare the only digit of a to b */
if (a->dp[0] > b) {
return FP_GT;
} else if (a->dp[0] < b) {
return FP_LT;
} else {
return FP_EQ;
}
}

31
fp_cmp_mag.c Normal file
View File

@ -0,0 +1,31 @@
/* TomsFastMath, a fast ISO C bignum library.
*
* This project is meant to fill in where LibTomMath
* falls short. That is speed ;-)
*
* This project is public domain and free for all purposes.
*
* Tom St Denis, tomstdenis@iahu.ca
*/
#include <tfm.h>
int fp_cmp_mag(fp_int *a, fp_int *b)
{
int x;
if (a->used > b->used) {
return FP_GT;
} else if (a->used < b->used) {
return FP_LT;
} else {
for (x = a->used - 1; x >= 0; x--) {
if (a->dp[x] > b->dp[x]) {
return FP_GT;
} else if (a->dp[x] < b->dp[x]) {
return FP_LT;
}
}
}
return FP_EQ;
}

42
fp_cnt_lsb.c Normal file
View File

@ -0,0 +1,42 @@
/* TomsFastMath, a fast ISO C bignum library.
*
* This project is meant to fill in where LibTomMath
* falls short. That is speed ;-)
*
* This project is public domain and free for all purposes.
*
* Tom St Denis, tomstdenis@iahu.ca
*/
#include <tfm.h>
static const int lnz[16] = {
4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0
};
/* Counts the number of lsbs which are zero before the first zero bit */
int fp_cnt_lsb(fp_int *a)
{
int x;
fp_digit q, qq;
/* easy out */
if (fp_iszero(a) == 1) {
return 0;
}
/* scan lower digits until non-zero */
for (x = 0; x < a->used && a->dp[x] == 0; x++);
q = a->dp[x];
x *= DIGIT_BIT;
/* now scan this digit until a 1 is found */
if ((q & 1) == 0) {
do {
qq = q & 15;
x += lnz[qq];
q >>= 4;
} while (qq == 0);
}
return x;
}

32
fp_count_bits.c Normal file
View File

@ -0,0 +1,32 @@
/* TomsFastMath, a fast ISO C bignum library.
*
* This project is meant to fill in where LibTomMath
* falls short. That is speed ;-)
*
* This project is public domain and free for all purposes.
*
* Tom St Denis, tomstdenis@iahu.ca
*/
#include <tfm.h>
int fp_count_bits (fp_int * a)
{
int r;
fp_digit q;
/* shortcut */
if (a->used == 0) {
return 0;
}
/* get number of digits and add that */
r = (a->used - 1) * DIGIT_BIT;
/* take the last digit and count the bits in it */
q = a->dp[a->used - 1];
while (q > ((fp_digit) 0)) {
++r;
q >>= ((fp_digit) 1);
}
return r;
}

153
fp_div.c Normal file
View File

@ -0,0 +1,153 @@
/* TomsFastMath, a fast ISO C bignum library.
*
* This project is meant to fill in where LibTomMath
* falls short. That is speed ;-)
*
* This project is public domain and free for all purposes.
*
* Tom St Denis, tomstdenis@iahu.ca
*/
#include <tfm.h>
/* a/b => cb + d == a */
int fp_div(fp_int *a, fp_int *b, fp_int *c, fp_int *d)
{
fp_int q, x, y, t1, t2;
int n, t, i, norm, neg;
/* is divisor zero ? */
if (fp_iszero (b) == 1) {
return FP_VAL;
}
/* if a < b then q=0, r = a */
if (fp_cmp_mag (a, b) == FP_LT) {
if (d != NULL) {
fp_copy (a, d);
}
if (c != NULL) {
fp_zero (c);
}
return FP_OKAY;
}
fp_init(&q);
q.used = a->used + 2;
fp_init(&t1);
fp_init(&t2);
fp_init_copy(&x, a);
fp_init_copy(&y, b);
/* fix the sign */
neg = (a->sign == b->sign) ? FP_ZPOS : FP_NEG;
x.sign = y.sign = FP_ZPOS;
/* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
norm = fp_count_bits(&y) % DIGIT_BIT;
if (norm < (int)(DIGIT_BIT-1)) {
norm = (DIGIT_BIT-1) - norm;
fp_mul_2d (&x, norm, &x);
fp_mul_2d (&y, norm, &y);
} else {
norm = 0;
}
/* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
n = x.used - 1;
t = y.used - 1;
/* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
fp_lshd (&y, n - t); /* y = y*b**{n-t} */
while (fp_cmp (&x, &y) != FP_LT) {
++(q.dp[n - t]);
fp_sub (&x, &y, &x);
}
/* reset y by shifting it back down */
fp_rshd (&y, n - t);
/* step 3. for i from n down to (t + 1) */
for (i = n; i >= (t + 1); i--) {
if (i > x.used) {
continue;
}
/* step 3.1 if xi == yt then set q{i-t-1} to b-1,
* otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
if (x.dp[i] == y.dp[t]) {
q.dp[i - t - 1] = ((((fp_word)1) << DIGIT_BIT) - 1);
} else {
fp_word tmp;
tmp = ((fp_word) x.dp[i]) << ((fp_word) DIGIT_BIT);
tmp |= ((fp_word) x.dp[i - 1]);
tmp /= ((fp_word) y.dp[t]);
q.dp[i - t - 1] = (fp_digit) (tmp);
}
/* while (q{i-t-1} * (yt * b + y{t-1})) >
xi * b**2 + xi-1 * b + xi-2
do q{i-t-1} -= 1;
*/
q.dp[i - t - 1] = (q.dp[i - t - 1] + 1);
do {
q.dp[i - t - 1] = (q.dp[i - t - 1] - 1);
/* find left hand */
fp_zero (&t1);
t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
t1.dp[1] = y.dp[t];
t1.used = 2;
fp_mul_d (&t1, q.dp[i - t - 1], &t1);
/* find right hand */
t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
t2.dp[2] = x.dp[i];
t2.used = 3;
} while (fp_cmp_mag(&t1, &t2) == FP_GT);
/* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
fp_mul_d (&y, q.dp[i - t - 1], &t1);
fp_lshd (&t1, i - t - 1);
fp_sub (&x, &t1, &x);
/* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
if (x.sign == FP_NEG) {
fp_copy (&y, &t1);
fp_lshd (&t1, i - t - 1);
fp_add (&x, &t1, &x);
q.dp[i - t - 1] = q.dp[i - t - 1] - 1;
}
}
/* now q is the quotient and x is the remainder
* [which we have to normalize]
*/
/* get sign before writing to c */
x.sign = x.used == 0 ? FP_ZPOS : a->sign;
if (c != NULL) {
fp_clamp (&q);
fp_copy (&q, c);
c->sign = neg;
}
if (d != NULL) {
fp_div_2d (&x, norm, &x, NULL);
/* the following is a kludge, essentially we were seeing the right remainder but
with excess digits that should have been zero
*/
for (i = b->used; i < x.used; i++) {
x.dp[i] = 0;
}
fp_clamp(&x);
fp_copy (&x, d);
}
return FP_OKAY;
}

49
fp_div_2.c Normal file
View File

@ -0,0 +1,49 @@
/* TomsFastMath, a fast ISO C bignum library.
*
* This project is meant to fill in where LibTomMath
* falls short. That is speed ;-)
*
* This project is public domain and free for all purposes.
*
* Tom St Denis, tomstdenis@iahu.ca
*/
#include <tfm.h>
/* b = a/2 */
void fp_div_2(fp_int * a, fp_int * b)
{
int x, oldused;
oldused = b->used;
b->used = a->used;
{
register fp_digit r, rr, *tmpa, *tmpb;
/* source alias */
tmpa = a->dp + b->used - 1;
/* dest alias */
tmpb = b->dp + b->used - 1;
/* carry */
r = 0;
for (x = b->used - 1; x >= 0; x--) {
/* get the carry for the next iteration */
rr = *tmpa & 1;
/* shift the current digit, add in carry and store */
*tmpb-- = (*tmpa-- >> 1) | (r << (DIGIT_BIT - 1));
/* forward carry to next iteration */
r = rr;
}
/* zero excess digits */
tmpb = b->dp + b->used;
for (x = b->used; x < oldused; x++) {
*tmpb++ = 0;
}
}
b->sign = a->sign;
fp_clamp (b);
}

75
fp_div_2d.c Normal file
View File

@ -0,0 +1,75 @@
/* TomsFastMath, a fast ISO C bignum library.
*
* This project is meant to fill in where LibTomMath
* falls short. That is speed ;-)
*
* This project is public domain and free for all purposes.
*
* Tom St Denis, tomstdenis@iahu.ca
*/
#include <tfm.h>
/* c = a / 2**b */
void fp_div_2d(fp_int *a, int b, fp_int *c, fp_int *d)
{
fp_digit D, r, rr;
int x;
fp_int t;
/* if the shift count is <= 0 then we do no work */
if (b <= 0) {
fp_copy (a, c);
if (d != NULL) {
fp_zero (d);
}
return;
}
fp_init(&t);
/* get the remainder */
if (d != NULL) {
fp_mod_2d (a, b, &t);
}
/* copy */
fp_copy(a, c);
/* shift by as many digits in the bit count */
if (b >= (int)DIGIT_BIT) {
fp_rshd (c, b / DIGIT_BIT);
}
/* shift any bit count < DIGIT_BIT */
D = (fp_digit) (b % DIGIT_BIT);
if (D != 0) {
register fp_digit *tmpc, mask, shift;
/* mask */
mask = (((fp_digit)1) << D) - 1;
/* shift for lsb */