GNU Octave  3.8.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-2013 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 = 0, yy = 1;
129  double lx = 1, ly = 0;
130 
131  while (bb != 0)
132  {
133  double qq = gnulib::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 (! xisinteger (a.real ()) || ! xisinteger (a.imag ())
160  || ! xisinteger (b.real ()) || ! xisinteger (b.imag ()))
161  (*current_liboctave_error_handler)
162  ("gcd: all complex parts must be integers");
163 
164  std::complex<FP> aa = a, bb = b;
165  bool swapped = false;
166  if (abs (aa) < abs (bb))
167  {
168  std::swap (aa, bb);
169  swapped = true;
170  }
171 
172  std::complex<FP> xx = 0, lx = 1;
173  std::complex<FP> yy = 1, ly = 0;
174 
175  while (abs(bb) != 0)
176  {
177  std::complex<FP> qq, rr;
178  divide (aa, bb, qq, rr);
179  aa = bb;
180  bb = rr;
181 
182  std::complex<FP> tx = lx - qq*xx;
183  lx = xx;
184  xx = tx;
185 
186  std::complex<FP> ty = ly - qq*yy;
187  ly = yy;
188  yy = ty;
189  }
190 
191  x = lx;
192  y = ly;
193 
194  if (swapped)
195  std::swap (x, y);
196 
197  return aa;
198 }
199 
200 template <class T>
201 static octave_int<T>
204 {
205  T aa = a.abs ().value ();
206  T bb = b.abs ().value ();
207  T xx = 0, lx = 1;
208  T yy = 1, ly = 0;
209 
210  while (bb != 0)
211  {
212  T qq = aa / bb;
213  T tt = aa % bb;
214  aa = bb;
215  bb = tt;
216 
217  T tx = lx - qq*xx;
218  lx = xx;
219  xx = tx;
220 
221  T ty = ly - qq*yy;
222  ly = yy;
223  yy = ty;
224  }
225 
226  x = octave_int<T> (lx) * a.signum ();
227  y = octave_int<T> (ly) * b.signum ();
228 
229  return aa;
230 }
231 
232 template<class NDA>
233 static octave_value
235 {
236  typedef typename NDA::element_type T;
237  octave_value retval;
238 
239  if (a.is_scalar_type () && b.is_scalar_type ())
240  {
241  // Optimize scalar case.
242  T aa = octave_value_extract<T> (a);
243  T bb = octave_value_extract<T> (b);
244  retval = simple_gcd (aa, bb);
245  }
246  else
247  {
248  NDA aa = octave_value_extract<NDA> (a);
249  NDA bb = octave_value_extract<NDA> (b);
250  retval = binmap<T> (aa, bb, simple_gcd, "gcd");
251  }
252 
253  return retval;
254 }
255 
256 // Dispatcher
257 static octave_value
258 do_simple_gcd (const octave_value& a, const octave_value& b)
259 {
260  octave_value retval;
262  b.builtin_type ());
263  switch (btyp)
264  {
265  case btyp_double:
266  if (a.is_sparse_type () && b.is_sparse_type ())
267  {
268  retval = do_simple_gcd<SparseMatrix> (a, b);
269  break;
270  }
271  // fall through!
272 
273  case btyp_float:
274  retval = do_simple_gcd<NDArray> (a, b);
275  break;
276 
277 #define MAKE_INT_BRANCH(X) \
278  case btyp_ ## X: \
279  retval = do_simple_gcd<X ## NDArray> (a, b); \
280  break
281 
282  MAKE_INT_BRANCH (int8);
283  MAKE_INT_BRANCH (int16);
284  MAKE_INT_BRANCH (int32);
285  MAKE_INT_BRANCH (int64);
286  MAKE_INT_BRANCH (uint8);
287  MAKE_INT_BRANCH (uint16);
288  MAKE_INT_BRANCH (uint32);
289  MAKE_INT_BRANCH (uint64);
290 
291 #undef MAKE_INT_BRANCH
292 
293  case btyp_complex:
294  retval = do_simple_gcd<ComplexNDArray> (a, b);
295  break;
296 
297  case btyp_float_complex:
298  retval = do_simple_gcd<FloatComplexNDArray> (a, b);
299  break;
300 
301  default:
302  error ("gcd: invalid class combination for gcd: %s and %s\n",
303  a.class_name ().c_str (), b.class_name ().c_str ());
304  }
305 
306  if (btyp == btyp_float)
307  retval = retval.float_array_value ();
308 
309  return retval;
310 }
311 
312 template<class NDA>
313 static octave_value
316 {
317  typedef typename NDA::element_type T;
318  octave_value retval;
319 
320  if (a.is_scalar_type () && b.is_scalar_type ())
321  {
322  // Optimize scalar case.
323  T aa = octave_value_extract<T> (a);
324  T bb = octave_value_extract<T> (b);
325  T xx, yy;
326  retval = extended_gcd (aa, bb, xx, yy);
327  x = xx;
328  y = yy;
329  }
330  else
331  {
332  NDA aa = octave_value_extract<NDA> (a);
333  NDA bb = octave_value_extract<NDA> (b);
334 
335  dim_vector dv = aa.dims ();
336  if (aa.numel () == 1)
337  dv = bb.dims ();
338  else if (bb.numel () != 1 && bb.dims () != dv)
339  gripe_nonconformant ("gcd", a.dims (), b.dims ());
340 
341  NDA gg (dv), xx (dv), yy (dv);
342 
343  const T *aptr = aa.fortran_vec ();
344  const T *bptr = bb.fortran_vec ();
345 
346  bool inca = aa.numel () != 1;
347  bool incb = bb.numel () != 1;
348 
349  T *gptr = gg.fortran_vec ();
350  T *xptr = xx.fortran_vec (), *yptr = yy.fortran_vec ();
351 
352  octave_idx_type n = gg.numel ();
353  for (octave_idx_type i = 0; i < n; i++)
354  {
355  octave_quit ();
356 
357  *gptr++ = extended_gcd (*aptr, *bptr, *xptr++, *yptr++);
358 
359  aptr += inca;
360  bptr += incb;
361  }
362 
363  x = xx;
364  y = yy;
365 
366  retval = gg;
367  }
368 
369  return retval;
370 }
371 
372 // Dispatcher
373 static octave_value
374 do_extended_gcd (const octave_value& a, const octave_value& b,
376 {
377  octave_value retval;
378 
380  b.builtin_type ());
381  switch (btyp)
382  {
383  case btyp_double:
384  case btyp_float:
385  retval = do_extended_gcd<NDArray> (a, b, x, y);
386  break;
387 
388 #define MAKE_INT_BRANCH(X) \
389  case btyp_ ## X: \
390  retval = do_extended_gcd<X ## NDArray> (a, b, x, y); \
391  break
392 
393  MAKE_INT_BRANCH (int8);
394  MAKE_INT_BRANCH (int16);
395  MAKE_INT_BRANCH (int32);
396  MAKE_INT_BRANCH (int64);
397  MAKE_INT_BRANCH (uint8);
398  MAKE_INT_BRANCH (uint16);
399  MAKE_INT_BRANCH (uint32);
400  MAKE_INT_BRANCH (uint64);
401 
402 #undef MAKE_INT_BRANCH
403 
404  case btyp_complex:
405  retval = do_extended_gcd<ComplexNDArray> (a, b, x, y);
406  break;
407 
408  case btyp_float_complex:
409  retval = do_extended_gcd<FloatComplexNDArray> (a, b, x, y);
410  break;
411 
412  default:
413  error ("gcd: invalid class combination for gcd: %s and %s\n",
414  a.class_name ().c_str (), b.class_name ().c_str ());
415  }
416 
417  // For consistency.
418  if (! error_state && a.is_sparse_type () && b.is_sparse_type ())
419  {
420  retval = retval.sparse_matrix_value ();
421  x = x.sparse_matrix_value ();
422  y = y.sparse_matrix_value ();
423  }
424 
425  if (btyp == btyp_float)
426  {
427  retval = retval.float_array_value ();
428  x = x.float_array_value ();
429  y = y.float_array_value ();
430  }
431 
432  return retval;
433 }
434 
435 DEFUN (gcd, args, nargout,
436  "-*- texinfo -*-\n\
437 @deftypefn {Built-in Function} {@var{g} =} gcd (@var{a1}, @var{a2}, @dots{})\n\
438 @deftypefnx {Built-in Function} {[@var{g}, @var{v1}, @dots{}] =} gcd (@var{a1}, @var{a2}, @dots{})\n\
439 \n\
440 Compute the greatest common divisor of @var{a1}, @var{a2}, @dots{}. If more\n\
441 than one argument is given all arguments must be the same size or scalar.\n\
442 In this case the greatest common divisor is calculated for each element\n\
443 individually. All elements must be ordinary or Gaussian (complex)\n\
444 integers. Note that for Gaussian integers, the gcd is not unique up to\n\
445 units (multiplication by 1, -1, @var{i} or -@var{i}), so an arbitrary\n\
446 greatest common divisor amongst four possible is returned.\n\
447 \n\
448 Example code:\n\
449 \n\
450 @example\n\
451 @group\n\
452 gcd ([15, 9], [20, 18])\n\
453  @result{} 5 9\n\
454 @end group\n\
455 @end example\n\
456 \n\
457 Optional return arguments @var{v1}, etc., contain integer vectors such\n\
458 that,\n\
459 \n\
460 @tex\n\
461 $g = v_1 a_1 + v_2 a_2 + \\cdots$\n\
462 @end tex\n\
463 @ifnottex\n\
464 \n\
465 @example\n\
466 @var{g} = @var{v1} .* @var{a1} + @var{v2} .* @var{a2} + @dots{}\n\
467 @end example\n\
468 \n\
469 @end ifnottex\n\
470 \n\
471 @seealso{lcm, factor}\n\
472 @end deftypefn")
473 {
474  octave_value_list retval;
475 
476  int nargin = args.length ();
477 
478  if (nargin > 1)
479  {
480  if (nargout > 1)
481  {
482  retval.resize (nargin + 1);
483 
484  retval(0) = do_extended_gcd (args(0), args(1), retval(1), retval(2));
485 
486  for (int j = 2; j < nargin; j++)
487  {
488  octave_value x;
489  retval(0) = do_extended_gcd (retval(0), args(j),
490  x, retval(j+1));
491  for (int i = 0; i < j; i++)
492  retval(i+1).assign (octave_value::op_el_mul_eq, x);
493 
494  if (error_state)
495  break;
496  }
497  }
498  else
499  {
500  retval(0) = do_simple_gcd (args(0), args(1));
501 
502  for (int j = 2; j < nargin; j++)
503  {
504  retval(0) = do_simple_gcd (retval(0), args(j));
505 
506  if (error_state)
507  break;
508  }
509  }
510  }
511  else
512  print_usage ();
513 
514  return retval;
515 }
516 
517 /*
518 %!assert (gcd (200, 300, 50, 35), 5)
519 %!assert (gcd (int16 (200), int16 (300), int16 (50), int16 (35)), int16 (5))
520 %!assert (gcd (uint64 (200), uint64 (300), uint64 (50), uint64 (35)), uint64 (5))
521 %!assert (gcd (18-i, -29+3i), -3-4i)
522 
523 %!error gcd ()
524 
525 %!test
526 %! s.a = 1;
527 %! fail ("gcd (s)");
528 */