diff --git a/demo/test.c b/demo/test.c index 2fa6e08d..20875da5 100644 --- a/demo/test.c +++ b/demo/test.c @@ -2142,8 +2142,71 @@ static int test_s_mp_sqr_toom(void) mp_clear_multi(&a, &b, &c, NULL); return EXIT_FAILURE; } +#if (MP_DIGIT_BIT > 16) +static int test_s_mp_mul_ntt(void) +{ + mp_int a, b, c, d; + int size; + +#if (MP_DIGIT_BIT == 60) + int tc_cutoff; +#endif + + DOR(mp_init_multi(&a, &b, &c, &d, NULL)); + /* This number construction is limb-size specific */ +#if (MP_DIGIT_BIT == 60) + DO(mp_rand(&a, 1196)); + DO(mp_mul_2d(&a,71787 - mp_count_bits(&a), &a)); + + DO(mp_rand(&b, 1338)); + DO(mp_mul_2d(&b, 80318 - mp_count_bits(&b), &b)); + DO(mp_mul_2d(&b, 6310, &b)); + DO(mp_2expt(&c, 99000 - 1000)); + DO(mp_add(&b, &c, &b)); + + tc_cutoff = MP_MUL_NTT_CUTOFF; + MP_MUL_NTT_CUTOFF = INT_MAX; + DO(mp_mul(&a, &b, &c)); + MP_MUL_NTT_CUTOFF = tc_cutoff; + DO(mp_mul(&a, &b, &d)); + EXPECT(mp_cmp(&c, &d) == MP_EQ); +#endif + + for (size = MP_MUL_NTT_CUTOFF; size < (MP_MUL_NTT_CUTOFF + 20); size++) { + DO(mp_rand(&a, size)); + DO(mp_rand(&b, size)); + DO(s_mp_mul_ntt(&a, &b, &c)); + DO(s_mp_mul_full(&a,&b,&d)); + EXPECT(mp_cmp(&c, &d) == MP_EQ); + } + + mp_clear_multi(&a, &b, &c, &d, NULL); + return EXIT_SUCCESS; +LBL_ERR: + mp_clear_multi(&a, &b, &c, &d, NULL); + return EXIT_FAILURE; +} +static int test_s_mp_sqr_ntt(void) +{ + mp_int a, b, c; + int size; + DOR(mp_init_multi(&a, &b, &c, NULL)); + for (size = MP_SQR_NTT_CUTOFF; size < (MP_SQR_NTT_CUTOFF + 20); size++) { + DO(mp_rand(&a, size)); + DO(s_mp_sqr_ntt(&a, &b)); + DO(s_mp_sqr(&a, &c)); + EXPECT(mp_cmp(&b, &c) == MP_EQ); + } + + mp_clear_multi(&a, &b, &c, NULL); + return EXIT_SUCCESS; +LBL_ERR: + mp_clear_multi(&a, &b, &c, NULL); + return EXIT_FAILURE; +} +#endif static int test_mp_radix_size(void) { mp_int a; @@ -2606,7 +2669,11 @@ static int unit_tests(int argc, char **argv) T2(s_mp_mul_karatsuba, ONLY_PUBLIC_API, S_MP_MUL_KARATSUBA), T2(s_mp_sqr_karatsuba, ONLY_PUBLIC_API, S_MP_SQR_KARATSUBA), T2(s_mp_mul_toom, ONLY_PUBLIC_API, S_MP_MUL_TOOM), - T2(s_mp_sqr_toom, ONLY_PUBLIC_API, S_MP_SQR_TOOM) + T2(s_mp_sqr_toom, ONLY_PUBLIC_API, S_MP_SQR_TOOM), +#if (MP_DIGIT_BIT > 16) + T2(s_mp_mul_ntt, ONLY_PUBLIC_API, S_MP_MUL_NTT), + T2(s_mp_sqr_ntt, ONLY_PUBLIC_API, S_MP_SQR_NTT) +#endif #undef T3 #undef T2 #undef T1 diff --git a/etc/tune.c b/etc/tune.c index 8ad20289..6b7afc57 100644 --- a/etc/tune.c +++ b/etc/tune.c @@ -20,8 +20,8 @@ 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 uint64_t s_time_mul(int size, int increment); +static uint64_t s_time_sqr(int size, int increment); static void s_usage(char *s); static uint64_t s_timer_function(void) @@ -61,7 +61,7 @@ static int s_stabilization_extra; static int s_offset = 1; #define s_mp_mul_full(a, b, c) s_mp_mul(a, b, c, (a)->used + (b)->used + 1) -static uint64_t s_time_mul(int size) +static uint64_t s_time_mul(int size, int increment) { int x; mp_err e; @@ -83,7 +83,7 @@ static uint64_t s_time_mul(int size) } s_timer_start(); - for (x = 0; x < s_number_of_test_loops; x++) { + for (x = 0; x < s_number_of_test_loops; x += increment) { if ((e = mp_mul(&a,&b,&c)) != MP_OKAY) { t1 = UINT64_MAX; goto LBL_ERR; @@ -107,7 +107,7 @@ static uint64_t s_time_mul(int size) return t1; } -static uint64_t s_time_sqr(int size) +static uint64_t s_time_sqr(int size, int increment) { int x; mp_err e; @@ -125,7 +125,7 @@ static uint64_t s_time_sqr(int size) } s_timer_start(); - for (x = 0; x < s_number_of_test_loops; x++) { + for (x = 0; x < s_number_of_test_loops; x+=increment) { if ((e = mp_sqr(&a,&b)) != MP_OKAY) { t1 = UINT64_MAX; goto LBL_ERR; @@ -158,7 +158,7 @@ struct tune_args { int increment_print; } args; -static void s_run(const char *name, uint64_t (*op)(int size), int *cutoff) +static void s_run(const char *name, uint64_t (*op)(int size, int increment), int *cutoff, int increment) { int x, count = 0; uint64_t t1, t2; @@ -167,14 +167,14 @@ static void s_run(const char *name, uint64_t (*op)(int size), int *cutoff) } for (x = 8; x < args.upper_limit_print; x += args.increment_print) { *cutoff = INT_MAX; - t1 = op(x); + t1 = op(x, increment); if ((t1 == 0u) || (t1 == UINT64_MAX)) { fprintf(stderr,"%s failed at x = INT_MAX (%s)\n", name, (t1 == 0u)?"wrong result":"internal error"); exit(EXIT_FAILURE); } *cutoff = x; - t2 = op(x); + t2 = op(x, increment); if ((t2 == 0u) || (t2 == UINT64_MAX)) { fprintf(stderr,"%s failed (%s)\n", name, (t2 == 0u)?"wrong result":"internal error"); @@ -243,6 +243,8 @@ static void s_usage(char *s) 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," nttm = NTT multiplication\n"); + fprintf(stderr," ntts = NTT squaring\n"); fprintf(stderr," Implies '-p'\n"); fprintf(stderr," -h this message\n"); exit(s_exit_code); @@ -251,10 +253,11 @@ static void s_usage(char *s) struct cutoffs { int MUL_KARATSUBA, SQR_KARATSUBA; int MUL_TOOM, SQR_TOOM; + int MUL_NTT, SQR_NTT; }; const struct cutoffs max_cutoffs = -{ INT_MAX, INT_MAX, INT_MAX, INT_MAX }; +{ INT_MAX, INT_MAX, INT_MAX, INT_MAX, INT_MAX, INT_MAX}; static void set_cutoffs(const struct cutoffs *c) { @@ -262,6 +265,8 @@ static void set_cutoffs(const struct cutoffs *c) MP_SQR_KARATSUBA_CUTOFF = c->SQR_KARATSUBA; MP_MUL_TOOM_CUTOFF = c->MUL_TOOM; MP_SQR_TOOM_CUTOFF = c->SQR_TOOM; + MP_MUL_NTT_CUTOFF = c->MUL_NTT; + MP_SQR_NTT_CUTOFF = c->SQR_NTT; } static void get_cutoffs(struct cutoffs *c) @@ -270,7 +275,8 @@ static void get_cutoffs(struct cutoffs *c) c->SQR_KARATSUBA = MP_SQR_KARATSUBA_CUTOFF; c->MUL_TOOM = MP_MUL_TOOM_CUTOFF; c->SQR_TOOM = MP_SQR_TOOM_CUTOFF; - + c->MUL_NTT = MP_MUL_NTT_CUTOFF; + c->SQR_NTT = MP_SQR_NTT_CUTOFF; } int main(int argc, char **argv) @@ -302,7 +308,7 @@ int main(int argc, char **argv) args.bncore = 0; args.terse = 0; - args.upper_limit_print = 3000; + args.upper_limit_print = 20000; args.increment_print = 1; /* Very simple option parser, please treat it nicely. */ @@ -423,6 +429,10 @@ int main(int argc, char **argv) MP_MUL_TOOM_CUTOFF = (int)s_strtol(str, &endptr, "[3/4] No value for MP_MUL_TOOM_CUTOFF given"); str = endptr + 1; MP_SQR_TOOM_CUTOFF = (int)s_strtol(str, &endptr, "[4/4] No value for MP_SQR_TOOM_CUTOFF given"); + str = endptr + 1; + MP_MUL_NTT_CUTOFF = (int)s_strtol(str, &endptr, "[3/4] No value for MP_MUL_NTT_CUTOFF given"); + str = endptr + 1; + MP_SQR_NTT_CUTOFF = (int)s_strtol(str, &endptr, "[4/4] No value for MP_SQR_NTT_CUTOFF given"); break; case 'h': s_exit_code = EXIT_SUCCESS; @@ -448,7 +458,7 @@ int main(int argc, char **argv) struct { const char *name; int *cutoff, *update; - uint64_t (*fn)(int size); + uint64_t (*fn)(int size, int increment); } test[] = { #define T_MUL_SQR(n, o, f) { #n, &MP_##o##_CUTOFF, &(updated.o), MP_HAS(S_MP_##o) ? f : NULL } /* @@ -457,33 +467,49 @@ int main(int argc, char **argv) of the macro MP_WPARRAY in tommath.h which needs to be changed manually (to 0 (zero)). */ - T_MUL_SQR("Karatsuba multiplication", MUL_KARATSUBA, s_time_mul), - T_MUL_SQR("Karatsuba squaring", SQR_KARATSUBA, s_time_sqr), - T_MUL_SQR("Toom-Cook 3-way multiplying", MUL_TOOM, s_time_mul), - T_MUL_SQR("Toom-Cook 3-way squaring", SQR_TOOM, s_time_sqr), + T_MUL_SQR("Karatsuba multiplication", MUL_KARATSUBA, s_time_mul), /* 0 */ + T_MUL_SQR("Karatsuba squaring", SQR_KARATSUBA, s_time_sqr), /* 1 */ + T_MUL_SQR("Toom-Cook 3-way multiplying", MUL_TOOM, s_time_mul), /* 2 */ + T_MUL_SQR("Toom-Cook 3-way squaring", SQR_TOOM, s_time_sqr), /* 3 */ + T_MUL_SQR("NTT multiplying", MUL_NTT, s_time_mul), /* 4 */ + T_MUL_SQR("NTT squaring", SQR_NTT, s_time_sqr), /* 5 */ #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) { + for (n = 0; n < 4; ++n) { + if (test[n].fn != NULL) { + s_run(test[n].name, test[n].fn, test[n].cutoff, 1); + *test[n].update = *test[n].cutoff; + /* *test[n].cutoff = INT_MAX; */ + } + } + for (n = 4; n < 6; ++n) { if (test[n].fn != NULL) { - s_run(test[n].name, test[n].fn, test[n].cutoff); + s_run(test[n].name, test[n].fn, test[n].cutoff, 53); *test[n].update = *test[n].cutoff; - *test[n].cutoff = INT_MAX; + /* Only test Karatsuba against school, + test the rest against the cutoff one lower */ + /* *test[n].cutoff = INT_MAX; */ } } } if (args.terse == 1) { - printf("%d %d %d %d\n", + printf("%d %d %d %d %d %d\n", updated.MUL_KARATSUBA, updated.SQR_KARATSUBA, updated.MUL_TOOM, - updated.SQR_TOOM); + updated.SQR_TOOM, + updated.MUL_NTT, + updated.SQR_NTT); + } else { printf("MUL_KARATSUBA_CUTOFF = %d\n", updated.MUL_KARATSUBA); printf("SQR_KARATSUBA_CUTOFF = %d\n", updated.SQR_KARATSUBA); printf("MUL_TOOM_CUTOFF = %d\n", updated.MUL_TOOM); printf("SQR_TOOM_CUTOFF = %d\n", updated.SQR_TOOM); + printf("MUL_NTT_CUTOFF = %d\n", updated.MUL_NTT); + printf("SQR_NTT_CUTOFF = %d\n", updated.SQR_NTT); } if (args.print == 1) { @@ -503,9 +529,9 @@ int main(int argc, char **argv) for (x = 8; x < args.upper_limit_print; x += args.increment_print) { set_cutoffs(&max_cutoffs); - t1 = s_time_mul(x); + t1 = s_time_mul(x,1); set_cutoffs(&orig); - t2 = s_time_mul(x); + t2 = s_time_mul(x,1); 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) { @@ -513,9 +539,9 @@ int main(int argc, char **argv) fflush(stdout); } set_cutoffs(&max_cutoffs); - t1 = s_time_sqr(x); + t1 = s_time_sqr(x,1); set_cutoffs(&orig); - t2 = s_time_sqr(x); + t2 = s_time_sqr(x,1); 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) { @@ -527,16 +553,20 @@ int main(int argc, char **argv) if (args.verbose == 1) { set_cutoffs(&orig); if (args.terse == 1) { - printf("%d %d %d %d\n", + printf("%d %d %d %d %d %d\n", MP_MUL_KARATSUBA_CUTOFF, MP_SQR_KARATSUBA_CUTOFF, MP_MUL_TOOM_CUTOFF, - MP_SQR_TOOM_CUTOFF); + MP_SQR_TOOM_CUTOFF, + MP_MUL_NTT_CUTOFF, + MP_SQR_NTT_CUTOFF); } else { printf("MUL_KARATSUBA_CUTOFF = %d\n", MP_MUL_KARATSUBA_CUTOFF); printf("SQR_KARATSUBA_CUTOFF = %d\n", MP_SQR_KARATSUBA_CUTOFF); printf("MUL_TOOM_CUTOFF = %d\n", MP_MUL_TOOM_CUTOFF); printf("SQR_TOOM_CUTOFF = %d\n", MP_SQR_TOOM_CUTOFF); + printf("MUL_NTT_CUTOFF = %d\n", MP_MUL_NTT_CUTOFF); + printf("SQR_NTT_CUTOFF = %d\n", MP_SQR_NTT_CUTOFF); } } } diff --git a/etc/tune_it.sh b/etc/tune_it.sh index dba5b696..2f03f133 100755 --- a/etc/tune_it.sh +++ b/etc/tune_it.sh @@ -33,10 +33,12 @@ RNUM=0 # test are running if you do not have enough memory to avoid page faults. # ############################################################################# +# TODO: Because of NTT the defaults of LIMIT is down from 100, RLOOPS down from 10 +# and LAG is up from 3 # Number of rounds overall. -LIMIT=100 +LIMIT=10 # Number of loops for each input. -RLOOPS=10 +RLOOPS=3 # 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 @@ -45,7 +47,7 @@ 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 +LAG=100 # 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 @@ -56,7 +58,7 @@ 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" $? +echo "km ks tc3m tc3s nttm ntts" > $FILE_NAME || die "Writing header to $FILE_NAME" $? i=1 while [ $i -le $LIMIT ]; do RNUM=$(LCG) @@ -104,3 +106,19 @@ echo "#define MP_DEFAULT_MUL_TOOM_CUTOFF $TMP" >> $TOMMATH_CUTOFFS_H || die TMP=$(median $FILE_NAME 4 $i) echo "#define MP_DEFAULT_SQR_TOOM_CUTOFF $TMP" echo "#define MP_DEFAULT_SQR_TOOM_CUTOFF $TMP" >> $TOMMATH_CUTOFFS_H || die "(tc3s) Appending to $TOMMATH_CUTOFFS_H" $? +TMP=$(median $FILE_NAME 5 $i) +echo "#define MP_DEFAULT_MUL_NTT_CUTOFF $TMP" +echo "#define MP_DEFAULT_MUL_NTT_CUTOFF $TMP" >> $TOMMATH_CUTOFFS_H || die "(nttm) Appending to $TOMMATH_CUTOFFS_H" $? +TMP=$(median $FILE_NAME 6 $i) +echo "#define MP_DEFAULT_SQR_NTT_CUTOFF $TMP" +echo "#define MP_DEFAULT_SQR_NTT_CUTOFF $TMP" >> $TOMMATH_CUTOFFS_H || die "(ntts) Appending to $TOMMATH_CUTOFFS_H" $? + + + + + + + + + + diff --git a/helper.pl b/helper.pl index 8a430f73..d024b30b 100755 --- a/helper.pl +++ b/helper.pl @@ -480,6 +480,8 @@ sub generate_def { MP_SQR_KARATSUBA_CUTOFF MP_MUL_TOOM_CUTOFF MP_SQR_TOOM_CUTOFF + MP_MUL_NTT_CUTOFF + MP_SQR_NTT_CUTOFF "; return 0; } diff --git a/libtommath_VS2008.vcproj b/libtommath_VS2008.vcproj index 080b7e0d..023ada29 100644 --- a/libtommath_VS2008.vcproj +++ b/libtommath_VS2008.vcproj @@ -888,10 +888,146 @@ RelativePath="s_mp_mul_karatsuba.c" > + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -928,6 +1064,10 @@ RelativePath="s_mp_sqr_karatsuba.c" > + + diff --git a/makefile b/makefile index 6daa522b..1f0e679f 100644 --- a/makefile +++ b/makefile @@ -47,11 +47,18 @@ mp_unpack.o mp_warray_free.o mp_xor.o mp_zero.o s_mp_add.o s_mp_copy_digs.o s_mp s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_fp_log.o \ s_mp_fp_log_d.o s_mp_get_bit.o s_mp_invmod.o s_mp_invmod_odd.o s_mp_log_2expt.o \ s_mp_montgomery_reduce_comba.o s_mp_mul.o s_mp_mul_balance.o s_mp_mul_comba.o s_mp_mul_high.o \ -s_mp_mul_high_comba.o s_mp_mul_karatsuba.o s_mp_mul_toom.o s_mp_prime_is_divisible.o s_mp_prime_tab.o \ -s_mp_radix_map.o s_mp_radix_size_overestimate.o s_mp_rand_platform.o s_mp_rand_source.o s_mp_sqr.o \ -s_mp_sqr_comba.o s_mp_sqr_karatsuba.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o s_mp_warray_get.o \ +s_mp_mul_high_comba.o s_mp_mul_karatsuba.o s_mp_mul_ntt.o s_mp_mul_toom.o s_mp_ntt_autoconvolve_single.o \ +s_mp_ntt_convolve_single.o s_mp_ntt_crt_assemble.o s_mp_ntt_crt_coeff.o s_mp_ntt_crt_init.o \ +s_mp_ntt_find_smooth_prime.o s_mp_ntt_from_mont.o s_mp_ntt_fwd2.o s_mp_ntt_fwd3.o s_mp_ntt_fwd5.o \ +s_mp_ntt_ilog2_ceil.o s_mp_ntt_inv2.o s_mp_ntt_inv3.o s_mp_ntt_inv5.o s_mp_ntt_invmod_d.o \ +s_mp_ntt_mont_add.o s_mp_ntt_mont_init.o s_mp_ntt_mont_mul.o s_mp_ntt_mont_pow.o s_mp_ntt_mont_redc.o \ +s_mp_ntt_mont_sub.o s_mp_ntt_small_is_prime.o s_mp_ntt_small_mulmod.o s_mp_ntt_small_powmod.o \ +s_mp_ntt_smooth_fwd.o s_mp_ntt_smooth_inv.o s_mp_ntt_smooth_length.o s_mp_ntt_smooth_length_weighted.o \ +s_mp_ntt_smooth_plan.o s_mp_ntt_sort_primes.o s_mp_ntt_sroots_free.o s_mp_ntt_sroots_init.o \ +s_mp_ntt_to_mont.o s_mp_prime_is_divisible.o s_mp_prime_tab.o s_mp_radix_map.o \ +s_mp_radix_size_overestimate.o s_mp_rand_platform.o s_mp_rand_source.o s_mp_sqr.o s_mp_sqr_comba.o \ +s_mp_sqr_karatsuba.o s_mp_sqr_ntt.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o s_mp_warray_get.o \ s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o - #END_INS $(LIBNAME): $(OBJECTS) diff --git a/makefile.mingw b/makefile.mingw index 10f2bd59..9de66275 100644 --- a/makefile.mingw +++ b/makefile.mingw @@ -49,11 +49,18 @@ mp_unpack.o mp_warray_free.o mp_xor.o mp_zero.o s_mp_add.o s_mp_copy_digs.o s_mp s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_fp_log.o \ s_mp_fp_log_d.o s_mp_get_bit.o s_mp_invmod.o s_mp_invmod_odd.o s_mp_log_2expt.o \ s_mp_montgomery_reduce_comba.o s_mp_mul.o s_mp_mul_balance.o s_mp_mul_comba.o s_mp_mul_high.o \ -s_mp_mul_high_comba.o s_mp_mul_karatsuba.o s_mp_mul_toom.o s_mp_prime_is_divisible.o s_mp_prime_tab.o \ -s_mp_radix_map.o s_mp_radix_size_overestimate.o s_mp_rand_platform.o s_mp_rand_source.o s_mp_sqr.o \ -s_mp_sqr_comba.o s_mp_sqr_karatsuba.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o s_mp_warray_get.o \ +s_mp_mul_high_comba.o s_mp_mul_karatsuba.o s_mp_mul_ntt.o s_mp_mul_toom.o s_mp_ntt_autoconvolve_single.o \ +s_mp_ntt_convolve_single.o s_mp_ntt_crt_assemble.o s_mp_ntt_crt_coeff.o s_mp_ntt_crt_init.o \ +s_mp_ntt_find_smooth_prime.o s_mp_ntt_from_mont.o s_mp_ntt_fwd2.o s_mp_ntt_fwd3.o s_mp_ntt_fwd5.o \ +s_mp_ntt_ilog2_ceil.o s_mp_ntt_inv2.o s_mp_ntt_inv3.o s_mp_ntt_inv5.o s_mp_ntt_invmod_d.o \ +s_mp_ntt_mont_add.o s_mp_ntt_mont_init.o s_mp_ntt_mont_mul.o s_mp_ntt_mont_pow.o s_mp_ntt_mont_redc.o \ +s_mp_ntt_mont_sub.o s_mp_ntt_small_is_prime.o s_mp_ntt_small_mulmod.o s_mp_ntt_small_powmod.o \ +s_mp_ntt_smooth_fwd.o s_mp_ntt_smooth_inv.o s_mp_ntt_smooth_length.o s_mp_ntt_smooth_length_weighted.o \ +s_mp_ntt_smooth_plan.o s_mp_ntt_sort_primes.o s_mp_ntt_sroots_free.o s_mp_ntt_sroots_init.o \ +s_mp_ntt_to_mont.o s_mp_prime_is_divisible.o s_mp_prime_tab.o s_mp_radix_map.o \ +s_mp_radix_size_overestimate.o s_mp_rand_platform.o s_mp_rand_source.o s_mp_sqr.o s_mp_sqr_comba.o \ +s_mp_sqr_karatsuba.o s_mp_sqr_ntt.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o s_mp_warray_get.o \ s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o - HEADERS_PUB=tommath.h HEADERS=tommath_private.h tommath_class.h tommath_superclass.h tommath_cutoffs.h $(HEADERS_PUB) diff --git a/makefile.msvc b/makefile.msvc index f160f2cc..1f015434 100644 --- a/makefile.msvc +++ b/makefile.msvc @@ -45,9 +45,17 @@ mp_unpack.obj mp_warray_free.obj mp_xor.obj mp_zero.obj s_mp_add.obj s_mp_copy_d s_mp_div_recursive.obj s_mp_div_school.obj s_mp_div_small.obj s_mp_exptmod.obj s_mp_exptmod_fast.obj s_mp_fp_log.obj \ s_mp_fp_log_d.obj s_mp_get_bit.obj s_mp_invmod.obj s_mp_invmod_odd.obj s_mp_log_2expt.obj \ s_mp_montgomery_reduce_comba.obj s_mp_mul.obj s_mp_mul_balance.obj s_mp_mul_comba.obj s_mp_mul_high.obj \ -s_mp_mul_high_comba.obj s_mp_mul_karatsuba.obj s_mp_mul_toom.obj s_mp_prime_is_divisible.obj s_mp_prime_tab.obj \ -s_mp_radix_map.obj s_mp_radix_size_overestimate.obj s_mp_rand_platform.obj s_mp_rand_source.obj s_mp_sqr.obj \ -s_mp_sqr_comba.obj s_mp_sqr_karatsuba.obj s_mp_sqr_toom.obj s_mp_sub.obj s_mp_warray.obj s_mp_warray_get.obj \ +s_mp_mul_high_comba.obj s_mp_mul_karatsuba.obj s_mp_mul_ntt.obj s_mp_mul_toom.obj s_mp_ntt_autoconvolve_single.obj \ +s_mp_ntt_convolve_single.obj s_mp_ntt_crt_assemble.obj s_mp_ntt_crt_coeff.obj s_mp_ntt_crt_init.obj \ +s_mp_ntt_find_smooth_prime.obj s_mp_ntt_from_mont.obj s_mp_ntt_fwd2.obj s_mp_ntt_fwd3.obj s_mp_ntt_fwd5.obj \ +s_mp_ntt_ilog2_ceil.obj s_mp_ntt_inv2.obj s_mp_ntt_inv3.obj s_mp_ntt_inv5.obj s_mp_ntt_invmod_d.obj \ +s_mp_ntt_mont_add.obj s_mp_ntt_mont_init.obj s_mp_ntt_mont_mul.obj s_mp_ntt_mont_pow.obj s_mp_ntt_mont_redc.obj \ +s_mp_ntt_mont_sub.obj s_mp_ntt_small_is_prime.obj s_mp_ntt_small_mulmod.obj s_mp_ntt_small_powmod.obj \ +s_mp_ntt_smooth_fwd.obj s_mp_ntt_smooth_inv.obj s_mp_ntt_smooth_length.obj s_mp_ntt_smooth_length_weighted.obj \ +s_mp_ntt_smooth_plan.obj s_mp_ntt_sort_primes.obj s_mp_ntt_sroots_free.obj s_mp_ntt_sroots_init.obj \ +s_mp_ntt_to_mont.obj s_mp_prime_is_divisible.obj s_mp_prime_tab.obj s_mp_radix_map.obj \ +s_mp_radix_size_overestimate.obj s_mp_rand_platform.obj s_mp_rand_source.obj s_mp_sqr.obj s_mp_sqr_comba.obj \ +s_mp_sqr_karatsuba.obj s_mp_sqr_ntt.obj s_mp_sqr_toom.obj s_mp_sub.obj s_mp_warray.obj s_mp_warray_get.obj \ s_mp_warray_put.obj s_mp_zero_buf.obj s_mp_zero_digs.obj HEADERS_PUB=tommath.h diff --git a/makefile.shared b/makefile.shared index 20bd8e03..e0cf9d78 100644 --- a/makefile.shared +++ b/makefile.shared @@ -44,13 +44,22 @@ mp_unpack.o mp_warray_free.o mp_xor.o mp_zero.o s_mp_add.o s_mp_copy_digs.o s_mp s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_fp_log.o \ s_mp_fp_log_d.o s_mp_get_bit.o s_mp_invmod.o s_mp_invmod_odd.o s_mp_log_2expt.o \ s_mp_montgomery_reduce_comba.o s_mp_mul.o s_mp_mul_balance.o s_mp_mul_comba.o s_mp_mul_high.o \ -s_mp_mul_high_comba.o s_mp_mul_karatsuba.o s_mp_mul_toom.o s_mp_prime_is_divisible.o s_mp_prime_tab.o \ -s_mp_radix_map.o s_mp_radix_size_overestimate.o s_mp_rand_platform.o s_mp_rand_source.o s_mp_sqr.o \ -s_mp_sqr_comba.o s_mp_sqr_karatsuba.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o s_mp_warray_get.o \ +s_mp_mul_high_comba.o s_mp_mul_karatsuba.o s_mp_mul_ntt.o s_mp_mul_toom.o s_mp_ntt_autoconvolve_single.o \ +s_mp_ntt_convolve_single.o s_mp_ntt_crt_assemble.o s_mp_ntt_crt_coeff.o s_mp_ntt_crt_init.o \ +s_mp_ntt_find_smooth_prime.o s_mp_ntt_from_mont.o s_mp_ntt_fwd2.o s_mp_ntt_fwd3.o s_mp_ntt_fwd5.o \ +s_mp_ntt_ilog2_ceil.o s_mp_ntt_inv2.o s_mp_ntt_inv3.o s_mp_ntt_inv5.o s_mp_ntt_invmod_d.o \ +s_mp_ntt_mont_add.o s_mp_ntt_mont_init.o s_mp_ntt_mont_mul.o s_mp_ntt_mont_pow.o s_mp_ntt_mont_redc.o \ +s_mp_ntt_mont_sub.o s_mp_ntt_small_is_prime.o s_mp_ntt_small_mulmod.o s_mp_ntt_small_powmod.o \ +s_mp_ntt_smooth_fwd.o s_mp_ntt_smooth_inv.o s_mp_ntt_smooth_length.o s_mp_ntt_smooth_length_weighted.o \ +s_mp_ntt_smooth_plan.o s_mp_ntt_sort_primes.o s_mp_ntt_sroots_free.o s_mp_ntt_sroots_init.o \ +s_mp_ntt_to_mont.o s_mp_prime_is_divisible.o s_mp_prime_tab.o s_mp_radix_map.o \ +s_mp_radix_size_overestimate.o s_mp_rand_platform.o s_mp_rand_source.o s_mp_sqr.o s_mp_sqr_comba.o \ +s_mp_sqr_karatsuba.o s_mp_sqr_ntt.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o s_mp_warray_get.o \ s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o #END_INS + objs: $(OBJECTS) .c.o: $(HEADERS) diff --git a/makefile.unix b/makefile.unix index 06dadddb..d4afe356 100644 --- a/makefile.unix +++ b/makefile.unix @@ -50,12 +50,19 @@ mp_unpack.o mp_warray_free.o mp_xor.o mp_zero.o s_mp_add.o s_mp_copy_digs.o s_mp s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_fp_log.o \ s_mp_fp_log_d.o s_mp_get_bit.o s_mp_invmod.o s_mp_invmod_odd.o s_mp_log_2expt.o \ s_mp_montgomery_reduce_comba.o s_mp_mul.o s_mp_mul_balance.o s_mp_mul_comba.o s_mp_mul_high.o \ -s_mp_mul_high_comba.o s_mp_mul_karatsuba.o s_mp_mul_toom.o s_mp_prime_is_divisible.o s_mp_prime_tab.o \ -s_mp_radix_map.o s_mp_radix_size_overestimate.o s_mp_rand_platform.o s_mp_rand_source.o s_mp_sqr.o \ -s_mp_sqr_comba.o s_mp_sqr_karatsuba.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o s_mp_warray_get.o \ +s_mp_mul_high_comba.o s_mp_mul_karatsuba.o s_mp_mul_ntt.o s_mp_mul_toom.o s_mp_ntt_autoconvolve_single.o \ +s_mp_ntt_convolve_single.o s_mp_ntt_crt_assemble.o s_mp_ntt_crt_coeff.o s_mp_ntt_crt_init.o \ +s_mp_ntt_find_smooth_prime.o s_mp_ntt_from_mont.o s_mp_ntt_fwd2.o s_mp_ntt_fwd3.o s_mp_ntt_fwd5.o \ +s_mp_ntt_ilog2_ceil.o s_mp_ntt_inv2.o s_mp_ntt_inv3.o s_mp_ntt_inv5.o s_mp_ntt_invmod_d.o \ +s_mp_ntt_mont_add.o s_mp_ntt_mont_init.o s_mp_ntt_mont_mul.o s_mp_ntt_mont_pow.o s_mp_ntt_mont_redc.o \ +s_mp_ntt_mont_sub.o s_mp_ntt_small_is_prime.o s_mp_ntt_small_mulmod.o s_mp_ntt_small_powmod.o \ +s_mp_ntt_smooth_fwd.o s_mp_ntt_smooth_inv.o s_mp_ntt_smooth_length.o s_mp_ntt_smooth_length_weighted.o \ +s_mp_ntt_smooth_plan.o s_mp_ntt_sort_primes.o s_mp_ntt_sroots_free.o s_mp_ntt_sroots_init.o \ +s_mp_ntt_to_mont.o s_mp_prime_is_divisible.o s_mp_prime_tab.o s_mp_radix_map.o \ +s_mp_radix_size_overestimate.o s_mp_rand_platform.o s_mp_rand_source.o s_mp_sqr.o s_mp_sqr_comba.o \ +s_mp_sqr_karatsuba.o s_mp_sqr_ntt.o s_mp_sqr_toom.o s_mp_sub.o s_mp_warray.o s_mp_warray_get.o \ s_mp_warray_put.o s_mp_zero_buf.o s_mp_zero_digs.o - HEADERS_PUB=tommath.h HEADERS=tommath_private.h tommath_class.h tommath_superclass.h tommath_cutoffs.h $(HEADERS_PUB) diff --git a/mp_cutoffs.c b/mp_cutoffs.c index 45b0beec..91547361 100644 --- a/mp_cutoffs.c +++ b/mp_cutoffs.c @@ -8,7 +8,9 @@ int MP_MUL_KARATSUBA_CUTOFF = MP_DEFAULT_MUL_KARATSUBA_CUTOFF, MP_SQR_KARATSUBA_CUTOFF = MP_DEFAULT_SQR_KARATSUBA_CUTOFF, MP_MUL_TOOM_CUTOFF = MP_DEFAULT_MUL_TOOM_CUTOFF, - MP_SQR_TOOM_CUTOFF = MP_DEFAULT_SQR_TOOM_CUTOFF; + MP_SQR_TOOM_CUTOFF = MP_DEFAULT_SQR_TOOM_CUTOFF, + MP_MUL_NTT_CUTOFF = MP_DEFAULT_MUL_NTT_CUTOFF, + MP_SQR_NTT_CUTOFF = MP_DEFAULT_SQR_NTT_CUTOFF; #endif #endif diff --git a/mp_mul.c b/mp_mul.c index 81807406..e7c5ed87 100644 --- a/mp_mul.c +++ b/mp_mul.c @@ -12,56 +12,85 @@ mp_err mp_mul(const mp_int *a, const mp_int *b, mp_int *c) digs = a->used + b->used + 1; bool neg = (a->sign != b->sign); +#if (MP_DIGIT_BIT > 16) if ((a == b) && - MP_HAS(S_MP_SQR_TOOM) && /* use Toom-Cook? */ - (a->used >= MP_SQR_TOOM_CUTOFF)) { - err = s_mp_sqr_toom(a, c); - } else if ((a == b) && - MP_HAS(S_MP_SQR_KARATSUBA) && /* Karatsuba? */ - (a->used >= MP_SQR_KARATSUBA_CUTOFF)) { - err = s_mp_sqr_karatsuba(a, c); - } else if ((a == b) && - MP_HAS(S_MP_SQR_COMBA) && /* can we use the fast comba multiplier? */ - (((a->used * 2) + 1) < MP_WARRAY) && - (a->used <= MP_MAX_COMBA)) { - err = s_mp_sqr_comba(a, c); - } else if ((a == b) && - MP_HAS(S_MP_SQR)) { - err = s_mp_sqr(a, c); - } else if (MP_HAS(S_MP_MUL_BALANCE) && - /* Check sizes. The smaller one needs to be larger than the Karatsuba cut-off. - * The bigger one needs to be at least about one MP_MUL_KARATSUBA_CUTOFF bigger - * to make some sense, but it depends on architecture, OS, position of the - * stars... so YMMV. - * Using it to cut the input into slices small enough for s_mp_mul_comba - * was actually slower on the author's machine, but YMMV. - */ - (min >= MP_MUL_KARATSUBA_CUTOFF) && - ((max / 2) >= MP_MUL_KARATSUBA_CUTOFF) && - /* Not much effect was observed below a ratio of 1:2, but again: YMMV. */ - (max >= (2 * min))) { - err = s_mp_mul_balance(a,b,c); - } else if (MP_HAS(S_MP_MUL_TOOM) && - (min >= MP_MUL_TOOM_CUTOFF)) { - err = s_mp_mul_toom(a, b, c); - } else if (MP_HAS(S_MP_MUL_KARATSUBA) && - (min >= MP_MUL_KARATSUBA_CUTOFF)) { - err = s_mp_mul_karatsuba(a, b, c); - } else if (MP_HAS(S_MP_MUL_COMBA) && - /* can we use the fast multiplier? - * - * The fast multiplier can be used if the output will - * have less than MP_WARRAY digits and the number of - * digits won't affect carry propagation - */ - (digs < MP_WARRAY) && - (min <= MP_MAX_COMBA)) { - err = s_mp_mul_comba(a, b, c, digs); - } else if (MP_HAS(S_MP_MUL)) { - err = s_mp_mul(a, b, c, digs); - } else { - err = MP_VAL; - } + MP_HAS(S_MP_SQR_NTT) && /* use NTT? */ + (a->used >= MP_SQR_NTT_CUTOFF)) { + err = s_mp_sqr_ntt(a, c); + } else +#endif + if ((a == b) && + MP_HAS(S_MP_SQR_TOOM) && /* use Toom-Cook? */ + (a->used >= MP_SQR_TOOM_CUTOFF)) { + err = s_mp_sqr_toom(a, c); + } else if ((a == b) && + MP_HAS(S_MP_SQR_KARATSUBA) && /* Karatsuba? */ + (a->used >= MP_SQR_KARATSUBA_CUTOFF)) { + err = s_mp_sqr_karatsuba(a, c); + } else if ((a == b) && + MP_HAS(S_MP_SQR_COMBA) && /* can we use the fast comba multiplier? */ + (((a->used * 2) + 1) < MP_WARRAY) && + (a->used <= MP_MAX_COMBA)) { + err = s_mp_sqr_comba(a, c); + } else if ((a == b) && + MP_HAS(S_MP_SQR)) { + err = s_mp_sqr(a, c); + /* NTT does not need balancing, it is even faster for unbalanced input */ + /* TODO: NTT is much faster when size(a) != size(b). + If one size divides the other it is the normal cutoff, e.g., + k * (nk). but when it is a bit off it can be much faster e.g., + k * (nk \pm k/epsilon) with epsilon < ~100. + + Much faster means about seven times faster on the author's machine. + + Acting according to the behaviour described above is not yet implemented. + + Would need to test the size of epsilon (both k+k/epsilon and k-k/epsilon) + at a couple of points. Maybe 100, 50, 10, 5, and 2. Including the symmetry + that would be 10 runs of the benchmark up to about 8,000 limbs. + That is a lot. + */ + } else if (MP_HAS(S_MP_MUL_BALANCE) && (max < MP_MUL_NTT_CUTOFF) && + /* Check sizes. The smaller one needs to be larger than the Karatsuba cut-off. + * The bigger one needs to be at least about one MP_MUL_KARATSUBA_CUTOFF bigger + * to make some sense, but it depends on architecture, OS, position of the + * stars... so YMMV. + * Using it to cut the input into slices small enough for s_mp_mul_comba + * was actually slower on the author's machine, but YMMV. + */ + (min >= MP_MUL_KARATSUBA_CUTOFF) && + ((max / 2) >= MP_MUL_KARATSUBA_CUTOFF) && + /* Not much effect was observed below a ratio of 1:2, but again: YMMV. */ + (max >= (2 * min))) { + err = s_mp_mul_balance(a,b,c); + } +#if (MP_DIGIT_BIT > 16) + else if (MP_HAS(S_MP_MUL_NTT) && + (min >= MP_MUL_NTT_CUTOFF)) { + err = s_mp_mul_ntt(a, b, c); + } +#endif + else if (MP_HAS(S_MP_MUL_TOOM) && + (min >= MP_MUL_TOOM_CUTOFF)) { + err = s_mp_mul_toom(a, b, c); + } else if (MP_HAS(S_MP_MUL_KARATSUBA) && + (min >= MP_MUL_KARATSUBA_CUTOFF)) { + err = s_mp_mul_karatsuba(a, b, c); + } else if (MP_HAS(S_MP_MUL_COMBA) && + /* can we use the fast multiplier? + * + * The fast multiplier can be used if the output will + * have less than MP_WARRAY digits and the number of + * digits won't affect carry propagation + */ + (digs < MP_WARRAY) && + (min <= MP_MAX_COMBA)) { + err = s_mp_mul_comba(a, b, c, digs); + } else if (MP_HAS(S_MP_MUL)) { + err = s_mp_mul(a, b, c, digs); + } else { + err = MP_VAL; + } c->sign = ((c->used > 0) && neg) ? MP_NEG : MP_ZPOS; return err; } diff --git a/s_mp_mul_ntt.c b/s_mp_mul_ntt.c new file mode 100644 index 00000000..98a4ebc4 --- /dev/null +++ b/s_mp_mul_ntt.c @@ -0,0 +1,132 @@ +#include "tommath_private.h" +#ifdef S_MP_MUL_NTT_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +/* The main function to multiply two mp_int's into another mp_int */ +mp_err s_mp_mul_ntt(const mp_int *a, const mp_int *b, mp_int *c) +{ + int ncoef, m, minl, bound_bits, nprimes, pi, i; + ntt_prime_ent *ents = NULL; + ntt_crt_ctx *cx = NULL; + mp_digit *primes = NULL, **residues = NULL, *ra = NULL, *rb = NULL; + mp_int P, pmp; + mp_err err = MP_OKAY; + + if (mp_iszero(a) || mp_iszero(b)) { + return mp_mul(a, b, c); + } + if (a->used < 8 || b->used < 8) { + return mp_mul(a, b, c); + } + + ncoef = a->used + b->used - 1; /* coefficients that carry the product */ + if ((ncoef > (1 << NTT_QUARTER_COEFF)) || (ncoef < 0)) { + return MP_VAL; + } + m = s_mp_ntt_smooth_length(ncoef); /* smallest 5-smooth transform length */ + if (m < ncoef) { + return MP_VAL; + } + + /* Convolution bound: each output coefficient < min(la,lb) * 2^120. */ + minl = (a->used < b->used) ? a->used : b->used; + bound_bits = (2 * NTT_LIMB_SIZE) + s_mp_ntt_ilog2_ceil(minl) + 1; + if (((bound_bits + NTT_LIMB_SIZE_M_ONE) / NTT_LIMB_SIZE) > NTT_MAXW) { + return MP_VAL; + } + + ents = (ntt_prime_ent *)MP_MALLOC((size_t)NTT_MAXP * sizeof(ntt_prime_ent)); + if (ents == NULL) { + return MP_MEM; + } + + if ((err = mp_init_multi(&P, &pmp, NULL)) != MP_OKAY) { + MP_FREE(ents, (size_t)NTT_MAXP * sizeof(ntt_prime_ent)); + return err; + } + /* TODO */ + mp_set(&P, 1); + + /* Find primes until their product > 2^bound_bits. should result in three primes */ + nprimes = 0; + while (mp_count_bits(&P) <= bound_bits) { + if (nprimes >= NTT_MAXP) { + err = MP_VAL; + goto LTM_ERR; + } + if (!s_mp_ntt_find_smooth_prime(&ents[nprimes], m, nprimes, 0)) { + err = MP_VAL; + goto LTM_ERR; + } + mp_set(&pmp, ents[nprimes].p); + if ((err = mp_mul(&P, &pmp, &P)) != MP_OKAY) goto LTM_ERR; + nprimes++; + } + + + s_mp_ntt_sort_primes(ents, nprimes); + + ra = (mp_digit *)MP_MALLOC((size_t)m * sizeof(mp_digit)); + if (ra == NULL) { + err = MP_MEM; + goto LTM_ERR; + } + rb = (mp_digit *)MP_MALLOC((size_t)m * sizeof(mp_digit)); + if (rb == NULL) { + err = MP_MEM; + goto LTM_ERR; + } + residues = (mp_digit **)MP_CALLOC((size_t)nprimes, sizeof(mp_digit *)); + if (residues == NULL) { + err = MP_MEM; + goto LTM_ERR; + } + primes = (mp_digit *)MP_MALLOC((size_t)nprimes * sizeof(mp_digit)); + if (primes == NULL) { + err = MP_MEM; + goto LTM_ERR; + } + + + /* convolve over each prime and keep the residues */ + for (pi = 0; pi < nprimes; pi++) { + residues[pi] = (mp_digit *)MP_MALLOC((size_t)ncoef * sizeof(mp_digit)); + if (residues[pi] == NULL) { + err = MP_MEM; + goto LTM_ERR; + } + if ((err = s_mp_ntt_convolve_single(a, b, &ents[pi], ra, rb, residues[pi], ncoef)) != MP_OKAY) goto LTM_ERR; + primes[pi] = ents[pi].p; + } + + /* CRT: residues -> binary, then one linear carry sweep */ + cx = (ntt_crt_ctx *)MP_MALLOC(sizeof(ntt_crt_ctx)); + if (cx == NULL) { + err = MP_MEM; + goto LTM_ERR; + } + s_mp_ntt_crt_init(cx, primes, nprimes, bound_bits); + err = s_mp_ntt_crt_assemble(cx, residues, nprimes, ncoef, c); + +LTM_ERR: + mp_clear_multi(&P, &pmp, NULL); + if (residues) { + for (i = 0; i < nprimes; i++) { + MP_FREE(residues[i], (size_t)ncoef * sizeof(mp_digit)); + } + MP_FREE(residues, (size_t)nprimes * sizeof(mp_digit *)); + } + MP_FREE(primes, (size_t)nprimes * sizeof(mp_digit)); + MP_FREE(ra, (size_t)m * sizeof(mp_digit)); + MP_FREE(rb, (size_t)m * sizeof(mp_digit)); + MP_FREE(cx, sizeof(ntt_crt_ctx)); + MP_FREE(ents, (size_t)NTT_MAXP * sizeof(ntt_prime_ent)); + return err; +} +#endif + +#endif diff --git a/s_mp_ntt_autoconvolve_single.c b/s_mp_ntt_autoconvolve_single.c new file mode 100644 index 00000000..740af535 --- /dev/null +++ b/s_mp_ntt_autoconvolve_single.c @@ -0,0 +1,47 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_AUTOCONVOLVE_SINGLE_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +mp_err s_mp_ntt_autoconvolve_single(const mp_int *a, const ntt_prime_ent *ent, mp_digit *ra, mp_digit *out, int ncoef) +{ + int m = ent->m, i; + ntt_mont mc; + ntt_sroots R; + mp_digit zeta; + mp_err err = MP_OKAY; + + s_mp_ntt_mont_init(&mc, ent->p); + zeta = s_mp_ntt_to_mont(&mc, ent->root); /* primitive m-th root (Montgomery) */ + + if ((err = s_mp_ntt_sroots_init(&R, &mc, m, zeta)) != MP_OKAY) return err; + + for (i = 0; i < m; i++) { + ra[i] = 0; + } + for (i = 0; i < a->used; i++) { + ra[i] = s_mp_ntt_to_mont(&mc, a->dp[i]); + } + + s_mp_ntt_smooth_fwd(&mc, ra, m, &R); + + for (i = 0; i < m; i++) { + ra[i] = s_mp_ntt_mont_mul(&mc, ra[i], ra[i]); + } + + s_mp_ntt_smooth_inv(&mc, ra, m, &R); /* includes the 1/m */ + + for (i = 0; i < ncoef; i++) { + out[i] = s_mp_ntt_from_mont(&mc, ra[i]); + } + + s_mp_ntt_sroots_free(&R); + return err; +} +#endif + +#endif diff --git a/s_mp_ntt_convolve_single.c b/s_mp_ntt_convolve_single.c new file mode 100644 index 00000000..5255d79e --- /dev/null +++ b/s_mp_ntt_convolve_single.c @@ -0,0 +1,63 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_CONVOLVE_SINGLE_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +/* + One prime's contribution: forward-transform both inputs, multiply pointwise, + inverse-transform, and return the convolution coefficients modulo that prime. + + The transform runs at m = ent->m (5-smooth, >= la + lb - 1), so its cyclic + convolution of length m agrees with the linear one on the first + ncoef = la + lb - 1 coefficients and is zero beyond - which is why only those + are read back. + */ +mp_err s_mp_ntt_convolve_single(const mp_int *a, const mp_int *b, const ntt_prime_ent *ent, mp_digit *ra, + mp_digit *rb, mp_digit *out, int ncoef) +{ + int m = ent->m, i; + ntt_mont mc; + ntt_sroots R; + mp_digit zeta; + mp_err err = MP_OKAY; + + s_mp_ntt_mont_init(&mc, ent->p); + zeta = s_mp_ntt_to_mont(&mc, ent->root); /* primitive m-th root (Montgomery) */ + + if ((err = s_mp_ntt_sroots_init(&R, &mc, m, zeta)) != MP_OKAY) return err; + + for (i = 0; i < m; i++) { + ra[i] = 0; + rb[i] = 0; + } + for (i = 0; i < a->used; i++) { + ra[i] = s_mp_ntt_to_mont(&mc, a->dp[i]); + } + for (i = 0; i < b->used; i++) { + rb[i] = s_mp_ntt_to_mont(&mc, b->dp[i]); + } + + s_mp_ntt_smooth_fwd(&mc, ra, m, &R); + s_mp_ntt_smooth_fwd(&mc, rb, m, &R); + + for (i = 0; i < m; i++) { + ra[i] = s_mp_ntt_mont_mul(&mc, ra[i], rb[i]); + } + + s_mp_ntt_smooth_inv(&mc, ra, m, &R); /* includes the 1/m */ + + for (i = 0; i < ncoef; i++) { + out[i] = s_mp_ntt_from_mont(&mc, ra[i]); + } + + s_mp_ntt_sroots_free(&R); + return err; +} +#endif + + +#endif diff --git a/s_mp_ntt_crt_assemble.c b/s_mp_ntt_crt_assemble.c new file mode 100644 index 00000000..ff0bb69f --- /dev/null +++ b/s_mp_ntt_crt_assemble.c @@ -0,0 +1,66 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_CRT_ASSEMBLE_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +/* + Residues -> the product. residues[i] holds the first ncoef convolution + coefficients modulo primes[i]; the result is written into c. + */ +mp_err s_mp_ntt_crt_assemble(ntt_crt_ctx *cx, mp_digit **residues, int nprimes, + int ncoef, mp_int *c) +{ + mp_err err = MP_OKAY; + mp_digit *r, *acc, limbs[NTT_MAXW]; + mp_digit carry; + int k, i, w, n; + mp_word v; + + r = (mp_digit *)MP_MALLOC((size_t)nprimes * sizeof(mp_digit)); + if (r == NULL) { + return MP_MEM; + } + /* lazy accumulator: one 60-bit slot per limb position; a slot receives at + * most one limb from each of the <= NTT_MAXW overlapping coefficients, so it + * stays below NTT_MAXW * 2^60 and cannot overflow */ + acc = (mp_digit *)MP_CALLOC((size_t)ncoef + NTT_MAXW, sizeof(mp_digit)); + if (acc == NULL) { + MP_FREE(r, (size_t)nprimes * sizeof(mp_digit)); + return MP_MEM; + } + + for (k = 0; k < ncoef; k++) { + for (i = 0; i < nprimes; i++) { + r[i] = residues[i][k]; + } + s_mp_ntt_crt_coeff(cx, r, limbs); + for (w = 0; w < cx->width; w++) { + acc[k + w] += limbs[w]; + } + } + + n = ncoef + NTT_MAXW; + if ((err = mp_grow(c, n)) != MP_OKAY) goto LTM_ERR; + carry = 0; + for (k = 0; k < n; k++) { + v = (mp_word)acc[k] + carry; + c->dp[k] = (mp_digit)(v & NTT_MASK); + carry = (mp_digit)(v >> NTT_LIMB_SIZE); + } + c->used = n; + c->sign = MP_ZPOS; + mp_clamp(c); + +LTM_ERR: + MP_FREE(r, (size_t)nprimes * sizeof(mp_digit)); + MP_FREE(acc, ((size_t)ncoef + NTT_MAXW) * sizeof(mp_digit)); + return err; +} +#endif + + +#endif diff --git a/s_mp_ntt_crt_coeff.c b/s_mp_ntt_crt_coeff.c new file mode 100644 index 00000000..697ebf54 --- /dev/null +++ b/s_mp_ntt_crt_coeff.c @@ -0,0 +1,56 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_CRT_COEFF_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +/* + Recover one coefficient (from its residue modulo every prime) as a short + array of 60-bit limbs. + + minv[i][j] = p_i^{-1} * 2^64 mod p_j, for i < j. Propagating the digit d_i to + the higher residues divides by p_i, so the update is + x_j <- (x_j - d_i) * p_i^{-1} mod p_j. + (M_i^{-1} = (p_0 ... p_{i-1})^{-1}, not p_i^{-1}, appears only when the partial + value rather than the digit is propagated.) Keeping p_i^{-1} in Montgomery + form makes that update one multiplication with no division. + */ +void s_mp_ntt_crt_coeff(ntt_crt_ctx *cx, const mp_digit *r, mp_digit *limbs) +{ + int t = cx->t, i, j, w; + mp_digit *x = cx->x; + mp_digit *d = cx->d; + mp_digit di, pj, diff; + mp_word carry, v; + + for (i = 0; i < t; i++) x[i] = r[i]; /* each < p_i */ + + for (i = 0; i < t; i++) { + di = x[i]; /* mixed-radix digit */ + d[i] = di; + for (j = i + 1; j < t; j++) { + pj = cx->p[j]; + diff = (x[j] >= di) ? (x[j] - di) : (x[j] + pj - di); + x[j] = s_mp_ntt_mont_mul(&cx->mc[j], diff, cx->minv[i][j]); + } + } + + /* Horner: c = d_0 + p_0 (d_1 + p_1 (d_2 + ... )) */ + for (w = 0; w < NTT_MAXW; w++) limbs[w] = 0; + limbs[0] = d[t - 1]; + for (i = t - 2; i >= 0; i--) { + carry = (mp_word)d[i]; + for (w = 0; w < NTT_MAXW; w++) { + v = (mp_word)limbs[w] * cx->p[i] + carry; + limbs[w] = (mp_digit)(v & NTT_MASK); + carry = v >> NTT_LIMB_SIZE; + } + } +} +#endif + + +#endif diff --git a/s_mp_ntt_crt_init.c b/s_mp_ntt_crt_init.c new file mode 100644 index 00000000..56c46d23 --- /dev/null +++ b/s_mp_ntt_crt_init.c @@ -0,0 +1,34 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_CRT_INIT_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +void s_mp_ntt_crt_init(ntt_crt_ctx *cx, const mp_digit *primes, int t, int bound_bits) +{ + int i, j; + + cx->t = t; + cx->width = (bound_bits + NTT_LIMB_SIZE_M_ONE) / NTT_LIMB_SIZE; + if (cx->width < 1) cx->width = 1; + + for (i = 0; i < t; i++) { + cx->p[i] = primes[i]; + s_mp_ntt_mont_init(&cx->mc[i], primes[i]); + } + for (i = 0; i < t; i++) { + for (j = i + 1; j < t; j++) { + mp_digit pj = cx->p[j]; + mp_digit inv = s_mp_ntt_small_powmod(cx->p[i] % pj, pj - 2, pj); /* p_i^{-1} mod p_j */ + cx->minv[i][j] = s_mp_ntt_to_mont(&cx->mc[j], inv); + } + } +} +#endif + + + +#endif diff --git a/s_mp_ntt_find_smooth_prime.c b/s_mp_ntt_find_smooth_prime.c new file mode 100644 index 00000000..b8e9dd44 --- /dev/null +++ b/s_mp_ntt_find_smooth_prime.c @@ -0,0 +1,83 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_FIND_SMOOTH_PRIME_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +/* + Search for the `index`-th largest usable prime of a 5-smooth transform length + m, i.e. p = k*m + 1 < 2^62 (index = 0 is the largest), requiring p > min_p. + + Largest rather than smallest for the same reason as ntt_find_large_prime: + three ~62-bit primes cover the whole convolution bound at every supported + size, so the prime count never grows with the input. + + The root search differs. For the truncated path a quadratic non-residue gives + a 2^s-th root for free (Euler's criterion), but m has 3s and 5s in it and no + single test covers those. So: w = a^((p-1)/m) has order dividing m, and it is + exactly m as soon as w^(m/r) != 1 for every prime r dividing m - because every + proper divisor of m divides one of those m/r. + */ +bool s_mp_ntt_find_smooth_prime(ntt_prime_ent *ent, int m, int index, mp_digit min_p) +{ + mp_digit mm, kmax, root; + int seen = 0; + mp_digit k, a, p, w; + + if (m < 2) { + return false; + } + mm = (mp_digit)m; + kmax = (((mp_digit)1 << NTT_PRIME_MAX_SIZE) - 1) / mm; + + for (k = kmax; k >= 1; k--) { + p = k * mm + 1; + if (p >= ((mp_digit)1 << NTT_PRIME_MAX_SIZE)) { + continue; + } + if (!s_mp_ntt_small_is_prime(p)) { + continue; + } + if (p <= min_p) { + return false; /* cannot cover the bound */ + } + if (seen++ != index) { + continue; /* want the index-th largest */ + } + + root = 0; + /* Max. "a" measured was 48 with 60-bit */ + for (a = 2; (a < 1000) && (root == 0); a++) { + w = s_mp_ntt_small_powmod(a, k, p); /* order divides m */ + if (w == 1) { + continue; + } + if (((m % 2) == 0) && s_mp_ntt_small_powmod(w, (mp_digit)m / 2, p) == 1) { + continue; + } + if (((m % 3) == 0) && s_mp_ntt_small_powmod(w, (mp_digit)m / 3, p) == 1) { + continue; + } + if (((m % 5) == 0) && s_mp_ntt_small_powmod(w, (mp_digit)m / 5, p) == 1) { + continue; + } + root = w; + } + if (root == 0) { + return false; /* no root of full order */ + } + ent->p = p; + ent->root = root; + ent->s = 0; + ent->m = m; + return true; + } + return false; +} +#endif + + +#endif diff --git a/s_mp_ntt_from_mont.c b/s_mp_ntt_from_mont.c new file mode 100644 index 00000000..217f9047 --- /dev/null +++ b/s_mp_ntt_from_mont.c @@ -0,0 +1,17 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_FROM_MONT_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +mp_digit s_mp_ntt_from_mont(const ntt_mont *m, mp_digit a) +{ + return s_mp_ntt_mont_redc(m, (mp_word)a); +} +#endif + + +#endif diff --git a/s_mp_ntt_fwd2.c b/s_mp_ntt_fwd2.c new file mode 100644 index 00000000..40f246bb --- /dev/null +++ b/s_mp_ntt_fwd2.c @@ -0,0 +1,35 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_FWD2_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +/* + Forward: transform first and do the twiddles after. + Inverse: undoes the twiddles first (reading the same table at m - e, + since omega^{-e} = omega^{m-e}) and inverse transform after. + The inverse is unscaled, because s_mp_ntt_smooth_inv divides by m once at + the end. Slot q carries index q in both directions. + */ +void s_mp_ntt_fwd2(const ntt_mont *mc, mp_digit *b, int L, int step, const ntt_sroots *R) +{ + const mp_digit *w = R->w; + int m = R->m, e = 0, j; + + for (j = 0; j < L; j++) { + mp_digit y0 = b[j], y1 = b[j + L]; + b[j] = s_mp_ntt_mont_add(mc, y0, y1); /* q = 0, W^0 = 1 */ + b[j + L] = s_mp_ntt_mont_mul(mc, s_mp_ntt_mont_sub(mc, y0, y1), w[e]); + e += step; + if (e >= m) { + e -= m; + } + } +} +#endif + + +#endif diff --git a/s_mp_ntt_fwd3.c b/s_mp_ntt_fwd3.c new file mode 100644 index 00000000..d2416233 --- /dev/null +++ b/s_mp_ntt_fwd3.c @@ -0,0 +1,41 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_FWD3_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +void s_mp_ntt_fwd3(const ntt_mont *mc, mp_digit *b, int L, int step, const ntt_sroots *R) +{ + const mp_digit *w = R->w; + int m = R->m, j; + int s1 = step, s2 = 2 * step; /* step = m/ell <= m/3, so 2*step < m */ + int e1 = 0, e2 = 0; + mp_digit t = R->c3; + + for (j = 0; j < L; j++) { + mp_digit y0 = b[j], y1 = b[j + L], y2 = b[j + 2 * L]; + mp_digit v = s_mp_ntt_mont_mul(mc, s_mp_ntt_mont_sub(mc, y1, y2), t); /* t*(y1-y2) */ + mp_digit z0 = s_mp_ntt_mont_add(mc, s_mp_ntt_mont_add(mc, y0, y1), y2); + mp_digit z1 = s_mp_ntt_mont_add(mc, s_mp_ntt_mont_sub(mc, y0, y2), v); + mp_digit z2 = s_mp_ntt_mont_sub(mc, s_mp_ntt_mont_sub(mc, y0, y1), v); + + b[j] = z0; + b[j + L] = s_mp_ntt_mont_mul(mc, z1, w[e1]); + b[j + 2 * L] = s_mp_ntt_mont_mul(mc, z2, w[e2]); + e1 += s1; + if (e1 >= m) { + e1 -= m; + } + e2 += s2; + if (e2 >= m) { + e2 -= m; + } + } +} +#endif + + +#endif diff --git a/s_mp_ntt_fwd5.c b/s_mp_ntt_fwd5.c new file mode 100644 index 00000000..333c8fc7 --- /dev/null +++ b/s_mp_ntt_fwd5.c @@ -0,0 +1,73 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_FWD5_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +void s_mp_ntt_fwd5(const ntt_mont *mc, mp_digit *b, int L, int step, const ntt_sroots *R) +{ + const mp_digit *w = R->w; + int m = R->m, j; + int s1 = step, s2 = 2 * step, s3 = 3 * step, s4 = 4 * step; /* all < m */ + int e1 = 0, e2 = 0, e3 = 0, e4 = 0; + mp_digit k0 = R->k5[0], k1 = R->k5[1], k2 = R->k5[2], k3 = R->k5[3]; + mp_digit y0, y1, y2, y3, y4, t1, t2, u1, u2, t1c1, t2c2, t1c2, t2c1, u1s1, u2s12, u1s12, u2s1, z1, z2, z3, z4; + + for (j = 0; j < L; j++) { + y0 = b[j]; + y1 = b[j + L]; + y2 = b[j + 2 * L]; + y3 = b[j + 3 * L]; + y4 = b[j + 4 * L]; + t1 = s_mp_ntt_mont_add(mc, y1, y4); + u1 = s_mp_ntt_mont_sub(mc, y1, y4); + t2 = s_mp_ntt_mont_add(mc, y2, y3); + u2 = s_mp_ntt_mont_sub(mc, y2, y3); + t1c1 = s_mp_ntt_mont_mul(mc, t1, k0); + u1s1 = s_mp_ntt_mont_mul(mc, u1, k2); + t2c2 = s_mp_ntt_mont_mul(mc, t2, k1); + u2s12 = s_mp_ntt_mont_mul(mc, u2, k3); + t1c2 = s_mp_ntt_mont_mul(mc, t1, k1); + u1s12 = s_mp_ntt_mont_mul(mc, u1, k3); + t2c1 = s_mp_ntt_mont_mul(mc, t2, k0); + u2s1 = s_mp_ntt_mont_mul(mc, u2, k2); + z1 = s_mp_ntt_mont_add(mc, s_mp_ntt_mont_add(mc, y0, s_mp_ntt_mont_add(mc, t1c1, u1s1)), s_mp_ntt_mont_add(mc, t2c2, + u2s12)); + z2 = s_mp_ntt_mont_add(mc, s_mp_ntt_mont_add(mc, y0, s_mp_ntt_mont_add(mc, t1c2, u1s12)), s_mp_ntt_mont_sub(mc, t2c1, + u2s1)); + z3 = s_mp_ntt_mont_add(mc, s_mp_ntt_mont_add(mc, y0, s_mp_ntt_mont_sub(mc, t1c2, u1s12)), s_mp_ntt_mont_add(mc, t2c1, + u2s1)); + z4 = s_mp_ntt_mont_add(mc, s_mp_ntt_mont_add(mc, y0, s_mp_ntt_mont_sub(mc, t1c1, u1s1)), s_mp_ntt_mont_sub(mc, t2c2, + u2s12)); + + b[j] = s_mp_ntt_mont_add(mc, s_mp_ntt_mont_add(mc, y0, t1), t2); + b[j + L] = s_mp_ntt_mont_mul(mc, z1, w[e1]); + b[j + 2 * L] = s_mp_ntt_mont_mul(mc, z2, w[e2]); + b[j + 3 * L] = s_mp_ntt_mont_mul(mc, z3, w[e3]); + b[j + 4 * L] = s_mp_ntt_mont_mul(mc, z4, w[e4]); + e1 += s1; + if (e1 >= m) { + e1 -= m; + } + e2 += s2; + if (e2 >= m) { + e2 -= m; + } + e3 += s3; + if (e3 >= m) { + e3 -= m; + } + e4 += s4; + if (e4 >= m) { + e4 -= m; + } + } +} +#endif + + + +#endif diff --git a/s_mp_ntt_ilog2_ceil.c b/s_mp_ntt_ilog2_ceil.c new file mode 100644 index 00000000..0146b282 --- /dev/null +++ b/s_mp_ntt_ilog2_ceil.c @@ -0,0 +1,21 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_ILOG2_CEIL_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +int s_mp_ntt_ilog2_ceil(int x) +{ + int s = 0; + while ((1 << s) < x) { + s++; + } + return s; +} +#endif + + +#endif diff --git a/s_mp_ntt_inv2.c b/s_mp_ntt_inv2.c new file mode 100644 index 00000000..fe1a0166 --- /dev/null +++ b/s_mp_ntt_inv2.c @@ -0,0 +1,30 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_INV2_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +void s_mp_ntt_inv2(const ntt_mont *mc, mp_digit *b, int L, int step, const ntt_sroots *R) +{ + const mp_digit *w = R->w; + int m = R->m, e = 0, j; + + for (j = 0; j < L; j++) { + mp_digit y0 = b[j]; + mp_digit y1 = s_mp_ntt_mont_mul(mc, b[j + L], w[e ? m - e : 0]); + b[j] = s_mp_ntt_mont_add(mc, y0, y1); + b[j + L] = s_mp_ntt_mont_sub(mc, y0, y1); + e += step; + if (e >= m) { + e -= m; + } + } +} +#endif + + + +#endif diff --git a/s_mp_ntt_inv3.c b/s_mp_ntt_inv3.c new file mode 100644 index 00000000..ec7c4f26 --- /dev/null +++ b/s_mp_ntt_inv3.c @@ -0,0 +1,41 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_INV3_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +void s_mp_ntt_inv3(const ntt_mont *mc, mp_digit *b, int L, int step, const ntt_sroots *R) +{ + const mp_digit *w = R->w; + int m = R->m, j; + int s1 = step, s2 = 2 * step; + int e1 = 0, e2 = 0; + mp_digit t = R->c3i; + + for (j = 0; j < L; j++) { + mp_digit y0 = b[j]; + mp_digit y1 = s_mp_ntt_mont_mul(mc, b[j + L], w[e1 ? m - e1 : 0]); + mp_digit y2 = s_mp_ntt_mont_mul(mc, b[j + 2 * L], w[e2 ? m - e2 : 0]); + mp_digit v = s_mp_ntt_mont_mul(mc, s_mp_ntt_mont_sub(mc, y1, y2), t); + + b[j] = s_mp_ntt_mont_add(mc, s_mp_ntt_mont_add(mc, y0, y1), y2); + b[j + L] = s_mp_ntt_mont_add(mc, s_mp_ntt_mont_sub(mc, y0, y2), v); + b[j + 2 * L] = s_mp_ntt_mont_sub(mc, s_mp_ntt_mont_sub(mc, y0, y1), v); + e1 += s1; + if (e1 >= m) { + e1 -= m; + } + e2 += s2; + if (e2 >= m) { + e2 -= m; + } + } +} +#endif + + + +#endif diff --git a/s_mp_ntt_inv5.c b/s_mp_ntt_inv5.c new file mode 100644 index 00000000..944f7ae2 --- /dev/null +++ b/s_mp_ntt_inv5.c @@ -0,0 +1,71 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_INV5_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +void s_mp_ntt_inv5(const ntt_mont *mc, mp_digit *b, int L, int step, const ntt_sroots *R) +{ + const mp_digit *w = R->w; + int m = R->m, j; + int s1 = step, s2 = 2 * step, s3 = 3 * step, s4 = 4 * step; + int e1 = 0, e2 = 0, e3 = 0, e4 = 0; + mp_digit k0 = R->k5i[0], k1 = R->k5i[1], k2 = R->k5i[2], k3 = R->k5i[3]; + mp_digit y0, y1, y2, y3, y4, t1, t2, u1, u2, t1c1, t2c2, t1c2, t2c1, u1s1, u2s12, u1s12, u2s1; + + + for (j = 0; j < L; j++) { + y0 = b[j]; + y1 = s_mp_ntt_mont_mul(mc, b[j + L], w[e1 ? m - e1 : 0]); + y2 = s_mp_ntt_mont_mul(mc, b[j + 2 * L], w[e2 ? m - e2 : 0]); + y3 = s_mp_ntt_mont_mul(mc, b[j + 3 * L], w[e3 ? m - e3 : 0]); + y4 = s_mp_ntt_mont_mul(mc, b[j + 4 * L], w[e4 ? m - e4 : 0]); + t1 = s_mp_ntt_mont_add(mc, y1, y4); + u1 = s_mp_ntt_mont_sub(mc, y1, y4); + t2 = s_mp_ntt_mont_add(mc, y2, y3); + u2 = s_mp_ntt_mont_sub(mc, y2, y3); + t1c1 = s_mp_ntt_mont_mul(mc, t1, k0); + u1s1 = s_mp_ntt_mont_mul(mc, u1, k2); + t2c2 = s_mp_ntt_mont_mul(mc, t2, k1); + u2s12 = s_mp_ntt_mont_mul(mc, u2, k3); + t1c2 = s_mp_ntt_mont_mul(mc, t1, k1); + u1s12 = s_mp_ntt_mont_mul(mc, u1, k3); + t2c1 = s_mp_ntt_mont_mul(mc, t2, k0); + u2s1 = s_mp_ntt_mont_mul(mc, u2, k2); + + b[j] = s_mp_ntt_mont_add(mc, s_mp_ntt_mont_add(mc, y0, t1), t2); + b[j + L] = s_mp_ntt_mont_add(mc, s_mp_ntt_mont_add(mc, y0, s_mp_ntt_mont_add(mc, t1c1, u1s1)), + s_mp_ntt_mont_add(mc, t2c2, u2s12)); + b[j + 2 * L] = s_mp_ntt_mont_add(mc, s_mp_ntt_mont_add(mc, y0, s_mp_ntt_mont_add(mc, t1c2, u1s12)), + s_mp_ntt_mont_sub(mc, t2c1, u2s1)); + b[j + 3 * L] = s_mp_ntt_mont_add(mc, s_mp_ntt_mont_add(mc, y0, s_mp_ntt_mont_sub(mc, t1c2, u1s12)), + s_mp_ntt_mont_add(mc, t2c1, u2s1)); + b[j + 4 * L] = s_mp_ntt_mont_add(mc, s_mp_ntt_mont_add(mc, y0, s_mp_ntt_mont_sub(mc, t1c1, u1s1)), + s_mp_ntt_mont_sub(mc, t2c2, u2s12)); + e1 += s1; + if (e1 >= m) { + e1 -= m; + } + e2 += s2; + if (e2 >= m) { + e2 -= m; + } + e3 += s3; + if (e3 >= m) { + e3 -= m; + } + e4 += s4; + if (e4 >= m) { + e4 -= m; + } + } +} +#endif + + + +#endif diff --git a/s_mp_ntt_invmod_d.c b/s_mp_ntt_invmod_d.c new file mode 100644 index 00000000..930a5dad --- /dev/null +++ b/s_mp_ntt_invmod_d.c @@ -0,0 +1,17 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_INVMOD_D_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + +mp_digit s_mp_ntt_invmod_d(mp_digit p) +{ + mp_digit x = 1; + int i; + for (i = 0; i < 6; i++) { + x = x * (mp_digit)(2 - p * x); + } + return x; +} + +#endif diff --git a/s_mp_ntt_mont_add.c b/s_mp_ntt_mont_add.c new file mode 100644 index 00000000..790b1274 --- /dev/null +++ b/s_mp_ntt_mont_add.c @@ -0,0 +1,19 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_MONT_ADD_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +/* Montgomery sub/add of residues in [0,p). */ +mp_digit s_mp_ntt_mont_add(const ntt_mont *m, mp_digit a, mp_digit b) +{ + mp_digit d = a + b - m->p; + return d + (m->p & ((mp_digit)0 - (d >> NTT_BASE_SIZE_M_ONE))); +} +#endif + + +#endif diff --git a/s_mp_ntt_mont_init.c b/s_mp_ntt_mont_init.c new file mode 100644 index 00000000..6d094b3a --- /dev/null +++ b/s_mp_ntt_mont_init.c @@ -0,0 +1,22 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_MONT_INIT_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +void s_mp_ntt_mont_init(ntt_mont *m, mp_digit p) +{ + mp_digit inv, r; + m->p = p; + inv = s_mp_ntt_invmod_d(p); + m->ninv = (mp_digit)(0 - inv); /* -p^{-1} mod 2^64 */ + r = (mp_digit)(((mp_word)1 << NTT_BASE_SIZE) % p); /* R mod p = 2^64 mod p */ + m->r2 = (mp_digit)((mp_word)r * r % p); /* R^2 mod p */ +} +#endif + + +#endif diff --git a/s_mp_ntt_mont_mul.c b/s_mp_ntt_mont_mul.c new file mode 100644 index 00000000..5cc004d5 --- /dev/null +++ b/s_mp_ntt_mont_mul.c @@ -0,0 +1,18 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_MONT_MUL_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +/* Montgomery multiply: a,b in [0,p) -> a*b*R^{-1} mod p. */ +mp_digit s_mp_ntt_mont_mul(const ntt_mont *m, mp_digit a, mp_digit b) +{ + return s_mp_ntt_mont_redc(m, (mp_word)a * b); +} +#endif + + +#endif diff --git a/s_mp_ntt_mont_pow.c b/s_mp_ntt_mont_pow.c new file mode 100644 index 00000000..403ae2d2 --- /dev/null +++ b/s_mp_ntt_mont_pow.c @@ -0,0 +1,29 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_MONT_POW_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +/* + Montgomery-form modular exponentiation: base^exp mod p, base in Montgomery + form, result in Montgomery form. + */ +mp_digit s_mp_ntt_mont_pow(const ntt_mont *m, mp_digit base, mp_digit exp) +{ + /* 1 in Montgomery form. No define because it is not a constant here */ + mp_digit result = s_mp_ntt_to_mont(m, 1); + mp_digit b = base; + while (exp) { + if (exp & 1) result = s_mp_ntt_mont_mul(m, result, b); + b = s_mp_ntt_mont_mul(m, b, b); + exp >>= 1; + } + return result; +} +#endif + + +#endif diff --git a/s_mp_ntt_mont_redc.c b/s_mp_ntt_mont_redc.c new file mode 100644 index 00000000..1e5b88e7 --- /dev/null +++ b/s_mp_ntt_mont_redc.c @@ -0,0 +1,36 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_MONT_REDC_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +/* We rely on p < 2^63, which the NTT primes (< 2^62) always satisfy: + - REDC's output r is in [0, 2p) and every input a, b, s is in [0, p); + - subtracting p from a value below p wraps to a value whose top bit is + set, because 2^64 - p > 2^63. The borrow is therefore just bit 63 of + the wrapped difference, and "add p back if we borrowed" is one mask and + one add. The same expression would NOT be valid for a modulus close to + 2^64 (e.g. 2^64 - 2^32 + 1): there 2^64 - p is small and the wrapped + difference does not reach bit 63. + */ +/* REDC: T in [0, p*R) -> T * R^{-1} mod p, in [0, p). */ +mp_digit s_mp_ntt_mont_redc(const ntt_mont *m, mp_word T) +{ + mp_digit q, r, d; + mp_word u; + q = (mp_digit)T * m->ninv; /* (T mod R) * (-p^{-1}) mod R */ + u = T + (mp_word)q * m->p; /* u < 2*p*R < 2^127 */ + r = (mp_digit)(u >> NTT_BASE_SIZE); /* r < 2*p */ + d = r - m->p; + d = d + (m->p & ((mp_digit)0u - (d >> NTT_BASE_SIZE_M_ONE))); + return d; +} +#endif + + + +#endif diff --git a/s_mp_ntt_mont_sub.c b/s_mp_ntt_mont_sub.c new file mode 100644 index 00000000..c965c42d --- /dev/null +++ b/s_mp_ntt_mont_sub.c @@ -0,0 +1,18 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_MONT_SUB_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +mp_digit s_mp_ntt_mont_sub(const ntt_mont *m, mp_digit a, mp_digit b) +{ + mp_digit d = a - b; + return d + (m->p & ((mp_digit)0 - (d >> NTT_BASE_SIZE_M_ONE))); +} +#endif + + +#endif diff --git a/s_mp_ntt_small_is_prime.c b/s_mp_ntt_small_is_prime.c new file mode 100644 index 00000000..e3c376f4 --- /dev/null +++ b/s_mp_ntt_small_is_prime.c @@ -0,0 +1,69 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_SMALL_IS_PRIME_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +bool s_mp_ntt_small_is_prime(mp_digit n) +{ + /* More, but all primes < 2^15 */ + /* const mp_digit bases[12] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}; */ + /* All primes < 2^31 */ + const mp_digit bases[7] = {2u, 325u, 9375u, 28178u, 450775u, 9780504u, 1795265022u}; + int i, j, r; + int composite; + mp_digit a, x, d, p; + + if (n < 2u) { + return false; + } + if (n == 2u) { + return true; + } + if ((n & 1) == 0) { + return false; + } + /* Yes, my dear linter, these are magic numbers, get used to it */ + for (p = 3u; (p < 1000u) && (p * p <= n); p+= 2u) { + if (n % p == 0) { + return (n == p); + } + } + + d = n - 1; + r = 0; + while (!(d & 1)) { + d >>= 1; + r++; + } + for (i = 0; i < 7; i++) { + a = bases[i] % n; + if (a == 0) { + continue; + } + x = s_mp_ntt_small_powmod(a, d, n); + if ((x == 1) || (x == n - 1)) { + continue; + } + composite = 1; + for (j = 1; j < r; j++) { + x = s_mp_ntt_small_mulmod(x, x, n); + if (x == (n - 1)) { + composite = 0; + break; + } + } + if (composite) { + return false; + } + } + return true; +} +#endif + + + +#endif diff --git a/s_mp_ntt_small_mulmod.c b/s_mp_ntt_small_mulmod.c new file mode 100644 index 00000000..e93eb915 --- /dev/null +++ b/s_mp_ntt_small_mulmod.c @@ -0,0 +1,19 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_SMALL_MULMOD_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +/* Non-montgomery for bootstrapping */ +mp_digit s_mp_ntt_small_mulmod(mp_digit a, mp_digit b, mp_digit m) +{ + return (mp_digit)((mp_word)a * b % m); +} +#endif + + + +#endif diff --git a/s_mp_ntt_small_powmod.c b/s_mp_ntt_small_powmod.c new file mode 100644 index 00000000..6bc53585 --- /dev/null +++ b/s_mp_ntt_small_powmod.c @@ -0,0 +1,28 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_SMALL_POWMOD_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +/* Non-montgomery for bootstrapping */ +mp_digit s_mp_ntt_small_powmod(mp_digit b, mp_digit e, mp_digit m) +{ + mp_digit r = 1u; + b %= m; + while (e) { + if (e & 1) { + r = s_mp_ntt_small_mulmod(r, b, m); + } + b = s_mp_ntt_small_mulmod(b, b, m); + e >>= 1; + } + return r; +} +#endif + + + +#endif diff --git a/s_mp_ntt_smooth_fwd.c b/s_mp_ntt_smooth_fwd.c new file mode 100644 index 00000000..a96f7fa9 --- /dev/null +++ b/s_mp_ntt_smooth_fwd.c @@ -0,0 +1,42 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_SMOOTH_FWD_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +void s_mp_ntt_smooth_fwd(ntt_mont *mc, mp_digit *a, int m, const ntt_sroots *R) +{ + int k, start, r, ell, L, step; + mp_digit *b; + + for (k = 0; k < R->nstages; k++) { + r = R->radix[k]; + ell = R->ell[k]; + L = ell / r; + step = m / ell; + for (start = 0; start < m; start += ell) { + b = a + start; + switch (r) { + case 5: + s_mp_ntt_fwd5(mc, b, L, step, R); + break; + case 3: + s_mp_ntt_fwd3(mc, b, L, step, R); + break; + default: + s_mp_ntt_fwd2(mc, b, L, step, R); + break; + } + } + } +} +#endif + + + + +#endif diff --git a/s_mp_ntt_smooth_inv.c b/s_mp_ntt_smooth_inv.c new file mode 100644 index 00000000..3ed49b07 --- /dev/null +++ b/s_mp_ntt_smooth_inv.c @@ -0,0 +1,45 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_SMOOTH_INV_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +void s_mp_ntt_smooth_inv(ntt_mont *mc, mp_digit *a, int m, const ntt_sroots *R) +{ + int k, start, i, r, ell, L, step; + mp_digit ninv; + mp_digit *b; + + for (k = R->nstages - 1; k >= 0; k--) { + r = R->radix[k]; + ell = R->ell[k]; + L = ell / r; + step = m / ell; + for (start = 0; start < m; start += ell) { + b = a + start; + switch (r) { + case 5: + s_mp_ntt_inv5(mc, b, L, step, R); + break; + case 3: + s_mp_ntt_inv3(mc, b, L, step, R); + break; + default: + s_mp_ntt_inv2(mc, b, L, step, R); + break; + } + } + } + ninv = s_mp_ntt_mont_pow(mc, s_mp_ntt_to_mont(mc, (mp_digit)m), mc->p - 2); /* 1/m, Montgomery */ + for (i = 0; i < m; i++) { + a[i] = s_mp_ntt_mont_mul(mc, a[i], ninv); + } +} +#endif + + + +#endif diff --git a/s_mp_ntt_smooth_length.c b/s_mp_ntt_smooth_length.c new file mode 100644 index 00000000..a082eb67 --- /dev/null +++ b/s_mp_ntt_smooth_length.c @@ -0,0 +1,19 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_SMOOTH_LENGTH_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +int s_mp_ntt_smooth_length(int need) +{ + int length; + length = s_mp_ntt_smooth_length_weighted(need, NTT_SMOOTH_W2, NTT_SMOOTH_W3, NTT_SMOOTH_W5); + return length; +} +#endif + + +#endif diff --git a/s_mp_ntt_smooth_length_weighted.c b/s_mp_ntt_smooth_length_weighted.c new file mode 100644 index 00000000..5168ec76 --- /dev/null +++ b/s_mp_ntt_smooth_length_weighted.c @@ -0,0 +1,69 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_SMOOTH_LENGTH_WEIGHTED_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +/* Gives the different radix sizes a different weight according to butterfly runtimes */ +int s_mp_ntt_smooth_length_weighted(int n_min, int w2, int w3, int w5) +{ + int64_t next_pow2 = 1; + int exp_pow2 = 0; + int a, b, c, p5, p5_p3; + int64_t n, current_n, size_limit, current_cost, cost; + + while (next_pow2 < n_min) { + next_pow2 *= 2; + exp_pow2++; + } + n = next_pow2; + a = exp_pow2; + b = 0; + c = 0; + cost = next_pow2 * (exp_pow2 * w2); + + size_limit = next_pow2 * 2; + + p5 = 1; + for (c = 0; p5 <= size_limit; c++, p5 *= 5) { + p5_p3 = p5; + for (b = 0; p5_p3 <= size_limit; b++, p5_p3 *= 3) { + /* Find the smallest 'a' to cross n_min */ + a = 0; + current_n = p5_p3; + while (current_n < n_min) { + current_n *= 2; + a++; + } + /* Just brute force */ + while (current_n <= size_limit) { + current_cost = current_n * (a * w2 + b * w3 + c * w5); + if (current_cost < cost) { + cost = current_cost; + n = current_n; + } + /* Pick the smaller in case of equal cost */ + else if (current_cost == cost && current_n < n) { + n = current_n; + } + current_n *= 2; + a++; + } + /* Check for overflow */ + if (p5_p3 > (INT_MAX / 3)) { + break; + } + } + if (p5 > (INT_MAX / 5)) { + break; + } + } + return (int)n; +} +#endif + + +#endif diff --git a/s_mp_ntt_smooth_plan.c b/s_mp_ntt_smooth_plan.c new file mode 100644 index 00000000..9152cdf8 --- /dev/null +++ b/s_mp_ntt_smooth_plan.c @@ -0,0 +1,50 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_SMOOTH_PLAN_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +/* + Radices of m in the order the stages apply them: all the 5s, then the 3s, + then the 2s. (Any order is correct; this one puts the long-radix work in the + early stages, where the block is large and the twiddles are few or trivial, + and leaves the lot of radix-2 stages for last, where a block is two adjacent + elements.) Returns the number of stages, or -1 if m is not 5-smooth or needs + more room than the caller has. + */ +/* TODO: output to given variable and only return errors */ +int s_mp_ntt_smooth_plan(int m, int *radix, int max) +{ + int n = 0, t = m; + + if (m < 2) return -1; + while ((t % 5) == 0) { + if (n >= max) { + return -1; + } + radix[n++] = 5; + t /= 5; + } + while ((t % 3) == 0) { + if (n >= max) { + return -1; + } + radix[n++] = 3; + t /= 3; + } + while ((t % 2) == 0) { + if (n >= max) { + return -1; + } + radix[n++] = 2; + t /= 2; + } + return (t == 1) ? n : -1; +} +#endif + + +#endif diff --git a/s_mp_ntt_sort_primes.c b/s_mp_ntt_sort_primes.c new file mode 100644 index 00000000..59798925 --- /dev/null +++ b/s_mp_ntt_sort_primes.c @@ -0,0 +1,32 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_SORT_PRIMES_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +/* TODO: There are at most four primes, no extra sort function necessary, do it inline */ +/* + Insertion sort by prime value. Both drivers select their primes largest + first and hand them to Garner ascending. + */ +void s_mp_ntt_sort_primes(ntt_prime_ent *ents, int n) +{ + int a, b; + + for (a = 1; a < n; a++) { + ntt_prime_ent key = ents[a]; + b = a - 1; + while (b >= 0 && ents[b].p > key.p) { + ents[b + 1] = ents[b]; + b--; + } + ents[b + 1] = key; + } +} +#endif + + +#endif diff --git a/s_mp_ntt_sroots_free.c b/s_mp_ntt_sroots_free.c new file mode 100644 index 00000000..1a184aa8 --- /dev/null +++ b/s_mp_ntt_sroots_free.c @@ -0,0 +1,22 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_SROOTS_FREE_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +void s_mp_ntt_sroots_free(ntt_sroots *R) +{ + if (!R) { + return; + } + MP_FREE(R->w, (size_t)R->m * sizeof(mp_digit)); + R->w = NULL; + R->nstages = 0; +} +#endif + + +#endif diff --git a/s_mp_ntt_sroots_init.c b/s_mp_ntt_sroots_init.c new file mode 100644 index 00000000..f4276711 --- /dev/null +++ b/s_mp_ntt_sroots_init.c @@ -0,0 +1,99 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_SROOTS_INIT_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +/* + Radix 5 needs more care than radix 2 or 3: the field has no imaginary part to + split the transform into. With t = zeta^{m/5} and the pairs (y1,y4), (y2,y3), + + t1 = y1+y4, u1 = y1-y4, t2 = y2+y3, u2 = y2-y3, + c1 = t + t^{-1}, c2 = t^2 + t^{-2}, + s1 = t - t^{-1}, s12 = t^2 - t^{-2} = s1*c1, + + the five outputs are + + Z_0 = y0 + t1 + t2 + Z_1 = y0 + (t1*c1 + u1*s1 + t2*c2 + u2*s12)/2 + Z_4 = y0 + (t1*c1 - u1*s1 + t2*c2 - u2*s12)/2 + Z_2 = y0 + (t1*c2 + u1*s12 + t2*c1 - u2*s1)/2 + Z_3 = y0 + (t1*c2 - u1*s12 + t2*c1 + u2*s1)/2 + + The 1/2 comes from writing each member of a pair through its sum and + difference; folding it into the constants leaves the butterfly with eight + multiplications per five elements - four distinct products per pair, each one + reused by a sum and a difference. The inverse uses the same constants with + s1 and s12 negated - those are the two that change sign under t -> t^{-1}; + c1 and c2 are symmetric. + */ +mp_err s_mp_ntt_sroots_init(ntt_sroots *R, ntt_mont *mc, int m, mp_digit zeta) +{ + int i, ell, n; + + R->m = m; + R->nstages = 0; + R->w = NULL; + R->has3 = ((m % 3) == 0); + R->has5 = ((m % 5) == 0); + R->c3 = 0; + R->c3i = 0; + + for (i = 0; i < 4; i++) { + R->k5[i] = 0; + R->k5i[i] = 0; + } + + if (m < 2) { + return MP_VAL; + } + n = s_mp_ntt_smooth_plan(m, R->radix, NTT_MR_MAXSTAGE); + if (n <= 0) { + return MP_VAL; + } + R->nstages = n; + for (i = 0, ell = m; i < n; i++) { + R->ell[i] = ell; + ell /= R->radix[i]; + } + + R->w = (mp_digit *)MP_MALLOC((size_t)m * sizeof(mp_digit)); + if (R->w == NULL) { + return MP_MEM; + } + R->w[0] = s_mp_ntt_to_mont(mc, 1); + for (i = 1; i < m; i++) { + R->w[i] = s_mp_ntt_mont_mul(mc, R->w[i - 1], zeta); + } + if (R->has3) { + R->c3 = R->w[m / 3]; /* primitive cube root */ + R->c3i = R->w[m - m / 3]; /* its inverse */ + } + if (R->has5) { + mp_digit t = R->w[m / 5]; /* primitive fifth root */ + mp_digit ti = R->w[m - m / 5]; /* its inverse */ + mp_digit h = s_mp_ntt_to_mont(mc, (mc->p + 1) >> 1); /* 2^{-1} mod p */ + mp_digit c1 = s_mp_ntt_mont_add(mc, t, ti); + /* TODO: is that readable or should distinct temporary variables be used? */ + mp_digit c2 = s_mp_ntt_mont_add(mc, s_mp_ntt_mont_mul(mc, t, t), s_mp_ntt_mont_mul(mc, ti, ti)); + mp_digit s1 = s_mp_ntt_mont_sub(mc, t, ti); + mp_digit s12 = s_mp_ntt_mont_mul(mc, s1, c1); + + R->k5[0] = s_mp_ntt_mont_mul(mc, c1, h); + R->k5[1] = s_mp_ntt_mont_mul(mc, c2, h); + R->k5[2] = s_mp_ntt_mont_mul(mc, s1, h); + R->k5[3] = s_mp_ntt_mont_mul(mc, s12, h); + R->k5i[0] = R->k5[0]; + R->k5i[1] = R->k5[1]; + R->k5i[2] = s_mp_ntt_mont_sub(mc, 0, R->k5[2]); /* -h*s1 */ + R->k5i[3] = s_mp_ntt_mont_sub(mc, 0, R->k5[3]); /* -h*s12 */ + } + return MP_OKAY; +} +#endif + + +#endif diff --git a/s_mp_ntt_to_mont.c b/s_mp_ntt_to_mont.c new file mode 100644 index 00000000..3eae9c50 --- /dev/null +++ b/s_mp_ntt_to_mont.c @@ -0,0 +1,15 @@ +#include "tommath_private.h" +#ifdef S_MP_NTT_TO_MONT_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +/* Convert x in [0,p) into Montgomery form (x * R mod p), and back. */ +mp_digit s_mp_ntt_to_mont(const ntt_mont *m, mp_digit x) +{ + return s_mp_ntt_mont_mul(m, x, m->r2); +} +#endif +#endif + diff --git a/s_mp_sqr_ntt.c b/s_mp_sqr_ntt.c new file mode 100644 index 00000000..01a81e5f --- /dev/null +++ b/s_mp_sqr_ntt.c @@ -0,0 +1,129 @@ +#include "tommath_private.h" +#ifdef S_MP_SQR_NTT_C +/* LibTomMath, multiple-precision integer library -- Tom St Denis */ +/* SPDX-License-Identifier: Unlicense */ + + + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) +/* The main function to square one mp_int into another mp_int */ +mp_err s_mp_sqr_ntt(const mp_int *a, mp_int *c) +{ + mp_err err = MP_OKAY; + /* TODO: nprimes is three, defined in NTT_MAXP, rewrite accordingly */ + int ncoef, m, minl, bound_bits, nprimes, pi, i; + ntt_prime_ent *ents = NULL; + ntt_crt_ctx *cx = NULL; + mp_digit *primes = NULL, **residues = NULL, *ra = NULL; + mp_int P, pmp; + + + if (mp_iszero(a)) { + return mp_sqr(a,c); + } + if (a->used < 8) { + return mp_sqr(a, c); + } + + ncoef = a->used + a->used - 1; /* coefficients that carry the product */ + if ((ncoef > (1 << NTT_QUARTER_COEFF)) || (ncoef < 0)) { + return MP_VAL; + } + m = s_mp_ntt_smooth_length(ncoef); /* smallest 5-smooth transform length */ + if (m < ncoef) { + return MP_VAL; + } + + /* Convolution bound: each output coefficient < min(la,lb) * 2^120. */ + minl = a->used; + bound_bits = (int)NTT_LARGE_SIZE + s_mp_ntt_ilog2_ceil(minl) + 1; + if (((bound_bits + NTT_LIMB_SIZE_M_ONE) / NTT_LIMB_SIZE) > NTT_MAXW) { + return MP_VAL; + } + + ents = (ntt_prime_ent *)MP_MALLOC((size_t)NTT_MAXP * sizeof(ntt_prime_ent)); + if (ents == NULL) { + return MP_MEM; + } + + if ((err = mp_init_multi(&P, &pmp, NULL)) != MP_OKAY) { + MP_FREE(ents, (size_t)NTT_MAXP * sizeof(ntt_prime_ent)); + return err; + } + /* TODO */ + mp_set(&P, 1); + + /* Find primes until their product > 2^bound_bits */ + nprimes = 0; + while (mp_count_bits(&P) <= bound_bits) { + if (nprimes >= NTT_MAXP) { + err = MP_VAL; + goto LTM_ERR; + } + if (!s_mp_ntt_find_smooth_prime(&ents[nprimes], m, nprimes, 0)) { + err = MP_VAL; + goto LTM_ERR; + } + mp_set(&pmp, ents[nprimes].p); + if ((err = mp_mul(&P, &pmp, &P)) != MP_OKAY) goto LTM_ERR; + nprimes++; + } + + s_mp_ntt_sort_primes(ents, nprimes); + + ra = (mp_digit *)MP_MALLOC((size_t)m * sizeof(mp_digit)); + if (ra == NULL) { + err = MP_MEM; + goto LTM_ERR; + } + residues = (mp_digit **)MP_CALLOC((size_t)nprimes, sizeof(mp_digit *)); + if (residues == NULL) { + err = MP_MEM; + goto LTM_ERR; + } + primes = (mp_digit *)MP_MALLOC((size_t)nprimes * sizeof(mp_digit)); + if (primes == NULL) { + err = MP_MEM; + goto LTM_ERR; + } + + + /* convolve over each prime and keep the residues */ + for (pi = 0; pi < nprimes; pi++) { + residues[pi] = (mp_digit *)MP_MALLOC((size_t)ncoef * sizeof(mp_digit)); + if (residues[pi] == NULL) { + err = MP_MEM; + goto LTM_ERR; + } + if ((err = s_mp_ntt_autoconvolve_single(a, &ents[pi], ra, residues[pi], ncoef)) != MP_OKAY) goto LTM_ERR; + primes[pi] = ents[pi].p; + } + + /* CRT: residues -> binary, then one linear carry sweep */ + cx = (ntt_crt_ctx *)MP_MALLOC(sizeof(ntt_crt_ctx)); + if (cx == NULL) { + err = MP_MEM; + goto LTM_ERR; + } + s_mp_ntt_crt_init(cx, primes, nprimes, bound_bits); + err = s_mp_ntt_crt_assemble(cx, residues, nprimes, ncoef, c); + +LTM_ERR: + mp_clear_multi(&P, &pmp, NULL); + if (residues) { + for (i = 0; i < nprimes; i++) { + MP_FREE(residues[i], (size_t)ncoef * sizeof(mp_digit)); + } + MP_FREE(residues, (size_t)nprimes * sizeof(mp_digit *)); + } + MP_FREE(primes, (size_t)nprimes * sizeof(mp_digit)); + MP_FREE(ra, (size_t)m * sizeof(mp_digit)); + MP_FREE(cx, sizeof(ntt_crt_ctx)); + MP_FREE(ents, (size_t)NTT_MAXP * sizeof(ntt_prime_ent)); + return err; +} +#endif + + +#endif diff --git a/sources.cmake b/sources.cmake index ba3012ef..3650b7d7 100644 --- a/sources.cmake +++ b/sources.cmake @@ -146,7 +146,41 @@ s_mp_mul_comba.c s_mp_mul_high.c s_mp_mul_high_comba.c s_mp_mul_karatsuba.c +s_mp_mul_ntt.c s_mp_mul_toom.c +s_mp_ntt_autoconvolve_single.c +s_mp_ntt_convolve_single.c +s_mp_ntt_crt_assemble.c +s_mp_ntt_crt_coeff.c +s_mp_ntt_crt_init.c +s_mp_ntt_find_smooth_prime.c +s_mp_ntt_from_mont.c +s_mp_ntt_fwd2.c +s_mp_ntt_fwd3.c +s_mp_ntt_fwd5.c +s_mp_ntt_ilog2_ceil.c +s_mp_ntt_inv2.c +s_mp_ntt_inv3.c +s_mp_ntt_inv5.c +s_mp_ntt_invmod_d.c +s_mp_ntt_mont_add.c +s_mp_ntt_mont_init.c +s_mp_ntt_mont_mul.c +s_mp_ntt_mont_pow.c +s_mp_ntt_mont_redc.c +s_mp_ntt_mont_sub.c +s_mp_ntt_small_is_prime.c +s_mp_ntt_small_mulmod.c +s_mp_ntt_small_powmod.c +s_mp_ntt_smooth_fwd.c +s_mp_ntt_smooth_inv.c +s_mp_ntt_smooth_length.c +s_mp_ntt_smooth_length_weighted.c +s_mp_ntt_smooth_plan.c +s_mp_ntt_sort_primes.c +s_mp_ntt_sroots_free.c +s_mp_ntt_sroots_init.c +s_mp_ntt_to_mont.c s_mp_prime_is_divisible.c s_mp_prime_tab.c s_mp_radix_map.c @@ -156,6 +190,7 @@ s_mp_rand_source.c s_mp_sqr.c s_mp_sqr_comba.c s_mp_sqr_karatsuba.c +s_mp_sqr_ntt.c s_mp_sqr_toom.c s_mp_sub.c s_mp_warray.c diff --git a/tommath.def b/tommath.def index ed5aa8b0..ad834099 100644 --- a/tommath.def +++ b/tommath.def @@ -132,3 +132,5 @@ EXPORTS MP_SQR_KARATSUBA_CUTOFF MP_MUL_TOOM_CUTOFF MP_SQR_TOOM_CUTOFF + MP_MUL_NTT_CUTOFF + MP_SQR_NTT_CUTOFF diff --git a/tommath.h b/tommath.h index 1820d243..f040bc1d 100644 --- a/tommath.h +++ b/tommath.h @@ -121,7 +121,9 @@ extern int MP_MUL_KARATSUBA_CUTOFF, MP_SQR_KARATSUBA_CUTOFF, MP_MUL_TOOM_CUTOFF, -MP_SQR_TOOM_CUTOFF; +MP_SQR_TOOM_CUTOFF, +MP_MUL_NTT_CUTOFF, +MP_SQR_NTT_CUTOFF; #endif /* define this to use lower memory usage routines (exptmods mostly) */ diff --git a/tommath_class.h b/tommath_class.h index 40ddfd4b..300559f6 100644 --- a/tommath_class.h +++ b/tommath_class.h @@ -155,7 +155,41 @@ # define S_MP_MUL_HIGH_C # define S_MP_MUL_HIGH_COMBA_C # define S_MP_MUL_KARATSUBA_C +# define S_MP_MUL_NTT_C # define S_MP_MUL_TOOM_C +# define S_MP_NTT_AUTOCONVOLVE_SINGLE_C +# define S_MP_NTT_CONVOLVE_SINGLE_C +# define S_MP_NTT_CRT_ASSEMBLE_C +# define S_MP_NTT_CRT_COEFF_C +# define S_MP_NTT_CRT_INIT_C +# define S_MP_NTT_FIND_SMOOTH_PRIME_C +# define S_MP_NTT_FROM_MONT_C +# define S_MP_NTT_FWD2_C +# define S_MP_NTT_FWD3_C +# define S_MP_NTT_FWD5_C +# define S_MP_NTT_ILOG2_CEIL_C +# define S_MP_NTT_INV2_C +# define S_MP_NTT_INV3_C +# define S_MP_NTT_INV5_C +# define S_MP_NTT_INVMOD_D_C +# define S_MP_NTT_MONT_ADD_C +# define S_MP_NTT_MONT_INIT_C +# define S_MP_NTT_MONT_MUL_C +# define S_MP_NTT_MONT_POW_C +# define S_MP_NTT_MONT_REDC_C +# define S_MP_NTT_MONT_SUB_C +# define S_MP_NTT_SMALL_IS_PRIME_C +# define S_MP_NTT_SMALL_MULMOD_C +# define S_MP_NTT_SMALL_POWMOD_C +# define S_MP_NTT_SMOOTH_FWD_C +# define S_MP_NTT_SMOOTH_INV_C +# define S_MP_NTT_SMOOTH_LENGTH_C +# define S_MP_NTT_SMOOTH_LENGTH_WEIGHTED_C +# define S_MP_NTT_SMOOTH_PLAN_C +# define S_MP_NTT_SORT_PRIMES_C +# define S_MP_NTT_SROOTS_FREE_C +# define S_MP_NTT_SROOTS_INIT_C +# define S_MP_NTT_TO_MONT_C # define S_MP_PRIME_IS_DIVISIBLE_C # define S_MP_PRIME_TAB_C # define S_MP_RADIX_MAP_C @@ -165,6 +199,7 @@ # define S_MP_SQR_C # define S_MP_SQR_COMBA_C # define S_MP_SQR_KARATSUBA_C +# define S_MP_SQR_NTT_C # define S_MP_SQR_TOOM_C # define S_MP_SUB_C # define S_MP_WARRAY_C @@ -557,10 +592,12 @@ # define S_MP_MUL_C # define S_MP_MUL_COMBA_C # define S_MP_MUL_KARATSUBA_C +# define S_MP_MUL_NTT_C # define S_MP_MUL_TOOM_C # define S_MP_SQR_C # define S_MP_SQR_COMBA_C # define S_MP_SQR_KARATSUBA_C +# define S_MP_SQR_NTT_C # define S_MP_SQR_TOOM_C #endif @@ -1209,6 +1246,21 @@ # define S_MP_SUB_C #endif +#if defined(S_MP_MUL_NTT_C) +# define MP_CLEAR_MULTI_C +# define MP_COUNT_BITS_C +# define MP_INIT_MULTI_C +# define MP_MUL_C +# define MP_SET_C +# define S_MP_NTT_CONVOLVE_SINGLE_C +# define S_MP_NTT_CRT_ASSEMBLE_C +# define S_MP_NTT_CRT_INIT_C +# define S_MP_NTT_FIND_SMOOTH_PRIME_C +# define S_MP_NTT_ILOG2_CEIL_C +# define S_MP_NTT_SMOOTH_LENGTH_C +# define S_MP_NTT_SORT_PRIMES_C +#endif + #if defined(S_MP_MUL_TOOM_C) # define MP_ADD_C # define MP_CLAMP_C @@ -1225,6 +1277,172 @@ # define S_MP_DIV_3_C #endif +#if defined(S_MP_NTT_AUTOCONVOLVE_SINGLE_C) +# define S_MP_NTT_FROM_MONT_C +# define S_MP_NTT_MONT_INIT_C +# define S_MP_NTT_MONT_MUL_C +# define S_MP_NTT_SMOOTH_FWD_C +# define S_MP_NTT_SMOOTH_INV_C +# define S_MP_NTT_SROOTS_FREE_C +# define S_MP_NTT_SROOTS_INIT_C +# define S_MP_NTT_TO_MONT_C +#endif + +#if defined(S_MP_NTT_CONVOLVE_SINGLE_C) +# define S_MP_NTT_FROM_MONT_C +# define S_MP_NTT_MONT_INIT_C +# define S_MP_NTT_MONT_MUL_C +# define S_MP_NTT_SMOOTH_FWD_C +# define S_MP_NTT_SMOOTH_INV_C +# define S_MP_NTT_SROOTS_FREE_C +# define S_MP_NTT_SROOTS_INIT_C +# define S_MP_NTT_TO_MONT_C +#endif + +#if defined(S_MP_NTT_CRT_ASSEMBLE_C) +# define MP_CLAMP_C +# define MP_GROW_C +# define S_MP_NTT_CRT_COEFF_C +#endif + +#if defined(S_MP_NTT_CRT_COEFF_C) +# define S_MP_NTT_MONT_MUL_C +#endif + +#if defined(S_MP_NTT_CRT_INIT_C) +# define S_MP_NTT_MONT_INIT_C +# define S_MP_NTT_SMALL_POWMOD_C +# define S_MP_NTT_TO_MONT_C +#endif + +#if defined(S_MP_NTT_FIND_SMOOTH_PRIME_C) +# define S_MP_NTT_SMALL_IS_PRIME_C +# define S_MP_NTT_SMALL_POWMOD_C +#endif + +#if defined(S_MP_NTT_FROM_MONT_C) +# define S_MP_NTT_MONT_REDC_C +#endif + +#if defined(S_MP_NTT_FWD2_C) +# define S_MP_NTT_MONT_ADD_C +# define S_MP_NTT_MONT_MUL_C +# define S_MP_NTT_MONT_SUB_C +#endif + +#if defined(S_MP_NTT_FWD3_C) +# define S_MP_NTT_MONT_ADD_C +# define S_MP_NTT_MONT_MUL_C +# define S_MP_NTT_MONT_SUB_C +#endif + +#if defined(S_MP_NTT_FWD5_C) +# define S_MP_NTT_MONT_ADD_C +# define S_MP_NTT_MONT_MUL_C +# define S_MP_NTT_MONT_SUB_C +#endif + +#if defined(S_MP_NTT_ILOG2_CEIL_C) +#endif + +#if defined(S_MP_NTT_INV2_C) +# define S_MP_NTT_MONT_ADD_C +# define S_MP_NTT_MONT_MUL_C +# define S_MP_NTT_MONT_SUB_C +#endif + +#if defined(S_MP_NTT_INV3_C) +# define S_MP_NTT_MONT_ADD_C +# define S_MP_NTT_MONT_MUL_C +# define S_MP_NTT_MONT_SUB_C +#endif + +#if defined(S_MP_NTT_INV5_C) +# define S_MP_NTT_MONT_ADD_C +# define S_MP_NTT_MONT_MUL_C +# define S_MP_NTT_MONT_SUB_C +#endif + +#if defined(S_MP_NTT_INVMOD_D_C) +#endif + +#if defined(S_MP_NTT_MONT_ADD_C) +#endif + +#if defined(S_MP_NTT_MONT_INIT_C) +# define S_MP_NTT_INVMOD_D_C +#endif + +#if defined(S_MP_NTT_MONT_MUL_C) +# define S_MP_NTT_MONT_REDC_C +#endif + +#if defined(S_MP_NTT_MONT_POW_C) +# define S_MP_NTT_MONT_MUL_C +# define S_MP_NTT_TO_MONT_C +#endif + +#if defined(S_MP_NTT_MONT_REDC_C) +#endif + +#if defined(S_MP_NTT_MONT_SUB_C) +#endif + +#if defined(S_MP_NTT_SMALL_IS_PRIME_C) +# define S_MP_NTT_SMALL_MULMOD_C +# define S_MP_NTT_SMALL_POWMOD_C +#endif + +#if defined(S_MP_NTT_SMALL_MULMOD_C) +#endif + +#if defined(S_MP_NTT_SMALL_POWMOD_C) +# define S_MP_NTT_SMALL_MULMOD_C +#endif + +#if defined(S_MP_NTT_SMOOTH_FWD_C) +# define S_MP_NTT_FWD2_C +# define S_MP_NTT_FWD3_C +# define S_MP_NTT_FWD5_C +#endif + +#if defined(S_MP_NTT_SMOOTH_INV_C) +# define S_MP_NTT_INV2_C +# define S_MP_NTT_INV3_C +# define S_MP_NTT_INV5_C +# define S_MP_NTT_MONT_MUL_C +# define S_MP_NTT_MONT_POW_C +# define S_MP_NTT_TO_MONT_C +#endif + +#if defined(S_MP_NTT_SMOOTH_LENGTH_C) +# define S_MP_NTT_SMOOTH_LENGTH_WEIGHTED_C +#endif + +#if defined(S_MP_NTT_SMOOTH_LENGTH_WEIGHTED_C) +#endif + +#if defined(S_MP_NTT_SMOOTH_PLAN_C) +#endif + +#if defined(S_MP_NTT_SORT_PRIMES_C) +#endif + +#if defined(S_MP_NTT_SROOTS_FREE_C) +#endif + +#if defined(S_MP_NTT_SROOTS_INIT_C) +# define S_MP_NTT_MONT_ADD_C +# define S_MP_NTT_MONT_MUL_C +# define S_MP_NTT_MONT_SUB_C +# define S_MP_NTT_SMOOTH_PLAN_C +# define S_MP_NTT_TO_MONT_C +#endif + +#if defined(S_MP_NTT_TO_MONT_C) +# define S_MP_NTT_MONT_MUL_C +#endif + #if defined(S_MP_PRIME_IS_DIVISIBLE_C) # define MP_DIV_D_C #endif @@ -1280,6 +1498,21 @@ # define S_MP_SUB_C #endif +#if defined(S_MP_SQR_NTT_C) +# define MP_CLEAR_MULTI_C +# define MP_COUNT_BITS_C +# define MP_INIT_MULTI_C +# define MP_MUL_C +# define MP_SET_C +# define S_MP_NTT_AUTOCONVOLVE_SINGLE_C +# define S_MP_NTT_CRT_ASSEMBLE_C +# define S_MP_NTT_CRT_INIT_C +# define S_MP_NTT_FIND_SMOOTH_PRIME_C +# define S_MP_NTT_ILOG2_CEIL_C +# define S_MP_NTT_SMOOTH_LENGTH_C +# define S_MP_NTT_SORT_PRIMES_C +#endif + #if defined(S_MP_SQR_TOOM_C) # define MP_ADD_C # define MP_CLAMP_C diff --git a/tommath_cutoffs.h b/tommath_cutoffs.h index fb841601..5459d26b 100644 --- a/tommath_cutoffs.h +++ b/tommath_cutoffs.h @@ -6,8 +6,9 @@ be aware that it may take a long time. It took 2:30 minutes on the aforementioned machine for example. */ - -#define MP_DEFAULT_MUL_KARATSUBA_CUTOFF 80 -#define MP_DEFAULT_SQR_KARATSUBA_CUTOFF 120 -#define MP_DEFAULT_MUL_TOOM_CUTOFF 350 -#define MP_DEFAULT_SQR_TOOM_CUTOFF 400 +#define MP_DEFAULT_MUL_KARATSUBA_CUTOFF 138 +#define MP_DEFAULT_SQR_KARATSUBA_CUTOFF 199 +#define MP_DEFAULT_MUL_TOOM_CUTOFF 248 +#define MP_DEFAULT_SQR_TOOM_CUTOFF 864 +#define MP_DEFAULT_MUL_NTT_CUTOFF 2900 +#define MP_DEFAULT_SQR_NTT_CUTOFF 2900 diff --git a/tommath_private.h b/tommath_private.h index be620dbc..a4036980 100644 --- a/tommath_private.h +++ b/tommath_private.h @@ -86,6 +86,8 @@ do { \ # define MP_SQR_KARATSUBA_CUTOFF MP_DEFAULT_SQR_KARATSUBA_CUTOFF # define MP_MUL_TOOM_CUTOFF MP_DEFAULT_MUL_TOOM_CUTOFF # define MP_SQR_TOOM_CUTOFF MP_DEFAULT_SQR_TOOM_CUTOFF +# define MP_MUL_NTT_CUTOFF MP_DEFAULT_MUL_NTT_CUTOFF +# define MP_SQR_NTT_CUTOFF MP_DEFAULT_SQR_NTT_CUTOFF #endif /* define heap macros */ @@ -217,6 +219,178 @@ MP_PRIVATE mp_err s_mp_mul_high(const mp_int *a, const mp_int *b, mp_int *c, int MP_PRIVATE mp_err s_mp_mul_high_comba(const mp_int *a, const mp_int *b, mp_int *c, int digs) MP_WUR; MP_PRIVATE mp_err s_mp_mul_karatsuba(const mp_int *a, const mp_int *b, mp_int *c) MP_WUR; MP_PRIVATE mp_err s_mp_mul_toom(const mp_int *a, const mp_int *b, mp_int *c) MP_WUR; + +/* NTT stuff, it's a lot */ + +/* MP_16BIT is just too small, not worth it */ +#if (MP_DIGIT_BIT > 16) + +#include + +/* Magic numbers */ +#define NTT_BASE_SIZE (sizeof(mp_digit) * CHAR_BIT) +#define NTT_BASE_SIZE_M_ONE (NTT_BASE_SIZE - 1) +#define NTT_LARGE_SIZE (2 * NTT_BASE_SIZE) +#define NTT_LIMB_SIZE MP_DIGIT_BIT +#define NTT_LIMB_SIZE_M_ONE (NTT_LIMB_SIZE - 1) + +#if (NTT_LIMB_SIZE == 60) +#define NTT_PRIME_MAX_SIZE (NTT_BASE_SIZE - 2) +#else +#define NTT_PRIME_MAX_SIZE (NTT_BASE_SIZE - 2) +#endif + +/* Indirectly influences the number of primes */ +#define NTT_BOUND_BIT_MULTIPLICATOR 2 + + +/* One coeff holds 2 limbs at most (120 bits for MP_64BIT) */ +#if (NTT_LIMB_SIZE == 60) +#define NTT_QUARTER_COEFF (NTT_LIMB_SIZE/2) +#else +#define NTT_QUARTER_COEFF (NTT_LIMB_SIZE) +#endif + +/* TODO: make tunable, although the tuning runs would take a a lot of time */ +/* + Weights for radix-n when computing optimal coefficient-array length. + Assuming radix 2 is the fastest, radix 5 the slowest. + Profiling would be a better idea than the general values below. + */ +#define NTT_SMOOTH_W2 1 +#define NTT_SMOOTH_W3 2 +#define NTT_SMOOTH_W5 6 + +#define NTT_MR_MAXSTAGE 48 /* stages of m: log2(m) plus the 3s and 5s */ + +/* Three primes for 60 bits, four for 28 bit */ +#define NTT_MAXP 6 /* max CRT primes */ + +#if (NTT_LIMB_SIZE == 60) +#define NTT_MAXW 4 /* max limbs per coefficient */ +#else +#define NTT_MAXW 6 +#endif + +/* TODO: Basically MP_MASK. Use directly? */ +#define NTT_MASK (((mp_digit)1 << NTT_LIMB_SIZE) - 1) + +/* Montgomery context for a prime p < 2^62, R = 2^64. */ +typedef struct { + mp_digit p; /* The prime p */ + mp_digit ninv; /* (-p^{-1}) mod R */ + mp_digit r2; /* R^2 mod p */ +} ntt_mont; + +typedef struct { + mp_digit p; /* The prime p */ + mp_digit root; /* primitive 2^s-th (or m-th) root mod p, ordinary residue */ + int s; /* root order is a power of two */ + int m; /* 5-smooth transform length of the root order */ +} ntt_prime_ent; + +/* Root table and butterfly constants of a single smooth transform. */ +typedef struct { + int m; + int nstages; + int radix[NTT_MR_MAXSTAGE]; + int ell[NTT_MR_MAXSTAGE]; /* block length at each stage */ + mp_digit *w; /* w[s] = zeta^s, Montgomery form, s < m */ + mp_digit c3, c3i; /* primitive cube root, and its inverse */ + mp_digit k5[4], k5i[4]; /* radix-5 constants, forward / inverse */ + bool has3; + bool has5; +} ntt_sroots; + +typedef struct { + int t; /* number of primes */ + int width; /* 60-bit limbs per coefficient */ + mp_digit p[NTT_MAXP]; /* the primes, ascending */ + ntt_mont mc[NTT_MAXP]; /* Montgomery context of each prime */ + mp_digit minv[NTT_MAXP][NTT_MAXP]; /* p_i^{-1} * 2^64 mod p_j, i < j */ + mp_digit x[NTT_MAXP]; /* scratch */ + mp_digit d[NTT_MAXP]; /* mixed-radix digits */ +} ntt_crt_ctx; + +/* The main function, the rest are helpers */ +MP_PRIVATE mp_err s_mp_mul_ntt(const mp_int *a, const mp_int *b, mp_int *c); +MP_PRIVATE mp_err s_mp_sqr_ntt(const mp_int *a, mp_int *c); +/* ceil(log_2(x)) */ +MP_PRIVATE int s_mp_ntt_ilog2_ceil(int x); +/* Convolution of a single radix */ +MP_PRIVATE mp_err s_mp_ntt_convolve_single(const mp_int *a, const mp_int *b, const ntt_prime_ent *ent, mp_digit *ra, + mp_digit *rb, mp_digit *out, int ncoef); +MP_PRIVATE mp_err s_mp_ntt_autoconvolve_single(const mp_int *a, const ntt_prime_ent *ent, mp_digit *ra, + mp_digit *out, int ncoef); + +/* Forward butterflies */ +MP_PRIVATE void s_mp_ntt_fwd2(const ntt_mont *mc, mp_digit *b, int L, int step, const ntt_sroots *R); +MP_PRIVATE void s_mp_ntt_fwd3(const ntt_mont *mc, mp_digit *b, int L, int step, const ntt_sroots *R); +MP_PRIVATE void s_mp_ntt_fwd5(const ntt_mont *mc, mp_digit *b, int L, int step, const ntt_sroots *R); + +/* Inverse butterflies */ +MP_PRIVATE void s_mp_ntt_inv2(const ntt_mont *mc, mp_digit *b, int L, int step, const ntt_sroots *R); +MP_PRIVATE void s_mp_ntt_inv3(const ntt_mont *mc, mp_digit *b, int L, int step, const ntt_sroots *R); +MP_PRIVATE void s_mp_ntt_inv5(const ntt_mont *mc, mp_digit *b, int L, int step, const ntt_sroots *R); + +/* Chinese remainder to get teh differrent radices together */ +MP_PRIVATE void s_mp_ntt_crt_init(ntt_crt_ctx *cx, const mp_digit *primes, int t, int bound_bits); +MP_PRIVATE void s_mp_ntt_crt_coeff(ntt_crt_ctx *cx, const mp_digit *r, mp_digit *limbs); +MP_PRIVATE mp_err s_mp_ntt_crt_assemble(ntt_crt_ctx *cx, mp_digit **residues, int nprimes,int ncoef, mp_int *c); + +/* Single mp_digit modular inverse. Was here but got out, div_exact* uses it, too */ +/* Initialize ntt_mont structure with prime p */ +MP_PRIVATE void s_mp_ntt_mont_init(ntt_mont *m, mp_digit p); +/* Montgomery reduction */ +MP_PRIVATE mp_digit s_mp_ntt_mont_redc(const ntt_mont *m, mp_word T); +/* Montgomery to plain */ +MP_PRIVATE mp_digit s_mp_ntt_from_mont(const ntt_mont *m, mp_digit a); +/* Plain to Montgomery */ +MP_PRIVATE mp_digit s_mp_ntt_to_mont(const ntt_mont *m, mp_digit x); + +/* Montgomery arithmetic */ +MP_PRIVATE mp_digit s_mp_ntt_mont_add(const ntt_mont *m, mp_digit a, mp_digit b); +MP_PRIVATE mp_digit s_mp_ntt_mont_mul(const ntt_mont *m, mp_digit a, mp_digit b); +MP_PRIVATE mp_digit s_mp_ntt_mont_pow(const ntt_mont *m, mp_digit base, mp_digit exp); +MP_PRIVATE mp_digit s_mp_ntt_mont_sub(const ntt_mont *m, mp_digit a, mp_digit b); +/* Plain arithmetic. Helpers for bootstrapping the primes */ +MP_PRIVATE mp_digit s_mp_ntt_small_mulmod(mp_digit a, mp_digit b, mp_digit m); +MP_PRIVATE mp_digit s_mp_ntt_small_powmod(mp_digit b, mp_digit e, mp_digit m); + +/* Selects forward butterfly according to radix */ +MP_PRIVATE void s_mp_ntt_smooth_fwd(ntt_mont *mc, mp_digit *a, int m, const ntt_sroots *R); +/* Selects forward butterfly according to radix */ +MP_PRIVATE void s_mp_ntt_smooth_inv(ntt_mont *mc, mp_digit *a, int m, const ntt_sroots *R); + +#ifdef USE_NAIVE +/* Choose length of the coefficient array */ +/* All weights 1 */ +MP_PRIVATE int s_mp_ntt_smooth_length_naive(int n_min); +#else +/* Takes different weights per radix*/ +MP_PRIVATE int s_mp_ntt_smooth_length_weighted(int n_min, int w2, int w3, int w5); +#endif +/* Choose length of the coefficient array, given minimum. Returns n >= need */ +MP_PRIVATE int s_mp_ntt_smooth_length(int need); +/* Organizes the number of radix calls */ +MP_PRIVATE int s_mp_ntt_smooth_plan(int m, int *radix, int max); + +/* Memory for root array */ +MP_PRIVATE mp_err s_mp_ntt_sroots_init(ntt_sroots *R, ntt_mont *mc, int m, mp_digit zeta); +MP_PRIVATE void s_mp_ntt_sroots_free(ntt_sroots *R); + +/* Gnerate prime for NTT round */ +MP_PRIVATE bool s_mp_ntt_find_smooth_prime(ntt_prime_ent *ent, int m, int index, mp_digit min_p); +/* Small number primetest */ +MP_PRIVATE bool s_mp_ntt_small_is_prime(mp_digit n); +/* The NTT primes are expected in a certain order */ +MP_PRIVATE void s_mp_ntt_sort_primes(ntt_prime_ent *ents, int n); + +#endif +/* NTT stuff ends here. Yes, it was a lot, as promised. */ + + + MP_PRIVATE mp_err s_mp_prime_is_divisible(const mp_int *a, bool *result) MP_WUR; MP_PRIVATE mp_err s_mp_rand_platform(void *p, size_t n) MP_WUR; MP_PRIVATE mp_err s_mp_sqr(const mp_int *a, mp_int *b) MP_WUR; @@ -231,6 +405,10 @@ MP_PRIVATE mp_err s_mp_radix_size_overestimate(const mp_int *a, const int radix, #define MP_PRECISION_FIXED_LOG ( (int) (((sizeof(mp_word) * CHAR_BIT) / 2) - 1)) #define MP_UPPER_LIMIT_FIXED_LOG ( (int) ( (sizeof(mp_word) * CHAR_BIT) - 1)) + +/* Compute the modular inverse b^(-1) mod 2^MP_DIGIT_BIT */ +MP_PRIVATE mp_digit s_mp_ntt_invmod_d(mp_digit p); + MP_PRIVATE mp_err s_mp_fp_log(const mp_int *a, mp_int *c) MP_WUR; MP_PRIVATE mp_err s_mp_fp_log_d(const mp_int *a, mp_word *c) MP_WUR;