From octave-maintainers-request at octave dot org Wed Apr 14 08:08:13 2004 Subject: GCD -> NDArrays, Matlab compatiable? From: David Bateman To: maintainers at octave dot org Date: Wed, 14 Apr 2004 15:04:15 +0200 --rJwd6BRFiFCcLxzm Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In my efforts to convert things to be compatiable with NDArrays, I picked on the gcd.m and lcm.m files. Basically not only are these candidates for conversion for use with NDArrays, but they are also not comaptiable with their equivalent matlab versions at the moment. Also when considering the conversion for NDArrays I saw no way to avoid for-loops in the dot-m files. Therefore I tried converting the gcd function to be an oct-file and will rewrite lcm.m to use this. I've attemptted to maintain compatiably with both the old octave version and the matlab code, but it goes beyond both. Like the old octave code any number of input args is valid, except now the args can be integer scalars or arrays. Likewise if all args are scalar or a single argument is given, then the previous behaviour of "[g,v] = gcd(...)" will be obtained. However, like the matlab version the output variable "v" can now be broken down into its components like "[g,v1,v2,...,vk] = gcd(a1,a2,...,ak)", this holds even in a case like "[g,v1,v2] = gcd(15,20)". In this manner compatiablity with both the old octave code and matlab is achieved. Does anyone see any objection to removing the old gcd.m code and using the attached oct-file instead? If not, John would you accept this if I rewrote it as a patch? Cheers David -- David Bateman David dot Bateman at motorola dot com Motorola CRM +33 1 69 35 48 04 (Ph) Parc Les Algorithmes, Commune de St Aubin +33 1 69 35 77 01 (Fax) 91193 Gif-Sur-Yvette FRANCE The information contained in this communication has been classified as: [x] General Business Information [ ] Motorola Internal Use Only [ ] Motorola Confidential Proprietary --rJwd6BRFiFCcLxzm Content-Type: text/plain; charset=us-ascii Content-Disposition: attachment; filename="gcd.cc" /* Copyright (C) 2004 David Bateman This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with Octave; see the file COPYING. If not, write to the Free Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. In addition to the terms of the GPL, you are permitted to link this program with any Open Source program, as defined by the Open Source Initiative (www.opensource.org) */ #include #include DEFUN_DLD (gcd, args, nargout, "-*- texinfo -*-\n\ at deftypefn {Loadable Function} {@var{g} =} gcd (@var{a1}, @code{...})\n\ at deftypefnx {Loadable Function} {[@var{g}, @var{v1}, @var{...}] =} gcd (@var{a1}, @code{...})\n\ \n\ If a single argument is given then compute the greatest common divisor of\n\ the elements of this argument. Otherwise if more than one argument is\n\ given all arguments must be the same size or scalar. In this case the\n\ greatest common divisor is calculated for element individually. All\n\ elements must be integers. For example,\n\ \n\ at example\n\ at group\n\ gcd ([15, 20])\n\ at result{} 5\n\ at end group\n\ at end example\n\ \n\ at noindent\n\ and\n\ \n\ at example\n\ at group\n\ gcd ([15, 9], [20 18])\n\ at result{} 5 9\n\ at end group\n\ at end example\n\ \n\ Optional return arguments at var{v1}, etc, contain integer vectors such\n\ that,\n\ \n\ at ifinfo\n\ at example\n\ at var{g} = @var{v1} .* @var{a1} + @var{v2} .* @var{a2} + @var{...}\n\ at end example\n\ at end ifinfo\n\ at iftex\n\ at tex\n\ $g = v_1 a_1 + v_2 a_2 + \cdots$\n\ at end tex\n\ at end iftex\n\ \n\ For backward compatiability with previous versions of this function, when\n\ all arguments are scalr, a single return argument at var{v1} containing\n\ all of the values of at var{v1}, @var{...} is acceptable.\n\ at end deftypefn\n\ at seealso{lcm, min, max, ceil, and floor}") { octave_value_list retval; int nargin = args.length (); bool all_args_scalar = true; dim_vector dv(1); for (int i = 0; i < nargin; i++) if (! args (i).is_scalar_type ()) { if (! args (i).is_matrix_type ()) { error ("gcd: invalid argument type"); return retval; } if (all_args_scalar) { all_args_scalar = false; dv = args (i).dims (); } else { if (dv != args (i).dims ()) { error ("gcd: all arguments must be the same size or scalar"); return retval; } } } if (nargin == 1) { NDArray gg = args (0). array_value (); int nel = dv.numel (); NDArray v (dv); RowVector x (3); RowVector y (3); double g = std::abs(gg(0)); if (g != (double)((int)g)) { error ("gcd: all arguments must be integer"); return retval; } v(0) = signum (gg(0)); for (int k = 1; k < nel; k++) { x(0) = g; x(1) = 1; x(2) = 0; y(0) = std::abs(gg(k)); y(1) = 0; y(2) = 1; if (y(0) != (double)((int) y(0))) { error ("gcd: all arguments must be integer"); return retval; } while (y(0) > 0) { RowVector r = x - y * ((int) ( x(0) / y(0))); x = y; y = r; } g = x(0); for (int i = 0; i < k; i++) v(i) *= x(1); v(k) = x(2) * signum (gg(k)); } retval (1) = v; retval (0) = g; } else if (all_args_scalar && nargout < 3) { double g = args (0).int_value (true); if (error_state) { error ("gcd: all arguments must be integer"); return retval; } RowVector v (nargin, 0); RowVector x (3); RowVector y (3); v(0) = signum (g); g = std::abs(g); for (int k = 1; k < nargin; k++) { x(0) = g; x(1) = 1; x(2) = 0; y(0) = args (k).int_value (true); y(1) = 0; y(2) = 1; double sgn = signum (y(0)); y(0) = std::abs(y(0)); if (error_state) { error ("gcd: all arguments must be integer"); return retval; } while (y(0) > 0) { RowVector r = x - y * ((int) ( x(0) / y(0))); x = y; y = r; } g = x(0); for (int i = 0; i < k; i++) v(i) *= x(1); v(k) = x(2) * sgn; } retval (1) = v; retval (0) = g; } else { NDArray g = args (0).array_value (); NDArray v[nargin]; int nel = dv.numel (); v[0].resize(dv); for (int i = 0; i < nel; i++) { v[0](i) = signum (g(i)); g(i) = std::abs (g(i)); if (g(i) != (double)((int) g(i))) { error ("gcd: all arguments must be integer"); return retval; } } RowVector x (3); RowVector y (3); for (int k = 1; k < nargin; k++) { NDArray gnew = args (k).array_value (); v[k].resize(dv); for (int n = 0; n < nel; n++) { x(0) = g(n); x(1) = 1; x(2) = 0; y(0) = std::abs(gnew(n)); y(1) = 0; y(2) = 1; if (y(0) != (double)((int) y(0))) { error ("gcd: all arguments must be integer"); return retval; } while (y(0) > 0) { RowVector r = x - y * ((int) ( x(0) / y(0))); x = y; y = r; } g(n) = x(0); for (int i = 0; i < k; i++) v[i](n) *= x(1); v[k](n) = x(2) * signum (gnew(n)); } } for (int k = 0; k < nargin; k++) retval (1+k) = v[k]; retval (0) = g; } return retval; } /* ;;; Local Variables: *** ;;; mode: C++ *** ;;; End: *** */ --rJwd6BRFiFCcLxzm--