|
4864
|
1 /* |
|
|
2 |
|
7017
|
3 Copyright (C) 2004, 2005, 2006, 2007 David Bateman |
|
4864
|
4 |
|
7016
|
5 This file is part of Octave. |
|
4864
|
6 |
|
7016
|
7 Octave is free software; you can redistribute it and/or modify it |
|
|
8 under the terms of the GNU General Public License as published by the |
|
|
9 Free Software Foundation; either version 3 of the License, or (at your |
|
|
10 option) any later version. |
|
|
11 |
|
|
12 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
15 for more details. |
|
4864
|
16 |
|
|
17 You should have received a copy of the GNU General Public License |
|
7016
|
18 along with Octave; see the file COPYING. If not, see |
|
|
19 <http://www.gnu.org/licenses/>. |
|
4864
|
20 |
|
|
21 */ |
|
|
22 |
|
|
23 #ifdef HAVE_CONFIG_H |
|
|
24 #include <config.h> |
|
|
25 #endif |
|
|
26 |
|
|
27 #include "dNDArray.h" |
|
|
28 #include "CNDArray.h" |
|
|
29 #include "lo-mappers.h" |
|
|
30 |
|
|
31 #include "defun-dld.h" |
|
|
32 #include "error.h" |
|
|
33 #include "oct-obj.h" |
|
|
34 |
|
5775
|
35 // FIXME -- should probably handle Inf, NaN. |
|
4864
|
36 |
|
|
37 static inline bool |
|
|
38 is_integer_value (double x) |
|
|
39 { |
|
5985
|
40 return x == std::floor (x); |
|
4864
|
41 } |
|
|
42 |
|
|
43 DEFUN_DLD (gcd, args, nargout, |
|
|
44 "-*- texinfo -*-\n\ |
|
6547
|
45 @deftypefn {Loadable Function} {@var{g} =} gcd (@var{a1}, @dots{})\n\ |
|
|
46 @deftypefnx {Loadable Function} {[@var{g}, @var{v1}, @dots{}] =} gcd (@var{a1}, @dots{})\n\ |
|
4864
|
47 \n\ |
|
|
48 If a single argument is given then compute the greatest common divisor of\n\ |
|
|
49 the elements of this argument. Otherwise if more than one argument is\n\ |
|
|
50 given all arguments must be the same size or scalar. In this case the\n\ |
|
|
51 greatest common divisor is calculated for element individually. All\n\ |
|
|
52 elements must be integers. For example,\n\ |
|
|
53 \n\ |
|
|
54 @example\n\ |
|
|
55 @group\n\ |
|
|
56 gcd ([15, 20])\n\ |
|
|
57 @result{} 5\n\ |
|
|
58 @end group\n\ |
|
|
59 @end example\n\ |
|
|
60 \n\ |
|
|
61 @noindent\n\ |
|
|
62 and\n\ |
|
|
63 \n\ |
|
|
64 @example\n\ |
|
|
65 @group\n\ |
|
|
66 gcd ([15, 9], [20 18])\n\ |
|
|
67 @result{} 5 9\n\ |
|
|
68 @end group\n\ |
|
|
69 @end example\n\ |
|
|
70 \n\ |
|
|
71 Optional return arguments @var{v1}, etc, contain integer vectors such\n\ |
|
|
72 that,\n\ |
|
|
73 \n\ |
|
|
74 @ifinfo\n\ |
|
|
75 @example\n\ |
|
6547
|
76 @var{g} = @var{v1} .* @var{a1} + @var{v2} .* @var{a2} + @dots{}\n\ |
|
4864
|
77 @end example\n\ |
|
|
78 @end ifinfo\n\ |
|
|
79 @iftex\n\ |
|
|
80 @tex\n\ |
|
|
81 $g = v_1 a_1 + v_2 a_2 + \\cdots$\n\ |
|
|
82 @end tex\n\ |
|
|
83 @end iftex\n\ |
|
|
84 \n\ |
|
7001
|
85 For backward compatibility with previous versions of this function, when\n\ |
|
|
86 all arguments are scalar, a single return argument @var{v1} containing\n\ |
|
6547
|
87 all of the values of @var{v1}, @dots{} is acceptable.\n\ |
|
5642
|
88 @seealso{lcm, min, max, ceil, floor}\n\ |
|
|
89 @end deftypefn") |
|
4864
|
90 { |
|
|
91 octave_value_list retval; |
|
|
92 |
|
|
93 int nargin = args.length (); |
|
|
94 |
|
|
95 if (nargin == 0) |
|
|
96 { |
|
5823
|
97 print_usage (); |
|
4864
|
98 return retval; |
|
|
99 } |
|
|
100 |
|
|
101 bool all_args_scalar = true; |
|
|
102 |
|
|
103 dim_vector dv(1); |
|
|
104 |
|
|
105 for (int i = 0; i < nargin; i++) |
|
|
106 { |
|
|
107 if (! args(i).is_scalar_type ()) |
|
|
108 { |
|
|
109 if (! args(i).is_matrix_type ()) |
|
|
110 { |
|
|
111 error ("gcd: invalid argument type"); |
|
|
112 return retval; |
|
|
113 } |
|
|
114 |
|
|
115 if (all_args_scalar) |
|
|
116 { |
|
|
117 all_args_scalar = false; |
|
|
118 dv = args(i).dims (); |
|
|
119 } |
|
|
120 else |
|
|
121 { |
|
|
122 if (dv != args(i).dims ()) |
|
|
123 { |
|
|
124 error ("gcd: all arguments must be the same size or scalar"); |
|
|
125 return retval; |
|
|
126 } |
|
|
127 } |
|
|
128 } |
|
|
129 } |
|
|
130 |
|
|
131 if (nargin == 1) |
|
|
132 { |
|
|
133 NDArray gg = args(0).array_value (); |
|
|
134 |
|
|
135 int nel = dv.numel (); |
|
|
136 |
|
|
137 NDArray v (dv); |
|
|
138 |
|
|
139 RowVector x (3); |
|
|
140 RowVector y (3); |
|
|
141 |
|
|
142 double g = std::abs (gg(0)); |
|
|
143 |
|
|
144 if (! is_integer_value (g)) |
|
|
145 { |
|
|
146 error ("gcd: all arguments must be integer"); |
|
|
147 return retval; |
|
|
148 } |
|
|
149 |
|
|
150 v(0) = signum (gg(0)); |
|
|
151 |
|
|
152 for (int k = 1; k < nel; k++) |
|
|
153 { |
|
|
154 x(0) = g; |
|
|
155 x(1) = 1; |
|
|
156 x(2) = 0; |
|
|
157 |
|
|
158 y(0) = std::abs (gg(k)); |
|
|
159 y(1) = 0; |
|
|
160 y(2) = 1; |
|
|
161 |
|
|
162 if (! is_integer_value (y(0))) |
|
|
163 { |
|
|
164 error ("gcd: all arguments must be integer"); |
|
|
165 return retval; |
|
|
166 } |
|
|
167 |
|
|
168 while (y(0) > 0) |
|
|
169 { |
|
5985
|
170 RowVector r = x - y * std::floor (x(0) / y(0)); |
|
4864
|
171 x = y; |
|
|
172 y = r; |
|
|
173 } |
|
|
174 |
|
|
175 g = x(0); |
|
|
176 |
|
|
177 for (int i = 0; i < k; i++) |
|
|
178 v(i) *= x(1); |
|
|
179 |
|
|
180 v(k) = x(2) * signum (gg(k)); |
|
|
181 } |
|
|
182 |
|
|
183 retval (1) = v; |
|
|
184 retval (0) = g; |
|
|
185 } |
|
|
186 else if (all_args_scalar && nargout < 3) |
|
|
187 { |
|
5985
|
188 double g = args(0).double_value (); |
|
4864
|
189 |
|
5985
|
190 if (error_state || ! is_integer_value (g)) |
|
4864
|
191 { |
|
|
192 error ("gcd: all arguments must be integer"); |
|
|
193 return retval; |
|
|
194 } |
|
|
195 |
|
|
196 RowVector v (nargin, 0); |
|
|
197 RowVector x (3); |
|
|
198 RowVector y (3); |
|
|
199 |
|
|
200 v(0) = signum (g); |
|
|
201 |
|
|
202 g = std::abs(g); |
|
|
203 |
|
|
204 for (int k = 1; k < nargin; k++) |
|
|
205 { |
|
|
206 x(0) = g; |
|
|
207 x(1) = 1; |
|
|
208 x(2) = 0; |
|
|
209 |
|
5985
|
210 y(0) = args(k).double_value (); |
|
4864
|
211 y(1) = 0; |
|
|
212 y(2) = 1; |
|
|
213 |
|
|
214 double sgn = signum (y(0)); |
|
|
215 |
|
|
216 y(0) = std::abs (y(0)); |
|
|
217 |
|
5985
|
218 if (error_state || ! is_integer_value (g)) |
|
4864
|
219 { |
|
|
220 error ("gcd: all arguments must be integer"); |
|
|
221 return retval; |
|
|
222 } |
|
|
223 |
|
|
224 while (y(0) > 0) |
|
|
225 { |
|
5985
|
226 RowVector r = x - y * std::floor (x(0) / y(0)); |
|
4864
|
227 x = y; |
|
|
228 y = r; |
|
|
229 } |
|
|
230 |
|
|
231 g = x(0); |
|
|
232 |
|
|
233 for (int i = 0; i < k; i++) |
|
|
234 v(i) *= x(1); |
|
|
235 |
|
|
236 v(k) = x(2) * sgn; |
|
|
237 } |
|
|
238 |
|
|
239 retval (1) = v; |
|
|
240 retval (0) = g; |
|
|
241 } |
|
|
242 else |
|
|
243 { |
|
5964
|
244 // FIXME -- we need to handle a possible mixture of scalar and |
|
|
245 // array values here. |
|
|
246 |
|
4864
|
247 NDArray g = args(0).array_value (); |
|
5275
|
248 |
|
|
249 OCTAVE_LOCAL_BUFFER (NDArray, v, nargin); |
|
4864
|
250 |
|
|
251 int nel = dv.numel (); |
|
|
252 |
|
|
253 v[0].resize(dv); |
|
|
254 |
|
|
255 for (int i = 0; i < nel; i++) |
|
|
256 { |
|
|
257 v[0](i) = signum (g(i)); |
|
|
258 g(i) = std::abs (g(i)); |
|
|
259 |
|
|
260 if (! is_integer_value (g(i))) |
|
|
261 { |
|
|
262 error ("gcd: all arguments must be integer"); |
|
|
263 return retval; |
|
|
264 } |
|
|
265 } |
|
|
266 |
|
|
267 RowVector x (3); |
|
|
268 RowVector y (3); |
|
|
269 |
|
|
270 for (int k = 1; k < nargin; k++) |
|
|
271 { |
|
|
272 NDArray gnew = args(k).array_value (); |
|
|
273 |
|
|
274 v[k].resize(dv); |
|
|
275 |
|
|
276 for (int n = 0; n < nel; n++) |
|
|
277 { |
|
|
278 x(0) = g(n); |
|
|
279 x(1) = 1; |
|
|
280 x(2) = 0; |
|
|
281 |
|
|
282 y(0) = std::abs (gnew(n)); |
|
|
283 y(1) = 0; |
|
|
284 y(2) = 1; |
|
|
285 |
|
|
286 if (! is_integer_value (y(0))) |
|
|
287 { |
|
|
288 error ("gcd: all arguments must be integer"); |
|
|
289 return retval; |
|
|
290 } |
|
|
291 |
|
|
292 while (y(0) > 0) |
|
|
293 { |
|
5985
|
294 RowVector r = x - y * std::floor (x(0) / y(0)); |
|
4864
|
295 x = y; |
|
|
296 y = r; |
|
|
297 } |
|
|
298 |
|
|
299 g(n) = x(0); |
|
|
300 |
|
|
301 for (int i = 0; i < k; i++) |
|
|
302 v[i](n) *= x(1); |
|
|
303 |
|
|
304 v[k](n) = x(2) * signum (gnew(n)); |
|
|
305 } |
|
|
306 } |
|
|
307 |
|
|
308 for (int k = 0; k < nargin; k++) |
|
|
309 retval(1+k) = v[k]; |
|
|
310 |
|
|
311 retval (0) = g; |
|
|
312 } |
|
|
313 |
|
|
314 return retval; |
|
|
315 } |
|
|
316 |
|
|
317 /* |
|
|
318 ;;; Local Variables: *** |
|
|
319 ;;; mode: C++ *** |
|
|
320 ;;; End: *** |
|
|
321 */ |
