1 /* Exercise mpz_fac_ui and mpz_bin_uiui. 2 3 Copyright 2000, 2001, 2002, 2012, 2013 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 http://www.gnu.org/licenses/. */ 19 20 #include <stdio.h> 21 #include <stdlib.h> 22 23 #include "testutils.h" 24 25 /* Usage: t-fac_ui [x|num] 26 27 With no arguments testing goes up to the initial value of "limit" below. 28 With a number argument tests are carried that far, or with a literal "x" 29 tests are continued without limit (this being meant only for development 30 purposes). */ 31 32 void 33 try_mpz_bin_uiui (mpz_srcptr want, unsigned long n, unsigned long k) 34 { 35 mpz_t got; 36 37 mpz_init (got); 38 mpz_bin_uiui (got, n, k); 39 if (mpz_cmp (got, want) != 0) 40 { 41 printf ("mpz_bin_uiui wrong\n"); 42 printf (" n=%lu\n", n); 43 printf (" k=%lu\n", k); 44 printf (" got="); mpz_out_str (stdout, 10, got); printf ("\n"); 45 printf (" want="); mpz_out_str (stdout, 10, want); printf ("\n"); 46 abort(); 47 } 48 mpz_clear (got); 49 } 50 51 /* Test all bin(n,k) cases, with 0 <= k <= n + 1 <= count. */ 52 void 53 bin_smallexaustive (unsigned int count) 54 { 55 mpz_t want; 56 unsigned long n, k; 57 58 mpz_init (want); 59 60 for (n = 0; n < count; n++) 61 { 62 mpz_set_ui (want, 1); 63 for (k = 0; k <= n; k++) 64 { 65 try_mpz_bin_uiui (want, n, k); 66 mpz_mul_ui (want, want, n - k); 67 mpz_fdiv_q_ui (want, want, k + 1); 68 } 69 try_mpz_bin_uiui (want, n, k); 70 } 71 72 mpz_clear (want); 73 } 74 75 /* Test all fac(n) cases, with 0 <= n <= limit. */ 76 void 77 fac_smallexaustive (unsigned int limit) 78 { 79 mpz_t f, r; 80 unsigned long n; 81 mpz_init_set_si (f, 1); /* 0! = 1 */ 82 mpz_init (r); 83 84 for (n = 0; n < limit; n++) 85 { 86 mpz_fac_ui (r, n); 87 88 if (mpz_cmp (f, r) != 0) 89 { 90 printf ("mpz_fac_ui(%lu) wrong\n", n); 91 printf (" got "); mpz_out_str (stdout, 10, r); printf("\n"); 92 printf (" want "); mpz_out_str (stdout, 10, f); printf("\n"); 93 abort (); 94 } 95 96 mpz_mul_ui (f, f, n+1); /* (n+1)! = n! * (n+1) */ 97 } 98 99 mpz_clear (f); 100 mpz_clear (r); 101 } 102 103 void checkWilson (mpz_t f, unsigned long n) 104 { 105 unsigned long m; 106 107 mpz_fac_ui (f, n - 1); 108 m = mpz_fdiv_ui (f, n); 109 if ( m != n - 1) 110 { 111 printf ("mpz_fac_ui(%lu) wrong\n", n - 1); 112 printf (" Wilson's theorem not verified: got %lu, expected %lu.\n",m ,n - 1); 113 abort (); 114 } 115 } 116 117 void 118 checkprimes (unsigned long p1, unsigned long p2, unsigned long p3) 119 { 120 mpz_t b, f; 121 122 if (p1 - 1 != p2 - 1 + p3 - 1) 123 { 124 printf ("checkprimes(%lu, %lu, %lu) wrong\n", p1, p2, p3); 125 printf (" %lu - 1 != %lu - 1 + %lu - 1 \n", p1, p2, p3); 126 abort (); 127 } 128 129 mpz_init (b); 130 mpz_init (f); 131 132 checkWilson (b, p1); /* b = (p1-1)! */ 133 checkWilson (f, p2); /* f = (p2-1)! */ 134 mpz_divexact (b, b, f); 135 checkWilson (f, p3); /* f = (p3-1)! */ 136 mpz_divexact (b, b, f); /* b = (p1-1)!/((p2-1)!(p3-1)!) */ 137 mpz_bin_uiui (f, p1 - 1, p2 - 1); 138 if (mpz_cmp (f, b) != 0) 139 { 140 printf ("checkprimes(%lu, %lu, %lu) wrong\n", p1, p2, p3); 141 printf (" got "); mpz_out_str (stdout, 10, b); printf("\n"); 142 printf (" want "); mpz_out_str (stdout, 10, f); printf("\n"); 143 abort (); 144 } 145 146 mpz_clear (b); 147 mpz_clear (f); 148 149 } 150 151 void 152 testmain (int argc, char *argv[]) 153 { 154 unsigned long limit = 128; 155 156 if (argc > 1 && argv[1][0] == 'x') 157 limit = ~ limit; 158 else if (argc > 1) 159 limit = atoi (argv[1]); 160 161 checkprimes(1009, 733, 277); 162 fac_smallexaustive (limit); 163 bin_smallexaustive (limit); 164 } 165