t-mul.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:6k
- /* Test mpz_cmp, mpz_mul.
- Copyright 1991, 1993, 1994, 1996, 1997, 2000, 2001, 2002, 2003, 2004 Free
- Software Foundation, Inc.
- This file is part of the GNU MP Library.
- The GNU MP Library is free software; you can redistribute it and/or modify
- it under the terms of the GNU Lesser General Public License as published by
- the Free Software Foundation; either version 3 of the License, or (at your
- option) any later version.
- The GNU MP Library is distributed in the hope that it will be useful, but
- WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
- or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
- License for more details.
- You should have received a copy of the GNU Lesser General Public License
- along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
- #include <stdio.h>
- #include <stdlib.h>
- #include "gmp.h"
- #include "gmp-impl.h"
- #include "longlong.h"
- #include "tests.h"
- void debug_mp __GMP_PROTO ((mpz_t));
- static void refmpz_mul __GMP_PROTO ((mpz_t, const mpz_t, const mpz_t));
- void dump_abort __GMP_PROTO ((int, char *, mpz_t, mpz_t, mpz_t, mpz_t));
- #define FFT_MIN_BITSIZE 100000
- char *extra_fft;
- void
- one (int i, mpz_t multiplicand, mpz_t multiplier)
- {
- mpz_t product, ref_product;
- mpz_init (product);
- mpz_init (ref_product);
- /* Test plain multiplication comparing results against reference code. */
- mpz_mul (product, multiplier, multiplicand);
- refmpz_mul (ref_product, multiplier, multiplicand);
- if (mpz_cmp (product, ref_product))
- dump_abort (i, "incorrect plain product",
- multiplier, multiplicand, product, ref_product);
- /* Test squaring, comparing results against plain multiplication */
- mpz_mul (product, multiplier, multiplier);
- mpz_set (multiplicand, multiplier);
- mpz_mul (ref_product, multiplier, multiplicand);
- if (mpz_cmp (product, ref_product))
- dump_abort (i, "incorrect square product",
- multiplier, multiplier, product, ref_product);
- mpz_clear (product);
- mpz_clear (ref_product);
- }
- int
- main (int argc, char **argv)
- {
- mpz_t op1, op2;
- int i;
- int fft_max_2exp;
- gmp_randstate_ptr rands;
- mpz_t bs;
- unsigned long bsi, size_range, fsize_range;
- tests_start ();
- rands = RANDS;
- extra_fft = getenv ("GMP_CHECK_FFT");
- fft_max_2exp = 0;
- if (extra_fft != NULL)
- fft_max_2exp = atoi (extra_fft);
- if (fft_max_2exp <= 1) /* compat with old use of GMP_CHECK_FFT */
- fft_max_2exp = 22; /* default limit, good for any machine */
- mpz_init (bs);
- mpz_init (op1);
- mpz_init (op2);
- fsize_range = 4 << 8; /* a fraction 1/256 of size_range */
- for (i = 0;; i++)
- {
- size_range = fsize_range >> 8;
- fsize_range = fsize_range * 33 / 32;
- if (size_range > fft_max_2exp)
- break;
- mpz_urandomb (bs, rands, size_range);
- mpz_rrandomb (op1, rands, mpz_get_ui (bs));
- if (i & 1)
- mpz_urandomb (bs, rands, size_range);
- mpz_rrandomb (op2, rands, mpz_get_ui (bs));
- mpz_urandomb (bs, rands, 4);
- bsi = mpz_get_ui (bs);
- if ((bsi & 0x3) == 0)
- mpz_neg (op1, op1);
- if ((bsi & 0xC) == 0)
- mpz_neg (op2, op2);
- /* printf ("%d %dn", SIZ (op1), SIZ (op2)); */
- one (i, op2, op1);
- }
- for (i = -50; i < 0; i++)
- {
- mpz_urandomb (bs, rands, 32);
- size_range = mpz_get_ui (bs) % fft_max_2exp;
- mpz_urandomb (bs, rands, size_range);
- mpz_rrandomb (op1, rands, mpz_get_ui (bs) + FFT_MIN_BITSIZE);
- mpz_urandomb (bs, rands, size_range);
- mpz_rrandomb (op2, rands, mpz_get_ui (bs) + FFT_MIN_BITSIZE);
- /* printf ("%d: %d %dn", i, SIZ (op1), SIZ (op2)); */
- fflush (stdout);
- one (-1, op2, op1);
- }
- mpz_clear (bs);
- mpz_clear (op1);
- mpz_clear (op2);
- tests_end ();
- exit (0);
- }
- static void
- refmpz_mul (mpz_t w, const mpz_t u, const mpz_t v)
- {
- mp_size_t usize = u->_mp_size;
- mp_size_t vsize = v->_mp_size;
- mp_size_t wsize;
- mp_size_t sign_product;
- mp_ptr up, vp;
- mp_ptr wp;
- mp_size_t talloc;
- sign_product = usize ^ vsize;
- usize = ABS (usize);
- vsize = ABS (vsize);
- if (usize == 0 || vsize == 0)
- {
- SIZ (w) = 0;
- return;
- }
- talloc = usize + vsize;
- up = u->_mp_d;
- vp = v->_mp_d;
- wp = __GMP_ALLOCATE_FUNC_LIMBS (talloc);
- if (usize > vsize)
- refmpn_mul (wp, up, usize, vp, vsize);
- else
- refmpn_mul (wp, vp, vsize, up, usize);
- wsize = usize + vsize;
- wsize -= wp[wsize - 1] == 0;
- MPZ_REALLOC (w, wsize);
- MPN_COPY (PTR(w), wp, wsize);
- SIZ(w) = sign_product < 0 ? -wsize : wsize;
- __GMP_FREE_FUNC_LIMBS (wp, talloc);
- }
- void
- dump_abort (int i, char *s,
- mpz_t op1, mpz_t op2, mpz_t product, mpz_t ref_product)
- {
- mp_size_t b, e;
- fprintf (stderr, "ERROR: %s in test %dn", s, i);
- fprintf (stderr, "op1 = "); debug_mp (op1);
- fprintf (stderr, "op2 = "); debug_mp (op2);
- fprintf (stderr, " product = "); debug_mp (product);
- fprintf (stderr, "ref_product = "); debug_mp (ref_product);
- for (b = 0; b < ABSIZ(ref_product); b++)
- if (PTR(ref_product)[b] != PTR(product)[b])
- break;
- for (e = ABSIZ(ref_product) - 1; e >= 0; e--)
- if (PTR(ref_product)[e] != PTR(product)[e])
- break;
- printf ("ERRORS in %ld--%ldn", b, e);
- abort();
- }
- void
- debug_mp (mpz_t x)
- {
- size_t siz = mpz_sizeinbase (x, 16);
- if (siz > 65)
- {
- mpz_t q;
- mpz_init (q);
- mpz_tdiv_q_2exp (q, x, 4 * (mpz_sizeinbase (x, 16) - 25));
- gmp_fprintf (stderr, "%ZX...", q);
- mpz_tdiv_r_2exp (q, x, 4 * 25);
- gmp_fprintf (stderr, "%025ZX [%d]n", q, (int) siz);
- mpz_clear (q);
- }
- else
- {
- gmp_fprintf (stderr, "%ZXn", x);
- }
- }