From e9b1837c8c7caf707f6a038456597568b39c3eeb Mon Sep 17 00:00:00 2001 From: Steffen Jaeckel Date: Thu, 13 Feb 2014 20:21:18 +0100 Subject: [PATCH 1/3] mp_expt_d: bring back pre 921be35779f7d71080ad85c27ed58671602d59b3 state The implementation of the expt_d functionality is now implemented in the mp_expt_d_ex() function. The user can now choose between the old (more timing resistant) version and the new version by modification of the parameter 'fast'. mp_expt_d() defaults to the old version --- bn_mp_expt_d.c | 35 ++------------------ bn_mp_expt_d_ex.c | 81 +++++++++++++++++++++++++++++++++++++++++++++++ makefile | 2 +- tommath.h | 1 + tommath_class.h | 5 +++ 5 files changed, 91 insertions(+), 33 deletions(-) create mode 100644 bn_mp_expt_d_ex.c diff --git a/bn_mp_expt_d.c b/bn_mp_expt_d.c index 2b0b095..ed5f28f 100644 --- a/bn_mp_expt_d.c +++ b/bn_mp_expt_d.c @@ -15,41 +15,12 @@ * Tom St Denis, tomstdenis@gmail.com, http://libtom.org */ -/* calculate c = a**b using a square-multiply algorithm */ +/* wrapper function for mp_expt_d_ex() */ int mp_expt_d (mp_int * a, mp_digit b, mp_int * c) { - int res; - mp_int g; - - if ((res = mp_init_copy (&g, a)) != MP_OKAY) { - return res; - } - - /* set initial result */ - mp_set (c, 1); - - while (b > 0) { - /* if the bit is set multiply */ - if (b & 1) { - if ((res = mp_mul (c, &g, c)) != MP_OKAY) { - mp_clear (&g); - return res; - } - } - - /* square */ - if (b > 1 && (res = mp_sqr (&g, &g)) != MP_OKAY) { - mp_clear (&g); - return res; - } - - /* shift to next bit */ - b >>= 1; - } - - mp_clear (&g); - return MP_OKAY; + return mp_expt_d_ex(a, b, c, 0); } + #endif /* $Source$ */ diff --git a/bn_mp_expt_d_ex.c b/bn_mp_expt_d_ex.c new file mode 100644 index 0000000..f8ff172 --- /dev/null +++ b/bn_mp_expt_d_ex.c @@ -0,0 +1,81 @@ +#include +#ifdef BN_MP_EXPT_D_EX_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis + * + * LibTomMath is a library that provides multiple-precision + * integer arithmetic as well as number theoretic functionality. + * + * The library was designed directly after the MPI library by + * Michael Fromberger but has been written from scratch with + * additional optimizations in place. + * + * The library is free for all purposes without any express + * guarantee it works. + * + * Tom St Denis, tomstdenis@gmail.com, http://libtom.org + */ + +/* calculate c = a**b using a square-multiply algorithm */ +int mp_expt_d_ex (mp_int * a, mp_digit b, mp_int * c, int fast) +{ + int res; + unsigned int x; + + mp_int g; + + if ((res = mp_init_copy (&g, a)) != MP_OKAY) { + return res; + } + + /* set initial result */ + mp_set (c, 1); + + if (fast) { + while (b > 0) { + /* if the bit is set multiply */ + if (b & 1) { + if ((res = mp_mul (c, &g, c)) != MP_OKAY) { + mp_clear (&g); + return res; + } + } + + /* square */ + if (b > 1 && (res = mp_sqr (&g, &g)) != MP_OKAY) { + mp_clear (&g); + return res; + } + + /* shift to next bit */ + b >>= 1; + } + } + else { + for (x = 0; x < DIGIT_BIT; x++) { + /* square */ + if ((res = mp_sqr (c, c)) != MP_OKAY) { + mp_clear (&g); + return res; + } + + /* if the bit is set multiply */ + if ((b & (mp_digit) (((mp_digit)1) << (DIGIT_BIT - 1))) != 0) { + if ((res = mp_mul (c, &g, c)) != MP_OKAY) { + mp_clear (&g); + return res; + } + } + + /* shift to next bit */ + b <<= 1; + } + } /* if ... else */ + + mp_clear (&g); + return MP_OKAY; +} +#endif + +/* $Source$ */ +/* $Revision$ */ +/* $Date$ */ diff --git a/makefile b/makefile index ed7bd91..e9b2eab 100644 --- a/makefile +++ b/makefile @@ -96,7 +96,7 @@ bn_mp_init_multi.o bn_mp_clear_multi.o bn_mp_exteuclid.o bn_mp_toradix_n.o \ bn_mp_prime_random_ex.o bn_mp_get_int.o bn_mp_sqrt.o bn_mp_is_square.o bn_mp_init_set.o \ bn_mp_init_set_int.o bn_mp_invmod_slow.o bn_mp_prime_rabin_miller_trials.o \ bn_mp_to_signed_bin_n.o bn_mp_to_unsigned_bin_n.o bn_mp_import.o bn_mp_export.o \ -bn_mp_balance_mul.o +bn_mp_balance_mul.o bn_mp_expt_d_ex.o $(LIBNAME): $(OBJECTS) $(AR) $(ARFLAGS) $@ $(OBJECTS) diff --git a/tommath.h b/tommath.h index 8a71d3c..7b5a703 100644 --- a/tommath.h +++ b/tommath.h @@ -364,6 +364,7 @@ int mp_div_3(mp_int *a, mp_int *c, mp_digit *d); /* c = a**b */ int mp_expt_d(mp_int *a, mp_digit b, mp_int *c); +int mp_expt_d_ex (mp_int * a, mp_digit b, mp_int * c, int fast); /* c = a mod b, 0 <= c < b */ int mp_mod_d(mp_int *a, mp_digit b, mp_digit *c); diff --git a/tommath_class.h b/tommath_class.h index 238e72d..6569f18 100644 --- a/tommath_class.h +++ b/tommath_class.h @@ -41,6 +41,7 @@ #define BN_MP_EXCH_C #define BN_MP_EXPORT_C #define BN_MP_EXPT_D_C +#define BN_MP_EXPT_D_EX_C #define BN_MP_EXPTMOD_C #define BN_MP_EXPTMOD_FAST_C #define BN_MP_EXTEUCLID_C @@ -333,6 +334,10 @@ #endif #if defined(BN_MP_EXPT_D_C) + #define BN_MP_EXPT_D_EX_C +#endif + +#if defined(BN_MP_EXPT_D_EX_C) #define BN_MP_INIT_COPY_C #define BN_MP_SET_C #define BN_MP_MUL_C From 52cfd5ff0aaf28d3b53bbc4c1308e5e4884214bb Mon Sep 17 00:00:00 2001 From: Steffen Jaeckel Date: Fri, 14 Feb 2014 11:26:07 +0100 Subject: [PATCH 2/3] mp_n_root: add mp_n_root_ex() with parameter 'fast' for mp_expt_d_ex() This change is introduced to be able to choose the underlying implementation of mp_expt_d_ex() The implementation of the root_n functionality is now implemented in the mp_n_root_ex() function. The parameter 'fast' is just passed over to mp_expt_d_ex(). mp_n_root() defaults to the pre 921be35779f7d71080ad85c27ed58671602d59b3 implementation --- bn_mp_n_root.c | 110 ++------------------------------------ bn_mp_n_root_ex.c | 132 ++++++++++++++++++++++++++++++++++++++++++++++ makefile | 2 +- tommath.h | 1 + tommath_class.h | 11 ++-- 5 files changed, 144 insertions(+), 112 deletions(-) create mode 100644 bn_mp_n_root_ex.c diff --git a/bn_mp_n_root.c b/bn_mp_n_root.c index 11b878b..56f1586 100644 --- a/bn_mp_n_root.c +++ b/bn_mp_n_root.c @@ -15,116 +15,14 @@ * Tom St Denis, tomstdenis@gmail.com, http://libtom.org */ -/* find the n'th root of an integer - * - * Result found such that (c)**b <= a and (c+1)**b > a - * - * This algorithm uses Newton's approximation - * x[i+1] = x[i] - f(x[i])/f'(x[i]) - * which will find the root in log(N) time where - * each step involves a fair bit. This is not meant to - * find huge roots [square and cube, etc]. +/* wrapper function for mp_n_root_ex() + * computes c = (a)**(1/b) such that (c)**b <= a and (c+1)**b > a */ int mp_n_root (mp_int * a, mp_digit b, mp_int * c) { - mp_int t1, t2, t3; - int res, neg; - - /* input must be positive if b is even */ - if ((b & 1) == 0 && a->sign == MP_NEG) { - return MP_VAL; - } - - if ((res = mp_init (&t1)) != MP_OKAY) { - return res; - } - - if ((res = mp_init (&t2)) != MP_OKAY) { - goto LBL_T1; - } - - if ((res = mp_init (&t3)) != MP_OKAY) { - goto LBL_T2; - } - - /* if a is negative fudge the sign but keep track */ - neg = a->sign; - a->sign = MP_ZPOS; - - /* t2 = 2 */ - mp_set (&t2, 2); - - do { - /* t1 = t2 */ - if ((res = mp_copy (&t2, &t1)) != MP_OKAY) { - goto LBL_T3; - } - - /* t2 = t1 - ((t1**b - a) / (b * t1**(b-1))) */ - - /* t3 = t1**(b-1) */ - if ((res = mp_expt_d (&t1, b - 1, &t3)) != MP_OKAY) { - goto LBL_T3; - } - - /* numerator */ - /* t2 = t1**b */ - if ((res = mp_mul (&t3, &t1, &t2)) != MP_OKAY) { - goto LBL_T3; - } - - /* t2 = t1**b - a */ - if ((res = mp_sub (&t2, a, &t2)) != MP_OKAY) { - goto LBL_T3; - } - - /* denominator */ - /* t3 = t1**(b-1) * b */ - if ((res = mp_mul_d (&t3, b, &t3)) != MP_OKAY) { - goto LBL_T3; - } - - /* t3 = (t1**b - a)/(b * t1**(b-1)) */ - if ((res = mp_div (&t2, &t3, &t3, NULL)) != MP_OKAY) { - goto LBL_T3; - } - - if ((res = mp_sub (&t1, &t3, &t2)) != MP_OKAY) { - goto LBL_T3; - } - } while (mp_cmp (&t1, &t2) != MP_EQ); - - /* result can be off by a few so check */ - for (;;) { - if ((res = mp_expt_d (&t1, b, &t2)) != MP_OKAY) { - goto LBL_T3; - } - - if (mp_cmp (&t2, a) == MP_GT) { - if ((res = mp_sub_d (&t1, 1, &t1)) != MP_OKAY) { - goto LBL_T3; - } - } else { - break; - } - } - - /* reset the sign of a first */ - a->sign = neg; - - /* set the result */ - mp_exch (&t1, c); - - /* set the sign of the result */ - c->sign = neg; - - res = MP_OKAY; - -LBL_T3:mp_clear (&t3); -LBL_T2:mp_clear (&t2); -LBL_T1:mp_clear (&t1); - return res; + return mp_n_root_ex(a, b, c, 0); } + #endif /* $Source$ */ diff --git a/bn_mp_n_root_ex.c b/bn_mp_n_root_ex.c new file mode 100644 index 0000000..6b42643 --- /dev/null +++ b/bn_mp_n_root_ex.c @@ -0,0 +1,132 @@ +#include +#ifdef BN_MP_N_ROOT_EX_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis + * + * LibTomMath is a library that provides multiple-precision + * integer arithmetic as well as number theoretic functionality. + * + * The library was designed directly after the MPI library by + * Michael Fromberger but has been written from scratch with + * additional optimizations in place. + * + * The library is free for all purposes without any express + * guarantee it works. + * + * Tom St Denis, tomstdenis@gmail.com, http://libtom.org + */ + +/* find the n'th root of an integer + * + * Result found such that (c)**b <= a and (c+1)**b > a + * + * This algorithm uses Newton's approximation + * x[i+1] = x[i] - f(x[i])/f'(x[i]) + * which will find the root in log(N) time where + * each step involves a fair bit. This is not meant to + * find huge roots [square and cube, etc]. + */ +int mp_n_root_ex (mp_int * a, mp_digit b, mp_int * c, int fast) +{ + mp_int t1, t2, t3; + int res, neg; + + /* input must be positive if b is even */ + if ((b & 1) == 0 && a->sign == MP_NEG) { + return MP_VAL; + } + + if ((res = mp_init (&t1)) != MP_OKAY) { + return res; + } + + if ((res = mp_init (&t2)) != MP_OKAY) { + goto LBL_T1; + } + + if ((res = mp_init (&t3)) != MP_OKAY) { + goto LBL_T2; + } + + /* if a is negative fudge the sign but keep track */ + neg = a->sign; + a->sign = MP_ZPOS; + + /* t2 = 2 */ + mp_set (&t2, 2); + + do { + /* t1 = t2 */ + if ((res = mp_copy (&t2, &t1)) != MP_OKAY) { + goto LBL_T3; + } + + /* t2 = t1 - ((t1**b - a) / (b * t1**(b-1))) */ + + /* t3 = t1**(b-1) */ + if ((res = mp_expt_d_ex (&t1, b - 1, &t3, fast)) != MP_OKAY) { + goto LBL_T3; + } + + /* numerator */ + /* t2 = t1**b */ + if ((res = mp_mul (&t3, &t1, &t2)) != MP_OKAY) { + goto LBL_T3; + } + + /* t2 = t1**b - a */ + if ((res = mp_sub (&t2, a, &t2)) != MP_OKAY) { + goto LBL_T3; + } + + /* denominator */ + /* t3 = t1**(b-1) * b */ + if ((res = mp_mul_d (&t3, b, &t3)) != MP_OKAY) { + goto LBL_T3; + } + + /* t3 = (t1**b - a)/(b * t1**(b-1)) */ + if ((res = mp_div (&t2, &t3, &t3, NULL)) != MP_OKAY) { + goto LBL_T3; + } + + if ((res = mp_sub (&t1, &t3, &t2)) != MP_OKAY) { + goto LBL_T3; + } + } while (mp_cmp (&t1, &t2) != MP_EQ); + + /* result can be off by a few so check */ + for (;;) { + if ((res = mp_expt_d_ex (&t1, b, &t2, fast)) != MP_OKAY) { + goto LBL_T3; + } + + if (mp_cmp (&t2, a) == MP_GT) { + if ((res = mp_sub_d (&t1, 1, &t1)) != MP_OKAY) { + goto LBL_T3; + } + } else { + break; + } + } + + /* reset the sign of a first */ + a->sign = neg; + + /* set the result */ + mp_exch (&t1, c); + + /* set the sign of the result */ + c->sign = neg; + + res = MP_OKAY; + +LBL_T3:mp_clear (&t3); +LBL_T2:mp_clear (&t2); +LBL_T1:mp_clear (&t1); + return res; +} +#endif + +/* $Source$ */ +/* $Revision$ */ +/* $Date$ */ diff --git a/makefile b/makefile index e9b2eab..6ef47b7 100644 --- a/makefile +++ b/makefile @@ -96,7 +96,7 @@ bn_mp_init_multi.o bn_mp_clear_multi.o bn_mp_exteuclid.o bn_mp_toradix_n.o \ bn_mp_prime_random_ex.o bn_mp_get_int.o bn_mp_sqrt.o bn_mp_is_square.o bn_mp_init_set.o \ bn_mp_init_set_int.o bn_mp_invmod_slow.o bn_mp_prime_rabin_miller_trials.o \ bn_mp_to_signed_bin_n.o bn_mp_to_unsigned_bin_n.o bn_mp_import.o bn_mp_export.o \ -bn_mp_balance_mul.o bn_mp_expt_d_ex.o +bn_mp_balance_mul.o bn_mp_expt_d_ex.o bn_mp_n_root_ex.o $(LIBNAME): $(OBJECTS) $(AR) $(ARFLAGS) $@ $(OBJECTS) diff --git a/tommath.h b/tommath.h index 7b5a703..13cc64a 100644 --- a/tommath.h +++ b/tommath.h @@ -400,6 +400,7 @@ int mp_lcm(mp_int *a, mp_int *b, mp_int *c); * returns error if a < 0 and b is even */ int mp_n_root(mp_int *a, mp_digit b, mp_int *c); +int mp_n_root_ex (mp_int * a, mp_digit b, mp_int * c, int fast); /* special sqrt algo */ int mp_sqrt(mp_int *arg, mp_int *ret); diff --git a/tommath_class.h b/tommath_class.h index 6569f18..afd08a4 100644 --- a/tommath_class.h +++ b/tommath_class.h @@ -77,6 +77,7 @@ #define BN_MP_MUL_D_C #define BN_MP_MULMOD_C #define BN_MP_N_ROOT_C +#define BN_MP_N_ROOT_EX_C #define BN_MP_NEG_C #define BN_MP_OR_C #define BN_MP_PRIME_FERMAT_C @@ -614,10 +615,14 @@ #endif #if defined(BN_MP_N_ROOT_C) + #define BN_MP_N_ROOT_EX_C +#endif + +#if defined(BN_MP_N_ROOT_EX_C) #define BN_MP_INIT_C #define BN_MP_SET_C #define BN_MP_COPY_C - #define BN_MP_EXPT_D_C + #define BN_MP_EXPT_D_EX_C #define BN_MP_MUL_C #define BN_MP_SUB_C #define BN_MP_MUL_D_C @@ -1023,7 +1028,3 @@ #else #define LTM_LAST #endif - -/* $Source$ */ -/* $Revision$ */ -/* $Date$ */ From 52bb535ff77285646af89faa5d7c9f893b8a909f Mon Sep 17 00:00:00 2001 From: Steffen Jaeckel Date: Fri, 14 Feb 2014 12:53:48 +0100 Subject: [PATCH 3/3] demo: test both mp_n_root() implementations --- demo/demo.c | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/demo/demo.c b/demo/demo.c index 37dba51..78e0903 100644 --- a/demo/demo.c +++ b/demo/demo.c @@ -179,8 +179,13 @@ printf("compare no compare!\n"); return EXIT_FAILURE; } printf("\nmp_sqrt() error!"); return EXIT_FAILURE; } - mp_n_root(&a, 2, &a); - if (mp_cmp_mag(&b, &a) != MP_EQ) { + mp_n_root_ex(&a, 2, &c, 0); + mp_n_root_ex(&a, 2, &d, 1); + if (mp_cmp_mag(&c, &d) != MP_EQ) { + printf("\nmp_n_root_ex() bad result!"); + return EXIT_FAILURE; + } + if (mp_cmp_mag(&b, &c) != MP_EQ) { printf("mp_sqrt() bad result!\n"); return 1; }