xref: /netbsd-src/external/lgpl3/gmp/dist/tests/mpz/t-primorial_ui.c (revision 72c7faa4dbb41dbb0238d6b4a109da0d4b236dd4)
1 /* Exercise mpz_primorial_ui.
2 
3 Copyright 2000-2002, 2012, 2015 Free Software Foundation, Inc.
4 
5 This file is part of the GNU MP Library test suite.
6 
7 The GNU MP Library test suite is free software; you can redistribute it
8 and/or modify it under the terms of the GNU General Public License as
9 published by the Free Software Foundation; either version 3 of the License,
10 or (at your option) any later version.
11 
12 The GNU MP Library test suite is distributed in the hope that it will be
13 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
15 Public License for more details.
16 
17 You should have received a copy of the GNU General Public License along with
18 the GNU MP Library test suite.  If not, see https://www.gnu.org/licenses/.  */
19 
20 #include <stdio.h>
21 #include <stdlib.h>
22 #include "gmp-impl.h"
23 #include "tests.h"
24 
25 
26 /* Usage: t-primorial_ui [x|num]
27 
28    With no arguments testing goes up to the initial value of "limit" below.
29    With a number argument tests are carried that far, or with a literal "x"
30    tests are continued without limit (this being meant only for development
31    purposes).  */
32 
33 static int isprime (unsigned long int t);
34 
35 int
main(int argc,char * argv[])36 main (int argc, char *argv[])
37 {
38   unsigned long  n;
39   unsigned long  limit = 2222;
40   gmp_randstate_ptr rands;
41   mpz_t          f, r, bs;
42 
43   tests_start ();
44   rands = RANDS;
45 
46   if (argc > 1 && argv[1][0] == 'x')
47     limit = ULONG_MAX;
48   else
49     TESTS_REPS (limit, argv, argc);
50 
51   /* for small limb testing */
52   limit = MIN (limit, MP_LIMB_T_MAX);
53 
54   mpz_init_set_ui (f, 1);  /* 0# = 1 */
55   mpz_init (r);
56 
57   n = 0;
58   do
59     {
60       mpz_primorial_ui (r, n);
61       MPZ_CHECK_FORMAT (r);
62 
63       if (mpz_cmp (f, r) != 0)
64 	{
65 	  printf ("mpz_primorial_ui(%lu) wrong\n", n);
66 	  printf ("  got  "); mpz_out_str (stdout, 10, r); printf("\n");
67 	  printf ("  want "); mpz_out_str (stdout, 10, f); printf("\n");
68 	  abort ();
69 	}
70 
71       if (isprime (++n))
72 	mpz_mul_ui (f, f, n);  /* p# = (p-1)# * (p) */
73       if (n%16 == 0) { mpz_clear (r); mpz_init (r); }
74     } while (n < limit);
75 
76   n = 0; limit =1;
77   mpz_init (bs);
78   do
79     {
80       unsigned long i, d;
81 
82       mpz_urandomb (bs, rands, 21);
83       i = mpz_get_ui (bs);
84       mpz_urandomb (bs, rands, 9);
85       d = mpz_get_ui (bs) + 3*64;
86       mpz_primorial_ui (f, i);
87       MPZ_CHECK_FORMAT (f);
88       mpz_primorial_ui (r, i+d);
89       MPZ_CHECK_FORMAT (r);
90 
91       do {
92 	if (isprime (++i))
93 	  mpz_mul_ui (f, f, i);
94       } while (--d != 0);
95 
96       if (mpz_cmp (f, r) != 0)
97 	{
98 	  printf ("mpz_primorial_ui(%lu) wrong\n", i);
99 	  printf ("  got  "); mpz_out_str (stdout, 10, r); printf("\n");
100 	  printf ("  want "); mpz_out_str (stdout, 10, f); printf("\n");
101 	  abort ();
102 	}
103     } while (++n < limit);
104   /* Chech a single "big" value, modulo a larger prime */
105   n = 2095637;
106   mpz_primorial_ui (r, n);
107   mpz_set_ui (f, 13);
108   mpz_setbit (f, 64); /* f = 2^64 + 13 */
109   mpz_tdiv_r (r, r, f);
110   mpz_set_str (f, "BAFCBF3C95B217D5", 16);
111 
112   if (mpz_cmp (f, r) != 0)
113     {
114       printf ("mpz_primorial_ui(%lu) wrong\n", n);
115       printf ("  got  "); mpz_out_str (stdout, 10, r); printf("\n");
116       printf ("  want "); mpz_out_str (stdout, 10, f); printf("\n");
117       abort ();
118     }
119 
120   mpz_clear (bs);
121   mpz_clear (f);
122   mpz_clear (r);
123 
124   tests_end ();
125 
126   exit (0);
127 }
128 
129 static int
isprime(unsigned long int t)130 isprime (unsigned long int t)
131 {
132   unsigned long int q, r, d;
133 
134   if (t < 3 || (t & 1) == 0)
135     return t == 2;
136 
137   for (d = 3, r = 1; r != 0; d += 2)
138     {
139       q = t / d;
140       r = t - q * d;
141       if (q < d)
142 	return 1;
143     }
144   return 0;
145 }
146