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