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