GNU Octave  4.2.1
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-2017 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 #if defined (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 "ovl.h"
38 
39 static double
40 simple_gcd (double a, double b)
41 {
43  error ("gcd: all values must be integers");
44 
45  double aa = fabs (a);
46  double bb = fabs (b);
47 
48  while (bb != 0)
49  {
50  double tt = fmod (aa, bb);
51  aa = bb;
52  bb = tt;
53  }
54 
55  return aa;
56 }
57 
58 // Don't use the Complex and FloatComplex typedefs because we need to
59 // refer to the actual float precision FP in the body (and when gcc
60 // implements template aliases from C++0x, can do a small fix here).
61 template <typename FP>
62 static void
63 divide (const std::complex<FP>& a, const std::complex<FP>& b,
64  std::complex<FP>& q, std::complex<FP>& r)
65 {
66  FP qr = std::floor ((a/b).real () + 0.5);
67  FP qi = std::floor ((a/b).imag () + 0.5);
68 
69  q = std::complex<FP> (qr, qi);
70 
71  r = a - q*b;
72 }
73 
74 template <typename FP>
75 static std::complex<FP>
76 simple_gcd (const std::complex<FP>& a, const std::complex<FP>& b)
77 {
78  if (! octave::math::isinteger (a.real ())
79  || ! octave::math::isinteger (a.imag ())
80  || ! octave::math::isinteger (b.real ())
81  || ! octave::math::isinteger (b.imag ()))
82  error ("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 <typename 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 {
122  error ("gcd: all values must be integers");
123 
124  double aa = fabs (a);
125  double bb = fabs (b);
126 
127  double xx, lx, yy, ly;
128  xx = 0, lx = 1;
129  yy = 1, ly = 0;
130 
131  while (bb != 0)
132  {
133  double qq = std::floor (aa / bb);
134  double tt = fmod (aa, bb);
135 
136  aa = bb;
137  bb = tt;
138 
139  double tx = lx - qq*xx;
140  lx = xx;
141  xx = tx;
142 
143  double ty = ly - qq*yy;
144  ly = yy;
145  yy = ty;
146  }
147 
148  x = a >= 0 ? lx : -lx;
149  y = b >= 0 ? ly : -ly;
150 
151  return aa;
152 }
153 
154 template <typename FP>
155 static std::complex<FP>
156 extended_gcd (const std::complex<FP>& a, const std::complex<FP>& b,
157  std::complex<FP>& x, std::complex<FP>& y)
158 {
159  if (! octave::math::isinteger (a.real ())
160  || ! octave::math::isinteger (a.imag ())
161  || ! octave::math::isinteger (b.real ())
162  || ! octave::math::isinteger (b.imag ()))
163  error ("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 <typename 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 <typename NDA>
237 static octave_value
239 {
240  typedef typename NDA::element_type T;
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 {
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);
290  MAKE_INT_BRANCH (uint8);
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 <typename NDA>
317 static octave_value
320 {
321  typedef typename NDA::element_type T;
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  octave::err_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 {
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);
402  MAKE_INT_BRANCH (uint8);
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 (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  doc: /* -*- texinfo -*-
442 @deftypefn {} {@var{g} =} gcd (@var{a1}, @var{a2}, @dots{})
443 @deftypefnx {} {[@var{g}, @var{v1}, @dots{}] =} gcd (@var{a1}, @var{a2}, @dots{})
444 Compute the greatest common divisor of @var{a1}, @var{a2}, @dots{}.
445 
446 If more than one argument is given then all arguments must be the same size
447 or scalar. In this case the greatest common divisor is calculated for each
448 element individually. All elements must be ordinary or Gaussian (complex)
449 integers. Note that for Gaussian integers, the gcd is only unique up to a
450 phase factor (multiplication by 1, -1, i, or -i), so an arbitrary greatest
451 common divisor among the four possible is returned.
452 
453 Optional return arguments @var{v1}, @dots{}, contain integer vectors such
454 that,
455 
456 @tex
457 $g = v_1 a_1 + v_2 a_2 + \cdots$
458 @end tex
459 @ifnottex
460 
461 @example
462 @var{g} = @var{v1} .* @var{a1} + @var{v2} .* @var{a2} + @dots{}
463 @end example
464 
465 @end ifnottex
466 
467 Example code:
468 
469 @example
470 @group
471 gcd ([15, 9], [20, 18])
472  @result{} 5 9
473 @end group
474 @end example
475 
476 @seealso{lcm, factor, isprime}
477 @end deftypefn */)
478 {
479  int nargin = args.length ();
480 
481  if (nargin < 2)
482  print_usage ();
483 
485 
486  if (nargout > 1)
487  {
488  retval.resize (nargin + 1);
489 
490  retval(0) = do_extended_gcd (args(0), args(1), retval(1), retval(2));
491 
492  for (int j = 2; j < nargin; j++)
493  {
494  octave_value x;
495  retval(0) = do_extended_gcd (retval(0), args(j), x, retval(j+1));
496  for (int i = 0; i < j; i++)
498  }
499  }
500  else
501  {
502  retval(0) = do_simple_gcd (args(0), args(1));
503 
504  for (int j = 2; j < nargin; j++)
505  retval(0) = do_simple_gcd (retval(0), args(j));
506  }
507 
508  return retval;
509 }
510 
511 /*
512 %!assert (gcd (200, 300, 50, 35), 5)
513 %!assert (gcd (int16 (200), int16 (300), int16 (50), int16 (35)), int16 (5))
514 %!assert (gcd (uint64 (200), uint64 (300), uint64 (50), uint64 (35)), uint64 (5))
515 %!assert (gcd (18-i, -29+3i), -3-4i)
516 
517 %!test
518 %! p = [953 967];
519 %! u = [953 + i*971, 967 + i*977];
520 %! [d, k(1), k(2)] = gcd (p(1), p(2));
521 %! [z, w(1), w(2)] = gcd (u(1), u(2));
522 %! assert (d, 1);
523 %! assert (sum (p.*k), d);
524 %! assert (abs (z), sqrt (2));
525 %! assert (abs (sum (u.*w)), sqrt (2));
526 
527 %!error <all values must be integers> gcd (1/2, 2)
528 %!error <all complex parts must be integers> gcd (e + i*pi, 1)
529 
530 %!error gcd ()
531 
532 %!test
533 %! s.a = 1;
534 %! fail ("gcd (s)");
535 */
bool isinteger(double x)
Definition: lo-mappers.h:255
OCTINTERP_API void print_usage(void)
Definition: defun.cc:52
OCTAVE_EXPORT octave_value_list uint16
Definition: ov.cc:1258
bool is_scalar_type(void) const
Definition: ov.h:673
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:46
void error(const char *fmt,...)
Definition: error.cc:570
octave_int< T > signum() const
Definition: oct-inttypes.h:919
#define MAKE_INT_BRANCH(X)
builtin_type_t
Definition: ov-base.h:61
static void divide(const std::complex< FP > &a, const std::complex< FP > &b, std::complex< FP > &q, std::complex< FP > &r)
Definition: gcd.cc:63
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:398
bool swap
Definition: load-save.cc:725
JNIEnv void * args
Definition: ov-java.cc:67
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:60
FloatNDArray float_array_value(bool frc_str_conv=false) const
Definition: ov.h:796
OCTAVE_EXPORT octave_value_list return the number of command line arguments passed to Octave If called with the optional argument the function xample nargout(@histc)
Definition: ov-usr-fcn.cc:935
bool is_sparse_type(void) const
Definition: ov.h:682
OCTAVE_EXPORT octave_value_list uint32
Definition: ov.cc:1258
octave_int< T > abs() const
Definition: oct-inttypes.h:918
int nargin
Definition: graphics.cc:10115
OCTAVE_EXPORT octave_value_list int16
Definition: ov.cc:1302
octave_value & assign(assign_op op, const std::string &type, const std::list< octave_value_list > &idx, const octave_value &rhs)
Definition: ov.cc:1556
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
octave_value retval
Definition: data.cc:6294
OCTAVE_EXPORT octave_value_list int32
Definition: ov.cc:1258
dim_vector dims(void) const
Definition: ov.h:486
OCTAVE_EXPORT octave_value_list int64
Definition: ov.cc:1258
Definition: mx-defs.h:77
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
OCTAVE_EXPORT octave_value_list return the value of the option it must match the dimension of the state and the relative tolerance must also be a vector of the same length tem it must match the dimension of the state and the absolute tolerance must also be a vector of the same length The local error test applied at each integration step is xample roup abs(local error in x(i))<
std::string class_name(void) const
Definition: ov.h:1234
the element is set to zero In other the statement xample y
Definition: data.cc:5342
b
Definition: cellfun.cc:398
void resize(octave_idx_type n, const octave_value &rfv=octave_value())
Definition: ovl.h:100
ColumnVector imag(const ComplexColumnVector &a)
Definition: dColVector.cc:142
double floor(double x)
Definition: lo-mappers.cc:330
SparseMatrix sparse_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:834
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:136
static double extended_gcd(double a, double b, double &x, double &y)
Definition: gcd.cc:119
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
builtin_type_t builtin_type(void) const
Definition: ov.h:619
dim_vector dv
Definition: sub2ind.cc:263
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE * x