# HG changeset patch
# User Jordi Gutiérrez Hermoso <address@hidden>
# Date 1285830331 18000
# Node ID 905492c1ed8cbb195456dea9c873ec09c8fede17
# Parent  a4d245b374aa79cd17c31f97690e9edfe805691a
Implement gcd for Gaussian (complex) integers

diff -r a4d245b374aa -r 905492c1ed8c src/ChangeLog
--- a/src/ChangeLog	Thu Sep 30 02:05:30 2010 -0500
+++ b/src/ChangeLog	Thu Sep 30 02:05:31 2010 -0500
@@ -1,3 +1,12 @@
+2010-09-30  Jordi Gutiérrez Hermoso  <address@hidden>
+	* DLD-FUNCTIONS/gcd.cc (divide): New function, complex integer
+	division with remainder.
+	(simple_gcd): Overload for complex values.
+	(extended_gcd): Ditto.
+	(do_simple_gcd): Dispatch for complex gcd.
+	(do_extended_gcd): Ditto.
+	(Fgcd): Mention that complex gcd is now also possible.
+
 2010-09-30  Jordi Gutiérrez Hermoso  <address@hidden>
 
 	* DLD-FUNCTIONS/gcd.cc (extended_gcd): Fix bug that didn't
diff -r a4d245b374aa -r 905492c1ed8c src/DLD-FUNCTIONS/gcd.cc
--- a/src/DLD-FUNCTIONS/gcd.cc	Thu Sep 30 02:05:30 2010 -0500
+++ b/src/DLD-FUNCTIONS/gcd.cc	Thu Sep 30 02:05:31 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 <typename FP>
+static void
+divide(const std::complex<FP>& a, const std::complex<FP>& b,
+       std::complex<FP>& q, std::complex<FP>& r)
+{
+  FP qr = floor((a/b).real() + 0.5);
+  FP qi = floor((a/b).imag() + 0.5);
+  q = std::complex<FP>(qr,qi);
+  r = a - q*b;
+}
+
+template <typename FP>
+static std::complex<FP>
+simple_gcd (const std::complex<FP>& a, const std::complex<FP>& 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<FP> aa = a, bb = b;
+
+  if ( abs(aa) < abs(bb) )
+    {
+      std::swap(aa,bb);
+    }
+
+  while ( abs(bb) != 0)
+    {
+      std::complex<FP> qq, rr;
+      divide (aa, bb, qq, rr);
+      aa = bb;
+      bb = rr;
+    }
+
+  return aa;
+}
+
 template <class T>
 static octave_int<T>
 simple_gcd (const octave_int<T>& a, const octave_int<T>& 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))
@@ -89,12 +130,54 @@
       ly = yy; yy = ty;
     }
 
-  x = a >= 0 ? lx : -lx; 
+  x = a >= 0 ? lx : -lx;
   y = b >= 0 ? ly : -ly;
 
   return aa;
 }
 
+template <typename FP>
+static std::complex<FP>
+extended_gcd (const std::complex<FP>& a, const std::complex<FP>& b,
+              std::complex<FP>& x, std::complex<FP>& 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<FP> aa = a, bb = b;
+  bool swapped = false;
+  if ( abs(aa) < abs(bb) )
+    {
+      std::swap(aa,bb);
+      swapped = true;
+    }
+
+  std::complex<FP> xx = 0, lx = 1, yy = 1, ly = 0;
+  while ( abs(bb) != 0)
+    {
+      std::complex<FP> qq, rr;
+      divide (aa, bb, qq, rr);
+      aa = bb; bb = rr;
+
+      std::complex<FP> tx = lx - qq*xx;
+      lx = xx; xx = tx;
+
+      std::complex<FP> ty = ly - qq*yy;
+      ly = yy; yy = ty;
+    }
+
+  x = lx;
+  y = ly;
+  if (swapped)
+    {
+      std::swap(x,y);
+    }
+
+  return aa;
+}
+
 template <class T>
 static octave_int<T>
 extended_gcd (const octave_int<T>& a, const octave_int<T>& b,
@@ -178,8 +261,11 @@
     MAKE_INT_BRANCH (uint64);
 #undef MAKE_INT_BRANCH
     case btyp_complex:
+      retval = do_simple_gcd<ComplexNDArray> (a,b);
+      break;
     case btyp_float_complex:
-      error ("gcd: complex numbers not allowed");
+      retval = do_simple_gcd<FloatComplexNDArray> (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<ComplexNDArray> (a, b, x, y);
+      break;
     case btyp_float_complex:
-      error ("gcd: complex numbers not allowed");
+      retval = do_extended_gcd<FloatComplexNDArray> (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\
@@ -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 <Invalid call to gcd.*> gcd ();
 
# HG changeset patch
# User Jordi Gutiérrez Hermoso <address@hidden>
# Date 1285830330 18000
# Node ID a4d245b374aa79cd17c31f97690e9edfe805691a
# Parent  b721e12140ccfb01d15eeeaf4cee8995ca80e996
Fix off-by-one and typo bugs in gcd

diff -r b721e12140cc -r a4d245b374aa src/ChangeLog
--- a/src/ChangeLog	Wed Sep 29 22:08:34 2010 +0200
+++ b/src/ChangeLog	Thu Sep 30 02:05:30 2010 -0500
@@ -1,3 +1,9 @@
+2010-09-30  Jordi Gutiérrez Hermoso  <address@hidden>
+
+	* DLD-FUNCTIONS/gcd.cc (extended_gcd): Fix bug that didn't
+	distinguish the two output coefficients.
+	(Fgcd): Fix off-by-one bug and typo from copy-pasted code.
+
 2010-09-29  Jaroslav Hajek  <address@hidden>
 
 	* oct-map.cc (octave_map::contents): Fix off-by-1 error.
diff -r b721e12140cc -r a4d245b374aa src/DLD-FUNCTIONS/gcd.cc
--- a/src/DLD-FUNCTIONS/gcd.cc	Wed Sep 29 22:08:34 2010 +0200
+++ b/src/DLD-FUNCTIONS/gcd.cc	Thu Sep 30 02:05:30 2010 -0500
@@ -75,7 +75,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);
@@ -101,7 +101,7 @@
               octave_int<T>& x, octave_int<T>& 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;
@@ -350,10 +350,10 @@
           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(j),
                                            x, retval(j+1));
               for (int i = 0; i < j; i++)
-                retval(i).assign (octave_value::op_el_mul_eq, x);
+                retval(i+1).assign (octave_value::op_el_mul_eq, x);
               if (error_state)
                 break;
             }