xref: /netbsd-src/external/lgpl3/gmp/dist/tests/mpz/t-root.c (revision 72c7faa4dbb41dbb0238d6b4a109da0d4b236dd4)
1 /* Test mpz_root, mpz_rootrem, and mpz_perfect_power_p.
2 
3 Copyright 1991, 1993, 1994, 1996, 2000, 2001, 2009, 2015 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 
24 #include "gmp-impl.h"
25 #include "tests.h"
26 
27 void debug_mp (mpz_t, int);
28 
29 void
check_one(mpz_t root1,mpz_t x2,unsigned long nth,int res,int i)30 check_one (mpz_t root1, mpz_t x2, unsigned long nth, int res, int i)
31 {
32   mpz_t temp, temp2;
33   mpz_t root2, rem2;
34 
35   mpz_init (root2);
36   mpz_init (rem2);
37   mpz_init (temp);
38   mpz_init (temp2);
39 
40   MPZ_CHECK_FORMAT (root1);
41 
42   mpz_rootrem (root2, rem2, x2, nth);
43   MPZ_CHECK_FORMAT (root2);
44   MPZ_CHECK_FORMAT (rem2);
45 
46   mpz_pow_ui (temp, root1, nth);
47   MPZ_CHECK_FORMAT (temp);
48 
49   mpz_add (temp2, temp, rem2);
50 
51   /* Is power of result > argument?  */
52   if (mpz_cmp (root1, root2) != 0 || mpz_cmp (x2, temp2) != 0 || mpz_cmpabs (temp, x2) > 0 || res == mpz_cmp_ui (rem2, 0))
53     {
54       fprintf (stderr, "ERROR after test %d\n", i);
55       debug_mp (x2, 10);
56       debug_mp (root1, 10);
57       debug_mp (root2, 10);
58       fprintf (stderr, "nth: %lu ,res: %i\n", nth, res);
59       abort ();
60     }
61 
62   if (nth > 1 && mpz_cmp_ui (temp, 1L) > 0 && ! mpz_perfect_power_p (temp))
63     {
64       fprintf (stderr, "ERROR in mpz_perfect_power_p after test %d\n", i);
65       debug_mp (temp, 10);
66       debug_mp (root1, 10);
67       fprintf (stderr, "nth: %lu\n", nth);
68       abort ();
69     }
70 
71   if (nth <= 10000 && mpz_sgn(x2) > 0)		/* skip too expensive test */
72     {
73       mpz_add_ui (temp2, root1, 1L);
74       mpz_pow_ui (temp2, temp2, nth);
75       MPZ_CHECK_FORMAT (temp2);
76 
77       /* Is square of (result + 1) <= argument?  */
78       if (mpz_cmp (temp2, x2) <= 0)
79 	{
80 	  fprintf (stderr, "ERROR after test %d\n", i);
81 	  debug_mp (x2, 10);
82 	  debug_mp (root1, 10);
83 	  fprintf (stderr, "nth: %lu\n", nth);
84 	  abort ();
85 	}
86     }
87 
88   mpz_clear (root2);
89   mpz_clear (rem2);
90   mpz_clear (temp);
91   mpz_clear (temp2);
92 }
93 
94 int
main(int argc,char ** argv)95 main (int argc, char **argv)
96 {
97   mpz_t x2;
98   mpz_t root1;
99   mp_size_t x2_size;
100   int i, res;
101   int reps = 500;
102   unsigned long nth;
103   gmp_randstate_ptr rands;
104   mpz_t bs;
105   unsigned long bsi, size_range;
106 
107   tests_start ();
108   TESTS_REPS (reps, argv, argc);
109 
110   rands = RANDS;
111 
112   mpz_init (bs);
113 
114   mpz_init (x2);
115   mpz_init (root1);
116 
117   /* This triggers a gcc 4.3.2 bug */
118   mpz_set_str (x2, "ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff80000000000000000000000000000000000000000000000000000000000000002", 16);
119   res = mpz_root (root1, x2, 2);
120   check_one (root1, x2, 2, res, -1);
121 
122   for (i = 0; i < reps; i++)
123     {
124       mpz_urandomb (bs, rands, 32);
125       size_range = mpz_get_ui (bs) % 17 + 2;
126 
127       mpz_urandomb (bs, rands, size_range);
128       x2_size = mpz_get_ui (bs) + 10;
129       mpz_rrandomb (x2, rands, x2_size);
130 
131       mpz_urandomb (bs, rands, 15);
132       nth = mpz_getlimbn (bs, 0) % mpz_sizeinbase (x2, 2) + 2;
133 
134       res = mpz_root (root1, x2, nth);
135 
136       mpz_urandomb (bs, rands, 4);
137       bsi = mpz_get_ui (bs);
138       if ((bsi & 1) != 0)
139 	{
140 	  /* With 50% probability, set x2 near a perfect power.  */
141 	  mpz_pow_ui (x2, root1, nth);
142 	  if ((bsi & 2) != 0)
143 	    {
144 	      mpz_sub_ui (x2, x2, bsi >> 2);
145 	      mpz_abs (x2, x2);
146 	    }
147 	  else
148 	    mpz_add_ui (x2, x2, bsi >> 2);
149 	  res = mpz_root (root1, x2, nth);
150 	}
151 
152       check_one (root1, x2, nth, res, i);
153 
154       if (((nth & 1) != 0) && ((bsi & 2) != 0))
155 	{
156 	  mpz_neg (x2, x2);
157 	  mpz_neg (root1, root1);
158 	  check_one (root1, x2, nth, res, i);
159 	}
160     }
161 
162   mpz_clear (bs);
163   mpz_clear (x2);
164   mpz_clear (root1);
165 
166   tests_end ();
167   exit (0);
168 }
169 
170 void
debug_mp(mpz_t x,int base)171 debug_mp (mpz_t x, int base)
172 {
173   mpz_out_str (stderr, base, x); fputc ('\n', stderr);
174 }
175