t-comb.c revision 1.1.1.1
1/* Exercise mpz_fac_ui and mpz_bin_uiui. 2 3Copyright 2000, 2001, 2002, 2012, 2013 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 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 32void 33try_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. */ 52void 53bin_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. */ 76void 77fac_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 103void 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 117void 118checkprimes (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 151void 152testmain (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