43 (*current_liboctave_error_handler)
44 (
"gcd: all values must be integers");
51 double tt = fmod (aa, bb);
62 template <
typename FP>
64 divide (
const std::complex<FP>& a,
const std::complex<FP>& b,
65 std::complex<FP>& q, std::complex<FP>& r)
70 q = std::complex<FP> (qr, qi);
75 template <
typename FP>
76 static std::complex<FP>
77 simple_gcd (
const std::complex<FP>& a,
const std::complex<FP>& b)
81 (*current_liboctave_error_handler)
82 (
"gcd: all complex parts must be integers");
84 std::complex<FP> aa = a;
85 std::complex<FP> bb = b;
92 std::complex<FP> qq, rr;
105 T aa = a.
abs ().value ();
106 T bb = b.
abs ().value ();
122 (*current_liboctave_error_handler)
123 (
"gcd: all values must be integers");
125 double aa = fabs (a);
126 double bb = fabs (b);
128 double xx = 0, yy = 1;
129 double lx = 1, ly = 0;
134 double tt = fmod (aa, bb);
139 double tx = lx - qq*xx;
143 double ty = ly - qq*yy;
148 x = a >= 0 ? lx : -lx;
149 y = b >= 0 ? ly : -ly;
154 template <
typename FP>
155 static std::complex<FP>
157 std::complex<FP>&
x, std::complex<FP>& y)
161 (*current_liboctave_error_handler)
162 (
"gcd: all complex parts must be integers");
164 std::complex<FP> aa = a, bb = b;
165 bool swapped =
false;
172 std::complex<FP> xx = 0, lx = 1;
173 std::complex<FP> yy = 1, ly = 0;
177 std::complex<FP> qq, rr;
182 std::complex<FP> tx = lx - qq*xx;
186 std::complex<FP> ty = ly - qq*yy;
205 T aa = a.
abs ().value ();
206 T bb = b.
abs ().value ();
236 typedef typename NDA::element_type T;
242 T aa = octave_value_extract<T> (a);
243 T bb = octave_value_extract<T> (b);
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");
268 retval = do_simple_gcd<SparseMatrix> (a, b);
274 retval = do_simple_gcd<NDArray> (a, b);
277 #define MAKE_INT_BRANCH(X) \
279 retval = do_simple_gcd<X ## NDArray> (a, b); \
291 #undef MAKE_INT_BRANCH
294 retval = do_simple_gcd<ComplexNDArray> (a, b);
298 retval = do_simple_gcd<FloatComplexNDArray> (a, b);
302 error (
"gcd: invalid class combination for gcd: %s and %s\n",
303 a.class_name ().c_str (), b.
class_name ().c_str ());
317 typedef typename NDA::element_type T;
323 T aa = octave_value_extract<T> (a);
324 T bb = octave_value_extract<T> (b);
332 NDA aa = octave_value_extract<NDA> (a);
333 NDA bb = octave_value_extract<NDA> (b);
336 if (aa.numel () == 1)
338 else if (bb.numel () != 1 && bb.dims () != dv)
341 NDA gg (dv), xx (dv), yy (dv);
343 const T *aptr = aa.fortran_vec ();
344 const T *bptr = bb.fortran_vec ();
346 bool inca = aa.numel () != 1;
347 bool incb = bb.numel () != 1;
349 T *gptr = gg.fortran_vec ();
350 T *xptr = xx.fortran_vec (), *yptr = yy.fortran_vec ();
357 *gptr++ =
extended_gcd (*aptr, *bptr, *xptr++, *yptr++);
385 retval = do_extended_gcd<NDArray> (a, b,
x, y);
388 #define MAKE_INT_BRANCH(X) \
390 retval = do_extended_gcd<X ## NDArray> (a, b, x, y); \
402 #undef MAKE_INT_BRANCH
405 retval = do_extended_gcd<ComplexNDArray> (a, b,
x, y);
409 retval = do_extended_gcd<FloatComplexNDArray> (a, b,
x, y);
413 error (
"gcd: invalid class combination for gcd: %s and %s\n",
414 a.class_name ().c_str (), b.class_name ().c_str ());
418 if (!
error_state && a.is_sparse_type () && b.is_sparse_type ())
421 x = x.sparse_matrix_value ();
428 x = x.float_array_value ();
435 DEFUN (gcd, args, nargout,
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\
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\
452 gcd ([15, 9], [20, 18])\n\
457 Optional return arguments @var{v1}, etc., contain integer vectors such\n\
461 $g = v_1 a_1 + v_2 a_2 + \\cdots$\n\
466 @var{g} = @var{v1} .* @var{a1} + @var{v2} .* @var{a2} + @dots{}\n\
471 @seealso{lcm, factor}\n\
476 int nargin = args.
length ();
482 retval.
resize (nargin + 1);
486 for (
int j = 2; j < nargin; j++)
491 for (
int i = 0; i < j; i++)
502 for (
int j = 2; j < nargin; j++)