1 /* 2 3 Copyright 2012, 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 <limits.h> 21 #include <stdlib.h> 22 #include <stdio.h> 23 24 #include "testutils.h" 25 26 #define MAXBITS 400 27 #define COUNT 10000 28 29 static void 30 dump (const char *label, const mpz_t x) 31 { 32 char *buf = mpz_get_str (NULL, 16, x); 33 fprintf (stderr, "%s: %s\n", label, buf); 34 testfree (buf); 35 } 36 37 /* Called when g is supposed to be gcd(a,b), and g = s a + t b. */ 38 static int 39 gcdext_valid_p (const mpz_t a, const mpz_t b, 40 const mpz_t g, const mpz_t s, const mpz_t t) 41 { 42 mpz_t ta, tb, r; 43 44 /* It's not clear that gcd(0,0) is well defined, but we allow it and 45 require that gcd(0,0) = 0. */ 46 if (mpz_sgn (g) < 0) 47 return 0; 48 49 if (mpz_sgn (a) == 0) 50 { 51 /* Must have g == abs (b). Any value for s is in some sense "correct", 52 but it makes sense to require that s == 0. */ 53 return mpz_cmpabs (g, b) == 0 && mpz_sgn (s) == 0; 54 } 55 else if (mpz_sgn (b) == 0) 56 { 57 /* Must have g == abs (a), s == sign (a) */ 58 return mpz_cmpabs (g, a) == 0 && mpz_cmp_si (s, mpz_sgn (a)) == 0; 59 } 60 61 if (mpz_sgn (g) <= 0) 62 return 0; 63 64 mpz_init (ta); 65 mpz_init (tb); 66 mpz_init (r); 67 68 mpz_mul (ta, s, a); 69 mpz_mul (tb, t, b); 70 mpz_add (ta, ta, tb); 71 72 if (mpz_cmp (ta, g) != 0) 73 { 74 fail: 75 mpz_clear (ta); 76 mpz_clear (tb); 77 mpz_clear (r); 78 return 0; 79 } 80 mpz_tdiv_qr (ta, r, a, g); 81 if (mpz_sgn (r) != 0) 82 goto fail; 83 84 mpz_tdiv_qr (tb, r, b, g); 85 if (mpz_sgn (r) != 0) 86 goto fail; 87 88 /* Require that 2 |s| < |b/g|, or |s| == 1. */ 89 if (mpz_cmpabs_ui (s, 1) > 0) 90 { 91 mpz_mul_2exp (r, s, 1); 92 if (mpz_cmpabs (r, tb) > 0) 93 goto fail; 94 } 95 96 /* Require that 2 |t| < |a/g| or |t| == 1*/ 97 if (mpz_cmpabs_ui (t, 1) > 0) 98 { 99 mpz_mul_2exp (r, t, 1); 100 if (mpz_cmpabs (r, ta) > 0) 101 return 0; 102 } 103 104 mpz_clear (ta); 105 mpz_clear (tb); 106 mpz_clear (r); 107 108 return 1; 109 } 110 111 void 112 testmain (int argc, char **argv) 113 { 114 unsigned i; 115 mpz_t a, b, g, s, t; 116 117 mpz_init (a); 118 mpz_init (b); 119 mpz_init (g); 120 mpz_init (s); 121 mpz_init (t); 122 123 for (i = 0; i < COUNT; i++) 124 { 125 mini_random_op3 (OP_GCD, MAXBITS, a, b, s); 126 mpz_gcd (g, a, b); 127 if (mpz_cmp (g, s)) 128 { 129 fprintf (stderr, "mpz_gcd failed:\n"); 130 dump ("a", a); 131 dump ("b", b); 132 dump ("r", g); 133 dump ("ref", s); 134 abort (); 135 } 136 } 137 138 for (i = 0; i < COUNT; i++) 139 { 140 unsigned flags; 141 mini_urandomb (a, 32); 142 flags = mpz_get_ui (a); 143 mini_rrandomb (a, MAXBITS); 144 mini_rrandomb (b, MAXBITS); 145 146 if (flags % 37 == 0) 147 mpz_mul (a, a, b); 148 if (flags % 37 == 1) 149 mpz_mul (b, a, b); 150 151 if (flags & 1) 152 mpz_neg (a, a); 153 if (flags & 2) 154 mpz_neg (b, b); 155 156 mpz_gcdext (g, s, t, a, b); 157 if (!gcdext_valid_p (a, b, g, s, t)) 158 { 159 fprintf (stderr, "mpz_gcdext failed:\n"); 160 dump ("a", a); 161 dump ("b", b); 162 dump ("g", g); 163 dump ("s", s); 164 dump ("t", t); 165 abort (); 166 } 167 168 mpz_gcd (s, a, b); 169 if (mpz_cmp (g, s)) 170 { 171 fprintf (stderr, "mpz_gcd failed:\n"); 172 dump ("a", a); 173 dump ("b", b); 174 dump ("r", g); 175 dump ("ref", s); 176 abort (); 177 } 178 } 179 mpz_clear (a); 180 mpz_clear (b); 181 mpz_clear (g); 182 mpz_clear (s); 183 mpz_clear (t); 184 } 185