gcd.cc

Go to the documentation of this file.
00001 /*
00002 
00003 Copyright (C) 2004-2012 David Bateman
00004 Copyright (C) 2010 Jaroslav Hajek, Jordi GutiƩrrez Hermoso
00005 
00006 This file is part of Octave.
00007 
00008 Octave is free software; you can redistribute it and/or modify it
00009 under the terms of the GNU General Public License as published by the
00010 Free Software Foundation; either version 3 of the License, or (at your
00011 option) any later version.
00012 
00013 Octave is distributed in the hope that it will be useful, but WITHOUT
00014 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00015 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00016 for more details.
00017 
00018 You should have received a copy of the GNU General Public License
00019 along with Octave; see the file COPYING.  If not, see
00020 <http://www.gnu.org/licenses/>.
00021 
00022 */
00023 
00024 #ifdef HAVE_CONFIG_H
00025 #include <config.h>
00026 #endif
00027 
00028 #include "dNDArray.h"
00029 #include "CNDArray.h"
00030 #include "fNDArray.h"
00031 #include "fCNDArray.h"
00032 #include "lo-mappers.h"
00033 #include "oct-binmap.h"
00034 
00035 #include "defun-dld.h"
00036 #include "error.h"
00037 #include "oct-obj.h"
00038 
00039 static double
00040 simple_gcd (double a, double b)
00041 {
00042   if (! xisinteger (a) || ! xisinteger (b))
00043     (*current_liboctave_error_handler)
00044       ("gcd: all values must be integers");
00045 
00046   double aa = fabs (a);
00047   double bb = fabs (b);
00048 
00049   while (bb != 0)
00050     {
00051       double tt = fmod (aa, bb);
00052       aa = bb;
00053       bb = tt;
00054     }
00055 
00056   return aa;
00057 }
00058 
00059 // Don't use the Complex and FloatComplex typedefs because we need to
00060 // refer to the actual float precision FP in the body (and when gcc
00061 // implements template aliases from C++0x, can do a small fix here).
00062 template <typename FP>
00063 static void
00064 divide (const std::complex<FP>& a, const std::complex<FP>& b,
00065         std::complex<FP>& q, std::complex<FP>& r)
00066 {
00067   FP qr = gnulib::floor ((a/b).real () + 0.5);
00068   FP qi = gnulib::floor ((a/b).imag () + 0.5);
00069 
00070   q = std::complex<FP> (qr, qi);
00071 
00072   r = a - q*b;
00073 }
00074 
00075 template <typename FP>
00076 static std::complex<FP>
00077 simple_gcd (const std::complex<FP>& a, const std::complex<FP>& b)
00078 {
00079   if (! xisinteger (a.real ()) || ! xisinteger(a.imag ())
00080       || ! xisinteger (b.real ()) || ! xisinteger(b.imag ()))
00081     (*current_liboctave_error_handler)
00082       ("gcd: all complex parts must be integers");
00083 
00084   std::complex<FP> aa = a;
00085   std::complex<FP> bb = b;
00086 
00087   if (abs (aa) < abs (bb))
00088     std::swap (aa, bb);
00089 
00090   while (abs (bb) != 0)
00091     {
00092       std::complex<FP> qq, rr;
00093       divide (aa, bb, qq, rr);
00094       aa = bb;
00095       bb = rr;
00096     }
00097 
00098   return aa;
00099 }
00100 
00101 template <class T>
00102 static octave_int<T>
00103 simple_gcd (const octave_int<T>& a, const octave_int<T>& b)
00104 {
00105   T aa = a.abs ().value ();
00106   T bb = b.abs ().value ();
00107 
00108   while (bb != 0)
00109     {
00110       T tt = aa % bb;
00111       aa = bb;
00112       bb = tt;
00113     }
00114 
00115   return aa;
00116 }
00117 
00118 static double
00119 extended_gcd (double a, double b, double& x, double& y)
00120 {
00121   if (! xisinteger (a) || ! xisinteger (b))
00122     (*current_liboctave_error_handler)
00123       ("gcd: all values must be integers");
00124 
00125   double aa = fabs (a);
00126   double bb = fabs (b);
00127 
00128   double xx = 0, yy = 1;
00129   double lx = 1, ly = 0;
00130 
00131   while (bb != 0)
00132     {
00133       double qq = gnulib::floor (aa / bb);
00134       double tt = fmod (aa, bb);
00135 
00136       aa = bb;
00137       bb = tt;
00138 
00139       double tx = lx - qq*xx;
00140       lx = xx;
00141       xx = tx;
00142 
00143       double ty = ly - qq*yy;
00144       ly = yy;
00145       yy = ty;
00146     }
00147 
00148   x = a >= 0 ? lx : -lx;
00149   y = b >= 0 ? ly : -ly;
00150 
00151   return aa;
00152 }
00153 
00154 template <typename FP>
00155 static std::complex<FP>
00156 extended_gcd (const std::complex<FP>& a, const std::complex<FP>& b,
00157               std::complex<FP>& x, std::complex<FP>& y)
00158 {
00159   if (! xisinteger (a.real ()) || ! xisinteger(a.imag ())
00160       || ! xisinteger (b.real ()) || ! xisinteger(b.imag ()))
00161     (*current_liboctave_error_handler)
00162       ("gcd: all complex parts must be integers");
00163 
00164   std::complex<FP> aa = a, bb = b;
00165   bool swapped = false;
00166   if (abs (aa) < abs (bb))
00167     {
00168       std::swap (aa, bb);
00169       swapped = true;
00170     }
00171 
00172   std::complex<FP> xx = 0, lx = 1;
00173   std::complex<FP> yy = 1, ly = 0;
00174 
00175   while (abs(bb) != 0)
00176     {
00177       std::complex<FP> qq, rr;
00178       divide (aa, bb, qq, rr);
00179       aa = bb;
00180       bb = rr;
00181 
00182       std::complex<FP> tx = lx - qq*xx;
00183       lx = xx;
00184       xx = tx;
00185 
00186       std::complex<FP> ty = ly - qq*yy;
00187       ly = yy;
00188       yy = ty;
00189     }
00190 
00191   x = lx;
00192   y = ly;
00193 
00194   if (swapped)
00195     std::swap (x, y);
00196 
00197   return aa;
00198 }
00199 
00200 template <class T>
00201 static octave_int<T>
00202 extended_gcd (const octave_int<T>& a, const octave_int<T>& b,
00203               octave_int<T>& x, octave_int<T>& y)
00204 {
00205   T aa = a.abs ().value ();
00206   T bb = b.abs ().value ();
00207   T xx = 0, lx = 1;
00208   T yy = 1, ly = 0;
00209 
00210   while (bb != 0)
00211     {
00212       T qq = aa / bb;
00213       T tt = aa % bb;
00214       aa = bb;
00215       bb = tt;
00216 
00217       T tx = lx - qq*xx;
00218       lx = xx;
00219       xx = tx;
00220 
00221       T ty = ly - qq*yy;
00222       ly = yy;
00223       yy = ty;
00224     }
00225 
00226   x = octave_int<T> (lx) * a.signum ();
00227   y = octave_int<T> (ly) * b.signum ();
00228 
00229   return aa;
00230 }
00231 
00232 template<class NDA>
00233 static octave_value
00234 do_simple_gcd (const octave_value& a, const octave_value& b)
00235 {
00236   typedef typename NDA::element_type T;
00237   octave_value retval;
00238 
00239   if (a.is_scalar_type () && b.is_scalar_type ())
00240     {
00241       // Optimize scalar case.
00242       T aa = octave_value_extract<T> (a);
00243       T bb = octave_value_extract<T> (b);
00244       retval = simple_gcd (aa, bb);
00245     }
00246   else
00247     {
00248       NDA aa = octave_value_extract<NDA> (a);
00249       NDA bb = octave_value_extract<NDA> (b);
00250       retval = binmap<T> (aa, bb, simple_gcd, "gcd");
00251     }
00252 
00253   return retval;
00254 }
00255 
00256 // Dispatcher
00257 static octave_value
00258 do_simple_gcd (const octave_value& a, const octave_value& b)
00259 {
00260   octave_value retval;
00261   builtin_type_t btyp = btyp_mixed_numeric (a.builtin_type (),
00262                                             b.builtin_type ());
00263   switch (btyp)
00264     {
00265     case btyp_double:
00266       if (a.is_sparse_type () && b.is_sparse_type ())
00267         {
00268           retval = do_simple_gcd<SparseMatrix> (a, b);
00269           break;
00270         }
00271       // fall through!
00272 
00273     case btyp_float:
00274       retval = do_simple_gcd<NDArray> (a, b);
00275       break;
00276 
00277 #define MAKE_INT_BRANCH(X) \
00278     case btyp_ ## X: \
00279       retval = do_simple_gcd<X ## NDArray> (a, b); \
00280       break
00281 
00282     MAKE_INT_BRANCH (int8);
00283     MAKE_INT_BRANCH (int16);
00284     MAKE_INT_BRANCH (int32);
00285     MAKE_INT_BRANCH (int64);
00286     MAKE_INT_BRANCH (uint8);
00287     MAKE_INT_BRANCH (uint16);
00288     MAKE_INT_BRANCH (uint32);
00289     MAKE_INT_BRANCH (uint64);
00290 
00291 #undef MAKE_INT_BRANCH
00292 
00293     case btyp_complex:
00294       retval = do_simple_gcd<ComplexNDArray> (a, b);
00295       break;
00296 
00297     case btyp_float_complex:
00298       retval = do_simple_gcd<FloatComplexNDArray> (a, b);
00299       break;
00300 
00301     default:
00302       error ("gcd: invalid class combination for gcd: %s and %s\n",
00303              a.class_name ().c_str (), b.class_name ().c_str ());
00304     }
00305 
00306   if (btyp == btyp_float)
00307     retval = retval.float_array_value ();
00308 
00309   return retval;
00310 }
00311 
00312 template<class NDA>
00313 static octave_value
00314 do_extended_gcd (const octave_value& a, const octave_value& b,
00315                  octave_value& x, octave_value& y)
00316 {
00317   typedef typename NDA::element_type T;
00318   octave_value retval;
00319 
00320   if (a.is_scalar_type () && b.is_scalar_type ())
00321     {
00322       // Optimize scalar case.
00323       T aa = octave_value_extract<T> (a);
00324       T bb = octave_value_extract<T> (b);
00325       T xx, yy;
00326       retval = extended_gcd (aa, bb, xx, yy);
00327       x = xx;
00328       y = yy;
00329     }
00330   else
00331     {
00332       NDA aa = octave_value_extract<NDA> (a);
00333       NDA bb = octave_value_extract<NDA> (b);
00334 
00335       dim_vector dv = aa.dims ();
00336       if (aa.numel () == 1)
00337         dv = bb.dims ();
00338       else if (bb.numel () != 1 && bb.dims () != dv)
00339         gripe_nonconformant ("gcd", a.dims (), b.dims ());
00340 
00341       NDA gg (dv), xx (dv), yy (dv);
00342 
00343       const T *aptr = aa.fortran_vec ();
00344       const T *bptr = bb.fortran_vec ();
00345 
00346       bool inca = aa.numel () != 1;
00347       bool incb = bb.numel () != 1;
00348 
00349       T *gptr = gg.fortran_vec ();
00350       T *xptr = xx.fortran_vec (), *yptr = yy.fortran_vec ();
00351 
00352       octave_idx_type n = gg.numel ();
00353       for (octave_idx_type i = 0; i < n; i++)
00354         {
00355           octave_quit ();
00356 
00357           *gptr++ = extended_gcd (*aptr, *bptr, *xptr++, *yptr++);
00358 
00359           aptr += inca;
00360           bptr += incb;
00361         }
00362 
00363       x = xx;
00364       y = yy;
00365 
00366       retval = gg;
00367     }
00368 
00369   return retval;
00370 }
00371 
00372 // Dispatcher
00373 static octave_value
00374 do_extended_gcd (const octave_value& a, const octave_value& b,
00375                  octave_value& x, octave_value& y)
00376 {
00377   octave_value retval;
00378 
00379   builtin_type_t btyp = btyp_mixed_numeric (a.builtin_type (),
00380                                             b.builtin_type ());
00381   switch (btyp)
00382     {
00383     case btyp_double:
00384     case btyp_float:
00385       retval = do_extended_gcd<NDArray> (a, b, x, y);
00386       break;
00387 
00388 #define MAKE_INT_BRANCH(X) \
00389     case btyp_ ## X: \
00390       retval = do_extended_gcd<X ## NDArray> (a, b, x, y); \
00391       break
00392 
00393     MAKE_INT_BRANCH (int8);
00394     MAKE_INT_BRANCH (int16);
00395     MAKE_INT_BRANCH (int32);
00396     MAKE_INT_BRANCH (int64);
00397     MAKE_INT_BRANCH (uint8);
00398     MAKE_INT_BRANCH (uint16);
00399     MAKE_INT_BRANCH (uint32);
00400     MAKE_INT_BRANCH (uint64);
00401 
00402 #undef MAKE_INT_BRANCH
00403 
00404     case btyp_complex:
00405       retval = do_extended_gcd<ComplexNDArray> (a, b, x, y);
00406       break;
00407 
00408     case btyp_float_complex:
00409       retval = do_extended_gcd<FloatComplexNDArray> (a, b, x, y);
00410       break;
00411 
00412     default:
00413       error ("gcd: invalid class combination for gcd: %s and %s\n",
00414              a.class_name ().c_str (), b.class_name ().c_str ());
00415     }
00416 
00417   // For consistency.
00418   if (! error_state && a.is_sparse_type () && b.is_sparse_type ())
00419     {
00420       retval = retval.sparse_matrix_value ();
00421       x = x.sparse_matrix_value ();
00422       y = y.sparse_matrix_value ();
00423     }
00424 
00425   if (btyp == btyp_float)
00426     {
00427       retval = retval.float_array_value ();
00428       x = x.float_array_value ();
00429       y = y.float_array_value ();
00430     }
00431 
00432   return retval;
00433 }
00434 
00435 DEFUN_DLD (gcd, args, nargout,
00436   "-*- texinfo -*-\n\
00437 @deftypefn  {Loadable Function} {@var{g} =} gcd (@var{a1}, @var{a2}, @dots{})\n\
00438 @deftypefnx {Loadable Function} {[@var{g}, @var{v1}, @dots{}] =} gcd (@var{a1}, @var{a2}, @dots{})\n\
00439 \n\
00440 Compute the greatest common divisor of @var{a1}, @var{a2}, @dots{}.  If more\n\
00441 than one argument is given all arguments must be the same size or scalar.\n\
00442 In this case the greatest common divisor is calculated for each element\n\
00443 individually.  All elements must be ordinary or Gaussian (complex)\n\
00444 integers.  Note that for Gaussian integers, the gcd is not unique up to\n\
00445 units (multiplication by 1, -1, @var{i} or -@var{i}), so an arbitrary\n\
00446 greatest common divisor amongst four possible is returned.  For example,\n\
00447 \n\
00448 @noindent\n\
00449 and\n\
00450 \n\
00451 @example\n\
00452 @group\n\
00453 gcd ([15, 9], [20, 18])\n\
00454     @result{}  5  9\n\
00455 @end group\n\
00456 @end example\n\
00457 \n\
00458 Optional return arguments @var{v1}, etc., contain integer vectors such\n\
00459 that,\n\
00460 \n\
00461 @tex\n\
00462 $g = v_1 a_1 + v_2 a_2 + \\cdots$\n\
00463 @end tex\n\
00464 @ifnottex\n\
00465 \n\
00466 @example\n\
00467 @var{g} = @var{v1} .* @var{a1} + @var{v2} .* @var{a2} + @dots{}\n\
00468 @end example\n\
00469 \n\
00470 @end ifnottex\n\
00471 \n\
00472 @seealso{lcm, factor}\n\
00473 @end deftypefn")
00474 {
00475   octave_value_list retval;
00476 
00477   int nargin = args.length ();
00478 
00479   if (nargin > 1)
00480     {
00481       if (nargout > 1)
00482         {
00483           retval.resize (nargin + 1);
00484 
00485           retval(0) = do_extended_gcd (args(0), args(1), retval(1), retval(2));
00486 
00487           for (int j = 2; j < nargin; j++)
00488             {
00489               octave_value x;
00490               retval(0) = do_extended_gcd (retval(0), args(j),
00491                                            x, retval(j+1));
00492               for (int i = 0; i < j; i++)
00493                 retval(i+1).assign (octave_value::op_el_mul_eq, x);
00494 
00495               if (error_state)
00496                 break;
00497             }
00498         }
00499       else
00500         {
00501           retval(0) = do_simple_gcd (args(0), args(1));
00502 
00503           for (int j = 2; j < nargin; j++)
00504             {
00505               retval(0) = do_simple_gcd (retval(0), args(j));
00506 
00507               if (error_state)
00508                 break;
00509             }
00510         }
00511     }
00512   else
00513     print_usage ();
00514 
00515   return retval;
00516 }
00517 
00518 /*
00519 
00520 %!assert(gcd (200, 300, 50, 35), 5)
00521 %!assert(gcd (int16(200), int16(300), int16(50), int16(35)), int16(5))
00522 %!assert(gcd (uint64(200), uint64(300), uint64(50), uint64(35)), uint64(5))
00523 %!assert(gcd (18-i, -29+3i), -3-4i)
00524 
00525 %!error <Invalid call to gcd> gcd ();
00526 
00527 %!test
00528 %! s.a = 1;
00529 %! fail("gcd (s)");
00530 
00531 */
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines