1/* Test mpn_fib2m.
2
3Copyright 2018 Free Software Foundation, Inc.
4
5This file is part of the GNU MP Library test suite.
6
7The GNU MP Library test suite is free software; you can redistribute it
8and/or modify it under the terms of the GNU General Public License as
9published by the Free Software Foundation; either version 3 of the License,
10or (at your option) any later version.
11
12The GNU MP Library test suite is distributed in the hope that it will be
13useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
15Public License for more details.
16
17You should have received a copy of the GNU General Public License along with
18the GNU MP Library test suite.  If not, see https://www.gnu.org/licenses/.  */
19
20#include <stdio.h>
21#include <stdlib.h>
22
23#include "gmp-impl.h"
24#include "tests.h"
25
26#define MAX_K_BITS 16
27#define MAX_K (1L << MAX_K_BITS)
28#define MIN_K 1
29
30#define MAX_MN 20
31#define MAX_KN 30
32
33#define COUNT 200
34
35static int
36test_fib2_fib2m (int count, gmp_randstate_ptr rands)
37{
38  int test;
39  mp_ptr fk, fks1, fkm, fks1m, mp, qp;
40  mp_size_t mn, fn, size, max_mn;
41  TMP_DECL;
42
43  size = MPN_FIB2_SIZE (MAX_K);
44  max_mn = size / 4 + 10;
45  ASSERT (max_mn < size);
46
47  TMP_MARK;
48  fk	 = TMP_ALLOC_LIMBS (size);
49  fks1	 = TMP_ALLOC_LIMBS (size);
50  qp	 = TMP_ALLOC_LIMBS (size);
51  mp	 = TMP_ALLOC_LIMBS (max_mn);
52  fkm	 = 1 + TMP_ALLOC_LIMBS (max_mn * 2 + 1 + 2);
53  fks1m	 = 1 + TMP_ALLOC_LIMBS (max_mn * 2 + 1 + 2);
54
55  for (test = 1; test <= count; ++test)
56    {
57      mp_limb_t fk_before, fk_after, fk1_before, fk1_after;
58      int signflip;
59      unsigned long k;
60
61      k = MIN_K +
62	gmp_urandomm_ui (rands, test < MAX_K_BITS ?
63			 MAX_K >> test : (MAX_K - MIN_K));
64
65      fn = mpn_fib2_ui (fk, fks1, k);
66      do {
67	mn = gmp_urandomm_ui (rands, MAX_K) % (fn / 4 + 10);
68      } while (mn == 0);
69      ASSERT (mn <= max_mn);
70      mpn_random2 (mp, mn);
71      ASSERT (mp [mn - 1] != 0);
72
73      if (fn >= mn)
74	{
75	  mpn_tdiv_qr (qp, fk, 0, fk, fn, mp, mn);
76	  mpn_tdiv_qr (qp, fks1, 0, fks1, fn, mp, mn);
77	}
78      else
79	{
80	  MPN_ZERO (fk + fn, mn - fn);
81	  MPN_ZERO (fks1 + fn, mn - fn);
82	}
83
84      mpn_random2 (fkm - 1, 2*mn+1+2);
85      fk_before = fkm [-1];
86      fk_after = fkm [2 * mn + 1];
87
88      mpn_random2 (fks1m - 1, 2*mn+1+2);
89      fk1_before = fks1m [-1];
90      fk1_after = fks1m [2 * mn + 1];
91
92      qp [0] = k;
93      signflip = mpn_fib2m (fkm, fks1m, qp, 1, mp, mn);
94      if (fkm [-1] != fk_before || fkm [2 * mn + 1] != fk_after
95	  || fks1m [-1] != fk1_before || fks1m [2 * mn + 1] != fk1_after)
96	{
97	  printf ("REDZONE violation in test %d, k = %lu, mn = %u\n",
98		  test, k, (unsigned) mn);
99	  if (fkm[-1] != fk_before)
100	    {
101	      printf ("before fkm:"); mpn_dump (fkm - 1, 1);
102	      printf ("keep:   "); mpn_dump (&fk_before, 1);
103	    }
104	  if (fkm[2 * mn + 1] != fk_after)
105	    {
106	      printf ("after fkm:"); mpn_dump (fkm + 2 * mn + 1, 1);
107	      printf ("keep:   "); mpn_dump (&fk_after, 1);
108	    }
109	  if (fks1m[-1] != fk1_before)
110	    {
111	      printf ("before fks1m:"); mpn_dump (fks1m - 1, 1);
112	      printf ("keep:   "); mpn_dump (&fk1_before, 1);
113	    }
114	  if (fks1m[2 * mn + 1] != fk1_after)
115	    {
116	      printf ("after fks1m:"); mpn_dump (fks1m + 2 * mn + 1, 1);
117	      printf ("keep:   "); mpn_dump (&fk1_after, 1);
118	    }
119	  abort();
120	}
121
122      if (mpn_cmp (fkm, fk, mn) != 0)
123	{
124	  if (mpn_sub_n (fk, mp, fk, mn) || mpn_cmp (fkm, fk, mn) != 0)
125	    {
126	      printf ("ERROR(k) in test %d, k = %lu, mn = %u\n",
127		      test, k, (unsigned) mn);
128	      mpn_dump (fk, mn);
129	      mpn_dump (fkm, mn);
130	      mpn_dump (mp, mn);
131	      abort();
132	    }
133	  signflip ^= 1;
134	}
135
136      if (mpn_cmp (fks1m, fks1, mn) != 0)
137	{
138	  if (mpn_sub_n (fks1, mp, fks1, mn) || mpn_cmp (fks1m, fks1, mn) != 0)
139	    {
140	      printf ("ERROR(k-1) in test %d, k = %lu, mn = %u\n",
141		      test, k, (unsigned) mn);
142	      mpn_dump (fks1, mn);
143	      mpn_dump (fks1m, mn);
144	      mpn_dump (mp, mn);
145	      abort();
146	    }
147	  signflip ^= 1;
148	}
149
150      if (signflip != 0 && ! mpn_zero_p (fks1m, mn) && ! mpn_zero_p (fkm, mn))
151	{
152	  if ((mp [0] & 1) == 0) /* Should we test only odd modulus-es? */
153	    {
154	      if (! mpn_lshift (fks1m, fks1m, mn, 1) &&
155		  mpn_cmp (mp, fks1m, mn) == 0)
156		continue;
157	      if (! mpn_lshift (fkm, fkm, mn, 1) &&
158		  mpn_cmp (mp, fkm, mn) == 0)
159		continue;
160	    }
161	  printf ("ERROR(sign) in test %d, k = %lu, mn = %u\n",
162		  test, k, (unsigned) mn);
163	  abort();
164	}
165    }
166  TMP_FREE;
167  return 0;
168}
169
170static int
171test_fib2m_2exp (int count, gmp_randstate_ptr rands)
172{
173  int test;
174  mp_ptr fka, fks1a, fkb, fks1b, mp, kp;
175  TMP_DECL;
176
177  TMP_MARK;
178  kp	 = TMP_ALLOC_LIMBS (MAX_KN);
179  mp	 = TMP_ALLOC_LIMBS (MAX_MN);
180  fka	 = 1 + TMP_ALLOC_LIMBS (MAX_MN * 2 + 1 + 2);
181  fks1a	 = 1 + TMP_ALLOC_LIMBS (MAX_MN * 2 + 1 + 2);
182  fkb	 = 1 + TMP_ALLOC_LIMBS (MAX_MN * 2 + 1 + 2);
183  fks1b	 = 1 + TMP_ALLOC_LIMBS (MAX_MN * 2 + 1 + 2);
184
185  for (test = 1; test <= count; ++test)
186    {
187      mp_limb_t fka_before, fka_after, fk1a_before, fk1a_after;
188      mp_limb_t fkb_before, fkb_after, fk1b_before, fk1b_after;
189      mp_size_t mn, kn;
190      int signflip;
191      mp_bitcnt_t exp2;
192
193      mn = gmp_urandomm_ui (rands, MAX_MN - 1) + 1;
194      mpn_random2 (mp, mn);
195
196      exp2 = MIN_K + 1 + gmp_urandomm_ui (rands, MAX_KN * GMP_NUMB_BITS - MIN_K - 1);
197
198      kn = BITS_TO_LIMBS (exp2);
199      MPN_ZERO (kp, kn - 1);
200      kp [kn - 1] = CNST_LIMB (1) << ((exp2 - 1) % GMP_NUMB_BITS);
201
202      mpn_random2 (fka - 1, 2*mn+1+2);
203      fka_before = fka [-1];
204      fka_after = fka [2 * mn + 1];
205
206      mpn_random2 (fks1a - 1, 2*mn+1+2);
207      fk1a_before = fks1a [-1];
208      fk1a_after = fks1a [2 * mn + 1];
209
210      signflip = mpn_fib2m (fka, fks1a, kp, kn, mp, mn);
211      if (fka [-1] != fka_before || fka [2 * mn + 1] != fka_after
212	  || fks1a [-1] != fk1a_before || fks1a [2 * mn + 1] != fk1a_after)
213	{
214	  printf ("REDZONE(a) violation in test %d, exp2 = %lu\n", test, exp2);
215	  if (fka[-1] != fka_before)
216	    {
217	      printf ("before fka:"); mpn_dump (fka - 1, 1);
218	      printf ("keep:   "); mpn_dump (&fka_before, 1);
219	    }
220	  if (fka[2 * mn + 1] != fka_after)
221	    {
222	      printf ("after fka:"); mpn_dump (fka + 2 * mn + 1, 1);
223	      printf ("keep:   "); mpn_dump (&fka_after, 1);
224	    }
225	  if (fks1a[-1] != fk1a_before)
226	    {
227	      printf ("before fks1a:"); mpn_dump (fks1a - 1, 1);
228	      printf ("keep:   "); mpn_dump (&fk1a_before, 1);
229	    }
230	  if (fks1a[2 * mn + 1] != fk1a_after)
231	    {
232	      printf ("after fks1a:"); mpn_dump (fks1a + 2 * mn + 1, 1);
233	      printf ("keep:   "); mpn_dump (&fk1a_after, 1);
234	    }
235	  abort();
236	}
237
238      if (signflip && ! mpn_zero_p (fks1a, mn))
239	mpn_sub_n (fks1a, mp, fks1a, mn);
240      if (mpn_sub_n (fka, fka, fks1a, mn))
241	ASSERT_CARRY (mpn_add_n (fka, fka, mp, mn));
242
243      mpn_sub_1 (kp, kp, kn, 1);
244      ASSERT (exp2 % GMP_NUMB_BITS == 1 || kp [kn - 1] != 0);
245      kn -= kp [kn - 1] == 0;
246
247      mpn_random2 (fkb - 1, 2*mn+1+2);
248      fkb_before = fkb [-1];
249      fkb_after = fkb [2 * mn + 1];
250
251      mpn_random2 (fks1b - 1, 2*mn+1+2);
252      fk1b_before = fks1b [-1];
253      fk1b_after = fks1b [2 * mn + 1];
254
255      signflip = mpn_fib2m (fkb, fks1b, kp, kn, mp, mn);
256      if (fkb [-1] != fkb_before || fkb [2 * mn + 1] != fkb_after
257	  || fks1b [-1] != fk1b_before || fks1b [2 * mn + 1] != fk1b_after)
258	{
259	  printf ("REDZONE(b) violation in test %d, exp2 = %lu\n", test, exp2);
260	  if (fkb[-1] != fkb_before)
261	    {
262	      printf ("before fkb:"); mpn_dump (fkb - 1, 1);
263	      printf ("keep:   "); mpn_dump (&fkb_before, 1);
264	    }
265	  if (fkb[2 * mn + 1] != fkb_after)
266	    {
267	      printf ("after fkb:"); mpn_dump (fkb + 2 * mn + 1, 1);
268	      printf ("keep:   "); mpn_dump (&fkb_after, 1);
269	    }
270	  if (fks1b[-1] != fk1b_before)
271	    {
272	      printf ("before fks1b:"); mpn_dump (fks1b - 1, 1);
273	      printf ("keep:   "); mpn_dump (&fk1b_before, 1);
274	    }
275	  if (fks1b[2 * mn + 1] != fk1b_after)
276	    {
277	      printf ("after fks1b:"); mpn_dump (fks1b + 2 * mn + 1, 1);
278	      printf ("keep:   "); mpn_dump (&fk1b_after, 1);
279	    }
280	  abort();
281	}
282
283      if (mpn_cmp (fks1a, fkb, mn) != 0)
284	{
285	  if (mpn_sub_n (fkb, mp, fkb, mn) || mpn_cmp (fks1a, fkb, mn) != 0)
286	    {
287	      printf ("ERROR(k) in test %d, exp2 = %lu\n", test, exp2);
288	      mpn_dump (fks1a, mn);
289	      mpn_dump (fkb, mn);
290	      mpn_dump (mp, mn);
291	      abort();
292	    }
293	  signflip ^= 1;
294	}
295
296      if (mpn_cmp (fka, fks1b, mn) != 0)
297	{
298	  if (mpn_sub_n (fks1b, mp, fks1b, mn) || mpn_cmp (fka, fks1b, mn) != 0)
299	    {
300	      printf ("ERROR(k-1) in test %d, exp2 = %lu\n", test, exp2);
301	      mpn_dump (fka, mn);
302	      mpn_dump (fks1b, mn);
303	      mpn_dump (mp, mn);
304	      abort();
305	    }
306	  signflip ^= 1;
307	}
308
309      if (signflip != 0 && ! mpn_zero_p (fks1b, mn) && ! mpn_zero_p (fkb, mn))
310	{
311	  if ((mp [0] & 1) == 0) /* Should we test only odd modulus-es? */
312	    {
313	      if (! mpn_lshift (fks1b, fks1b, mn, 1) &&
314		  mpn_cmp (mp, fks1b, mn) == 0)
315		continue;
316	      if (! mpn_lshift (fkb, fkb, mn, 1) &&
317		  mpn_cmp (mp, fkb, mn) == 0)
318		continue;
319	    }
320	  printf ("ERROR(sign) in test %d, exp2 = %lu\n",
321		  test, exp2);
322	  abort();
323	}
324    }
325  TMP_FREE;
326  return 0;
327}
328
329int
330main (int argc, char **argv)
331{
332  int count = COUNT;
333  gmp_randstate_ptr rands;
334
335  tests_start ();
336  TESTS_REPS (count, argv, argc);
337  rands = RANDS;
338
339  test_fib2_fib2m (count / 2, rands);
340  test_fib2m_2exp (count / 2, rands);
341
342  tests_end ();
343  exit (0);
344}
345