1 // ========================================================================== 2 // $Source: /var/lib/cvs/Givaro/src/kernel/gmp++/gmp++_int_gcd.C,v $ 3 // Copyright(c)'1994-2009 by The Givaro group 4 // This file is part of Givaro. 5 // Givaro is governed by the CeCILL-B license under French law 6 // and abiding by the rules of distribution of free software. 7 // see the COPYRIGHT file for more details. 8 // Authors: M. Samama, T. Gautier 9 // $Id: gmp++_int_gcd.C,v 1.4 2011-01-18 17:49:06 jgdumas Exp $ 10 // ========================================================================== 11 // Description: 12 // 13 /** @file gmp++/gmp++_int_gcd.C 14 * gcding stuff. 15 */ 16 17 #ifndef __GIVARO_gmpxx_gmpxx_int_gcd_C 18 #define __GIVARO_gmpxx_gmpxx_int_gcd_C 19 20 #ifndef __GIVARO_INLINE_ALL 21 #include "gmp++/gmp++.h" 22 #endif 23 #include "givaro/giverror.h" 24 namespace Givaro { 25 // ========================================================================== 26 // Computes and returns the lcm of the two integers a and b. lcm(const Integer & a,const Integer & b)27 Integer lcm(const Integer& a, const Integer& b) { 28 Integer Res(Integer::one); 29 mpz_lcm( (mpz_ptr)&(Res.gmp_rep), (mpz_srcptr)&(a.gmp_rep), (mpz_srcptr)&(b.gmp_rep) ) ; 30 if (Res.priv_sign() <0) return -Res; 31 else return Res ; 32 } 33 lcm(Integer & g,const Integer & a,const Integer & b)34 Integer& lcm(Integer& g, const Integer& a, const Integer& b) { 35 mpz_lcm( (mpz_ptr)&(g.gmp_rep), (mpz_srcptr)&(a.gmp_rep), (mpz_srcptr)&(b.gmp_rep) ) ; 36 if (g.priv_sign() <0) return Integer::negin(g); 37 else return g ; 38 } 39 40 41 // ========================================================================== 42 // Computes and returns the gcd of the two integers a and b. gcd(const Integer & a,const Integer & b)43 Integer gcd(const Integer& a, const Integer& b) 44 { 45 Integer Res(Integer::one); 46 mpz_gcd( (mpz_ptr)&(Res.gmp_rep), (mpz_srcptr)&(a.gmp_rep), (mpz_srcptr)&(b.gmp_rep) ) ; 47 if (Res.priv_sign() <0) return -Res; 48 return Res ; 49 } 50 gcd(Integer & g,const Integer & a,const Integer & b)51 Integer& gcd(Integer& g, const Integer& a, const Integer& b) 52 { 53 mpz_gcd( (mpz_ptr)&(g.gmp_rep), (mpz_srcptr)&(a.gmp_rep), (mpz_srcptr)&(b.gmp_rep) ) ; 54 if (g.priv_sign() <0) return Integer::negin(g); 55 return g ; 56 } 57 inv(Integer & u,const Integer & a,const Integer & b)58 Integer& inv(Integer& u, const Integer& a, const Integer& b) { 59 #ifdef GIVARO_DEBUG 60 const int res = 61 #endif 62 mpz_invert( (mpz_ptr)&(u.gmp_rep), (mpz_srcptr)&(a.gmp_rep), (mpz_srcptr)&(b.gmp_rep) ) ; 63 #ifdef GIVARO_DEBUG 64 if(! res) { 65 throw GivMathDivZero("*** Error: division by zero, in operator Integer::inv in gmp++_int_gcd.C") ; 66 } 67 #endif 68 return u ; 69 } 70 invin(Integer & u,const Integer & b)71 Integer& invin(Integer& u, const Integer& b) { 72 return inv(u,u,b); 73 } 74 75 // ========================================================================== 76 // Computes and returns the gcd g of the two integers a and b such that 77 // g = a*u + b*v . 78 // The algorithm used is this of Gmp. gcd(Integer & u,Integer & v,const Integer & a,const Integer & b)79 Integer gcd (Integer& u, Integer& v, 80 const Integer& a, const Integer& b ) 81 { 82 v = 1; // v must not be 0 to be computed. 83 Integer Res(Integer::one); 84 mpz_gcdext( (mpz_ptr)&(Res.gmp_rep), (mpz_ptr)&(u.gmp_rep), (mpz_ptr)&(v.gmp_rep), 85 (mpz_srcptr)&(a.gmp_rep), (mpz_srcptr)&(b.gmp_rep) ) ; 86 if (Res.priv_sign() < 0) { 87 Integer::negin(u); 88 Integer::negin(v); 89 return Integer::negin(Res); 90 } 91 // { u = -u ; v = -v ; return -Res;} 92 return Res; 93 } 94 gcd(Integer & g,Integer & u,Integer & v,const Integer & a,const Integer & b)95 Integer& gcd (Integer& g, Integer& u, Integer& v, 96 const Integer& a, const Integer& b) 97 { 98 v = 1; // v must not be 0 to be computed. 99 mpz_gcdext( (mpz_ptr)&(g.gmp_rep), (mpz_ptr)&(u.gmp_rep), (mpz_ptr)&(v.gmp_rep), 100 (mpz_srcptr)&(a.gmp_rep), (mpz_srcptr)&(b.gmp_rep) ) ; 101 if (g.priv_sign() < 0) { Integer::negin(u); Integer::negin(v); return Integer::negin(g);} 102 return g; 103 } 104 105 pp(const Integer & P,const Integer & Q)106 Integer pp( const Integer& P, const Integer& Q ) 107 { 108 Integer U = P ; 109 Integer V = gcd(P,Q) ; 110 // -- computes the prime part U of g relatively to U 111 while ( V != Integer::one ) 112 { 113 U = U / V ; 114 V = gcd( U,V) ; 115 } 116 return U ; 117 } 118 119 } 120 121 #endif // __GIVARO_gmpxx_gmpxx_int_gcd_C 122 123 /* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 124 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 125