xref: /netbsd/external/lgpl3/gmp/dist/tests/mpz/t-lucm.c (revision 671ea119)
1 /* Test mpz_powm, mpz_lucas_mod.
2 
3 Copyright 1991, 1993, 1994, 1996, 1999-2001, 2009, 2012, 2018 Free Software
4 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 8
31 
32 /* FIXME: Should we implement another sequence to test lucas mod?	*/
33 /* Eg: a generalisation of what we use for Fibonacci:	*/
34 /* U_{2n-1} = U_n^2 - Q*U_{n-1}^2	*/
35 /* U_{2n+1} = D*U_n^2  + Q*U_{2n-1} + 2*Q^n ; whith D = (P^2-4*Q)	*/
36 /* P*U_{2n} = U_{2n+1} + Q*U_{2n-1}	*/
37 
38 int
main(int argc,char ** argv)39 main (int argc, char **argv)
40 {
41   mpz_t base, exp, mod;
42   mpz_t r1, r2, t1, t2;
43   mp_size_t base_size, exp_size, mod_size;
44   int i, res;
45   int reps = 1000;
46   long Q;
47   gmp_randstate_ptr rands;
48   mpz_t bs;
49   unsigned long bsi, size_range;
50 
51   tests_start ();
52   TESTS_REPS (reps, argv, argc);
53 
54   rands = RANDS;
55 
56   mpz_init (bs);
57 
58   mpz_init (base);
59   mpz_init (exp);
60   mpz_init (mod);
61   mpz_init (r1);
62   mpz_init (r2);
63   mpz_init (t1);
64   mpz_init (t2);
65 
66   for (i = 0; i < reps; i++)
67     {
68       mpz_urandomb (bs, rands, 32);
69       size_range = mpz_get_ui (bs) % SIZEM + 1;
70 
71       do  /* Loop until base >= 2 and fits in a long.  */
72 	{
73 	  mpz_urandomb (base, rands, BITS_PER_ULONG - 2);
74 	}
75       while (mpz_cmp_ui (base, 2) < 0 || mpz_fits_slong_p (base) == 0);
76 
77       Q = mpz_get_ui (base);
78 
79       do
80         {
81 	  ++size_range;
82 	  size_range = MIN (size_range, SIZEM);
83 	  mpz_urandomb (bs, rands, size_range);
84 	  mod_size = mpz_get_ui (bs);
85 	  mpz_rrandomb (mod, rands, mod_size);
86 	  mpz_add_ui (mod, mod, 16);
87 	}
88       while (mpz_gcd_ui (NULL, mod, Q) != 1);
89 
90       mod_size = mpz_sizeinbase (mod, 2) - 3;
91       mpz_urandomb (bs, rands, 32);
92       exp_size = mpz_get_ui (bs) % mod_size + 2;
93 
94       mpz_tdiv_q_2exp (exp, mod, exp_size);
95       mpz_add_ui (exp, exp, 1);
96 
97       mpz_urandomb (bs, rands, 2);
98       bsi = mpz_get_ui (bs);
99       if ((bsi & 1) != 0)
100 	{
101 	  mpz_neg (base, base);
102 	  Q = -Q;
103 	}
104 
105       res = mpz_lucas_mod (t1, r2, Q, exp_size, mod, t2, r1);
106       if (res && ++reps)
107 	continue;
108       MPZ_CHECK_FORMAT (r2);
109       if (mpz_cmp_ui (r2, 0) < 0)
110 	mpz_add (r2, r2, mod);
111       mpz_powm (r1, base, exp, mod);
112 
113       if (mpz_cmp (r1, r2) != 0)
114 	{
115 	  fprintf (stderr, "\nIncorrect results in test %d for operands:\n", i);
116 	  debug_mp (base, -16);
117 	  debug_mp (exp, -16);
118 	  debug_mp (mod, -16);
119 	  fprintf (stderr, "mpz_powm result:\n");
120 	  debug_mp (r1, -16);
121 	  fprintf (stderr, "mpz_lucas_mod result (%d) Q=%ld:\n", res, Q);
122 	  debug_mp (r2, -16);
123 	  abort ();
124 	}
125     }
126 
127   mpz_clear (bs);
128   mpz_clear (base);
129   mpz_clear (exp);
130   mpz_clear (mod);
131   mpz_clear (r1);
132   mpz_clear (r2);
133   mpz_clear (t1);
134   mpz_clear (t2);
135 
136   tests_end ();
137   exit (0);
138 }
139 
140 void
debug_mp(mpz_t x,int base)141 debug_mp (mpz_t x, int base)
142 {
143   mpz_out_str (stderr, base, x); fputc ('\n', stderr);
144 }
145