1/* Test for sqrmod_bnm1 function.
2
3   Contributed to the GNU project by Marco Bodrato.
4
5Copyright 2009 Free Software Foundation, Inc.
6
7This file is part of the GNU MP Library.
8
9The GNU MP Library is free software; you can redistribute it and/or modify
10it under the terms of the GNU Lesser General Public License as published by
11the Free Software Foundation; either version 3 of the License, or (at your
12option) any later version.
13
14The GNU MP Library is distributed in the hope that it will be useful, but
15WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
17License for more details.
18
19You should have received a copy of the GNU Lesser General Public License
20along 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
52static void
53ref_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
74int
75main (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}
192