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