diff -r b099acf06b55 src/ChangeLog
--- a/src/ChangeLog Wed Sep 29 04:25:57 2010 -0400
+++ b/src/ChangeLog Wed Sep 29 03:56:56 2010 -0500
@@ -1,3 +1,14 @@
+2010-09-29 Jordi GutiƩrrez Hermoso
+
+ * DLD-FUNCTIONS/gcd.cc (divide): New function, complex integer
+ division with remainder.
+ (simple_gcd): Overload for complex values.
+ (extended_gcd): Ditto. Also fix small bug that didn't distinguish
+ the two output coefficients.
+ (do_simple_gcd): Dispatch for complex gcd.
+ (do_extended_gcd): Ditto.
+ (Fgcd): Mention that complex gcd is now also possible.
+
2010-09-29 John W. Eaton
* DLD-FUNCTIONS/urlwrite.cc (F__ftp_dir__): Use octave_scalar_map
diff -r b099acf06b55 src/DLD-FUNCTIONS/gcd.cc
--- a/src/DLD-FUNCTIONS/gcd.cc Wed Sep 29 04:25:57 2010 -0400
+++ b/src/DLD-FUNCTIONS/gcd.cc Wed Sep 29 03:56:56 2010 -0500
@@ -1,7 +1,7 @@
/*
Copyright (C) 2004, 2005, 2006, 2007, 2008, 2009 David Bateman
-Copyirght (C) 2010 Jaroslav Hajek
+Copyright (C) 2010 Jaroslav Hajek, Jordi GutiƩrrez Hermoso
This file is part of Octave.
@@ -36,7 +36,7 @@
#include "error.h"
#include "oct-obj.h"
-static double
+static double
simple_gcd (double a, double b)
{
if (! xisinteger (a) || ! xisinteger (b))
@@ -53,6 +53,47 @@
return aa;
}
+// Don't use the Complex and FloatComplex typedefs because we need to
+// refer to the actual float precision FP in the body (and when gcc
+// implements template aliases from C++0x, can do a small fix here).
+template
+static void
+divide(const std::complex& a, const std::complex& b,
+ std::complex& q, std::complex& r)
+{
+ FP qr = floor((a/b).real() + 0.5);
+ FP qi = floor((a/b).imag() + 0.5);
+ q = std::complex(qr,qi);
+ r = a - q*b;
+}
+
+template
+static std::complex
+simple_gcd (const std::complex& a, const std::complex& b)
+{
+ if ( ! xisinteger (a.real ()) || ! xisinteger(a.imag ()) ||
+ ! xisinteger (b.real ()) || ! xisinteger(b.imag ()) )
+ (*current_liboctave_error_handler)
+ ("gcd: all complex parts must be integers");
+
+ std::complex aa = a, bb = b;
+
+ if ( abs(aa) < abs(bb) )
+ {
+ std::swap(aa,bb);
+ }
+
+ while ( abs(bb) != 0)
+ {
+ std::complex qq, rr;
+ divide (aa, bb, qq, rr);
+ aa = bb;
+ bb = rr;
+ }
+
+ return aa;
+}
+
template
static octave_int
simple_gcd (const octave_int& a, const octave_int& b)
@@ -67,7 +108,7 @@
return aa;
}
-static double
+static double
extended_gcd (double a, double b, double& x, double& y)
{
if (! xisinteger (a) || ! xisinteger (b))
@@ -75,7 +116,7 @@
("gcd: all values must be integers");
double aa = fabs (a), bb = fabs (b);
- double xx = 0, lx = 1, yy = 0, ly = 1;
+ double xx = 0, lx = 1, yy = 1, ly = 0;
while (bb != 0)
{
double qq = floor (aa / bb);
@@ -89,19 +130,61 @@
ly = yy; yy = ty;
}
- x = a >= 0 ? lx : -lx;
+ x = a >= 0 ? lx : -lx;
y = b >= 0 ? ly : -ly;
return aa;
}
+template
+static std::complex
+extended_gcd (const std::complex& a, const std::complex& b,
+ std::complex& x, std::complex& y)
+{
+ if ( ! xisinteger (a.real ()) || ! xisinteger(a.imag ()) ||
+ ! xisinteger (b.real ()) || ! xisinteger(b.imag ()) )
+ (*current_liboctave_error_handler)
+ ("gcd: all complex parts must be integers");
+
+ std::complex aa = a, bb = b;
+ bool swapped = false;
+ if ( abs(aa) < abs(bb) )
+ {
+ std::swap(aa,bb);
+ swapped = true;
+ }
+
+ std::complex xx = 0, lx = 1, yy = 1, ly = 0;
+ while ( abs(bb) != 0)
+ {
+ std::complex qq, rr;
+ divide (aa, bb, qq, rr);
+ aa = bb; bb = rr;
+
+ std::complex tx = lx - qq*xx;
+ lx = xx; xx = tx;
+
+ std::complex ty = ly - qq*yy;
+ ly = yy; yy = ty;
+ }
+
+ x = lx;
+ y = ly;
+ if (swapped)
+ {
+ std::swap(x,y);
+ }
+
+ return aa;
+}
+
template
static octave_int
extended_gcd (const octave_int& a, const octave_int& b,
octave_int& x, octave_int& y)
{
T aa = a.abs ().value (), bb = b.abs ().value ();
- T xx = 0, lx = 1, yy = 0, ly = 1;
+ T xx = 0, lx = 1, yy = 1, ly = 0;
while (bb != 0)
{
T qq = aa / bb;
@@ -178,8 +261,11 @@
MAKE_INT_BRANCH (uint64);
#undef MAKE_INT_BRANCH
case btyp_complex:
+ retval = do_simple_gcd (a,b);
+ break;
case btyp_float_complex:
- error ("gcd: complex numbers not allowed");
+ retval = do_simple_gcd (a,b);
+ break;
default:
error ("gcd: invalid class combination for gcd: %s and %s\n",
a.class_name ().c_str (), b.class_name ().c_str ());
@@ -275,8 +361,11 @@
MAKE_INT_BRANCH (uint64);
#undef MAKE_INT_BRANCH
case btyp_complex:
+ retval = do_extended_gcd (a, b, x, y);
+ break;
case btyp_float_complex:
- error ("gcd: complex numbers not allowed");
+ retval = do_extended_gcd (a, b, x, y);
+ break;
default:
error ("gcd: invalid class combination for gcd: %s and %s\n",
a.class_name ().c_str (), b.class_name ().c_str ());
@@ -309,7 +398,10 @@
Compute the greatest common divisor of @var{a1}, @var{a2}, @dots{}. If more\n\
than one argument is given all arguments must be the same size or scalar.\n\
In this case the greatest common divisor is calculated for each element\n\
-individually. All elements must be integers. For example,\n\
+individually. All elements must be ordinary or Gaussian (complex)\n\
+integers. Note that for Gaussian integers, the gcd is not unique up to\n\
+units (multiplication by 1, -1, @var{i} or address@hidden), so an arbitrary\n\
+greatest common divisor amongst four possible is returned. For example,\n\
\n\
@noindent\n\
and\n\
@@ -350,7 +442,7 @@
for (int j = 2; j < nargin; j++)
{
octave_value x;
- retval(0) = do_extended_gcd (retval(0), args(1),
+ retval(0) = do_extended_gcd (retval(0), args(1),
x, retval(j+1));
for (int i = 0; i < j; i++)
retval(i).assign (octave_value::op_el_mul_eq, x);
@@ -380,6 +472,7 @@
%!assert(gcd (200, 300, 50, 35), 5)
%!assert(gcd (int16(200), int16(300), int16(50), int16(35)), int16(5))
%!assert(gcd (uint64(200), uint64(300), uint64(50), uint64(35)), uint64(5))
+%!assert(gcd (18-i, -29+3i), -3-4i)
%!error gcd ();