tests/mpn/t-mulmod_bnm1.c

   1 /* Test for mulmod_bnm1 function.
   2
   3    Contributed to the GNU project by Marco Bodrato.
   4
   5 Copyright 2009 Free Software Foundation, Inc.
   6
   7 This file is part of the GNU MP Library.
   8
   9 The GNU MP Library is free software; you can redistribute it and/or modify
  10 it under the terms of the GNU Lesser General Public License as published by
  11 the Free Software Foundation; either version 3 of the License, or (at your
  12 option) any later version.
  13
  14 The GNU MP Library is distributed in the hope that it will be useful, but
  15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  16 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  17 License for more details.
  18
  19 You should have received a copy of the GNU Lesser General Public License
  20 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  21
  22
  23 #include "gmp.h"
  24 #include "gmp-impl.h"
  25 #include "tests.h"
  26
  27 #include <stdlib.h>
  28 #include <stdio.h>
  29
  30 /* Sizes are up to 2^SIZE_LOG limbs */
  31 #ifndef SIZE_LOG
  32 #define SIZE_LOG 11
  33 #endif
  34
  35 #ifndef COUNT
  36 #define COUNT 5000
  37 #endif
  38
  39 #define MAX_N (1L << SIZE_LOG)
  40 #define MIN_N 1
  41
  42 /*
  43   Reference function for multiplication modulo B^rn-1.
  44
  45   The result is expected to be ZERO if and only if one of the operand
  46   already is. Otherwise the class [0] Mod(B^rn-1) is represented by
  47   B^rn-1. This should not be a problem if mulmod_bnm1 is used to
  48   combine results and obtain a natural number when one knows in
  49   advance that the final value is less than (B^rn-1).
  50 */
  51
  52 static void
  53 ref_mulmod_bnm1 (mp_ptr rp, mp_size_t rn, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
  54 {
  55   mp_limb_t cy;
  56
  57   ASSERT (0 < an && an <= rn);
  58   ASSERT (0 < bn && bn <= rn);
  59
  60   if (an >= bn)
  61     refmpn_mul (rp, ap, an, bp, bn);
  62   else
  63     refmpn_mul (rp, bp, bn, ap, an);
  64   an += bn;
  65   if (an > rn) {
  66     cy = mpn_add (rp, rp, rn, rp + rn, an - rn);
  67     /* If cy == 1, then the value of rp is at most B^rn - 2, so there can
  68      * be no overflow when adding in the carry. */
  69     MPN_INCR_U (rp, rn, cy);
  70   }
  71 }
  72
  73 /*
  74   Compare the result of the mpn_mulmod_bnm1 function in the library
  75   with the reference function above.
  76 */
  77
  78 int
  79 main (int argc, char **argv)
  80 {
  81   mp_ptr ap, bp, refp, pp, scratch;
  82   int count = COUNT;
  83   int test;
  84   gmp_randstate_ptr rands;
  85   TMP_DECL;
  86   TMP_MARK;
  87
  88   if (argc > 1)
  89     {
  90       char *end;
  91       count = strtol (argv[1], &end, 0);
  92       if (*end || count <= 0)
  93         {
  94           fprintf (stderr, "Invalid test count: %s.\n", argv[1]);
  95           return 1;
  96         }
  97     }
  98
  99   tests_start ();
 100   rands = RANDS;
 101
 102   ASSERT_ALWAYS (mpn_mulmod_bnm1_next_size (MAX_N) == MAX_N);
 103
 104   ap = TMP_ALLOC_LIMBS (MAX_N);
 105   bp = TMP_ALLOC_LIMBS (MAX_N);
 106   refp = TMP_ALLOC_LIMBS (MAX_N * 4);
 107   pp = 1+TMP_ALLOC_LIMBS (MAX_N + 2);
 108   scratch
 109     = 1+TMP_ALLOC_LIMBS (mpn_mulmod_bnm1_itch (MAX_N, MAX_N, MAX_N) + 2);
 110
 111   for (test = 0; test < count; test++)
 112     {
 113       unsigned size_min;
 114       unsigned size_range;
 115       mp_size_t an,bn,rn,n;
 116       mp_size_t itch;
 117       mp_limb_t p_before, p_after, s_before, s_after;
 118
 119       for (size_min = 1; (1L << size_min) < MIN_N; size_min++)
 120         ;
 121
 122       /* We generate an in the MIN_N <= n <= (1 << size_range). */
 123       size_range = size_min
 124         + gmp_urandomm_ui (rands, SIZE_LOG + 1 - size_min);
 125
 126       n = MIN_N
 127         + gmp_urandomm_ui (rands, (1L << size_range) + 1 - MIN_N);
 128       n = mpn_mulmod_bnm1_next_size (n);
 129
 130       if ( (test & 1) || n == 1) {
 131         /* Half of the tests are done with the main scenario in mind:
 132            both an and bn >= rn/2 */
 133         an = ((n+1) >> 1) + gmp_urandomm_ui (rands, (n+1) >> 1);
 134         bn = ((n+1) >> 1) + gmp_urandomm_ui (rands, (n+1) >> 1);
 135       } else {
 136         /* Second half of the tests are done using mulmod to compute a
 137            full product with n/2 < an+bn <= n. */
 138         an = 1 + gmp_urandomm_ui (rands, n - 1);
 139         if (an >= n/2)
 140           bn = 1 + gmp_urandomm_ui (rands, n - an);
 141         else
 142           bn = n/2 + 1 - an + gmp_urandomm_ui (rands, (n+1)/2);
 143       }
 144
 145       /* Make sure an >= bn */
 146       if (an < bn)
 147         MP_SIZE_T_SWAP (an, bn);
 148
 149       mpn_random2 (ap, an);
 150       mpn_random2 (bp, bn);
 151
 152       /* Sometime trigger the borderline conditions
 153          A = -1,0,+1 or B = -1,0,+1 or A*B == -1,0,1 Mod(B^{n/2}+1).
 154          This only makes sense if there is at least a split, i.e. n is even. */
 155       if ((test & 0x1f) == 1 && (n & 1) == 0) {
 156         mp_size_t x;
 157         MPN_COPY (ap, ap + (n >> 1), an - (n >> 1));
 158         MPN_ZERO (ap + an - (n >> 1) , n - an);
 159         MPN_COPY (bp, bp + (n >> 1), bn - (n >> 1));
 160         MPN_ZERO (bp + bn - (n >> 1) , n - bn);
 161         x = (n == an) ? 0 : gmp_urandomm_ui (rands, n - an);
 162         ap[x] += gmp_urandomm_ui (rands, 3) - 1;
 163         x = (n >> 1) - x % (n >> 1);
 164         bp[x] += gmp_urandomm_ui (rands, 3) - 1;
 165         /* We don't propagate carry, this means that the desired condition
 166            is not triggered all the times. A few times are enough anyway. */
 167       }
 168       rn = MIN(n, an + bn);
 169       mpn_random2 (pp-1, rn + 2);
 170       p_before = pp[-1];
 171       p_after = pp[rn];
 172
 173       itch = mpn_mulmod_bnm1_itch (n, an, bn);
 174       ASSERT_ALWAYS (itch <= mpn_mulmod_bnm1_itch (MAX_N, MAX_N, MAX_N));
 175       mpn_random2 (scratch-1, itch+2);
 176       s_before = scratch[-1];
 177       s_after = scratch[itch];
 178
 179       mpn_mulmod_bnm1 (  pp, n, ap, an, bp, bn, scratch);
 180       ref_mulmod_bnm1 (refp, n, ap, an, bp, bn);
 181       if (pp[-1] != p_before || pp[rn] != p_after
 182           || scratch[-1] != s_before || scratch[itch] != s_after
 183           || mpn_cmp (refp, pp, rn) != 0)
 184         {
 185           printf ("ERROR in test %d, an = %d, bn = %d, n = %d\n",
 186                   test, (int) an, (int) bn, (int) n);
 187           if (pp[-1] != p_before)
 188             {
 189               printf ("before pp:"); mpn_dump (pp -1, 1);
 190               printf ("keep:   "); mpn_dump (&p_before, 1);
 191             }
 192           if (pp[rn] != p_after)
 193             {
 194               printf ("after pp:"); mpn_dump (pp + rn, 1);
 195               printf ("keep:   "); mpn_dump (&p_after, 1);
 196             }
 197           if (scratch[-1] != s_before)
 198             {
 199               printf ("before scratch:"); mpn_dump (scratch-1, 1);
 200               printf ("keep:   "); mpn_dump (&s_before, 1);
 201             }
 202           if (scratch[itch] != s_after)
 203             {
 204               printf ("after scratch:"); mpn_dump (scratch + itch, 1);
 205               printf ("keep:   "); mpn_dump (&s_after, 1);
 206             }
 207           mpn_dump (ap, an);
 208           mpn_dump (bp, bn);
 209           mpn_dump (pp, rn);
 210           mpn_dump (refp, rn);
 211
 212           abort();
 213         }
 214     }
 215   TMP_FREE;
 216   tests_end ();
 217   return 0;
 218 }