1 // ========================================================================== 2 // $Source: /var/lib/cvs/Givaro/src/kernel/rational/givratreconstruct.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: Jean-Guillaume Dumas 9 // $Id: givratreconstruct.C,v 1.4 2010-04-12 15:54:39 jgdumas Exp $ 10 // ========================================================================== 11 // Description: 12 13 #include "givaro/givrational.h" 14 #include <iostream> 15 16 namespace Givaro { 17 18 Rational(const Integer & f,const Integer & m,const Integer & k,bool recurs)19 Rational::Rational(const Integer& f, const Integer& m, const Integer& k, bool recurs ) 20 { 21 bool res = this->ratrecon(f,m,k,recurs); 22 if (recurs) 23 for( Integer newk = k + 1; (!res) && (newk<f) ; ++newk) 24 res = this->ratrecon(f,m,newk,true); 25 } 26 27 28 ratrecon(const Integer & f,const Integer & m,const Integer & k,bool recurs)29 bool Rational::ratrecon(const Integer& f, const Integer& m, const Integer& k, bool recurs ) 30 { 31 32 // #ifdef GIVARO_DEBUG 33 // std::cerr << "RatRecon : " << f << " " << m << " " << k << std::endl; 34 // #endif 35 36 37 Integer r0, t0, r1, t1, q, u; 38 r0=m; 39 t0=0; 40 r1=f; 41 if (f<0) r1+= m; 42 t1=1; 43 while(r1>=k) 44 { 45 // #ifdef GIVARO_DEBUG 46 // std::cerr << "r0: " << r0 47 // << ", r1: " << r1 48 // << ", q: " << q 49 // << ", t0: " << t0 50 // << ", t1: " << t1 51 // << std::endl; 52 // #endif 53 q = r0; 54 q /= r1; // r0/r1 55 56 u = r1; 57 r1 = r0; // r1 <-- r0 58 r0 = u; // r0 <-- r1 59 u *= q; 60 r1 -= u; // r1 <-- r0-q*r1 61 if (r1 == 0) break; 62 63 u = t1; 64 t1 = t0; // r1 <-- r0 65 t0 = u; // r0 <-- r1 66 u *= q; 67 t1 -= u; // r1 <-- r0-q*r1 68 } 69 70 // [GG, MCA, 1999] Theorem 5.26 71 // (i) 72 if (t1 < 0) { 73 num = -r1; 74 den = -t1; 75 } else { 76 num = r1; 77 den = t1; 78 } 79 80 // (ii) 81 if (gcd(num,den) != 1) { 82 // #ifdef GIVARO_DEBUG 83 // std::cerr << "num: " << num 84 // << ", den: " << den 85 // << ", g: " << gcd(num,den) 86 // << std::endl; 87 // std::cerr << "r0: " << r0 88 // << ", r1: " << r1 89 // << ", k: " << k 90 // << std::endl; 91 // #endif 92 q = r0; 93 q += r1; 94 q -= k; 95 q /= r1; 96 97 // #ifdef GIVARO_DEBUG 98 // std::cerr << "q: " << q << std::endl; 99 // #endif 100 101 r0 -= q * r1; 102 t0 -= q * t1; 103 104 if (t0 < 0) { 105 num = -r0; 106 den = -t0; 107 } else { 108 num = r0; 109 den = t0; 110 } 111 112 if (t0 > m/k) { 113 if (!recurs) 114 std::cerr 115 << "*** Error *** No rational reconstruction of " 116 << f 117 << " modulo " 118 << m 119 << " with denominator <= " 120 << (m/k) 121 << std::endl; 122 } 123 if (gcd(num,den) != 1) { 124 if (!recurs) 125 std::cerr 126 << "*** Error *** There exists no rational reconstruction of " 127 << f 128 << " modulo " 129 << m 130 << " with |numerator| < " 131 << k 132 << std::endl 133 << "*** Error *** But " 134 << num 135 << " = " 136 << den 137 << " * " 138 << f 139 << " modulo " 140 << m 141 << std::endl; 142 return false; 143 } 144 } 145 // #ifdef GIVARO_DEBUG 146 // std::cerr << "RatRecon End " << std::endl; 147 // #endif 148 return true; 149 } 150 151 152 } // namespace Givaro 153 154 /* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 155 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 156