24 #if defined (HAVE_CONFIG_H) 43 error (
"gcd: all values must be integers");
50 double tt = fmod (aa, bb);
61 template <
typename FP>
63 divide (
const std::complex<FP>&
a,
const std::complex<FP>&
b,
64 std::complex<FP>& q, std::complex<FP>& r)
69 q = std::complex<FP> (
qr, qi);
74 template <
typename FP>
75 static std::complex<FP>
82 error (
"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;
101 template <
typename T>
105 T aa =
a.abs ().value ();
106 T bb =
b.abs ().value ();
122 error (
"gcd: all values must be integers");
124 double aa = fabs (
a);
125 double bb = fabs (
b);
127 double xx, lx, yy, ly;
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)
163 error (
"gcd: all complex parts must be integers");
165 std::complex<FP> aa =
a;
166 std::complex<FP> bb =
b;
167 bool swapped =
false;
174 std::complex<FP> xx, lx, yy, ly;
180 std::complex<FP> qq, rr;
185 std::complex<FP> tx = lx - qq*xx;
189 std::complex<FP> ty = ly - qq*yy;
203 template <
typename T>
208 T aa =
a.abs ().value ();
209 T bb =
b.abs ().value ();
236 template <
typename NDA>
240 typedef typename NDA::element_type T;
243 if (
a.is_scalar_type () &&
b.is_scalar_type ())
246 T aa = octave_value_extract<T> (
a);
247 T bb = octave_value_extract<T> (
b);
252 NDA aa = octave_value_extract<NDA> (
a);
253 NDA bb = octave_value_extract<NDA> (
b);
270 if (
a.issparse () &&
b.issparse ())
272 retval = do_simple_gcd<SparseMatrix> (
a,
b);
278 retval = do_simple_gcd<NDArray> (
a,
b);
281 #define MAKE_INT_BRANCH(X) \ 283 retval = do_simple_gcd<X ## NDArray> (a, b); \ 295 #undef MAKE_INT_BRANCH 298 retval = do_simple_gcd<ComplexNDArray> (
a,
b);
302 retval = do_simple_gcd<FloatComplexNDArray> (
a,
b);
306 error (
"gcd: invalid class combination for gcd: %s and %s\n",
307 a.class_name ().c_str (),
b.class_name ().c_str ());
316 template <
typename NDA>
321 typedef typename NDA::element_type T;
324 if (
a.is_scalar_type () &&
b.is_scalar_type ())
327 T aa = octave_value_extract<T> (
a);
328 T bb = octave_value_extract<T> (
b);
336 NDA aa = octave_value_extract<NDA> (
a);
337 NDA bb = octave_value_extract<NDA> (
b);
340 if (aa.numel () == 1)
342 else if (bb.numel () != 1 && bb.dims () !=
dv)
345 NDA gg (
dv), xx (
dv), yy (
dv);
347 const T *aptr = aa.fortran_vec ();
348 const T *bptr = bb.fortran_vec ();
350 bool inca = aa.numel () != 1;
351 bool incb = bb.numel () != 1;
353 T *gptr = gg.fortran_vec ();
354 T *xptr = xx.fortran_vec ();
355 T *yptr = yy.fortran_vec ();
362 *gptr++ =
extended_gcd (*aptr, *bptr, *xptr++, *yptr++);
390 retval = do_extended_gcd<NDArray> (
a,
b,
x,
y);
393 #define MAKE_INT_BRANCH(X) \ 395 retval = do_extended_gcd<X ## NDArray> (a, b, x, y); \ 407 #undef MAKE_INT_BRANCH 410 retval = do_extended_gcd<ComplexNDArray> (
a,
b,
x,
y);
414 retval = do_extended_gcd<FloatComplexNDArray> (
a,
b,
x,
y);
418 error (
"gcd: invalid class combination for gcd: %s and %s\n",
419 a.class_name ().c_str (),
b.class_name ().c_str ());
423 if (
a.issparse () &&
b.issparse ())
426 x =
x.sparse_matrix_value ();
427 y =
y.sparse_matrix_value ();
433 x =
x.float_array_value ();
434 y =
y.float_array_value ();
479 int nargin = args.length ();
492 for (
int j = 2; j <
nargin; j++)
496 for (
int i = 0;
i < j;
i++)
504 for (
int j = 2; j <
nargin; j++)
OCTINTERP_API void print_usage(void)
static octave_value do_extended_gcd(const octave_value &a, const octave_value &b, octave_value &x, octave_value &y)
#define DEFUN(name, args_name, nargout_name, doc)
Macro to define a builtin function.
void error(const char *fmt,...)
std::complex< T > floor(const std::complex< T > &x)
FloatNDArray float_array_value(bool frc_str_conv=false) const
#define MAKE_INT_BRANCH(X)
static void divide(const std::complex< FP > &a, const std::complex< FP > &b, std::complex< FP > &q, std::complex< FP > &r)
octave_value resize(const dim_vector &dv, bool fill=false) const
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
static double simple_gcd(double a, double b)
builtin_type_t btyp_mixed_numeric(builtin_type_t x, builtin_type_t y)
Determine the resulting type for a possible mixed-type operation.
octave_value & assign(assign_op op, const std::string &type, const std::list< octave_value_list > &idx, const octave_value &rhs)
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)
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
SparseMatrix sparse_matrix_value(bool frc_str_conv=false) const
the element is set to zero In other the statement xample y
ColumnVector imag(const ComplexColumnVector &a)
static octave_value do_simple_gcd(const octave_value &a, const octave_value &b)
ColumnVector real(const ComplexColumnVector &a)
static double extended_gcd(double a, double b, double &x, double &y)
Vector representing the dimensions (size) of an Array.
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 const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T const F77_DBLE * x