|
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 |
|
|
17 Software Foundation, 59 Temple Place - Suite 330, Boston, MA |
|
|
18 02111-1307, USA. |
|
|
19 |
|
|
20 In addition to the terms of the GPL, you are permitted to link |
|
|
21 this program with any Open Source program, as defined by the |
|
|
22 Open Source Initiative (www.opensource.org) |
|
|
23 |
|
|
24 */ |
|
|
25 |
|
|
26 #ifdef HAVE_CONFIG_H |
|
|
27 #include <config.h> |
|
|
28 #endif |
|
|
29 |
|
|
30 #include "dNDArray.h" |
|
|
31 #include "CNDArray.h" |
|
|
32 #include "lo-mappers.h" |
|
|
33 |
|
|
34 #include "defun-dld.h" |
|
|
35 #include "error.h" |
|
|
36 #include "oct-obj.h" |
|
|
37 |
|
|
38 // XXX FIXME XXX -- should probably handle Inf, NaN. |
|
|
39 |
|
|
40 static inline bool |
|
|
41 is_integer_value (double x) |
|
|
42 { |
|
|
43 return x == static_cast<double> (static_cast<long> (x)); |
|
|
44 } |
|
|
45 |
|
|
46 DEFUN_DLD (gcd, args, nargout, |
|
|
47 "-*- texinfo -*-\n\ |
|
|
48 @deftypefn {Loadable Function} {@var{g} =} gcd (@var{a1}, @code{...})\n\ |
|
|
49 @deftypefnx {Loadable Function} {[@var{g}, @var{v1}, @var{...}] =} gcd (@var{a1}, @code{...})\n\ |
|
|
50 \n\ |
|
|
51 If a single argument is given then compute the greatest common divisor of\n\ |
|
|
52 the elements of this argument. Otherwise if more than one argument is\n\ |
|
|
53 given all arguments must be the same size or scalar. In this case the\n\ |
|
|
54 greatest common divisor is calculated for element individually. All\n\ |
|
|
55 elements must be integers. For example,\n\ |
|
|
56 \n\ |
|
|
57 @example\n\ |
|
|
58 @group\n\ |
|
|
59 gcd ([15, 20])\n\ |
|
|
60 @result{} 5\n\ |
|
|
61 @end group\n\ |
|
|
62 @end example\n\ |
|
|
63 \n\ |
|
|
64 @noindent\n\ |
|
|
65 and\n\ |
|
|
66 \n\ |
|
|
67 @example\n\ |
|
|
68 @group\n\ |
|
|
69 gcd ([15, 9], [20 18])\n\ |
|
|
70 @result{} 5 9\n\ |
|
|
71 @end group\n\ |
|
|
72 @end example\n\ |
|
|
73 \n\ |
|
|
74 Optional return arguments @var{v1}, etc, contain integer vectors such\n\ |
|
|
75 that,\n\ |
|
|
76 \n\ |
|
|
77 @ifinfo\n\ |
|
|
78 @example\n\ |
|
|
79 @var{g} = @var{v1} .* @var{a1} + @var{v2} .* @var{a2} + @var{...}\n\ |
|
|
80 @end example\n\ |
|
|
81 @end ifinfo\n\ |
|
|
82 @iftex\n\ |
|
|
83 @tex\n\ |
|
|
84 $g = v_1 a_1 + v_2 a_2 + \\cdots$\n\ |
|
|
85 @end tex\n\ |
|
|
86 @end iftex\n\ |
|
|
87 \n\ |
|
|
88 For backward compatiability with previous versions of this function, when\n\ |
|
|
89 all arguments are scalr, a single return argument @var{v1} containing\n\ |
|
|
90 all of the values of @var{v1}, @var{...} is acceptable.\n\ |
|
|
91 @end deftypefn\n\ |
|
|
92 @seealso{lcm, min, max, ceil, and floor}") |
|
|
93 { |
|
|
94 octave_value_list retval; |
|
|
95 |
|
|
96 int nargin = args.length (); |
|
|
97 |
|
|
98 if (nargin == 0) |
|
|
99 { |
|
|
100 print_usage ("gcd"); |
|
|
101 return retval; |
|
|
102 } |
|
|
103 |
|
|
104 bool all_args_scalar = true; |
|
|
105 |
|
|
106 dim_vector dv(1); |
|
|
107 |
|
|
108 for (int i = 0; i < nargin; i++) |
|
|
109 { |
|
|
110 if (! args(i).is_scalar_type ()) |
|
|
111 { |
|
|
112 if (! args(i).is_matrix_type ()) |
|
|
113 { |
|
|
114 error ("gcd: invalid argument type"); |
|
|
115 return retval; |
|
|
116 } |
|
|
117 |
|
|
118 if (all_args_scalar) |
|
|
119 { |
|
|
120 all_args_scalar = false; |
|
|
121 dv = args(i).dims (); |
|
|
122 } |
|
|
123 else |
|
|
124 { |
|
|
125 if (dv != args(i).dims ()) |
|
|
126 { |
|
|
127 error ("gcd: all arguments must be the same size or scalar"); |
|
|
128 return retval; |
|
|
129 } |
|
|
130 } |
|
|
131 } |
|
|
132 } |
|
|
133 |
|
|
134 if (nargin == 1) |
|
|
135 { |
|
|
136 NDArray gg = args(0).array_value (); |
|
|
137 |
|
|
138 int nel = dv.numel (); |
|
|
139 |
|
|
140 NDArray v (dv); |
|
|
141 |
|
|
142 RowVector x (3); |
|
|
143 RowVector y (3); |
|
|
144 |
|
|
145 double g = std::abs (gg(0)); |
|
|
146 |
|
|
147 if (! is_integer_value (g)) |
|
|
148 { |
|
|
149 error ("gcd: all arguments must be integer"); |
|
|
150 return retval; |
|
|
151 } |
|
|
152 |
|
|
153 v(0) = signum (gg(0)); |
|
|
154 |
|
|
155 for (int k = 1; k < nel; k++) |
|
|
156 { |
|
|
157 x(0) = g; |
|
|
158 x(1) = 1; |
|
|
159 x(2) = 0; |
|
|
160 |
|
|
161 y(0) = std::abs (gg(k)); |
|
|
162 y(1) = 0; |
|
|
163 y(2) = 1; |
|
|
164 |
|
|
165 if (! is_integer_value (y(0))) |
|
|
166 { |
|
|
167 error ("gcd: all arguments must be integer"); |
|
|
168 return retval; |
|
|
169 } |
|
|
170 |
|
|
171 while (y(0) > 0) |
|
|
172 { |
|
|
173 RowVector r = x - y * (static_cast<int> ( x(0) / y(0))); |
|
|
174 x = y; |
|
|
175 y = r; |
|
|
176 } |
|
|
177 |
|
|
178 g = x(0); |
|
|
179 |
|
|
180 for (int i = 0; i < k; i++) |
|
|
181 v(i) *= x(1); |
|
|
182 |
|
|
183 v(k) = x(2) * signum (gg(k)); |
|
|
184 } |
|
|
185 |
|
|
186 retval (1) = v; |
|
|
187 retval (0) = g; |
|
|
188 } |
|
|
189 else if (all_args_scalar && nargout < 3) |
|
|
190 { |
|
|
191 double g = args(0).int_value (true); |
|
|
192 |
|
|
193 if (error_state) |
|
|
194 { |
|
|
195 error ("gcd: all arguments must be integer"); |
|
|
196 return retval; |
|
|
197 } |
|
|
198 |
|
|
199 RowVector v (nargin, 0); |
|
|
200 RowVector x (3); |
|
|
201 RowVector y (3); |
|
|
202 |
|
|
203 v(0) = signum (g); |
|
|
204 |
|
|
205 g = std::abs(g); |
|
|
206 |
|
|
207 for (int k = 1; k < nargin; k++) |
|
|
208 { |
|
|
209 x(0) = g; |
|
|
210 x(1) = 1; |
|
|
211 x(2) = 0; |
|
|
212 |
|
|
213 y(0) = args(k).int_value (true); |
|
|
214 y(1) = 0; |
|
|
215 y(2) = 1; |
|
|
216 |
|
|
217 double sgn = signum (y(0)); |
|
|
218 |
|
|
219 y(0) = std::abs (y(0)); |
|
|
220 |
|
|
221 if (error_state) |
|
|
222 { |
|
|
223 error ("gcd: all arguments must be integer"); |
|
|
224 return retval; |
|
|
225 } |
|
|
226 |
|
|
227 while (y(0) > 0) |
|
|
228 { |
|
|
229 RowVector r = x - y * (static_cast<int> (x(0) / y(0))); |
|
|
230 x = y; |
|
|
231 y = r; |
|
|
232 } |
|
|
233 |
|
|
234 g = x(0); |
|
|
235 |
|
|
236 for (int i = 0; i < k; i++) |
|
|
237 v(i) *= x(1); |
|
|
238 |
|
|
239 v(k) = x(2) * sgn; |
|
|
240 } |
|
|
241 |
|
|
242 retval (1) = v; |
|
|
243 retval (0) = g; |
|
|
244 } |
|
|
245 else |
|
|
246 { |
|
|
247 NDArray g = args(0).array_value (); |
|
|
248 NDArray v[nargin]; |
|
|
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 { |
|
|
293 RowVector r = x - y * (static_cast<int> (x(0) / y(0))); |
|
|
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 */ |
