tests/mpz/t-div_2exp.c

   1 /* Test mpz_[cft]div_[qr]_2exp.
   2
   3 Copyright 2001 Free Software Foundation, Inc.
   4
   5 This file is part of the GNU MP Library.
   6
   7 The GNU MP Library is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU Lesser General Public License as published by
   9 the Free Software Foundation; either version 3 of the License, or (at your
  10 option) any later version.
  11
  12 The GNU MP Library is distributed in the hope that it will be useful, but
  13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  15 License for more details.
  16
  17 You should have received a copy of the GNU Lesser General Public License
  18 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  19
  20 #include <stdio.h>
  21 #include <stdlib.h>
  22
  23 #include "gmp.h"
  24 #include "gmp-impl.h"
  25 #include "tests.h"
  26
  27
  28 /* If the remainder is in the correct range and q*d+r is correct, then q
  29    must have rounded correctly.  */
  30
  31 void
  32 check_one (mpz_srcptr a, unsigned long d)
  33 {
  34   mpz_t  q, r, p, d2exp;
  35   int    inplace;
  36
  37   mpz_init (d2exp);
  38   mpz_init (q);
  39   mpz_init (r);
  40   mpz_init (p);
  41
  42   mpz_set_ui (d2exp, 1L);
  43   mpz_mul_2exp (d2exp, d2exp, d);
  44
  45 #define INPLACE(fun,dst,src,d)  \
  46   if (inplace)                  \
  47     {                           \
  48       mpz_set (dst, src);       \
  49       fun (dst, dst, d);        \
  50     }                           \
  51   else                          \
  52     fun (dst, src, d);
  53
  54   for (inplace = 0; inplace <= 1; inplace++)
  55     {
  56       INPLACE (mpz_fdiv_q_2exp, q, a, d);
  57       INPLACE (mpz_fdiv_r_2exp, r, a, d);
  58
  59       mpz_mul_2exp (p, q, d);
  60       mpz_add (p, p, r);
  61       if (mpz_sgn (r) < 0 || mpz_cmp (r, d2exp) >= 0)
  62         {
  63           printf ("mpz_fdiv_r_2exp result out of range\n");
  64           goto error;
  65         }
  66       if (mpz_cmp (p, a) != 0)
  67         {
  68           printf ("mpz_fdiv_[qr]_2exp doesn't multiply back\n");
  69           goto error;
  70         }
  71
  72
  73       INPLACE (mpz_cdiv_q_2exp, q, a, d);
  74       INPLACE (mpz_cdiv_r_2exp, r, a, d);
  75
  76       mpz_mul_2exp (p, q, d);
  77       mpz_add (p, p, r);
  78       if (mpz_sgn (r) > 0 || mpz_cmpabs (r, d2exp) >= 0)
  79         {
  80           printf ("mpz_cdiv_r_2exp result out of range\n");
  81           goto error;
  82         }
  83       if (mpz_cmp (p, a) != 0)
  84         {
  85           printf ("mpz_cdiv_[qr]_2exp doesn't multiply back\n");
  86           goto error;
  87         }
  88
  89
  90       INPLACE (mpz_tdiv_q_2exp, q, a, d);
  91       INPLACE (mpz_tdiv_r_2exp, r, a, d);
  92
  93       mpz_mul_2exp (p, q, d);
  94       mpz_add (p, p, r);
  95       if (mpz_sgn (r) != 0 && mpz_sgn (r) != mpz_sgn (a))
  96         {
  97           printf ("mpz_tdiv_r_2exp result wrong sign\n");
  98           goto error;
  99         }
 100       if (mpz_cmpabs (r, d2exp) >= 0)
 101         {
 102           printf ("mpz_tdiv_r_2exp result out of range\n");
 103           goto error;
 104         }
 105       if (mpz_cmp (p, a) != 0)
 106         {
 107           printf ("mpz_tdiv_[qr]_2exp doesn't multiply back\n");
 108           goto error;
 109         }
 110     }
 111
 112   mpz_clear (d2exp);
 113   mpz_clear (q);
 114   mpz_clear (r);
 115   mpz_clear (p);
 116   return;
 117
 118
 119  error:
 120   mpz_trace ("a", a);
 121   printf    ("d=%lu\n", d);
 122   mpz_trace ("q", q);
 123   mpz_trace ("r", r);
 124   mpz_trace ("p", p);
 125
 126   mp_trace_base = -16;
 127   mpz_trace ("a", a);
 128   printf    ("d=0x%lX\n", d);
 129   mpz_trace ("q", q);
 130   mpz_trace ("r", r);
 131   mpz_trace ("p", p);
 132
 133   abort ();
 134 }
 135
 136
 137 void
 138 check_all (mpz_ptr a, unsigned long d)
 139 {
 140   check_one (a, d);
 141   mpz_neg (a, a);
 142   check_one (a, d);
 143 }
 144
 145
 146 void
 147 check_various (void)
 148 {
 149   static const unsigned long  table[] = {
 150     0, 1, 2, 3, 4, 5,
 151     GMP_NUMB_BITS-1, GMP_NUMB_BITS, GMP_NUMB_BITS+1,
 152     2*GMP_NUMB_BITS-1, 2*GMP_NUMB_BITS, 2*GMP_NUMB_BITS+1,
 153     3*GMP_NUMB_BITS-1, 3*GMP_NUMB_BITS, 3*GMP_NUMB_BITS+1,
 154     4*GMP_NUMB_BITS-1, 4*GMP_NUMB_BITS, 4*GMP_NUMB_BITS+1
 155   };
 156
 157   int            i, j;
 158   unsigned long  n, d;
 159   mpz_t          a;
 160
 161   mpz_init (a);
 162
 163   /* a==0, and various d */
 164   mpz_set_ui (a, 0L);
 165   for (i = 0; i < numberof (table); i++)
 166     check_one (a, table[i]);
 167
 168   /* a==2^n, and various d */
 169   for (i = 0; i < numberof (table); i++)
 170     {
 171       n = table[i];
 172       mpz_set_ui (a, 1L);
 173       mpz_mul_2exp (a, a, n);
 174
 175       for (j = 0; j < numberof (table); j++)
 176         {
 177           d = table[j];
 178           check_all (a, d);
 179         }
 180     }
 181
 182   mpz_clear (a);
 183 }
 184
 185
 186 void
 187 check_random (int argc, char *argv[])
 188 {
 189   gmp_randstate_ptr  rands = RANDS;
 190   int            reps = 100;
 191   mpz_t          a;
 192   unsigned long  d;
 193   int            i;
 194
 195   if (argc == 2)
 196     reps = atoi (argv[1]);
 197
 198   mpz_init (a);
 199
 200   for (i = 0; i < reps; i++)
 201     {
 202       /* exponentially within 2 to 257 bits */
 203       mpz_erandomb (a, rands, urandom () % 8 + 2);
 204
 205       d = urandom () % 256;
 206
 207       check_all (a, d);
 208     }
 209
 210   mpz_clear (a);
 211 }
 212
 213
 214 int
 215 main (int argc, char *argv[])
 216 {
 217   tests_start ();
 218
 219   check_various ();
 220   check_random (argc, argv);
 221
 222   tests_end ();
 223   exit (0);
 224 }