xref: /netbsd-src/external/lgpl3/gmp/dist/tests/mpz/t-powm.c (revision 70f7362772ba52b749c976fb5e86e39a8b2c9afc)
1 /* Test mpz_powm, mpz_mul, mpz_mod, mpz_mod_ui, mpz_div_ui.
2 
3 Copyright 1991, 1993, 1994, 1996, 1999-2001, 2009, 2012, 2019 Free
4 Software Foundation, Inc.
5 
6 This file is part of the GNU MP Library test suite.
7 
8 The GNU MP Library test suite is free software; you can redistribute it
9 and/or modify it under the terms of the GNU General Public License as
10 published by the Free Software Foundation; either version 3 of the License,
11 or (at your option) any later version.
12 
13 The GNU MP Library test suite is distributed in the hope that it will be
14 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
16 Public License for more details.
17 
18 You should have received a copy of the GNU General Public License along with
19 the GNU MP Library test suite.  If not, see https://www.gnu.org/licenses/.  */
20 
21 #include <stdio.h>
22 #include <stdlib.h>
23 #include <string.h>
24 
25 #include "gmp-impl.h"
26 #include "tests.h"
27 
28 void debug_mp (mpz_t, int);
29 
30 #define SIZEM 13
31 
32 /* Check that all sizes up to just above MUL_TOOM22_THRESHOLD have been tested
33    a few times.  FIXME: If SIZEM is set too low, this will never happen.  */
34 int
35 allsizes_seen (unsigned int *allsizes)
36 {
37   mp_size_t i;
38 
39   for (i = 1; i < MUL_TOOM22_THRESHOLD + 4; i++)
40     if (allsizes[i] < 4)
41       return 0;
42   return 1;
43 }
44 
45 int
46 main (int argc, char **argv)
47 {
48   mpz_t base, exp, mod;
49   mpz_t r1, r2, t1, exp2, base2;
50   mp_size_t base_size, exp_size, mod_size;
51   int i;
52   int reps = 1000;
53   gmp_randstate_ptr rands;
54   mpz_t bs;
55   unsigned long bsi, size_range;
56   unsigned int allsizes[1 << (SIZEM + 2 - 1)];
57 
58   tests_start ();
59   TESTS_REPS (reps, argv, argc);
60 
61   rands = RANDS;
62 
63   mpz_init (bs);
64 
65   mpz_init (base);
66   mpz_init (exp);
67   mpz_init (mod);
68   mpz_init (r1);
69   mpz_init (r2);
70   mpz_init (t1);
71   mpz_init (exp2);
72   mpz_init (base2);
73 
74   memset (allsizes, 0, (1 << (SIZEM + 2 - 1)) * sizeof (int));
75 
76   reps += reps >> 3;
77   for (i = 0; i < reps || ! allsizes_seen (allsizes); i++)
78     {
79       mpz_urandomb (bs, rands, 32);
80       size_range = mpz_get_ui (bs) % SIZEM + 2;
81 
82       if ((i & 7) == 0)
83 	{
84 	  mpz_set_ui (exp, 1);
85 
86 	  do  /* Loop until mathematically well-defined.  */
87 	    {
88 	      mpz_urandomb (bs, rands, size_range / 2 + 2);
89 	      base_size = mpz_get_ui (bs);
90 	      mpz_rrandomb (base, rands, base_size);
91 	    }
92 	  while (mpz_cmp_ui (base, 0) == 0);
93 
94 	  mpz_urandomb (bs, rands, size_range / 2);
95 	  mod_size = mpz_get_ui (bs);
96 	  mod_size = MIN (mod_size, base_size);
97 	  mpz_rrandomb (mod, rands, mod_size);
98 
99 	  mpz_urandomb (bs, rands, size_range);
100 	  mod_size = mpz_get_ui (bs) + base_size + 2;
101 	  if ((i & 8) == 0)
102 	    mod_size += (GMP_NUMB_BITS - mod_size) % GMP_NUMB_BITS;
103 	  mpz_setbit (mod, mod_size);
104 
105 	  mpz_sub (base, base, mod);
106 	}
107       else
108 	{
109       do  /* Loop until mathematically well-defined.  */
110 	{
111 	  mpz_urandomb (bs, rands, size_range);
112 	  base_size = mpz_get_ui (bs);
113 	  mpz_rrandomb (base, rands, base_size);
114 
115 	  mpz_urandomb (bs, rands, 7L);
116 	  exp_size = mpz_get_ui (bs);
117 	  mpz_rrandomb (exp, rands, exp_size);
118 	}
119       while (mpz_cmp_ui (base, 0) == 0 && mpz_cmp_ui (exp, 0) == 0);
120 
121       do
122         {
123 	  mpz_urandomb (bs, rands, size_range);
124 	  mod_size = mpz_get_ui (bs);
125 	  mpz_rrandomb (mod, rands, mod_size);
126 	}
127       while (mpz_cmp_ui (mod, 0) == 0);
128 
129       allsizes[SIZ(mod)] += 1;
130 
131       mpz_urandomb (bs, rands, 2);
132       bsi = mpz_get_ui (bs);
133       if ((bsi & 1) != 0)
134 	mpz_neg (base, base);
135 
136       /* printf ("%ld %ld %ld\n", SIZ (base), SIZ (exp), SIZ (mod)); */
137 	}
138 
139       mpz_set_ui (r2, 1);
140       mpz_mod (base2, base, mod);
141       mpz_set (exp2, exp);
142       mpz_mod (r2, r2, mod);
143 
144       for (;;)
145 	{
146 	  if (mpz_tstbit (exp2, 0))
147 	    {
148 	      mpz_mul (r2, r2, base2);
149 	      mpz_mod (r2, r2, mod);
150 	    }
151 	  if  (mpz_cmp_ui (exp2, 1) <= 0)
152 	    break;
153 	  mpz_mul (base2, base2, base2);
154 	  mpz_mod (base2, base2, mod);
155 	  mpz_tdiv_q_2exp (exp2, exp2, 1);
156 	}
157 
158       mpz_powm (r1, base, exp, mod);
159       MPZ_CHECK_FORMAT (r1);
160 
161       if (mpz_cmp (r1, r2) != 0)
162 	{
163 	  fprintf (stderr, "\nIncorrect results in test %d for operands:\n", i);
164 	  debug_mp (base, -16);
165 	  debug_mp (exp, -16);
166 	  debug_mp (mod, -16);
167 	  fprintf (stderr, "mpz_powm result:\n");
168 	  debug_mp (r1, -16);
169 	  fprintf (stderr, "reference result:\n");
170 	  debug_mp (r2, -16);
171 	  abort ();
172 	}
173 
174       if (mpz_tdiv_ui (mod, 2) == 0)
175 	continue;
176 
177       mpz_powm_sec (r1, base, exp, mod);
178       MPZ_CHECK_FORMAT (r1);
179 
180       if (mpz_cmp (r1, r2) != 0)
181 	{
182 	  fprintf (stderr, "\nIncorrect results in test %d for operands:\n", i);
183 	  debug_mp (base, -16);
184 	  debug_mp (exp, -16);
185 	  debug_mp (mod, -16);
186 	  fprintf (stderr, "mpz_powm_sec result:\n");
187 	  debug_mp (r1, -16);
188 	  fprintf (stderr, "reference result:\n");
189 	  debug_mp (r2, -16);
190 	  abort ();
191 	}
192     }
193 
194   mpz_clear (bs);
195   mpz_clear (base);
196   mpz_clear (exp);
197   mpz_clear (mod);
198   mpz_clear (r1);
199   mpz_clear (r2);
200   mpz_clear (t1);
201   mpz_clear (exp2);
202   mpz_clear (base2);
203 
204   tests_end ();
205   exit (0);
206 }
207 
208 void
209 debug_mp (mpz_t x, int base)
210 {
211   mpz_out_str (stderr, base, x); fputc ('\n', stderr);
212 }
213