tests/mpn/t-sqrmod_bnm1.c

   1 /* Test for sqrmod_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 12
  33 #endif
  34
  35 #ifndef COUNT
  36 #define COUNT 3000
  37 #endif
  38
  39 #define MAX_N (1L << SIZE_LOG)
  40 #define MIN_N 1
  41
  42 /*
  43   Reference function for squaring 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 sqrmod_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_sqrmod_bnm1 (mp_ptr rp, mp_size_t rn, mp_srcptr ap, mp_size_t an)
  54 {
  55   mp_limb_t cy;
  56
  57   ASSERT (0 < an && an <= rn);
  58
  59   refmpn_mul (rp, ap, an, ap, an);
  60   an *= 2;
  61   if (an > rn) {
  62     cy = mpn_add (rp, rp, rn, rp + rn, an - rn);
  63     /* If cy == 1, then the value of rp is at most B^rn - 2, so there can
  64      * be no overflow when adding in the carry. */
  65     MPN_INCR_U (rp, rn, cy);
  66   }
  67 }
  68
  69 /*
  70   Compare the result of the mpn_sqrmod_bnm1 function in the library
  71   with the reference function above.
  72 */
  73
  74 int
  75 main (int argc, char **argv)
  76 {
  77   mp_ptr ap, refp, pp, scratch;
  78   int count = COUNT;
  79   int test;
  80   gmp_randstate_ptr rands;
  81   TMP_DECL;
  82   TMP_MARK;
  83
  84   if (argc > 1)
  85     {
  86       char *end;
  87       count = strtol (argv[1], &end, 0);
  88       if (*end || count <= 0)
  89         {
  90           fprintf (stderr, "Invalid test count: %s.\n", argv[1]);
  91           return 1;
  92         }
  93     }
  94
  95   tests_start ();
  96   rands = RANDS;
  97
  98   ASSERT_ALWAYS (mpn_sqrmod_bnm1_next_size (MAX_N) == MAX_N);
  99
 100   ap = TMP_ALLOC_LIMBS (MAX_N);
 101   refp = TMP_ALLOC_LIMBS (MAX_N * 4);
 102   pp = 1+TMP_ALLOC_LIMBS (MAX_N + 2);
 103   scratch
 104     = 1+TMP_ALLOC_LIMBS (mpn_sqrmod_bnm1_itch (MAX_N, MAX_N) + 2);
 105
 106   for (test = 0; test < count; test++)
 107     {
 108       unsigned size_min;
 109       unsigned size_range;
 110       mp_size_t an,rn,n;
 111       mp_size_t itch;
 112       mp_limb_t p_before, p_after, s_before, s_after;
 113
 114       for (size_min = 1; (1L << size_min) < MIN_N; size_min++)
 115         ;
 116
 117       /* We generate an in the MIN_N <= n <= (1 << size_range). */
 118       size_range = size_min
 119         + gmp_urandomm_ui (rands, SIZE_LOG + 1 - size_min);
 120
 121       n = MIN_N
 122         + gmp_urandomm_ui (rands, (1L << size_range) + 1 - MIN_N);
 123       n = mpn_sqrmod_bnm1_next_size (n);
 124
 125       if (n == 1)
 126         an = 1;
 127       else
 128         an = ((n+1) >> 1) + gmp_urandomm_ui (rands, (n+1) >> 1);
 129
 130       mpn_random2 (ap, an);
 131
 132       /* Sometime trigger the borderline conditions
 133          A = -1,0,+1 Mod(B^{n/2}+1).
 134          This only makes sense if there is at least a split, i.e. n is even. */
 135       if ((test & 0x1f) == 1 && (n & 1) == 0) {
 136         mp_size_t x;
 137         MPN_COPY (ap, ap + (n >> 1), an - (n >> 1));
 138         MPN_ZERO (ap + an - (n >> 1) , n - an);
 139         x = (n == an) ? 0 : gmp_urandomm_ui (rands, n - an);
 140         ap[x] += gmp_urandomm_ui (rands, 3) - 1;
 141       }
 142       rn = MIN(n, 2*an);
 143       mpn_random2 (pp-1, rn + 2);
 144       p_before = pp[-1];
 145       p_after = pp[rn];
 146
 147       itch = mpn_sqrmod_bnm1_itch (n, an);
 148       ASSERT_ALWAYS (itch <= mpn_sqrmod_bnm1_itch (MAX_N, MAX_N));
 149       mpn_random2 (scratch-1, itch+2);
 150       s_before = scratch[-1];
 151       s_after = scratch[itch];
 152
 153       mpn_sqrmod_bnm1 (  pp, n, ap, an, scratch);
 154       ref_sqrmod_bnm1 (refp, n, ap, an);
 155       if (pp[-1] != p_before || pp[rn] != p_after
 156           || scratch[-1] != s_before || scratch[itch] != s_after
 157           || mpn_cmp (refp, pp, rn) != 0)
 158         {
 159           printf ("ERROR in test %d, an = %d, n = %d\n",
 160                   test, (int) an, (int) n);
 161           if (pp[-1] != p_before)
 162             {
 163               printf ("before pp:"); mpn_dump (pp -1, 1);
 164               printf ("keep:   "); mpn_dump (&p_before, 1);
 165             }
 166           if (pp[rn] != p_after)
 167             {
 168               printf ("after pp:"); mpn_dump (pp + rn, 1);
 169               printf ("keep:   "); mpn_dump (&p_after, 1);
 170             }
 171           if (scratch[-1] != s_before)
 172             {
 173               printf ("before scratch:"); mpn_dump (scratch-1, 1);
 174               printf ("keep:   "); mpn_dump (&s_before, 1);
 175             }
 176           if (scratch[itch] != s_after)
 177             {
 178               printf ("after scratch:"); mpn_dump (scratch + itch, 1);
 179               printf ("keep:   "); mpn_dump (&s_after, 1);
 180             }
 181           mpn_dump (ap, an);
 182           mpn_dump (pp, rn);
 183           mpn_dump (refp, rn);
 184
 185           abort();
 186         }
 187     }
 188   TMP_FREE;
 189   tests_end ();
 190   return 0;
 191 }