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