misc.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:13k
源码类别:

数学计算

开发平台:

Unix_Linux

  1. /* Miscellaneous test program support routines.
  2. Copyright 2000, 2001, 2002, 2003, 2005 Free Software Foundation, Inc.
  3. This file is part of the GNU MP Library.
  4. The GNU MP Library is free software; you can redistribute it and/or modify
  5. it under the terms of the GNU Lesser General Public License as published by
  6. the Free Software Foundation; either version 3 of the License, or (at your
  7. option) any later version.
  8. The GNU MP Library is distributed in the hope that it will be useful, but
  9. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  10. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  11. License for more details.
  12. You should have received a copy of the GNU Lesser General Public License
  13. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  14. #include "config.h"
  15. #include <ctype.h>
  16. #include <signal.h>
  17. #include <stdio.h>
  18. #include <stdlib.h>     /* for getenv */
  19. #include <string.h>
  20. #if HAVE_FLOAT_H
  21. #include <float.h>      /* for DBL_MANT_DIG */
  22. #endif
  23. #if TIME_WITH_SYS_TIME
  24. # include <sys/time.h>  /* for struct timeval */
  25. # include <time.h>
  26. #else
  27. # if HAVE_SYS_TIME_H
  28. #  include <sys/time.h>
  29. # else
  30. #  include <time.h>
  31. # endif
  32. #endif
  33. #include "gmp.h"
  34. #include "gmp-impl.h"
  35. #include "tests.h"
  36. /* The various tests setups and final checks, collected up together. */
  37. void
  38. tests_start (void)
  39. {
  40.   /* don't buffer, so output is not lost if a test causes a segv etc */
  41.   setbuf (stdout, NULL);
  42.   setbuf (stderr, NULL);
  43.   tests_memory_start ();
  44.   tests_rand_start ();
  45. }
  46. void
  47. tests_end (void)
  48. {
  49.   tests_rand_end ();
  50.   tests_memory_end ();
  51. }
  52. void
  53. tests_rand_start (void)
  54. {
  55.   gmp_randstate_ptr  rands;
  56.   char           *perform_seed;
  57.   unsigned long  seed;
  58.   if (__gmp_rands_initialized)
  59.     {
  60.       printf ("Please let tests_start() initialize the global __gmp_rands.n");
  61.       printf ("ie. ensure that function is called before the first use of RANDS.n");
  62.       abort ();
  63.     }
  64.   gmp_randinit_default (__gmp_rands);
  65.   __gmp_rands_initialized = 1;
  66.   rands = __gmp_rands;
  67.   perform_seed = getenv ("GMP_CHECK_RANDOMIZE");
  68.   if (perform_seed != NULL)
  69.     {
  70. #ifdef HAVE_STRTOUL
  71.       seed = strtoul (perform_seed, 0, 0);
  72. #else
  73.       /* This will not work right for seeds >= 2^31 on 64-bit machines.
  74.  Perhaps use atol unconditionally?  Is that ubiquitous?  */
  75.       seed = atoi (perform_seed);
  76. #endif
  77.       if (! (seed == 0 || seed == 1))
  78.         {
  79.           printf ("Re-seeding with GMP_CHECK_RANDOMIZE=%lun", seed);
  80.           gmp_randseed_ui (rands, seed);
  81.         }
  82.       else
  83.         {
  84. #if HAVE_GETTIMEOFDAY
  85.           struct timeval  tv;
  86.           gettimeofday (&tv, NULL);
  87.           seed = tv.tv_sec ^ (tv.tv_usec << 12);
  88.   seed &= 0xffffffff;
  89. #else
  90.           time_t  tv;
  91.           time (&tv);
  92.           seed = tv;
  93. #endif
  94.           gmp_randseed_ui (rands, seed);
  95.           printf ("Seed GMP_CHECK_RANDOMIZE=%lu (include this in bug reports)n", seed);
  96.         }
  97.       fflush (stdout);
  98.     }
  99. }
  100. void
  101. tests_rand_end (void)
  102. {
  103.   RANDS_CLEAR ();
  104. }
  105. /* Only used if CPU calling conventions checking is available. */
  106. mp_limb_t (*calling_conventions_function) __GMP_PROTO ((ANYARGS));
  107. /* Return p advanced to the next multiple of "align" bytes.  "align" must be
  108.    a power of 2.  Care is taken not to assume sizeof(int)==sizeof(pointer).
  109.    Using "unsigned long" avoids a warning on hpux.  */
  110. void *
  111. align_pointer (void *p, size_t align)
  112. {
  113.   gmp_intptr_t d;
  114.   d = ((gmp_intptr_t) p) & (align-1);
  115.   d = (d != 0 ? align-d : 0);
  116.   return (void *) (((char *) p) + d);
  117. }
  118. /* Note that memory allocated with this function can never be freed, because
  119.    the start address of the block allocated is lost. */
  120. void *
  121. __gmp_allocate_func_aligned (size_t bytes, size_t align)
  122. {
  123.   return align_pointer ((*__gmp_allocate_func) (bytes + align-1), align);
  124. }
  125. void *
  126. __gmp_allocate_or_reallocate (void *ptr, size_t oldsize, size_t newsize)
  127. {
  128.   if (ptr == NULL)
  129.     return (*__gmp_allocate_func) (newsize);
  130.   else
  131.     return (*__gmp_reallocate_func) (ptr, oldsize, newsize);
  132. }
  133. char *
  134. __gmp_allocate_strdup (const char *s)
  135. {
  136.   size_t  len;
  137.   char    *t;
  138.   len = strlen (s);
  139.   t = (*__gmp_allocate_func) (len+1);
  140.   memcpy (t, s, len+1);
  141.   return t;
  142. }
  143. char *
  144. strtoupper (char *s_orig)
  145. {
  146.   char  *s;
  147.   for (s = s_orig; *s != ''; s++)
  148.     if (isascii (*s))
  149.       *s = toupper (*s);
  150.   return s_orig;
  151. }
  152. void
  153. mpz_set_n (mpz_ptr z, mp_srcptr p, mp_size_t size)
  154. {
  155.   ASSERT (size >= 0);
  156.   MPN_NORMALIZE (p, size);
  157.   MPZ_REALLOC (z, size);
  158.   MPN_COPY (PTR(z), p, size);
  159.   SIZ(z) = size;
  160. }
  161. void
  162. mpz_init_set_n (mpz_ptr z, mp_srcptr p, mp_size_t size)
  163. {
  164.   ASSERT (size >= 0);
  165.   MPN_NORMALIZE (p, size);
  166.   ALLOC(z) = MAX (size, 1);
  167.   PTR(z) = __GMP_ALLOCATE_FUNC_LIMBS (ALLOC(z));
  168.   SIZ(z) = size;
  169.   MPN_COPY (PTR(z), p, size);
  170. }
  171. /* Find least significant limb position where p1,size and p2,size differ.  */
  172. mp_size_t
  173. mpn_diff_lowest (mp_srcptr p1, mp_srcptr p2, mp_size_t size)
  174. {
  175.   mp_size_t  i;
  176.   for (i = 0; i < size; i++)
  177.     if (p1[i] != p2[i])
  178.       return i;
  179.   /* no differences */
  180.   return -1;
  181. }
  182. /* Find most significant limb position where p1,size and p2,size differ.  */
  183. mp_size_t
  184. mpn_diff_highest (mp_srcptr p1, mp_srcptr p2, mp_size_t size)
  185. {
  186.   mp_size_t  i;
  187.   for (i = size-1; i >= 0; i--)
  188.     if (p1[i] != p2[i])
  189.       return i;
  190.   /* no differences */
  191.   return -1;
  192. }
  193. /* Find least significant byte position where p1,size and p2,size differ.  */
  194. mp_size_t
  195. byte_diff_lowest (const void *p1, const void *p2, mp_size_t size)
  196. {
  197.   mp_size_t  i;
  198.   for (i = 0; i < size; i++)
  199.     if (((const char *) p1)[i] != ((const char *) p2)[i])
  200.       return i;
  201.   /* no differences */
  202.   return -1;
  203. }
  204. /* Find most significant limb position where p1,size and p2,size differ.  */
  205. mp_size_t
  206. byte_diff_highest (const void *p1, const void *p2, mp_size_t size)
  207. {
  208.   mp_size_t  i;
  209.   for (i = size-1; i >= 0; i--)
  210.     if (((const char *) p1)[i] != ((const char *) p2)[i])
  211.       return i;
  212.   /* no differences */
  213.   return -1;
  214. }
  215. void
  216. mpz_set_str_or_abort (mpz_ptr z, const char *str, int base)
  217. {
  218.   if (mpz_set_str (z, str, base) != 0)
  219.     {
  220.       fprintf (stderr, "ERROR: mpz_set_str failedn");
  221.       fprintf (stderr, "   str  = "%s"n", str);
  222.       fprintf (stderr, "   base = %dn", base);
  223.       abort();
  224.     }
  225. }
  226. void
  227. mpq_set_str_or_abort (mpq_ptr q, const char *str, int base)
  228. {
  229.   if (mpq_set_str (q, str, base) != 0)
  230.     {
  231.       fprintf (stderr, "ERROR: mpq_set_str failedn");
  232.       fprintf (stderr, "   str  = "%s"n", str);
  233.       fprintf (stderr, "   base = %dn", base);
  234.       abort();
  235.     }
  236. }
  237. void
  238. mpf_set_str_or_abort (mpf_ptr f, const char *str, int base)
  239. {
  240.   if (mpf_set_str (f, str, base) != 0)
  241.     {
  242.       fprintf (stderr, "ERROR mpf_set_str failedn");
  243.       fprintf (stderr, "   str  = "%s"n", str);
  244.       fprintf (stderr, "   base = %dn", base);
  245.       abort();
  246.     }
  247. }
  248. /* Whether the absolute value of z is a power of 2. */
  249. int
  250. mpz_pow2abs_p (mpz_srcptr z)
  251. {
  252.   mp_size_t  size, i;
  253.   mp_srcptr  ptr;
  254.   size = SIZ (z);
  255.   if (size == 0)
  256.     return 0;  /* zero is not a power of 2 */
  257.   size = ABS (size);
  258.   ptr = PTR (z);
  259.   for (i = 0; i < size-1; i++)
  260.     if (ptr[i] != 0)
  261.       return 0;  /* non-zero low limb means not a power of 2 */
  262.   return POW2_P (ptr[i]);  /* high limb power of 2 */
  263. }
  264. /* Exponentially distributed between 0 and 2^nbits-1, meaning the number of
  265.    bits in the result is uniformly distributed between 0 and nbits-1.
  266.    FIXME: This is not a proper exponential distribution, since the
  267.    probability function will have a stepped shape due to using a uniform
  268.    distribution after choosing how many bits.  */
  269. void
  270. mpz_erandomb (mpz_ptr rop, gmp_randstate_t rstate, unsigned long nbits)
  271. {
  272.   mpz_urandomb (rop, rstate, gmp_urandomm_ui (rstate, nbits));
  273. }
  274. void
  275. mpz_erandomb_nonzero (mpz_ptr rop, gmp_randstate_t rstate, unsigned long nbits)
  276. {
  277.   mpz_erandomb (rop, rstate, nbits);
  278.   if (mpz_sgn (rop) == 0)
  279.     mpz_set_ui (rop, 1L);
  280. }
  281. void
  282. mpz_errandomb (mpz_ptr rop, gmp_randstate_t rstate, unsigned long nbits)
  283. {
  284.   mpz_rrandomb (rop, rstate, gmp_urandomm_ui (rstate, nbits));
  285. }
  286. void
  287. mpz_errandomb_nonzero (mpz_ptr rop, gmp_randstate_t rstate, unsigned long nbits)
  288. {
  289.   mpz_errandomb (rop, rstate, nbits);
  290.   if (mpz_sgn (rop) == 0)
  291.     mpz_set_ui (rop, 1L);
  292. }
  293. void
  294. mpz_negrandom (mpz_ptr rop, gmp_randstate_t rstate)
  295. {
  296.   mp_limb_t  n;
  297.   _gmp_rand (&n, rstate, 1);
  298.   if (n != 0)
  299.     mpz_neg (rop, rop);
  300. }
  301. mp_limb_t
  302. urandom (void)
  303. {
  304. #if GMP_NAIL_BITS == 0
  305.   mp_limb_t  n;
  306.   _gmp_rand (&n, RANDS, GMP_LIMB_BITS);
  307.   return n;
  308. #else
  309.   mp_limb_t n[2];
  310.   _gmp_rand (n, RANDS, GMP_LIMB_BITS);
  311.   return n[0] + (n[1] << GMP_NUMB_BITS);
  312. #endif
  313. }
  314. /* Call (*func)() with various random number generators. */
  315. void
  316. call_rand_algs (void (*func) __GMP_PROTO ((const char *, gmp_randstate_ptr)))
  317. {
  318.   gmp_randstate_t  rstate;
  319.   mpz_t            a;
  320.   mpz_init (a);
  321.   gmp_randinit_default (rstate);
  322.   (*func) ("gmp_randinit_default", rstate);
  323.   gmp_randclear (rstate);
  324.   gmp_randinit_mt (rstate);
  325.   (*func) ("gmp_randinit_mt", rstate);
  326.   gmp_randclear (rstate);
  327.   gmp_randinit_lc_2exp_size (rstate, 8L);
  328.   (*func) ("gmp_randinit_lc_2exp_size 8", rstate);
  329.   gmp_randclear (rstate);
  330.   gmp_randinit_lc_2exp_size (rstate, 16L);
  331.   (*func) ("gmp_randinit_lc_2exp_size 16", rstate);
  332.   gmp_randclear (rstate);
  333.   gmp_randinit_lc_2exp_size (rstate, 128L);
  334.   (*func) ("gmp_randinit_lc_2exp_size 128", rstate);
  335.   gmp_randclear (rstate);
  336.   /* degenerate always zeros */
  337.   mpz_set_ui (a, 0L);
  338.   gmp_randinit_lc_2exp (rstate, a, 0L, 8L);
  339.   (*func) ("gmp_randinit_lc_2exp a=0 c=0 m=8", rstate);
  340.   gmp_randclear (rstate);
  341.   /* degenerate always FFs */
  342.   mpz_set_ui (a, 0L);
  343.   gmp_randinit_lc_2exp (rstate, a, 0xFFL, 8L);
  344.   (*func) ("gmp_randinit_lc_2exp a=0 c=0xFF m=8", rstate);
  345.   gmp_randclear (rstate);
  346.   mpz_clear (a);
  347. }
  348. /* Return +infinity if available, or 0 if not.
  349.    We don't want to use libm, so INFINITY or other system values are not
  350.    used here.  */
  351. double
  352. tests_infinity_d (void)
  353. {
  354. #if _GMP_IEEE_FLOATS
  355.   union ieee_double_extract x;
  356.   x.s.exp = 2047;
  357.   x.s.manl = 0;
  358.   x.s.manh = 0;
  359.   x.s.sig = 0;
  360.   return x.d;
  361. #else
  362.   return 0;
  363. #endif
  364. }
  365. /* Return non-zero if d is an infinity (either positive or negative).
  366.    Don't want libm, so don't use isinf() or other system tests.  */
  367. int
  368. tests_isinf (double d)
  369. {
  370. #if _GMP_IEEE_FLOATS
  371.   union ieee_double_extract x;
  372.   x.d = d;
  373.   return (x.s.exp == 2047 && x.s.manl == 0 && x.s.manh == 0);
  374. #else
  375.   return 0;
  376. #endif
  377. }
  378. /* Set the hardware floating point rounding mode.  Same mode values as mpfr,
  379.    namely 0=nearest, 1=tozero, 2=up, 3=down.  Return 1 if successful, 0 if
  380.    not.  */
  381. int
  382. tests_hardware_setround (int mode)
  383. {
  384. #if HAVE_HOST_CPU_FAMILY_x86
  385.   int  rc;
  386.   switch (mode) {
  387.   case 0: rc = 0; break;  /* nearest */
  388.   case 1: rc = 3; break;  /* tozero  */
  389.   case 2: rc = 2; break;  /* up      */
  390.   case 3: rc = 1; break;  /* down    */
  391.   default:
  392.     return 0;
  393.   }
  394.   x86_fldcw ((x86_fstcw () & ~0xC00) | (rc << 10));
  395.   return 1;
  396. #endif
  397.   return 0;
  398. }
  399. /* Return the hardware floating point rounding mode, or -1 if unknown. */
  400. int
  401. tests_hardware_getround (void)
  402. {
  403. #if HAVE_HOST_CPU_FAMILY_x86
  404.   switch ((x86_fstcw () & ~0xC00) >> 10) {
  405.   case 0: return 0; break;  /* nearest */
  406.   case 1: return 3; break;  /* down    */
  407.   case 2: return 2; break;  /* up      */
  408.   case 3: return 1; break;  /* tozero  */
  409.   }
  410. #endif
  411.   return -1;
  412. }
  413. /* tests_dbl_mant_bits() determines by experiment the number of bits in the
  414.    mantissa of a "double".  If it's not possible to find a value (perhaps
  415.    due to the compiler optimizing too aggressively), then return 0.
  416.    This code is used rather than DBL_MANT_DIG from <float.h> since ancient
  417.    systems like SunOS don't have that file, and since one GNU/Linux ARM
  418.    system was seen where the float emulation seemed to have only 32 working
  419.    bits, not the 53 float.h claimed.  */
  420. int
  421. tests_dbl_mant_bits (void)
  422. {
  423.   static int n = -1;
  424.   volatile double x, y, d;
  425.   if (n != -1)
  426.     return n;
  427.   n = 1;
  428.   x = 2.0;
  429.   for (;;)
  430.     {
  431.       /* see if 2^(n+1)+1 can be formed without rounding, if so then
  432.          continue, if not then "n" is the answer */
  433.       y = x + 1.0;
  434.       d = y - x;
  435.       if (d != 1.0)
  436.         {
  437. #if defined (DBL_MANT_DIG) && DBL_RADIX == 2
  438.           if (n != DBL_MANT_DIG)
  439.             printf ("Warning, tests_dbl_mant_bits got %d but DBL_MANT_DIG says %dn", n, DBL_MANT_DIG);
  440. #endif
  441.           break;
  442.         }
  443.       x *= 2;
  444.       n++;
  445.       if (n > 1000)
  446.         {
  447.           printf ("Oops, tests_dbl_mant_bits can't determine mantissa sizen");
  448.           n = 0;
  449.           break;
  450.         }
  451.     }
  452.   return n;
  453. }
  454. /* See tests_setjmp_sigfpe in tests.h. */
  455. jmp_buf    tests_sigfpe_target;
  456. RETSIGTYPE
  457. tests_sigfpe_handler (int sig)
  458. {
  459.   longjmp (tests_sigfpe_target, 1);
  460. }
  461. void
  462. tests_sigfpe_done (void)
  463. {
  464.   signal (SIGFPE, SIG_DFL);
  465. }