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