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

数学计算

开发平台:

Unix_Linux

  1. /* Test mpz_[cft]div_[qr]_2exp.
  2. Copyright 2001 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 <stdio.h>
  15. #include <stdlib.h>
  16. #include "gmp.h"
  17. #include "gmp-impl.h"
  18. #include "tests.h"
  19. /* If the remainder is in the correct range and q*d+r is correct, then q
  20.    must have rounded correctly.  */
  21. void
  22. check_one (mpz_srcptr a, unsigned long d)
  23. {
  24.   mpz_t  q, r, p, d2exp;
  25.   int    inplace;
  26.   mpz_init (d2exp);
  27.   mpz_init (q);
  28.   mpz_init (r);
  29.   mpz_init (p);
  30.   mpz_set_ui (d2exp, 1L);
  31.   mpz_mul_2exp (d2exp, d2exp, d);
  32. #define INPLACE(fun,dst,src,d)  
  33.   if (inplace)                  
  34.     {                           
  35.       mpz_set (dst, src);       
  36.       fun (dst, dst, d);        
  37.     }                           
  38.   else                          
  39.     fun (dst, src, d);
  40.   for (inplace = 0; inplace <= 1; inplace++)
  41.     {
  42.       INPLACE (mpz_fdiv_q_2exp, q, a, d);
  43.       INPLACE (mpz_fdiv_r_2exp, r, a, d);
  44.       mpz_mul_2exp (p, q, d);
  45.       mpz_add (p, p, r);
  46.       if (mpz_sgn (r) < 0 || mpz_cmp (r, d2exp) >= 0)
  47. {
  48.   printf ("mpz_fdiv_r_2exp result out of rangen");
  49.   goto error;
  50. }
  51.       if (mpz_cmp (p, a) != 0)
  52. {
  53.   printf ("mpz_fdiv_[qr]_2exp doesn't multiply backn");
  54.   goto error;
  55. }
  56.       INPLACE (mpz_cdiv_q_2exp, q, a, d);
  57.       INPLACE (mpz_cdiv_r_2exp, r, a, d);
  58.       mpz_mul_2exp (p, q, d);
  59.       mpz_add (p, p, r);
  60.       if (mpz_sgn (r) > 0 || mpz_cmpabs (r, d2exp) >= 0)
  61. {
  62.   printf ("mpz_cdiv_r_2exp result out of rangen");
  63.   goto error;
  64. }
  65.       if (mpz_cmp (p, a) != 0)
  66. {
  67.   printf ("mpz_cdiv_[qr]_2exp doesn't multiply backn");
  68.   goto error;
  69. }
  70.       INPLACE (mpz_tdiv_q_2exp, q, a, d);
  71.       INPLACE (mpz_tdiv_r_2exp, r, a, d);
  72.       mpz_mul_2exp (p, q, d);
  73.       mpz_add (p, p, r);
  74.       if (mpz_sgn (r) != 0 && mpz_sgn (r) != mpz_sgn (a))
  75. {
  76.   printf ("mpz_tdiv_r_2exp result wrong signn");
  77.   goto error;
  78. }
  79.       if (mpz_cmpabs (r, d2exp) >= 0)
  80. {
  81.   printf ("mpz_tdiv_r_2exp result out of rangen");
  82.   goto error;
  83. }
  84.       if (mpz_cmp (p, a) != 0)
  85. {
  86.   printf ("mpz_tdiv_[qr]_2exp doesn't multiply backn");
  87.   goto error;
  88. }
  89.     }
  90.   mpz_clear (d2exp);
  91.   mpz_clear (q);
  92.   mpz_clear (r);
  93.   mpz_clear (p);
  94.   return;
  95.  error:
  96.   mpz_trace ("a", a);
  97.   printf    ("d=%lun", d);
  98.   mpz_trace ("q", q);
  99.   mpz_trace ("r", r);
  100.   mpz_trace ("p", p);
  101.   mp_trace_base = -16;
  102.   mpz_trace ("a", a);
  103.   printf    ("d=0x%lXn", d);
  104.   mpz_trace ("q", q);
  105.   mpz_trace ("r", r);
  106.   mpz_trace ("p", p);
  107.   abort ();
  108. }
  109. void
  110. check_all (mpz_ptr a, unsigned long d)
  111. {
  112.   check_one (a, d);
  113.   mpz_neg (a, a);
  114.   check_one (a, d);
  115. }
  116. void
  117. check_various (void)
  118. {
  119.   static const unsigned long  table[] = {
  120.     0, 1, 2, 3, 4, 5,
  121.     GMP_NUMB_BITS-1, GMP_NUMB_BITS, GMP_NUMB_BITS+1,
  122.     2*GMP_NUMB_BITS-1, 2*GMP_NUMB_BITS, 2*GMP_NUMB_BITS+1,
  123.     3*GMP_NUMB_BITS-1, 3*GMP_NUMB_BITS, 3*GMP_NUMB_BITS+1,
  124.     4*GMP_NUMB_BITS-1, 4*GMP_NUMB_BITS, 4*GMP_NUMB_BITS+1
  125.   };
  126.   int            i, j;
  127.   unsigned long  n, d;
  128.   mpz_t          a;
  129.   mpz_init (a);
  130.   /* a==0, and various d */
  131.   mpz_set_ui (a, 0L);
  132.   for (i = 0; i < numberof (table); i++)
  133.     check_one (a, table[i]);
  134.   /* a==2^n, and various d */
  135.   for (i = 0; i < numberof (table); i++)
  136.     {
  137.       n = table[i];
  138.       mpz_set_ui (a, 1L);
  139.       mpz_mul_2exp (a, a, n);
  140.       for (j = 0; j < numberof (table); j++)
  141. {
  142.   d = table[j];
  143.   check_all (a, d);
  144. }
  145.     }
  146.   mpz_clear (a);
  147. }
  148. void
  149. check_random (int argc, char *argv[])
  150. {
  151.   gmp_randstate_ptr  rands = RANDS;
  152.   int            reps = 100;
  153.   mpz_t          a;
  154.   unsigned long  d;
  155.   int            i;
  156.   if (argc == 2)
  157.     reps = atoi (argv[1]);
  158.   mpz_init (a);
  159.   for (i = 0; i < reps; i++)
  160.     {
  161.       /* exponentially within 2 to 257 bits */
  162.       mpz_erandomb (a, rands, urandom () % 8 + 2);
  163.       d = urandom () % 256;
  164.       check_all (a, d);
  165.     }
  166.   mpz_clear (a);
  167. }
  168. int
  169. main (int argc, char *argv[])
  170. {
  171.   tests_start ();
  172.   check_various ();
  173.   check_random (argc, argv);
  174.   tests_end ();
  175.   exit (0);
  176. }