1 /* mpz_gcdext(g, s, t, a, b) -- Set G to gcd(a, b), and S and T such that 2 g = as + bt. 3 4 Copyright 1991, 1993-1997, 2000, 2001, 2005, 2011, 2012 Free Software 5 Foundation, Inc. 6 7 This file is part of the GNU MP Library. 8 9 The GNU MP Library is free software; you can redistribute it and/or modify 10 it under the terms of either: 11 12 * the GNU Lesser General Public License as published by the Free 13 Software Foundation; either version 3 of the License, or (at your 14 option) any later version. 15 16 or 17 18 * the GNU General Public License as published by the Free Software 19 Foundation; either version 2 of the License, or (at your option) any 20 later version. 21 22 or both in parallel, as here. 23 24 The GNU MP Library is distributed in the hope that it will be useful, but 25 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 26 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 27 for more details. 28 29 You should have received copies of the GNU General Public License and the 30 GNU Lesser General Public License along with the GNU MP Library. If not, 31 see https://www.gnu.org/licenses/. */ 32 33 #include <stdio.h> /* for NULL */ 34 #include "gmp.h" 35 #include "gmp-impl.h" 36 37 void 38 mpz_gcdext (mpz_ptr g, mpz_ptr s, mpz_ptr t, mpz_srcptr a, mpz_srcptr b) 39 { 40 mp_size_t asize, bsize; 41 mp_ptr tmp_ap, tmp_bp; 42 mp_size_t gsize, ssize, tmp_ssize; 43 mp_ptr gp, tmp_gp, tmp_sp; 44 TMP_DECL; 45 46 /* mpn_gcdext requires that Usize >= Vsize. Therefore, we often 47 have to swap U and V. The computed cofactor will be the 48 "smallest" one, which is faster to produce. The wanted one will 49 be computed here; this is needed anyway when both are requested. */ 50 51 asize = ABSIZ (a); 52 bsize = ABSIZ (b); 53 54 if (asize < bsize) 55 { 56 MPZ_SRCPTR_SWAP (a, b); 57 MP_SIZE_T_SWAP (asize, bsize); 58 MPZ_PTR_SWAP (s, t); 59 } 60 61 if (bsize == 0) 62 { 63 /* g = |a|, s = sgn(a), t = 0. */ 64 ssize = SIZ (a) >= 0 ? (asize != 0) : -1; 65 66 gp = MPZ_REALLOC (g, asize); 67 MPN_COPY (gp, PTR (a), asize); 68 SIZ (g) = asize; 69 70 if (t != NULL) 71 SIZ (t) = 0; 72 if (s != NULL) 73 { 74 SIZ (s) = ssize; 75 PTR (s)[0] = 1; 76 } 77 return; 78 } 79 80 TMP_MARK; 81 82 TMP_ALLOC_LIMBS_2 (tmp_ap, asize, tmp_bp, bsize); 83 MPN_COPY (tmp_ap, PTR (a), asize); 84 MPN_COPY (tmp_bp, PTR (b), bsize); 85 86 TMP_ALLOC_LIMBS_2 (tmp_gp, bsize, tmp_sp, bsize + 1); 87 88 gsize = mpn_gcdext (tmp_gp, tmp_sp, &tmp_ssize, tmp_ap, asize, tmp_bp, bsize); 89 90 ssize = ABS (tmp_ssize); 91 tmp_ssize = SIZ (a) >= 0 ? tmp_ssize : -tmp_ssize; 92 93 if (t != NULL) 94 { 95 mpz_t x; 96 __mpz_struct gtmp, stmp; 97 98 PTR (>mp) = tmp_gp; 99 SIZ (>mp) = gsize; 100 101 PTR (&stmp) = tmp_sp; 102 SIZ (&stmp) = tmp_ssize; 103 104 MPZ_TMP_INIT (x, ssize + asize + 1); 105 mpz_mul (x, &stmp, a); 106 mpz_sub (x, >mp, x); 107 mpz_divexact (t, x, b); 108 } 109 110 if (s != NULL) 111 { 112 mp_ptr sp; 113 114 sp = MPZ_REALLOC (s, ssize); 115 MPN_COPY (sp, tmp_sp, ssize); 116 SIZ (s) = tmp_ssize; 117 } 118 119 gp = MPZ_REALLOC (g, gsize); 120 MPN_COPY (gp, tmp_gp, gsize); 121 SIZ (g) = gsize; 122 123 TMP_FREE; 124 } 125