diff --git a/tests/fft/t-adjust.c b/tests/fft/t-adjust.c new file mode 100644 index 00000000..0e4cfff5 --- /dev/null +++ b/tests/fft/t-adjust.c @@ -0,0 +1,125 @@ +/* + +Copyright 2009, 2011 William Hart. All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, are +permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of + conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright notice, this list + of conditions and the following disclaimer in the documentation and/or other materials + provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY William Hart ``AS IS'' AND ANY EXPRESS OR IMPLIED +WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL William Hart OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +The views and conclusions contained in the software and documentation are those of the +authors and should not be interpreted as representing official policies, either expressed +or implied, of William Hart. + +*/ + +#include +#include +#include +#include "gmp-impl.h" +#include "tests.h" + +/* set p = 2^wn + 1 */ +void set_p(mpz_t p, mp_size_t n, mp_bitcnt_t w) +{ + mpz_set_ui(p, 1); + mpz_mul_2exp(p, p, n*w); + mpz_add_ui(p, p, 1); +} + +void ref_adjust(mpz_t r, mpz_t i1, mpz_t p, mp_size_t i, mp_size_t w) +{ + mpz_mul_2exp(r, i1, w*i); + mpz_mod(r, r, p); +} + +int +main(void) +{ + mp_size_t c, bits, j, k, n, w, limbs; + mpz_t p, m2a, m2b, mn1; + mp_limb_t * nn1, * r1; + + gmp_randstate_t state; + + tests_start(); + fflush(stdout); + + gmp_randinit_default(state); + + mpz_init(p); + mpz_init(m2a); + mpz_init(m2b); + mpz_init(mn1); + + for (bits = GMP_LIMB_BITS; bits < 20*GMP_LIMB_BITS; bits += GMP_LIMB_BITS) + { + for (j = 1; j < 10; j++) + { + for (k = 1; k <= GMP_LIMB_BITS; k <<= 1) + { + n = bits/k; + w = j*k; + + limbs = (n*w)/GMP_LIMB_BITS; + + for (c = 0; c < n; c++) + { + set_p(p, n, w); + + nn1 = malloc((limbs+1)*sizeof(mp_limb_t)); + r1 = malloc((limbs+1)*sizeof(mp_limb_t)); + + random_fermat(nn1, state, limbs); + fermat_to_mpz(mn1, nn1, limbs); + ref_adjust(m2a, mn1, p, c, w); + + fft_adjust(r1, nn1, c, limbs, w); + fermat_to_mpz(m2b, r1, limbs); + + mpz_mod(m2a, m2a, p); + mpz_mod(m2b, m2b, p); + + if (mpz_cmp(m2a, m2a) != 0) + { + printf("FAIL:\n"); + printf("adjust error a\n"); + printf("limbs = %ld\n", limbs); + printf("n = %ld, w = %ld, c = %ld\n", n, w, c); + gmp_printf("want %Zx\n\n", m2a); + gmp_printf("got %Zx\n", m2b); + abort(); + } + + free(nn1); + free(r1); + } + } + } + } + + mpz_clear(p); + mpz_clear(m2a); + mpz_clear(m2b); + mpz_clear(mn1); + + gmp_randclear(state); + + tests_end(); + return 0; +} diff --git a/tests/fft/t-adjust_sqrt2.c b/tests/fft/t-adjust_sqrt2.c new file mode 100644 index 00000000..e609cbef --- /dev/null +++ b/tests/fft/t-adjust_sqrt2.c @@ -0,0 +1,133 @@ +/* + +Copyright 2009, 2011 William Hart. All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, are +permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of + conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright notice, this list + of conditions and the following disclaimer in the documentation and/or other materials + provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY William Hart ``AS IS'' AND ANY EXPRESS OR IMPLIED +WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL William Hart OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +The views and conclusions contained in the software and documentation are those of the +authors and should not be interpreted as representing official policies, either expressed +or implied, of William Hart. + +*/ + +#include +#include +#include +#include "gmp-impl.h" +#include "tests.h" + +/* set p = 2^wn + 1 */ +void set_p(mpz_t p, mp_size_t n, mp_bitcnt_t w) +{ + mpz_set_ui(p, 1); + mpz_mul_2exp(p, p, n*w); + mpz_add_ui(p, p, 1); +} + +void ref_adjust_sqrt2(mpz_t r, mpz_t i1, mpz_t p, mp_size_t i, mp_size_t limbs, mp_size_t w) +{ + mpz_mul_2exp(r, i1, (w/2)*i + i/2); + if (i & 1) + { + mpz_mul_2exp(i1, r, 3*limbs*GMP_LIMB_BITS/4); + mpz_mul_2exp(r, r, limbs*GMP_LIMB_BITS/4); + mpz_sub(r, i1, r); + } + mpz_mod(r, r, p); +} + +int +main(void) +{ + mp_size_t c, bits, j, k, n, w, limbs; + mpz_t p, m2a, m2b, mn1; + mp_limb_t * nn1, * r1, * temp; + + gmp_randstate_t state; + + tests_start(); + fflush(stdout); + + gmp_randinit_default(state); + + mpz_init(p); + mpz_init(m2a); + mpz_init(m2b); + mpz_init(mn1); + + for (bits = GMP_LIMB_BITS; bits < 20*GMP_LIMB_BITS; bits += GMP_LIMB_BITS) + { + for (j = 1; j < 10; j++) + { + for (k = 1; k <= GMP_LIMB_BITS; k <<= 1) + { + n = bits/k; + w = j*k; + + limbs = (n*w)/GMP_LIMB_BITS; + + for (c = 1; c < 2*n; c+=2) + { + set_p(p, n, w); + + nn1 = malloc((limbs+1)*sizeof(mp_limb_t)); + r1 = malloc((limbs+1)*sizeof(mp_limb_t)); + temp = malloc((limbs+1)*sizeof(mp_limb_t)); + + random_fermat(nn1, state, limbs); + fermat_to_mpz(mn1, nn1, limbs); + ref_adjust_sqrt2(m2a, mn1, p, c, limbs, w); + + fft_adjust_sqrt2(r1, nn1, c, limbs, w, temp); + fermat_to_mpz(m2b, r1, limbs); + + mpz_mod(m2a, m2a, p); + mpz_mod(m2b, m2b, p); + + if (mpz_cmp(m2a, m2b) != 0) + { + printf("FAIL:\n"); + printf("adjust_sqrt2 error a\n"); + printf("limbs = %ld\n", limbs); + printf("n = %ld, w = %ld, c = %ld\n", n, w, c); + gmp_printf("want %Zx\n\n", m2a); + gmp_printf("got %Zx\n", m2b); + abort(); + } + + free(temp); + free(nn1); + free(r1); + } + } + } + } + + mpz_clear(p); + mpz_clear(m2a); + mpz_clear(m2b); + mpz_clear(mn1); + + gmp_randclear(state); + + tests_end(); + return 0; +} diff --git a/tests/fft/t-butterfly.c b/tests/fft/t-butterfly.c new file mode 100644 index 00000000..aa0d6aba --- /dev/null +++ b/tests/fft/t-butterfly.c @@ -0,0 +1,223 @@ +/* + +Copyright 2009, 2011 William Hart. All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, are +permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of + conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright notice, this list + of conditions and the following disclaimer in the documentation and/or other materials + provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY William Hart ``AS IS'' AND ANY EXPRESS OR IMPLIED +WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL William Hart OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +The views and conclusions contained in the software and documentation are those of the +authors and should not be interpreted as representing official policies, either expressed +or implied, of William Hart. + +*/ + +#include +#include +#include +#include "gmp-impl.h" +#include "tests.h" + +/* set p = 2^wn + 1 */ +void set_p(mpz_t p, mp_size_t n, mp_bitcnt_t w) +{ + mpz_set_ui(p, 1); + mpz_mul_2exp(p, p, n*w); + mpz_add_ui(p, p, 1); +} + +void ref_fft_butterfly(mpz_t s, mpz_t t, mpz_t i1, mpz_t i2, + mpz_t p, mp_size_t i, mp_size_t w) +{ + mpz_add(s, i1, i2); + mpz_sub(t, i1, i2); + mpz_mul_2exp(t, t, i*w); + mpz_mod(s, s, p); + mpz_mod(t, t, p); +} + +void ref_ifft_butterfly(mpz_t s, mpz_t t, mpz_t i1, mpz_t i2, + mpz_t p, mp_size_t i, mp_size_t n, mp_size_t w) +{ + mpz_mul_2exp(i2, i2, 2*n*w - i*w); + mpz_add(s, i1, i2); + mpz_sub(t, i1, i2); + mpz_mod(s, s, p); + mpz_mod(t, t, p); +} + +int +main(void) +{ + mp_size_t c, bits, j, k, n, w, limbs; + mpz_t p, ma, mb, m2a, m2b, mn1, mn2; + mp_limb_t * nn1, * nn2, * r1, * r2; + + gmp_randstate_t state; + + tests_start(); + fflush(stdout); + + gmp_randinit_default(state); + + mpz_init(p); + mpz_init(ma); + mpz_init(mb); + mpz_init(m2a); + mpz_init(m2b); + mpz_init(mn1); + mpz_init(mn2); + + for (bits = GMP_LIMB_BITS; bits < 20*GMP_LIMB_BITS; bits += GMP_LIMB_BITS) + { + for (j = 1; j < 10; j++) + { + for (k = 1; k <= GMP_LIMB_BITS; k <<= 1) + { + n = bits/k; + w = j*k; + + limbs = (n*w)/GMP_LIMB_BITS; + + for (c = 0; c < n; c++) + { + nn1 = malloc((limbs + 1)*sizeof(mp_limb_t)); + nn2 = malloc((limbs + 1)*sizeof(mp_limb_t)); + r1 = malloc((limbs + 1)*sizeof(mp_limb_t)); + r2 = malloc((limbs + 1)*sizeof(mp_limb_t)); + random_fermat(nn1, state, limbs); + random_fermat(nn2, state, limbs); + + fermat_to_mpz(mn1, nn1, limbs); + fermat_to_mpz(mn2, nn2, limbs); + set_p(p, n, w); + + fft_butterfly(r1, r2, nn1, nn2, c, limbs, w); + fermat_to_mpz(m2a, r1, limbs); + fermat_to_mpz(m2b, r2, limbs); + + mpz_mod(m2a, m2a, p); + mpz_mod(m2b, m2b, p); + ref_fft_butterfly(ma, mb, mn1, mn2, p, c, w); + + if (mpz_cmp(ma, m2a) != 0) + { + printf("FAIL:\n"); + printf("fft_butterfly error a\n"); + printf("limbs = %ld\n", limbs); + printf("n = %ld, w = %ld, k = %ld, c = %ld\n", n, w, k, c); + gmp_printf("want %Zx\n\n", ma); + gmp_printf("got %Zx\n", m2a); + abort(); + } + if (mpz_cmp(mb, m2b) != 0) + { + printf("FAIL:\n"); + printf("fft_butterfly error b\n"); + printf("limbs = %ld\n", limbs); + printf("n = %ld, w = %ld, k = %ld, c = %ld\n", n, w, k, c); + gmp_printf("want %Zx\n\n", mb); + gmp_printf("got %Zx\n", m2b); + abort(); + } + + free(nn1); + free(nn2); + free(r1); + free(r2); + } + } + } + } + + for (bits = GMP_LIMB_BITS; bits < 20*GMP_LIMB_BITS; bits += GMP_LIMB_BITS) + { + for (j = 1; j < 10; j++) + { + for (k = 1; k <= GMP_LIMB_BITS; k <<= 1) + { + n = bits/k; + w = j*k; + + limbs = (n*w)/GMP_LIMB_BITS; + + for (c = 0; c < n; c++) + { + nn1 = malloc((limbs + 1)*sizeof(mp_limb_t)); + nn2 = malloc((limbs + 1)*sizeof(mp_limb_t)); + r1 = malloc((limbs + 1)*sizeof(mp_limb_t)); + r2 = malloc((limbs + 1)*sizeof(mp_limb_t)); + random_fermat(nn1, state, limbs); + random_fermat(nn2, state, limbs); + + fermat_to_mpz(mn1, nn1, limbs); + fermat_to_mpz(mn2, nn2, limbs); + set_p(p, n, w); + + ifft_butterfly(r1, r2, nn1, nn2, c, limbs, w); + fermat_to_mpz(m2a, r1, limbs); + fermat_to_mpz(m2b, r2, limbs); + + mpz_mod(m2a, m2a, p); + mpz_mod(m2b, m2b, p); + ref_ifft_butterfly(ma, mb, mn1, mn2, p, c, n, w); + + if (mpz_cmp(ma, m2a) != 0) + { + printf("FAIL:\n"); + printf("ifft_butterfly error a\n"); + printf("limbs = %ld\n", limbs); + printf("n = %ld, w = %ld, k = %ld, c = %ld\n", n, w, k, c); + gmp_printf("want %Zx\n\n", ma); + gmp_printf("got %Zx\n", m2a); + abort(); + } + if (mpz_cmp(mb, m2b) != 0) + { + printf("FAIL:\n"); + printf("ifft_butterfly error b\n"); + printf("limbs = %ld\n", limbs); + printf("n = %ld, w = %ld, k = %ld, c = %ld\n", n, w, k, c); + gmp_printf("want %Zx\n\n", mb); + gmp_printf("got %Zx\n", m2b); + abort(); + } + + free(nn1); + free(nn2); + free(r1); + free(r2); + } + } + } + } + + mpz_clear(p); + mpz_clear(ma); + mpz_clear(mb); + mpz_clear(m2a); + mpz_clear(m2b); + mpz_clear(mn1); + mpz_clear(mn2); + + gmp_randclear(state); + + tests_end(); + return 0; +} diff --git a/tests/fft/t-butterfly_lshB.c b/tests/fft/t-butterfly_lshB.c new file mode 100644 index 00000000..98764e2a --- /dev/null +++ b/tests/fft/t-butterfly_lshB.c @@ -0,0 +1,159 @@ +/* + +Copyright 2009, 2011 William Hart. All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, are +permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of + conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright notice, this list + of conditions and the following disclaimer in the documentation and/or other materials + provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY William Hart ``AS IS'' AND ANY EXPRESS OR IMPLIED +WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL William Hart OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +The views and conclusions contained in the software and documentation are those of the +authors and should not be interpreted as representing official policies, either expressed +or implied, of William Hart. + +*/ + +#include +#include +#include +#include "gmp-impl.h" +#include "tests.h" + +/* set p = 2^wn + 1 */ +void set_p(mpz_t p, mp_size_t n, mp_bitcnt_t w) +{ + mpz_set_ui(p, 1); + mpz_mul_2exp(p, p, n*w); + mpz_add_ui(p, p, 1); +} + +void ref_butterfly_lshB(mpz_t t, mpz_t u, mpz_t i1, mpz_t i2, + mpz_t p, mp_size_t x, mp_size_t y) +{ + mpz_add(t, i1, i2); + mpz_sub(u, i1, i2); + mpz_mul_2exp(t, t, x*GMP_LIMB_BITS); + mpz_mul_2exp(u, u, y*GMP_LIMB_BITS); + mpz_mod(t, t, p); + mpz_mod(u, u, p); +} + +int +main(void) +{ + mp_size_t c, bits, j, k, n, w, limbs; + mp_limb_t x, y; + mpz_t p, ma, mb, m2a, m2b, mn1, mn2; + mp_limb_t * nn1, * nn2, * r1, * r2; + + gmp_randstate_t state; + + tests_start(); + fflush(stdout); + + gmp_randinit_default(state); + + mpz_init(p); + mpz_init(ma); + mpz_init(mb); + mpz_init(m2a); + mpz_init(m2b); + mpz_init(mn1); + mpz_init(mn2); + + for (bits = GMP_LIMB_BITS; bits < 20*GMP_LIMB_BITS; bits += GMP_LIMB_BITS) + { + for (j = 1; j < 10; j++) + { + for (k = 1; k <= GMP_LIMB_BITS; k <<= 1) + { + n = bits/k; + w = j*k; + + limbs = (n*w)/GMP_LIMB_BITS; + + for (c = 0; c < limbs; c++) + { + mpn_rrandom(&x, state, 1); + mpn_rrandom(&y, state, 1); + x %= limbs; + y %= limbs; + + nn1 = malloc((limbs + 1)*sizeof(mp_limb_t)); + nn2 = malloc((limbs + 1)*sizeof(mp_limb_t)); + r1 = malloc((limbs + 1)*sizeof(mp_limb_t)); + r2 = malloc((limbs + 1)*sizeof(mp_limb_t)); + mpn_rrandom(nn1, state, limbs); + random_fermat(nn1, state, limbs); + random_fermat(nn2, state, limbs); + + fermat_to_mpz(mn1, nn1, limbs); + fermat_to_mpz(mn2, nn2, limbs); + set_p(p, n, w); + + butterfly_lshB(r1, r2, nn1, nn2, limbs, x, y); + fermat_to_mpz(m2a, r1, limbs); + fermat_to_mpz(m2b, r2, limbs); + + mpz_mod(m2a, m2a, p); + mpz_mod(m2b, m2b, p); + ref_butterfly_lshB(ma, mb, mn1, mn2, p, x, y); + + if (mpz_cmp(ma, m2a) != 0) + { + printf("FAIL:\n"); + printf("butterfly_lshB error a\n"); + printf("x = %ld, y = %ld, limbs = %ld\n", x, y, limbs); + printf("n = %ld, w = %ld, k = %ld, c = %ld\n", n, w, k, c); + gmp_printf("want %Zx\n\n", ma); + gmp_printf("got %Zx\n", m2a); + abort(); + } + if (mpz_cmp(mb, m2b) != 0) + { + printf("FAIL:\n"); + printf("butterfly_lshB error b\n"); + printf("x = %ld, y = %ld, limbs = %ld\n", x, y, limbs); + printf("n = %ld, w = %ld, k = %ld, c = %ld\n", n, w, k, c); + gmp_printf("want %Zx\n\n", mb); + gmp_printf("got %Zx\n", m2b); + abort(); + } + + free(nn1); + free(nn2); + free(r1); + free(r2); + } + } + } + } + + mpz_clear(p); + mpz_clear(ma); + mpz_clear(mb); + mpz_clear(m2a); + mpz_clear(m2b); + mpz_clear(mn1); + mpz_clear(mn2); + + gmp_randclear(state); + + tests_end(); + return 0; +} diff --git a/tests/fft/t-butterfly_rshB.c b/tests/fft/t-butterfly_rshB.c new file mode 100644 index 00000000..4cea6543 --- /dev/null +++ b/tests/fft/t-butterfly_rshB.c @@ -0,0 +1,173 @@ +/* + +Copyright 2009, 2011 William Hart. All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, are +permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of + conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright notice, this list + of conditions and the following disclaimer in the documentation and/or other materials + provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY William Hart ``AS IS'' AND ANY EXPRESS OR IMPLIED +WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL William Hart OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +The views and conclusions contained in the software and documentation are those of the +authors and should not be interpreted as representing official policies, either expressed +or implied, of William Hart. + +*/ + +#include +#include +#include +#include "gmp-impl.h" +#include "tests.h" + +/* set p = 2^wn + 1 */ +void set_p(mpz_t p, mp_size_t n, mp_bitcnt_t w) +{ + mpz_set_ui(p, 1); + mpz_mul_2exp(p, p, n*w); + mpz_add_ui(p, p, 1); +} + +void ref_butterfly_rshB(mpz_t t, mpz_t u, mpz_t i1, mpz_t i2, + mpz_t p, mp_size_t x, mp_size_t y) +{ + mpz_t mult1, mult2; + + mpz_init(mult1); + mpz_init(mult2); + + mpz_set_ui(mult1, 1); + mpz_mul_2exp(mult1, mult1, x*GMP_LIMB_BITS); + mpz_invert(mult1, mult1, p); + mpz_set_ui(mult2, 1); + mpz_mul_2exp(mult2, mult2, y*GMP_LIMB_BITS); + mpz_invert(mult2, mult2, p); + mpz_mul(mult1, mult1, i1); + mpz_mul(mult2, mult2, i2); + mpz_add(t, mult1, mult2); + mpz_sub(u, mult1, mult2); + mpz_mod(t, t, p); + mpz_mod(u, u, p); + + mpz_clear(mult1); + mpz_clear(mult2); +} + +int +main(void) +{ + mp_size_t c, bits, j, k, n, w, limbs; + mp_limb_t x, y; + mpz_t p, ma, mb, m2a, m2b, mn1, mn2; + mp_limb_t * nn1, * nn2, * r1, * r2; + + gmp_randstate_t state; + + tests_start(); + fflush(stdout); + + gmp_randinit_default(state); + + mpz_init(p); + mpz_init(ma); + mpz_init(mb); + mpz_init(m2a); + mpz_init(m2b); + mpz_init(mn1); + mpz_init(mn2); + + for (bits = GMP_LIMB_BITS; bits < 20*GMP_LIMB_BITS; bits += GMP_LIMB_BITS) + { + for (j = 1; j < 10; j++) + { + for (k = 1; k <= GMP_LIMB_BITS; k <<= 1) + { + n = bits/k; + w = j*k; + + limbs = (n*w)/GMP_LIMB_BITS; + + for (c = 0; c < limbs; c++) + { + mpn_rrandom(&x, state, 1); + x %= limbs; + mpn_rrandom(&y, state, 1); + y %= limbs; + + nn1 = malloc((limbs + 1)*sizeof(mp_limb_t)); + nn2 = malloc((limbs + 1)*sizeof(mp_limb_t)); + r1 = malloc((limbs + 1)*sizeof(mp_limb_t)); + r2 = malloc((limbs + 1)*sizeof(mp_limb_t)); + mpn_rrandom(nn1, state, limbs); + random_fermat(nn1, state, limbs); + random_fermat(nn2, state, limbs); + + fermat_to_mpz(mn1, nn1, limbs); + fermat_to_mpz(mn2, nn2, limbs); + set_p(p, n, w); + + butterfly_rshB(r1, r2, nn1, nn2, limbs, x, y); + fermat_to_mpz(m2a, r1, limbs); + fermat_to_mpz(m2b, r2, limbs); + + mpz_mod(m2a, m2a, p); + mpz_mod(m2b, m2b, p); + ref_butterfly_rshB(ma, mb, mn1, mn2, p, x, y); + + if (mpz_cmp(ma, m2a) != 0) + { + printf("FAIL:\n"); + printf("butterfly_rshB error a\n"); + printf("x = %ld, y = %ld, limbs = %ld\n", x, y, limbs); + printf("n = %ld, w = %ld, k = %ld, c = %ld\n", n, w, k, c); + gmp_printf("want %Zx\n\n", ma); + gmp_printf("got %Zx\n", m2a); + abort(); + } + if (mpz_cmp(mb, m2b) != 0) + { + printf("FAIL:\n"); + printf("butterfly_rshB error b\n"); + printf("x = %ld, y = %ld, limbs = %ld\n", x, y, limbs); + printf("n = %ld, w = %ld, k = %ld, c = %ld\n", n, w, k, c); + gmp_printf("want %Zx\n\n", mb); + gmp_printf("got %Zx\n", m2b); + abort(); + } + + free(nn1); + free(nn2); + free(r1); + free(r2); + } + } + } + } + + mpz_clear(p); + mpz_clear(ma); + mpz_clear(mb); + mpz_clear(m2a); + mpz_clear(m2b); + mpz_clear(mn1); + mpz_clear(mn2); + + gmp_randclear(state); + + tests_end(); + return 0; +} diff --git a/tests/fft/t-butterfly_sqrt2.c b/tests/fft/t-butterfly_sqrt2.c new file mode 100644 index 00000000..d988cfa8 --- /dev/null +++ b/tests/fft/t-butterfly_sqrt2.c @@ -0,0 +1,235 @@ +/* + +Copyright 2009, 2011 William Hart. All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, are +permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of + conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright notice, this list + of conditions and the following disclaimer in the documentation and/or other materials + provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY William Hart ``AS IS'' AND ANY EXPRESS OR IMPLIED +WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL William Hart OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +The views and conclusions contained in the software and documentation are those of the +authors and should not be interpreted as representing official policies, either expressed +or implied, of William Hart. + +*/ + +#include +#include +#include +#include "gmp-impl.h" +#include "tests.h" + +/* set p = 2^wn + 1 */ +void set_p(mpz_t p, mp_size_t n, mp_bitcnt_t w) +{ + mpz_set_ui(p, 1); + mpz_mul_2exp(p, p, n*w); + mpz_add_ui(p, p, 1); +} + +void ref_fft_butterfly_sqrt2(mpz_t s, mpz_t t, mpz_t i1, mpz_t i2, + mpz_t p, mp_size_t i, mp_size_t limbs, mp_size_t w) +{ + mpz_sub(t, i1, i2); + mpz_mul_2exp(t, t, i*(w/2) + i/2); + mpz_mul_2exp(s, t, 3*limbs*GMP_LIMB_BITS/4); + mpz_mul_2exp(t, t, limbs*GMP_LIMB_BITS/4); + mpz_sub(t, s, t); + mpz_add(s, i1, i2); + mpz_mod(s, s, p); + mpz_mod(t, t, p); +} + +void ref_ifft_butterfly_sqrt2(mpz_t s, mpz_t t, mpz_t i1, mpz_t i2, + mpz_t p, mp_size_t i, mp_size_t n, mp_size_t limbs, mp_size_t w) +{ + mpz_mul_2exp(s, i2, 2*n*w - i*(w/2) - 1 - i/2); + mpz_mul_2exp(t, s, 3*limbs*GMP_LIMB_BITS/4); + mpz_mul_2exp(s, s, limbs*GMP_LIMB_BITS/4); + mpz_sub(i2, t, s); + mpz_add(s, i1, i2); + mpz_sub(t, i1, i2); + mpz_mod(s, s, p); + mpz_mod(t, t, p); +} + +int +main(void) +{ + mp_size_t c, bits, j, k, n, w, limbs; + mpz_t p, ma, mb, m2a, m2b, mn1, mn2; + mp_limb_t * nn1, * nn2, * r1, * r2, * temp; + + gmp_randstate_t state; + + tests_start(); + fflush(stdout); + + gmp_randinit_default(state); + + mpz_init(p); + mpz_init(ma); + mpz_init(mb); + mpz_init(m2a); + mpz_init(m2b); + mpz_init(mn1); + mpz_init(mn2); + + for (bits = GMP_LIMB_BITS; bits < 20*GMP_LIMB_BITS; bits += GMP_LIMB_BITS) + { + for (j = 1; j < 10; j++) + { + for (k = 1; k <= GMP_LIMB_BITS; k <<= 1) + { + n = bits/k; + w = j*k; + + if ((w & 1) == 0) continue; /* w must be odd here */ + + limbs = (n*w)/GMP_LIMB_BITS; + + for (c = 1; c < 2*n; c+=2) + { + nn1 = malloc((limbs + 1)*sizeof(mp_limb_t)); + nn2 = malloc((limbs + 1)*sizeof(mp_limb_t)); + temp = malloc((limbs + 1)*sizeof(mp_limb_t)); + r1 = malloc((limbs + 1)*sizeof(mp_limb_t)); + r2 = malloc((limbs + 1)*sizeof(mp_limb_t)); + random_fermat(nn1, state, limbs); + random_fermat(nn2, state, limbs); + + fermat_to_mpz(mn1, nn1, limbs); + fermat_to_mpz(mn2, nn2, limbs); + set_p(p, n, w); + + fft_butterfly_sqrt2(r1, r2, nn1, nn2, c, limbs, w, temp); + fermat_to_mpz(m2a, r1, limbs); + fermat_to_mpz(m2b, r2, limbs); + + mpz_mod(m2a, m2a, p); + mpz_mod(m2b, m2b, p); + ref_fft_butterfly_sqrt2(ma, mb, mn1, mn2, p, c, limbs, w); + + if (mpz_cmp(ma, m2a) != 0) + { + printf("FAIL:\n"); + printf("fft_butterfly_sqrt2 error a\n"); + printf("limbs = %ld\n", limbs); + printf("n = %ld, w = %ld, k = %ld, c = %ld\n", n, w, k, c); + gmp_printf("want %Zx\n\n", ma); + gmp_printf("got %Zx\n", m2a); + abort(); + } + if (mpz_cmp(mb, m2b) != 0) + { + printf("FAIL:\n"); + printf("fft_butterfly_sqrt2 error b\n"); + printf("limbs = %ld\n", limbs); + printf("n = %ld, w = %ld, k = %ld, c = %ld\n", n, w, k, c); + gmp_printf("want %Zx\n\n", mb); + gmp_printf("got %Zx\n", m2b); + abort(); + } + + free(temp); + free(nn1); + free(nn2); + free(r1); + free(r2); + } + } + } + } + + for (bits = GMP_LIMB_BITS; bits < 20*GMP_LIMB_BITS; bits += GMP_LIMB_BITS) + { + for (j = 1; j < 10; j++) + { + for (k = 1; k <= GMP_LIMB_BITS; k <<= 1) + { + n = bits/k; + w = j*k; + if ((w & 1) == 0) continue; /* w must be odd here */ + + limbs = (n*w)/GMP_LIMB_BITS; + + for (c = 1; c < 2*n; c+=2) + { + nn1 = malloc((limbs + 1)*sizeof(mp_limb_t)); + nn2 = malloc((limbs + 1)*sizeof(mp_limb_t)); + temp = malloc((limbs + 1)*sizeof(mp_limb_t)); + r1 = malloc((limbs + 1)*sizeof(mp_limb_t)); + r2 = malloc((limbs + 1)*sizeof(mp_limb_t)); + random_fermat(nn1, state, limbs); + random_fermat(nn2, state, limbs); + + fermat_to_mpz(mn1, nn1, limbs); + fermat_to_mpz(mn2, nn2, limbs); + set_p(p, n, w); + + ifft_butterfly_sqrt2(r1, r2, nn1, nn2, c, limbs, w, temp); + fermat_to_mpz(m2a, r1, limbs); + fermat_to_mpz(m2b, r2, limbs); + + mpz_mod(m2a, m2a, p); + mpz_mod(m2b, m2b, p); + ref_ifft_butterfly_sqrt2(ma, mb, mn1, mn2, p, c, n, limbs, w); + + if (mpz_cmp(ma, m2a) != 0) + { + printf("FAIL:\n"); + printf("ifft_butterfly_sqrt2 error a\n"); + printf("limbs = %ld\n", limbs); + printf("n = %ld, w = %ld, k = %ld, c = %ld\n", n, w, k, c); + gmp_printf("want %Zx\n\n", ma); + gmp_printf("got %Zx\n", m2a); + abort(); + } + if (mpz_cmp(mb, m2b) != 0) + { + printf("FAIL:\n"); + printf("ifft_butterfly_sqrt2 error b\n"); + printf("limbs = %ld\n", limbs); + printf("n = %ld, w = %ld, k = %ld, c = %ld\n", n, w, k, c); + gmp_printf("want %Zx\n\n", mb); + gmp_printf("got %Zx\n", m2b); + abort(); + } + free(temp); + free(nn1); + free(nn2); + free(r1); + free(r2); + } + } + } + } + + mpz_clear(p); + mpz_clear(ma); + mpz_clear(mb); + mpz_clear(m2a); + mpz_clear(m2b); + mpz_clear(mn1); + mpz_clear(mn2); + + gmp_randclear(state); + + tests_end(); + return 0; +} diff --git a/tests/fft/t-butterfly_twiddle.c b/tests/fft/t-butterfly_twiddle.c new file mode 100644 index 00000000..4d016349 --- /dev/null +++ b/tests/fft/t-butterfly_twiddle.c @@ -0,0 +1,234 @@ +/* + +Copyright 2009, 2011 William Hart. All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, are +permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of + conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright notice, this list + of conditions and the following disclaimer in the documentation and/or other materials + provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY William Hart ``AS IS'' AND ANY EXPRESS OR IMPLIED +WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL William Hart OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +The views and conclusions contained in the software and documentation are those of the +authors and should not be interpreted as representing official policies, either expressed +or implied, of William Hart. + +*/ + +#include +#include +#include +#include "gmp-impl.h" +#include "tests.h" + +/* set p = 2^wn + 1 */ +void set_p(mpz_t p, mp_size_t n, mp_bitcnt_t w) +{ + mpz_set_ui(p, 1); + mpz_mul_2exp(p, p, n*w); + mpz_add_ui(p, p, 1); +} + +void ref_fft_butterfly_twiddle(mpz_t s, mpz_t t, mpz_t i1, mpz_t i2, + mpz_t p, mp_size_t i, mp_size_t w, mp_bitcnt_t b1, mp_bitcnt_t b2) +{ + mpz_add(s, i1, i2); + mpz_sub(t, i1, i2); + mpz_mul_2exp(s, s, b1); + mpz_mul_2exp(t, t, b2); + mpz_mod(s, s, p); + mpz_mod(t, t, p); +} + +void ref_ifft_butterfly_twiddle(mpz_t s, mpz_t t, mpz_t i1, mpz_t i2, + mpz_t p, mp_size_t i, mp_size_t n, mp_size_t w, mp_bitcnt_t b1, mp_bitcnt_t b2) +{ + mpz_mul_2exp(i1, i1, 2*n*w - b1); + mpz_mul_2exp(i2, i2, 2*n*w - b2); + mpz_add(s, i1, i2); + mpz_sub(t, i1, i2); + mpz_mod(s, s, p); + mpz_mod(t, t, p); +} + +int +main(void) +{ + mp_size_t c, bits, j, k, n, w, limbs; + mpz_t p, ma, mb, m2a, m2b, mn1, mn2; + mp_limb_t * nn1, * nn2, * r1, * r2; + mp_bitcnt_t b1, b2; + + gmp_randstate_t state; + + tests_start(); + fflush(stdout); + + gmp_randinit_default(state); + + mpz_init(p); + mpz_init(ma); + mpz_init(mb); + mpz_init(m2a); + mpz_init(m2b); + mpz_init(mn1); + mpz_init(mn2); + + for (bits = GMP_LIMB_BITS; bits < 20*GMP_LIMB_BITS; bits += GMP_LIMB_BITS) + { + for (j = 1; j < 10; j++) + { + for (k = 1; k <= GMP_LIMB_BITS; k <<= 1) + { + n = bits/k; + w = j*k; + + limbs = (n*w)/GMP_LIMB_BITS; + + for (c = 0; c < n; c++) + { + mpn_rrandom(&b1, state, 1); + b1 %= n * w; + mpn_rrandom(&b2, state, 1); + b2 %= n * w; + nn1 = malloc((limbs + 1)*sizeof(mp_limb_t)); + nn2 = malloc((limbs + 1)*sizeof(mp_limb_t)); + r1 = malloc((limbs + 1)*sizeof(mp_limb_t)); + r2 = malloc((limbs + 1)*sizeof(mp_limb_t)); + random_fermat(nn1, state, limbs); + random_fermat(nn2, state, limbs); + + fermat_to_mpz(mn1, nn1, limbs); + fermat_to_mpz(mn2, nn2, limbs); + set_p(p, n, w); + + fft_butterfly_twiddle(r1, r2, nn1, nn2, limbs, b1, b2); + fermat_to_mpz(m2a, r1, limbs); + fermat_to_mpz(m2b, r2, limbs); + + mpz_mod(m2a, m2a, p); + mpz_mod(m2b, m2b, p); + ref_fft_butterfly_twiddle(ma, mb, mn1, mn2, p, c, w, b1, b2); + + if (mpz_cmp(ma, m2a) != 0) + { + printf("FAIL:\n"); + printf("fft_butterfly_twiddle error a\n"); + printf("limbs = %ld\n", limbs); + printf("n = %ld, w = %ld, k = %ld, c = %ld\n", n, w, k, c); + gmp_printf("want %Zx\n\n", ma); + gmp_printf("got %Zx\n", m2a); + abort(); + } + if (mpz_cmp(mb, m2b) != 0) + { + printf("FAIL:\n"); + printf("fft_butterfly_twiddle error b\n"); + printf("limbs = %ld\n", limbs); + printf("n = %ld, w = %ld, k = %ld, c = %ld\n", n, w, k, c); + gmp_printf("want %Zx\n\n", mb); + gmp_printf("got %Zx\n", m2b); + abort(); + } + + free(nn1); + free(nn2); + free(r1); + free(r2); + } + } + } + } + + for (bits = GMP_LIMB_BITS; bits < 20*GMP_LIMB_BITS; bits += GMP_LIMB_BITS) + { + for (j = 1; j < 10; j++) + { + for (k = 1; k <= GMP_LIMB_BITS; k <<= 1) + { + n = bits/k; + w = j*k; + + limbs = (n*w)/GMP_LIMB_BITS; + + for (c = 0; c < n; c++) + { + mpn_rrandom(&b1, state, 1); + b1 %= n * w; + mpn_rrandom(&b2, state, 1); + b2 %= n * w; + nn1 = malloc((limbs + 1)*sizeof(mp_limb_t)); + nn2 = malloc((limbs + 1)*sizeof(mp_limb_t)); + r1 = malloc((limbs + 1)*sizeof(mp_limb_t)); + r2 = malloc((limbs + 1)*sizeof(mp_limb_t)); + random_fermat(nn1, state, limbs); + random_fermat(nn2, state, limbs); + + fermat_to_mpz(mn1, nn1, limbs); + fermat_to_mpz(mn2, nn2, limbs); + set_p(p, n, w); + + ifft_butterfly_twiddle(r1, r2, nn1, nn2, limbs, b1, b2); + fermat_to_mpz(m2a, r1, limbs); + fermat_to_mpz(m2b, r2, limbs); + + mpz_mod(m2a, m2a, p); + mpz_mod(m2b, m2b, p); + ref_ifft_butterfly_twiddle(ma, mb, mn1, mn2, p, c, n, w, b1, b2); + + if (mpz_cmp(ma, m2a) != 0) + { + printf("FAIL:\n"); + printf("ifft_butterfly_twiddle error a\n"); + printf("limbs = %ld\n", limbs); + printf("n = %ld, w = %ld, k = %ld, c = %ld\n", n, w, k, c); + gmp_printf("want %Zx\n\n", ma); + gmp_printf("got %Zx\n", m2a); + abort(); + } + if (mpz_cmp(mb, m2b) != 0) + { + printf("FAIL:\n"); + printf("ifft_butterfly_twiddle error b\n"); + printf("limbs = %ld\n", limbs); + printf("n = %ld, w = %ld, k = %ld, c = %ld\n", n, w, k, c); + gmp_printf("want %Zx\n\n", mb); + gmp_printf("got %Zx\n", m2b); + abort(); + } + + free(nn1); + free(nn2); + free(r1); + free(r2); + } + } + } + } + + mpz_clear(p); + mpz_clear(ma); + mpz_clear(mb); + mpz_clear(m2a); + mpz_clear(m2b); + mpz_clear(mn1); + mpz_clear(mn2); + + gmp_randclear(state); + + tests_end(); + return 0; +} diff --git a/tests/fft/t-div_2expmod_2expp1.c b/tests/fft/t-div_2expmod_2expp1.c new file mode 100644 index 00000000..2423fefa --- /dev/null +++ b/tests/fft/t-div_2expmod_2expp1.c @@ -0,0 +1,119 @@ +/* + +Copyright 2009, 2011 William Hart. All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, are +permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of + conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright notice, this list + of conditions and the following disclaimer in the documentation and/or other materials + provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY William Hart ``AS IS'' AND ANY EXPRESS OR IMPLIED +WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL William Hart OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +The views and conclusions contained in the software and documentation are those of the +authors and should not be interpreted as representing official policies, either expressed +or implied, of William Hart. + +*/ + +#include +#include +#include +#include "gmp-impl.h" +#include "tests.h" + +/* set p = 2^wn + 1 */ +void set_p(mpz_t p, mp_size_t n, mp_bitcnt_t w) +{ + mpz_set_ui(p, 1); + mpz_mul_2exp(p, p, n*w); + mpz_add_ui(p, p, 1); +} + +int +main(void) +{ + mp_bitcnt_t bits; + mp_size_t j, k, n, w, limbs, d; + mp_limb_t * nn, * r; + mpz_t p, m1, m2, mn1, mn2; + + gmp_randstate_t state; + + tests_start(); + fflush(stdout); + + gmp_randinit_default(state); + + mpz_init(m1); + mpz_init(m2); + mpz_init(mn1); + mpz_init(mn2); + mpz_init(p); + + /* normalisation mod p = 2^wn + 1 where B divides nw and n is a power of 2 */ + for (bits = GMP_LIMB_BITS; bits < 16*GMP_LIMB_BITS; bits += GMP_LIMB_BITS) + { + for (j = 1; j < 32; j++) + { + for (k = 1; k <= GMP_LIMB_BITS; k <<= 1) + { + for (d = 0; d < GMP_LIMB_BITS; d++) + { + n = bits/k; + w = j*k; + limbs = (n*w)/GMP_LIMB_BITS; + + nn = malloc((limbs + 1)*sizeof(mp_limb_t)); + r = malloc((limbs + 1)*sizeof(mp_limb_t)); + random_fermat(nn, state, limbs); + fermat_to_mpz(mn1, nn, limbs); + set_p(p, n, w); + + mpn_div_2expmod_2expp1(r, nn, limbs, d); + fermat_to_mpz(m2, r, limbs); + mpz_mod(m2, m2, p); + + mpz_mod(m1, mn1, p); + mpz_mul_2exp(m2, m2, d); + mpz_mod(m2, m2, p); + + if (mpz_cmp(m1, m2) != 0) + { + printf("FAIL:\n"); + printf("mpn_div_2expmod_2expp1 error\n"); + gmp_printf("want %Zx\n\n", m1); + gmp_printf("got %Zx\n", m2); + abort(); + } + } + + free(nn); + free(r); + } + } + } + + mpz_clear(mn2); + mpz_clear(mn1); + mpz_clear(m2); + mpz_clear(m1); + mpz_clear(p); + + gmp_randclear(state); + + tests_end(); + return 0; +} diff --git a/tests/fft/t-fft_ifft_mfa_truncate_sqrt2.c b/tests/fft/t-fft_ifft_mfa_truncate_sqrt2.c new file mode 100644 index 00000000..70502862 --- /dev/null +++ b/tests/fft/t-fft_ifft_mfa_truncate_sqrt2.c @@ -0,0 +1,114 @@ +/* + +Copyright 2009, 2011 William Hart. All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, are +permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of + conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright notice, this list + of conditions and the following disclaimer in the documentation and/or other materials + provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY William Hart ``AS IS'' AND ANY EXPRESS OR IMPLIED +WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL William Hart OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +The views and conclusions contained in the software and documentation are those of the +authors and should not be interpreted as representing official policies, either expressed +or implied, of William Hart. + +*/ + +#include +#include +#include +#include "gmp-impl.h" +#include "tests.h" + +int +main(void) +{ + mp_bitcnt_t depth, w; + + gmp_randstate_t state; + + tests_start(); + fflush(stdout); + + gmp_randinit_default(state); + + for (depth = 6; depth <= 13; depth++) + { + for (w = 1; w <= 5; w++) + { + mp_size_t n = (((mp_limb_t)1)< +#include +#include +#include "gmp-impl.h" +#include "tests.h" + +int +main(void) +{ + mp_bitcnt_t depth, w; + + gmp_randstate_t state; + + tests_start(); + fflush(stdout); + + gmp_randinit_default(state); + + for (depth = 6; depth <= 12; depth++) + { + for (w = 1; w <= 5; w++) + { + mp_size_t n = (((mp_limb_t)1)< +#include +#include +#include "gmp-impl.h" +#include "tests.h" + +int +main(void) +{ + mp_bitcnt_t depth, w; + + gmp_randstate_t state; + + tests_start(); + fflush(stdout); + + gmp_randinit_default(state); + + for (depth = 6; depth <= 12; depth++) + { + for (w = 1; w <= 5; w++) + { + mp_size_t n = (((mp_limb_t)1)< +#include +#include +#include "gmp-impl.h" +#include "tests.h" + +int +main(void) +{ + mp_bitcnt_t depth, w; + + gmp_randstate_t state; + + tests_start(); + fflush(stdout); + + gmp_randinit_default(state); + + for (depth = 6; depth <= 12; depth++) + { + for (w = 1; w <= 5; w++) + { + mp_size_t n = (((mp_limb_t)1)< +#include +#include +#include "gmp-impl.h" +#include "tests.h" + +int +main(void) +{ + mp_bitcnt_t depth, w; + + gmp_randstate_t state; + + tests_start(); + fflush(stdout); + + gmp_randinit_default(state); + + for (depth = 6; depth <= 12; depth++) + { + for (w = 1; w <= 5; w++) + { + mp_size_t n = (((mp_limb_t)1)< +#include +#include +#include "gmp-impl.h" +#include "tests.h" + +/* set p = 2^wn + 1 */ +void set_p(mpz_t p, mp_size_t n, mp_bitcnt_t w) +{ + mpz_set_ui(p, 1); + mpz_mul_2exp(p, p, n*w); + mpz_add_ui(p, p, 1); +} + +int +main(void) +{ + mp_bitcnt_t bits; + mp_size_t j, k, n, w, limbs, d; + mp_limb_t * nn, * r; + mpz_t p, m1, m2, mn1, mn2; + + gmp_randstate_t state; + + tests_start(); + fflush(stdout); + + gmp_randinit_default(state); + + mpz_init(m1); + mpz_init(m2); + mpz_init(mn1); + mpz_init(mn2); + mpz_init(p); + + /* normalisation mod p = 2^wn + 1 where B divides nw and n is a power of 2 */ + for (bits = GMP_LIMB_BITS; bits < 16*GMP_LIMB_BITS; bits += GMP_LIMB_BITS) + { + for (j = 1; j < 32; j++) + { + for (k = 1; k <= GMP_LIMB_BITS; k <<= 1) + { + for (d = 0; d < GMP_LIMB_BITS; d++) + { + n = bits/k; + w = j*k; + limbs = (n*w)/GMP_LIMB_BITS; + + nn = malloc((limbs + 1)*sizeof(mp_limb_t)); + r = malloc((limbs + 1)*sizeof(mp_limb_t)); + random_fermat(nn, state, limbs); + fermat_to_mpz(mn1, nn, limbs); + set_p(p, n, w); + + mpn_mul_2expmod_2expp1(r, nn, limbs, d); + fermat_to_mpz(m2, r, limbs); + mpz_mod(m2, m2, p); + + mpz_mod(m1, mn1, p); + mpz_mul_2exp(m1, m1, d); + mpz_mod(m1, m1, p); + + if (mpz_cmp(m1, m2) != 0) + { + printf("FAIL:\n"); + printf("mpn_mul_2expmod_2expp1 error\n"); + gmp_printf("want %Zx\n\n", m1); + gmp_printf("got %Zx\n", m2); + abort(); + } + } + + free(nn); + free(r); + } + } + } + + mpz_clear(mn2); + mpz_clear(mn1); + mpz_clear(m2); + mpz_clear(m1); + mpz_clear(p); + + gmp_randclear(state); + + tests_end(); + return 0; +} diff --git a/tests/fft/t-mul_fft_main.c b/tests/fft/t-mul_fft_main.c new file mode 100644 index 00000000..4ce7914c --- /dev/null +++ b/tests/fft/t-mul_fft_main.c @@ -0,0 +1,122 @@ +/* + +Copyright 2009, 2011 William Hart. All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, are +permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of + conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright notice, this list + of conditions and the following disclaimer in the documentation and/or other materials + provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY William Hart ``AS IS'' AND ANY EXPRESS OR IMPLIED +WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL William Hart OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +The views and conclusions contained in the software and documentation are those of the +authors and should not be interpreted as representing official policies, either expressed +or implied, of William Hart. + +*/ + +#include +#include +#include +#include +#include "gmp-impl.h" +#include "tests.h" + +int +main(void) +{ + mp_bitcnt_t depth, w; + + gmp_randstate_t state; + + tests_start(); + fflush(stdout); + + gmp_randinit_default(state); + + for (depth = 6; depth <= 13; depth++) + { + for (w = 1; w <= 3 - (depth >= 12); w++) + { + int iter = 1 + 200*(depth <= 8) + 80*(depth <= 9) + 10*(depth <= 10), i; + + for (i = 0; i < iter; i++) + { + mp_size_t n = (((mp_limb_t)1)<= b2 */ + { + mp_size_t t = n1; + mp_bitcnt_t tb = b1; + n1 = n2; + b1 = b2; + n2 = t; + b2 = tb; + } + + i1 = malloc(3*(n1 + n2)*sizeof(mp_limb_t)); + i2 = i1 + n1; + r1 = i2 + n2; + r2 = r1 + n1 + n2; + + mpn_urandomb(i1, state, b1); + mpn_urandomb(i2, state, b2); + + mpn_mul(r2, i1, n1, i2, n2); + mpn_mul_fft_main(r1, i1, n1, i2, n2); + + for (j = 0; j < n1 + n2; j++) + { + if (r1[j] != r2[j]) + { + printf("error in limb %ld, %lx != %lx\n", j, r1[j], r2[j]); + abort(); + } + } + + free(i1); + } + } + } + + gmp_randclear(state); + + tests_end(); + return 0; +} diff --git a/tests/fft/t-mul_mfa_truncate_sqrt2.c b/tests/fft/t-mul_mfa_truncate_sqrt2.c new file mode 100644 index 00000000..46daa39c --- /dev/null +++ b/tests/fft/t-mul_mfa_truncate_sqrt2.c @@ -0,0 +1,93 @@ +/* + +Copyright 2009, 2011 William Hart. All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, are +permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of + conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright notice, this list + of conditions and the following disclaimer in the documentation and/or other materials + provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY William Hart ``AS IS'' AND ANY EXPRESS OR IMPLIED +WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL William Hart OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +The views and conclusions contained in the software and documentation are those of the +authors and should not be interpreted as representing official policies, either expressed +or implied, of William Hart. + +*/ + +#include +#include +#include +#include "gmp-impl.h" +#include "tests.h" + +int +main(void) +{ + mp_bitcnt_t depth, w; + + gmp_randstate_t state; + + tests_start(); + fflush(stdout); + + gmp_randinit_default(state); + + for (depth = 6; depth <= 14; depth++) + { + for (w = 1; w <= 3 - (depth >= 12); w++) + { + mp_size_t n = (((mp_limb_t)1)< +#include +#include +#include "gmp-impl.h" +#include "tests.h" + +int +main(void) +{ + mp_bitcnt_t depth, w; + + gmp_randstate_t state; + + tests_start(); + fflush(stdout); + + gmp_randinit_default(state); + + for (depth = 6; depth <= 12; depth++) + { + for (w = 1; w <= 5; w++) + { + mp_size_t n = (((mp_limb_t)1)< +#include +#include +#include "gmp-impl.h" +#include "tests.h" + +int +main(void) +{ + mp_bitcnt_t depth, w; + int iters; + + gmp_randstate_t state; + + tests_start(); + fflush(stdout); + + gmp_randinit_default(state); + + for (iters = 0; iters < 100; iters++) + { + for (depth = 6; depth <= 18; depth++) + { + for (w = 1; w <= 2; w++) + { + mp_size_t n = (((mp_limb_t)1)< +#include +#include +#include "gmp-impl.h" +#include "tests.h" + +/* set p = 2^wn + 1 */ +void set_p(mpz_t p, mp_size_t n, mp_bitcnt_t w) +{ + mpz_set_ui(p, 1); + mpz_mul_2exp(p, p, n*w); + mpz_add_ui(p, p, 1); +} + +int +main(void) +{ + mp_bitcnt_t bits; + mp_size_t j, k, n, w, limbs; + mp_limb_t * nn; + mpz_t p, m1, m2; + + gmp_randstate_t state; + + tests_start(); + fflush(stdout); + + gmp_randinit_default(state); + + mpz_init(m1); + mpz_init(m2); + mpz_init(p); + + /* normalisation mod p = 2^wn + 1 where B divides nw and n is a power of 2 */ + for (bits = GMP_LIMB_BITS; bits < 32*GMP_LIMB_BITS; bits += GMP_LIMB_BITS) + { + for (j = 1; j < 32; j++) + { + for (k = 1; k <= GMP_NUMB_BITS; k <<= 1) + { + n = bits/k; + w = j*k; + limbs = (n*w)/GMP_LIMB_BITS; + + nn = malloc((limbs + 1)*sizeof(mp_limb_t)); + mpn_rrandom(nn, state, limbs + 1); + fermat_to_mpz(m1, nn, limbs); + set_p(p, n, w); + + mpn_normmod_2expp1(nn, limbs); + fermat_to_mpz(m2, nn, limbs); + mpz_mod(m1, m1, p); + + if (mpz_cmp(m1, m2) != 0) + { + printf("FAIL:\n"); + printf("mpn_normmod_2expp1 error\n"); + gmp_printf("want %Zx\n\n", m1); + gmp_printf("got %Zx\n", m2); + abort(); + } + + free(nn); + } + } + } + + mpz_clear(m2); + mpz_clear(m1); + mpz_clear(p); + + gmp_randclear(state); + + tests_end(); + return 0; +} diff --git a/tests/fft/t-split_combine_bits.c b/tests/fft/t-split_combine_bits.c new file mode 100644 index 00000000..82e9073a --- /dev/null +++ b/tests/fft/t-split_combine_bits.c @@ -0,0 +1,101 @@ +/* + +Copyright 2009, 2011 William Hart. All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, are +permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this list of + conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright notice, this list + of conditions and the following disclaimer in the documentation and/or other materials + provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY William Hart ``AS IS'' AND ANY EXPRESS OR IMPLIED +WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL William Hart OR +CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +The views and conclusions contained in the software and documentation are those of the +authors and should not be interpreted as representing official policies, either expressed +or implied, of William Hart. + +*/ + +#include +#include +#include +#include "gmp-impl.h" +#include "tests.h" + +int +main(void) +{ + int i; + mp_size_t j; + + gmp_randstate_t state; + + tests_start(); + fflush(stdout); + + gmp_randinit_default(state); + + for (i = 0; i < 10000; i++) + { + mp_limb_t total_limbs; + mp_limb_t * in; + mp_limb_t * out; + mp_bitcnt_t bits; + mp_size_t limbs; + long length; + mp_limb_t ** poly; + + mpn_rrandom(&total_limbs, state, 1); + total_limbs = total_limbs % 1000 + 1; + in = malloc(total_limbs*sizeof(mp_limb_t)); + out = calloc(total_limbs, sizeof(mp_limb_t)); + mpn_rrandom(&bits, state, 1); + bits = bits % 200 + 1; + limbs = (2*bits - 1)/GMP_LIMB_BITS + 1; + length = (total_limbs*GMP_LIMB_BITS - 1)/bits + 1; + + poly = malloc(length*sizeof(mp_limb_t *)); + for (j = 0; j < length; j++) + poly[j] = malloc((limbs + 1)*sizeof(mp_limb_t)); + + mpn_urandomb(in, state, total_limbs*GMP_LIMB_BITS); + + fft_split_bits(poly, in, total_limbs, bits, limbs); + fft_combine_bits(out, poly, length, bits, limbs, total_limbs); + + for (j = 0; j < total_limbs; j++) + { + if (in[j] != out[j]) + { + printf("FAIL:\n"); + printf("Error in limb %ld, %lu != %lu\n", j, in[j], out[j]); + abort(); + } + } + + free(in); + free(out); + + for (j = 0; j < length; j++) + free(poly[j]); + + free(poly); + } + + gmp_randclear(state); + + tests_end(); + return 0; +}