# HG changeset patch # User John W. Eaton # Date 1285838111 14400 # Node ID 61a6c99d8ff791a4ba9fd767b70833d78b222e8b # Parent e378c0176a38203b1babb2c792a18e2565376d1a gcd.cc: style fixes diff -r e378c0176a38 -r 61a6c99d8ff7 src/ChangeLog --- a/src/ChangeLog Thu Sep 30 05:11:59 2010 -0400 +++ b/src/ChangeLog Thu Sep 30 05:15:11 2010 -0400 @@ -1,9 +1,14 @@ +2010-09-30 John W. Eaton + + * DLD-FUNCTIONS/gcd.cc: Style fixes. + 2010-09-30 John W. Eaton * oct-errno.cc.in (octave_errno::do_list, octave_errno::list): Use octave_scalar_map instead of Octave_map. 2010-09-30 Jordi GutiƩrrez Hermoso + * DLD-FUNCTIONS/gcd.cc (divide): New function, complex integer division with remainder. (simple_gcd): Overload for complex values. diff -r e378c0176a38 -r 61a6c99d8ff7 src/DLD-FUNCTIONS/gcd.cc --- a/src/DLD-FUNCTIONS/gcd.cc Thu Sep 30 05:11:59 2010 -0400 +++ b/src/DLD-FUNCTIONS/gcd.cc Thu Sep 30 05:15:11 2010 -0400 @@ -43,11 +43,14 @@ (*current_liboctave_error_handler) ("gcd: all values must be integers"); - double aa = fabs (a), bb = fabs (b); + double aa = fabs (a); + double bb = fabs (b); + while (bb != 0) { double tt = fmod (aa, bb); - aa = bb; bb = tt; + aa = bb; + bb = tt; } return aa; @@ -61,9 +64,11 @@ divide(const std::complex& a, const std::complex& b, std::complex& q, std::complex& r) { - FP qr = floor((a/b).real() + 0.5); - FP qi = floor((a/b).imag() + 0.5); - q = std::complex(qr,qi); + FP qr = floor ((a/b).real () + 0.5); + FP qi = floor ((a/b).imag () + 0.5); + + q = std::complex (qr, qi); + r = a - q*b; } @@ -71,19 +76,18 @@ static std::complex simple_gcd (const std::complex& a, const std::complex& b) { - if ( ! xisinteger (a.real ()) || ! xisinteger(a.imag ()) || - ! xisinteger (b.real ()) || ! xisinteger(b.imag ()) ) + if (! xisinteger (a.real ()) || ! xisinteger(a.imag ()) + || ! xisinteger (b.real ()) || ! xisinteger(b.imag ())) (*current_liboctave_error_handler) ("gcd: all complex parts must be integers"); - std::complex aa = a, bb = b; + std::complex aa = a; + std::complex bb = b; - if ( abs(aa) < abs(bb) ) - { - std::swap(aa,bb); - } + if (abs (aa) < abs (bb)) + std::swap (aa, bb); - while ( abs(bb) != 0) + while (abs (bb) != 0) { std::complex qq, rr; divide (aa, bb, qq, rr); @@ -98,11 +102,14 @@ static octave_int simple_gcd (const octave_int& a, const octave_int& b) { - T aa = a.abs ().value (), bb = b.abs ().value (); + T aa = a.abs ().value (); + T bb = b.abs ().value (); + while (bb != 0) { T tt = aa % bb; - aa = bb; bb = tt; + aa = bb; + bb = tt; } return aa; @@ -115,19 +122,27 @@ (*current_liboctave_error_handler) ("gcd: all values must be integers"); - double aa = fabs (a), bb = fabs (b); - double xx = 0, lx = 1, yy = 1, ly = 0; + double aa = fabs (a); + double bb = fabs (b); + + double xx = 0, yy = 1; + double lx = 1, ly = 0; + while (bb != 0) { double qq = floor (aa / bb); double tt = fmod (aa, bb); - aa = bb; bb = tt; + + aa = bb; + bb = tt; double tx = lx - qq*xx; - lx = xx; xx = tx; + lx = xx; + xx = tx; double ty = ly - qq*yy; - ly = yy; yy = ty; + ly = yy; + yy = ty; } x = a >= 0 ? lx : -lx; @@ -141,39 +156,43 @@ extended_gcd (const std::complex& a, const std::complex& b, std::complex& x, std::complex& y) { - if ( ! xisinteger (a.real ()) || ! xisinteger(a.imag ()) || - ! xisinteger (b.real ()) || ! xisinteger(b.imag ()) ) + if (! xisinteger (a.real ()) || ! xisinteger(a.imag ()) + || ! xisinteger (b.real ()) || ! xisinteger(b.imag ())) (*current_liboctave_error_handler) ("gcd: all complex parts must be integers"); std::complex aa = a, bb = b; bool swapped = false; - if ( abs(aa) < abs(bb) ) + if (abs (aa) < abs (bb)) { - std::swap(aa,bb); + std::swap (aa, bb); swapped = true; } - std::complex xx = 0, lx = 1, yy = 1, ly = 0; - while ( abs(bb) != 0) + std::complex xx = 0, lx = 1; + std::complex yy = 1, ly = 0; + + while (abs(bb) != 0) { std::complex qq, rr; divide (aa, bb, qq, rr); - aa = bb; bb = rr; + aa = bb; + bb = rr; std::complex tx = lx - qq*xx; - lx = xx; xx = tx; + lx = xx; + xx = tx; std::complex ty = ly - qq*yy; - ly = yy; yy = ty; + ly = yy; + yy = ty; } x = lx; y = ly; + if (swapped) - { - std::swap(x,y); - } + std::swap (x, y); return aa; } @@ -183,19 +202,25 @@ extended_gcd (const octave_int& a, const octave_int& b, octave_int& x, octave_int& y) { - T aa = a.abs ().value (), bb = b.abs ().value (); - T xx = 0, lx = 1, yy = 1, ly = 0; + T aa = a.abs ().value (); + T bb = b.abs ().value (); + T xx = 0, lx = 1; + T yy = 1, ly = 0; + while (bb != 0) { T qq = aa / bb; T tt = aa % bb; - aa = bb; bb = tt; + aa = bb; + bb = tt; T tx = lx - qq*xx; - lx = xx; xx = tx; + lx = xx; + xx = tx; T ty = ly - qq*yy; - ly = yy; yy = ty; + ly = yy; + yy = ty; } x = octave_int (lx) * a.signum (); @@ -243,14 +268,17 @@ retval = do_simple_gcd (a, b); break; } - else { } // fall through! + // fall through! + case btyp_float: retval = do_simple_gcd (a, b); break; + #define MAKE_INT_BRANCH(X) \ case btyp_ ## X: \ retval = do_simple_gcd (a, b); \ break + MAKE_INT_BRANCH (int8); MAKE_INT_BRANCH (int16); MAKE_INT_BRANCH (int32); @@ -259,13 +287,17 @@ MAKE_INT_BRANCH (uint16); MAKE_INT_BRANCH (uint32); MAKE_INT_BRANCH (uint64); + #undef MAKE_INT_BRANCH + case btyp_complex: - retval = do_simple_gcd (a,b); + retval = do_simple_gcd (a, b); break; + case btyp_float_complex: - retval = do_simple_gcd (a,b); + retval = do_simple_gcd (a, b); break; + default: error ("gcd: invalid class combination for gcd: %s and %s\n", a.class_name ().c_str (), b.class_name ().c_str ()); @@ -308,9 +340,11 @@ NDA gg (dv), xx (dv), yy (dv); - const T *aptr = aa.fortran_vec (), *bptr = bb.fortran_vec (); + const T *aptr = aa.fortran_vec (); + const T *bptr = bb.fortran_vec (); - bool inca = aa.numel () != 1, incb = bb.numel () != 1; + bool inca = aa.numel () != 1; + bool incb = bb.numel () != 1; T *gptr = gg.fortran_vec (); T *xptr = xx.fortran_vec (), *yptr = yy.fortran_vec (); @@ -321,12 +355,14 @@ octave_quit (); *gptr++ = extended_gcd (*aptr, *bptr, *xptr++, *yptr++); + aptr += inca; bptr += incb; } x = xx; y = yy; + retval = gg; } @@ -339,6 +375,7 @@ octave_value& x, octave_value& y) { octave_value retval; + builtin_type_t btyp = btyp_mixed_numeric (a.builtin_type (), b.builtin_type ()); switch (btyp) @@ -347,10 +384,12 @@ case btyp_float: retval = do_extended_gcd (a, b, x, y); break; + #define MAKE_INT_BRANCH(X) \ case btyp_ ## X: \ retval = do_extended_gcd (a, b, x, y); \ break + MAKE_INT_BRANCH (int8); MAKE_INT_BRANCH (int16); MAKE_INT_BRANCH (int32); @@ -359,13 +398,17 @@ MAKE_INT_BRANCH (uint16); MAKE_INT_BRANCH (uint32); MAKE_INT_BRANCH (uint64); + #undef MAKE_INT_BRANCH + case btyp_complex: retval = do_extended_gcd (a, b, x, y); break; + case btyp_float_complex: retval = do_extended_gcd (a, b, x, y); break; + default: error ("gcd: invalid class combination for gcd: %s and %s\n", a.class_name ().c_str (), b.class_name ().c_str ()); @@ -389,7 +432,6 @@ return retval; } - DEFUN_DLD (gcd, args, nargout, "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {@var{g} =} gcd (@var{a1}, @var{a2}, @dots{})\n\ @@ -431,6 +473,7 @@ @end deftypefn") { octave_value_list retval; + int nargin = args.length (); if (nargin > 1) @@ -438,7 +481,9 @@ if (nargout > 1) { retval.resize (nargin + 1); + retval(0) = do_extended_gcd (args(0), args(1), retval(1), retval(2)); + for (int j = 2; j < nargin; j++) { octave_value x; @@ -446,6 +491,7 @@ x, retval(j+1)); for (int i = 0; i < j; i++) retval(i+1).assign (octave_value::op_el_mul_eq, x); + if (error_state) break; } @@ -453,9 +499,11 @@ else { retval(0) = do_simple_gcd (args(0), args(1)); + for (int j = 2; j < nargin; j++) { retval(0) = do_simple_gcd (retval(0), args(j)); + if (error_state) break; }