GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
gcd.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2004-2015 David Bateman
4 Copyright (C) 2010 Jaroslav Hajek, Jordi GutiĆ©rrez Hermoso
5 
6 This file is part of Octave.
7 
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21 
22 */
23 
24 #ifdef HAVE_CONFIG_H
25 #include <config.h>
26 #endif
27 
28 #include "dNDArray.h"
29 #include "CNDArray.h"
30 #include "fNDArray.h"
31 #include "fCNDArray.h"
32 #include "lo-mappers.h"
33 #include "oct-binmap.h"
34 
35 #include "defun.h"
36 #include "error.h"
37 #include "oct-obj.h"
38 
39 static double
40 simple_gcd (double a, double b)
41 {
42  if (! xisinteger (a) || ! xisinteger (b))
43  (*current_liboctave_error_handler)
44  ("gcd: all values must be integers");
45 
46  double aa = fabs (a);
47  double bb = fabs (b);
48 
49  while (bb != 0)
50  {
51  double tt = fmod (aa, bb);
52  aa = bb;
53  bb = tt;
54  }
55 
56  return aa;
57 }
58 
59 // Don't use the Complex and FloatComplex typedefs because we need to
60 // refer to the actual float precision FP in the body (and when gcc
61 // implements template aliases from C++0x, can do a small fix here).
62 template <typename FP>
63 static void
64 divide (const std::complex<FP>& a, const std::complex<FP>& b,
65  std::complex<FP>& q, std::complex<FP>& r)
66 {
67  FP qr = gnulib::floor ((a/b).real () + 0.5);
68  FP qi = gnulib::floor ((a/b).imag () + 0.5);
69 
70  q = std::complex<FP> (qr, qi);
71 
72  r = a - q*b;
73 }
74 
75 template <typename FP>
76 static std::complex<FP>
77 simple_gcd (const std::complex<FP>& a, const std::complex<FP>& b)
78 {
79  if (! xisinteger (a.real ()) || ! xisinteger (a.imag ())
80  || ! xisinteger (b.real ()) || ! xisinteger (b.imag ()))
81  (*current_liboctave_error_handler)
82  ("gcd: all complex parts must be integers");
83 
84  std::complex<FP> aa = a;
85  std::complex<FP> bb = b;
86 
87  if (abs (aa) < abs (bb))
88  std::swap (aa, bb);
89 
90  while (abs (bb) != 0)
91  {
92  std::complex<FP> qq, rr;
93  divide (aa, bb, qq, rr);
94  aa = bb;
95  bb = rr;
96  }
97 
98  return aa;
99 }
100 
101 template <class T>
102 static octave_int<T>
104 {
105  T aa = a.abs ().value ();
106  T bb = b.abs ().value ();
107 
108  while (bb != 0)
109  {
110  T tt = aa % bb;
111  aa = bb;
112  bb = tt;
113  }
114 
115  return aa;
116 }
117 
118 static double
119 extended_gcd (double a, double b, double& x, double& y)
120 {
121  if (! xisinteger (a) || ! xisinteger (b))
122  (*current_liboctave_error_handler)
123  ("gcd: all values must be integers");
124 
125  double aa = fabs (a);
126  double bb = fabs (b);
127 
128  double xx, lx, yy, ly;
129  xx = 0, lx = 1;
130  yy = 1, ly = 0;
131 
132  while (bb != 0)
133  {
134  double qq = gnulib::floor (aa / bb);
135  double tt = fmod (aa, bb);
136 
137  aa = bb;
138  bb = tt;
139 
140  double tx = lx - qq*xx;
141  lx = xx;
142  xx = tx;
143 
144  double ty = ly - qq*yy;
145  ly = yy;
146  yy = ty;
147  }
148 
149  x = a >= 0 ? lx : -lx;
150  y = b >= 0 ? ly : -ly;
151 
152  return aa;
153 }
154 
155 template <typename FP>
156 static std::complex<FP>
157 extended_gcd (const std::complex<FP>& a, const std::complex<FP>& b,
158  std::complex<FP>& x, std::complex<FP>& y)
159 {
160  if (! xisinteger (a.real ()) || ! xisinteger (a.imag ())
161  || ! xisinteger (b.real ()) || ! xisinteger (b.imag ()))
162  (*current_liboctave_error_handler)
163  ("gcd: all complex parts must be integers");
164 
165  std::complex<FP> aa = a;
166  std::complex<FP> bb = b;
167  bool swapped = false;
168  if (abs (aa) < abs (bb))
169  {
170  std::swap (aa, bb);
171  swapped = true;
172  }
173 
174  std::complex<FP> xx, lx, yy, ly;
175  xx = 0, lx = 1;
176  yy = 1, ly = 0;
177 
178  while (abs(bb) != 0)
179  {
180  std::complex<FP> qq, rr;
181  divide (aa, bb, qq, rr);
182  aa = bb;
183  bb = rr;
184 
185  std::complex<FP> tx = lx - qq*xx;
186  lx = xx;
187  xx = tx;
188 
189  std::complex<FP> ty = ly - qq*yy;
190  ly = yy;
191  yy = ty;
192  }
193 
194  x = lx;
195  y = ly;
196 
197  if (swapped)
198  std::swap (x, y);
199 
200  return aa;
201 }
202 
203 template <class T>
204 static octave_int<T>
207 {
208  T aa = a.abs ().value ();
209  T bb = b.abs ().value ();
210  T xx, lx, yy, ly;
211  xx = 0, lx = 1;
212  yy = 1, ly = 0;
213 
214  while (bb != 0)
215  {
216  T qq = aa / bb;
217  T tt = aa % bb;
218  aa = bb;
219  bb = tt;
220 
221  T tx = lx - qq*xx;
222  lx = xx;
223  xx = tx;
224 
225  T ty = ly - qq*yy;
226  ly = yy;
227  yy = ty;
228  }
229 
230  x = octave_int<T> (lx) * a.signum ();
231  y = octave_int<T> (ly) * b.signum ();
232 
233  return aa;
234 }
235 
236 template<class NDA>
237 static octave_value
239 {
240  typedef typename NDA::element_type T;
241  octave_value retval;
242 
243  if (a.is_scalar_type () && b.is_scalar_type ())
244  {
245  // Optimize scalar case.
246  T aa = octave_value_extract<T> (a);
247  T bb = octave_value_extract<T> (b);
248  retval = simple_gcd (aa, bb);
249  }
250  else
251  {
252  NDA aa = octave_value_extract<NDA> (a);
253  NDA bb = octave_value_extract<NDA> (b);
254  retval = binmap<T> (aa, bb, simple_gcd, "gcd");
255  }
256 
257  return retval;
258 }
259 
260 // Dispatcher
261 static octave_value
262 do_simple_gcd (const octave_value& a, const octave_value& b)
263 {
264  octave_value retval;
266  b.builtin_type ());
267  switch (btyp)
268  {
269  case btyp_double:
270  if (a.is_sparse_type () && b.is_sparse_type ())
271  {
272  retval = do_simple_gcd<SparseMatrix> (a, b);
273  break;
274  }
275  // fall through!
276 
277  case btyp_float:
278  retval = do_simple_gcd<NDArray> (a, b);
279  break;
280 
281 #define MAKE_INT_BRANCH(X) \
282  case btyp_ ## X: \
283  retval = do_simple_gcd<X ## NDArray> (a, b); \
284  break
285 
286  MAKE_INT_BRANCH (int8);
287  MAKE_INT_BRANCH (int16);
288  MAKE_INT_BRANCH (int32);
289  MAKE_INT_BRANCH (int64);
290  MAKE_INT_BRANCH (uint8);
291  MAKE_INT_BRANCH (uint16);
292  MAKE_INT_BRANCH (uint32);
293  MAKE_INT_BRANCH (uint64);
294 
295 #undef MAKE_INT_BRANCH
296 
297  case btyp_complex:
298  retval = do_simple_gcd<ComplexNDArray> (a, b);
299  break;
300 
301  case btyp_float_complex:
302  retval = do_simple_gcd<FloatComplexNDArray> (a, b);
303  break;
304 
305  default:
306  error ("gcd: invalid class combination for gcd: %s and %s\n",
307  a.class_name ().c_str (), b.class_name ().c_str ());
308  }
309 
310  if (btyp == btyp_float)
311  retval = retval.float_array_value ();
312 
313  return retval;
314 }
315 
316 template<class NDA>
317 static octave_value
320 {
321  typedef typename NDA::element_type T;
322  octave_value retval;
323 
324  if (a.is_scalar_type () && b.is_scalar_type ())
325  {
326  // Optimize scalar case.
327  T aa = octave_value_extract<T> (a);
328  T bb = octave_value_extract<T> (b);
329  T xx, yy;
330  retval = extended_gcd (aa, bb, xx, yy);
331  x = xx;
332  y = yy;
333  }
334  else
335  {
336  NDA aa = octave_value_extract<NDA> (a);
337  NDA bb = octave_value_extract<NDA> (b);
338 
339  dim_vector dv = aa.dims ();
340  if (aa.numel () == 1)
341  dv = bb.dims ();
342  else if (bb.numel () != 1 && bb.dims () != dv)
343  gripe_nonconformant ("gcd", a.dims (), b.dims ());
344 
345  NDA gg (dv), xx (dv), yy (dv);
346 
347  const T *aptr = aa.fortran_vec ();
348  const T *bptr = bb.fortran_vec ();
349 
350  bool inca = aa.numel () != 1;
351  bool incb = bb.numel () != 1;
352 
353  T *gptr = gg.fortran_vec ();
354  T *xptr = xx.fortran_vec ();
355  T *yptr = yy.fortran_vec ();
356 
357  octave_idx_type n = gg.numel ();
358  for (octave_idx_type i = 0; i < n; i++)
359  {
360  octave_quit ();
361 
362  *gptr++ = extended_gcd (*aptr, *bptr, *xptr++, *yptr++);
363 
364  aptr += inca;
365  bptr += incb;
366  }
367 
368  x = xx;
369  y = yy;
370 
371  retval = gg;
372  }
373 
374  return retval;
375 }
376 
377 // Dispatcher
378 static octave_value
379 do_extended_gcd (const octave_value& a, const octave_value& b,
381 {
382  octave_value retval;
383 
385  b.builtin_type ());
386  switch (btyp)
387  {
388  case btyp_double:
389  case btyp_float:
390  retval = do_extended_gcd<NDArray> (a, b, x, y);
391  break;
392 
393 #define MAKE_INT_BRANCH(X) \
394  case btyp_ ## X: \
395  retval = do_extended_gcd<X ## NDArray> (a, b, x, y); \
396  break
397 
398  MAKE_INT_BRANCH (int8);
399  MAKE_INT_BRANCH (int16);
400  MAKE_INT_BRANCH (int32);
401  MAKE_INT_BRANCH (int64);
402  MAKE_INT_BRANCH (uint8);
403  MAKE_INT_BRANCH (uint16);
404  MAKE_INT_BRANCH (uint32);
405  MAKE_INT_BRANCH (uint64);
406 
407 #undef MAKE_INT_BRANCH
408 
409  case btyp_complex:
410  retval = do_extended_gcd<ComplexNDArray> (a, b, x, y);
411  break;
412 
413  case btyp_float_complex:
414  retval = do_extended_gcd<FloatComplexNDArray> (a, b, x, y);
415  break;
416 
417  default:
418  error ("gcd: invalid class combination for gcd: %s and %s\n",
419  a.class_name ().c_str (), b.class_name ().c_str ());
420  }
421 
422  // For consistency.
423  if (! error_state && a.is_sparse_type () && b.is_sparse_type ())
424  {
425  retval = retval.sparse_matrix_value ();
426  x = x.sparse_matrix_value ();
427  y = y.sparse_matrix_value ();
428  }
429 
430  if (btyp == btyp_float)
431  {
432  retval = retval.float_array_value ();
433  x = x.float_array_value ();
434  y = y.float_array_value ();
435  }
436 
437  return retval;
438 }
439 
440 DEFUN (gcd, args, nargout,
441  "-*- texinfo -*-\n\
442 @deftypefn {Built-in Function} {@var{g} =} gcd (@var{a1}, @var{a2}, @dots{})\n\
443 @deftypefnx {Built-in Function} {[@var{g}, @var{v1}, @dots{}] =} gcd (@var{a1}, @var{a2}, @dots{})\n\
444 Compute the greatest common divisor of @var{a1}, @var{a2}, @dots{}.\n\
445 \n\
446 If more than one argument is given then all arguments must be the same size\n\
447 or scalar. In this case the greatest common divisor is calculated for each\n\
448 element individually. All elements must be ordinary or Gaussian (complex)\n\
449 integers. Note that for Gaussian integers, the gcd is only unique up to a\n\
450 phase factor (multiplication by 1, -1, i, or -i), so an arbitrary greatest\n\
451 common divisor among the four possible is returned.\n\
452 \n\
453 Optional return arguments @var{v1}, @dots{}, contain integer vectors such\n\
454 that,\n\
455 \n\
456 @tex\n\
457 $g = v_1 a_1 + v_2 a_2 + \\cdots$\n\
458 @end tex\n\
459 @ifnottex\n\
460 \n\
461 @example\n\
462 @var{g} = @var{v1} .* @var{a1} + @var{v2} .* @var{a2} + @dots{}\n\
463 @end example\n\
464 \n\
465 @end ifnottex\n\
466 \n\
467 Example code:\n\
468 \n\
469 @example\n\
470 @group\n\
471 gcd ([15, 9], [20, 18])\n\
472  @result{} 5 9\n\
473 @end group\n\
474 @end example\n\
475 \n\
476 @seealso{lcm, factor, isprime}\n\
477 @end deftypefn")
478 {
479  octave_value_list retval;
480 
481  int nargin = args.length ();
482 
483  if (nargin > 1)
484  {
485  if (nargout > 1)
486  {
487  retval.resize (nargin + 1);
488 
489  retval(0) = do_extended_gcd (args(0), args(1), retval(1), retval(2));
490 
491  for (int j = 2; j < nargin; j++)
492  {
493  octave_value x;
494  retval(0) = do_extended_gcd (retval(0), args(j),
495  x, retval(j+1));
496  for (int i = 0; i < j; i++)
497  retval(i+1).assign (octave_value::op_el_mul_eq, x);
498 
499  if (error_state)
500  break;
501  }
502  }
503  else
504  {
505  retval(0) = do_simple_gcd (args(0), args(1));
506 
507  for (int j = 2; j < nargin; j++)
508  {
509  retval(0) = do_simple_gcd (retval(0), args(j));
510 
511  if (error_state)
512  break;
513  }
514  }
515  }
516  else
517  print_usage ();
518 
519  return retval;
520 }
521 
522 /*
523 %!assert (gcd (200, 300, 50, 35), 5)
524 %!assert (gcd (int16 (200), int16 (300), int16 (50), int16 (35)), int16 (5))
525 %!assert (gcd (uint64 (200), uint64 (300), uint64 (50), uint64 (35)), uint64 (5))
526 %!assert (gcd (18-i, -29+3i), -3-4i)
527 
528 %!test
529 %! p = [953 967];
530 %! u = [953 + i*971, 967 + i*977];
531 %! [d, k(1), k(2)] = gcd (p(1), p(2));
532 %! [z, w(1), w(2)] = gcd (u(1), u(2));
533 %! assert (d, 1)
534 %! assert (sum (p.*k), d)
535 %! assert (abs (z), sqrt (2))
536 %! assert (abs (sum (u.*w)), sqrt (2))
537 
538 %!error <all values must be integers> gcd (1/2, 2);
539 %!error <all complex parts must be integers> gcd (e + i*pi, 1);
540 
541 %!error gcd ()
542 
543 %!test
544 %! s.a = 1;
545 %! fail ("gcd (s)");
546 */
void gripe_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
bool is_scalar_type(void) const
Definition: ov.h:657
static octave_value do_extended_gcd(const octave_value &a, const octave_value &b, octave_value &x, octave_value &y)
Definition: gcd.cc:318
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
void error(const char *fmt,...)
Definition: error.cc:476
octave_int< T > abs() const
Definition: oct-inttypes.h:904
#define MAKE_INT_BRANCH(X)
builtin_type_t
Definition: ov-base.h:59
static void divide(const std::complex< FP > &a, const std::complex< FP > &b, std::complex< FP > &q, std::complex< FP > &r)
Definition: gcd.cc:64
static double simple_gcd(double a, double b)
Definition: gcd.cc:40
builtin_type_t btyp_mixed_numeric(builtin_type_t x, builtin_type_t y)
Definition: ov-base.cc:59
FloatNDArray float_array_value(bool frc_str_conv=false) const
Definition: ov.h:782
bool is_sparse_type(void) const
Definition: ov.h:666
int error_state
Definition: error.cc:101
octave_idx_type length(void) const
Definition: ov.cc:1525
dim_vector dims(void) const
Definition: ov.h:470
bool xisinteger(double x)
Definition: lo-mappers.h:206
octave_int< T > signum() const
Definition: oct-inttypes.h:905
std::string class_name(void) const
Definition: ov.h:1049
void resize(octave_idx_type n, const octave_value &rfv=octave_value())
Definition: oct-obj.h:93
ColumnVector imag(const ComplexColumnVector &a)
Definition: dColVector.cc:162
SparseMatrix sparse_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:820
std::complex< T > floor(const std::complex< T > &x)
Definition: lo-mappers.h:282
static octave_value do_simple_gcd(const octave_value &a, const octave_value &b)
Definition: gcd.cc:238
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:156
static double extended_gcd(double a, double b, double &x, double &y)
Definition: gcd.cc:119
T abs(T x)
Definition: pr-output.cc:3062
builtin_type_t builtin_type(void) const
Definition: ov.h:603
F77_RET_T const double * x