songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

new / more robust version .-)

This commit is contained in:
Michael Hemmer 2008-03-31 12:43:38 +00:00
parent c05f1298fa
commit 7dcb9b1258
1 changed files with 40 additions and 33 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

73
Algebraic_foundations/include/CGAL/extended_euclidean_algorithm.h
View File

@ -23,46 +23,53 @@
#define CGAL_EXTENDED_EUCLIDEAN_ALGORITHM_H 1
#include <CGAL/basic.h>
#include <vector>
CGAL_BEGIN_NAMESPACE
template< class NT >
NT extended_euclidean_algorithm(const NT& a_, const NT& b_, NT& u, NT& v){
// EEA computing the normalized gcd
// Modern Computer Algebra (Hardcover)
// by Joachim von zur Gathen (Author), Jürgen Gerhard (Author)
// Publisher: Cambridge University Press; 2 edition (September 1, 2003)
// Language: English
// ISBN-10: 0521826462
// ISBN-13: 978-0521826464
// pp.: 55
template< class AS >
AS extended_euclidean_algorithm(const AS& f, const AS& g, AS& s_, AS& t_){
typename Algebraic_structure_traits<AS>::Integral_division idiv;
typename Algebraic_structure_traits<AS>::Div div;
typename Algebraic_structure_traits<AS>::Unit_part unit_part;
typedef Algebraic_structure_traits<NT> AST;
typename AST::Div_mod div_mod;
typename AST::Unit_part unit_part;
typename AST::Integral_division idiv;
std::vector<AS> p,r,s,t,q;
p.push_back(unit_part(f));
r.push_back(idiv(f,p[0]));
s.push_back(idiv(AS(1),p[0]));
t.push_back(AS(0));
q.push_back(AS(0));
NT unit_part_a(unit_part(a_));
NT unit_part_b(unit_part(b_));
p.push_back(unit_part(g));
r.push_back(idiv(g,p[1]));
s.push_back(AS(0));
t.push_back(idiv(AS(1),p[1]));
NT a(idiv(a_,unit_part_a));
NT b(idiv(b_,unit_part_b));
NT x(0),y(1),last_x(1),last_y(0);
NT temp, quotient;
//TODO: unroll to avoid swapping
while (b != 0){
temp = b;
div_mod(a,b,quotient,b);
a = temp;
temp = x;
x = last_x-quotient*x;
last_x = temp;
temp = y;
y = last_y-quotient*y;
last_y = temp;
int i = 1;
while(!is_zero(r[i])){
q.push_back(div(r[i-1],r[i]));
r.push_back(r[i-1]-q[i]*r[i]);
p.push_back(unit_part(r[i+1]));
r[i+1] = idiv(r[i+1],p[i+1]);
s.push_back(idiv(s[i-1]-q[i]*s[i],p[i+1]));
t.push_back(idiv(t[i-1]-q[i]*t[i],p[i+1]));
i++;
}
u = last_x * unit_part_a;
v = last_y * unit_part_b;
CGAL_precondition(unit_part(a) == NT(1));
CGAL_precondition(a == a_*u + b_*v);
return a;
s_=s[i-1];
t_=t[i-1];
AS h = r[i-1];
CGAL_precondition( h == f*s_ + g*t_);
return h;
}
CGAL_END_NAMESPACE