|
1
|
1 /* |
|
|
2 |
|
2847
|
3 Copyright (C) 1996, 1997 John W. Eaton |
|
1
|
4 |
|
|
5 This file is part of Octave. |
|
|
6 |
|
|
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 2, or (at your option) any |
|
|
10 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. |
|
|
16 |
|
|
17 You should have received a copy of the GNU General Public License |
|
|
18 along with Octave; see the file COPYING. If not, write to the Free |
|
1315
|
19 Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
|
1
|
20 |
|
|
21 */ |
|
|
22 |
|
240
|
23 #ifdef HAVE_CONFIG_H |
|
1192
|
24 #include <config.h> |
|
1
|
25 #endif |
|
|
26 |
|
1343
|
27 #include <cassert> |
|
1580
|
28 #include <climits> |
|
1343
|
29 |
|
4669
|
30 #include "Array-util.h" |
|
1352
|
31 #include "CColVector.h" |
|
453
|
32 #include "CDiagMatrix.h" |
|
1352
|
33 #include "CMatrix.h" |
|
453
|
34 #include "EIG.h" |
|
1352
|
35 #include "dDiagMatrix.h" |
|
|
36 #include "dMatrix.h" |
|
3585
|
37 #include "mx-cm-cdm.h" |
|
1651
|
38 #include "oct-cmplx.h" |
|
4153
|
39 #include "quit.h" |
|
1352
|
40 |
|
|
41 #include "error.h" |
|
4055
|
42 #include "oct-obj.h" |
|
1567
|
43 #include "utils.h" |
|
1352
|
44 #include "xpow.h" |
|
1
|
45 |
|
1567
|
46 static inline int |
|
|
47 xisint (double x) |
|
|
48 { |
|
|
49 return (D_NINT (x) == x |
|
|
50 && ((x >= 0 && x < INT_MAX) |
|
|
51 || (x <= 0 && x > INT_MIN))); |
|
|
52 } |
|
|
53 |
|
767
|
54 // Safer pow functions. |
|
|
55 // |
|
|
56 // op2 \ op1: s m cs cm |
|
|
57 // +-- +---+---+----+----+ |
|
|
58 // scalar | | 1 | 5 | 7 | 11 | |
|
|
59 // +---+---+----+----+ |
|
2365
|
60 // matrix | 2 | * | 8 | * | |
|
767
|
61 // +---+---+----+----+ |
|
|
62 // complex_scalar | 3 | 6 | 9 | 12 | |
|
|
63 // +---+---+----+----+ |
|
2365
|
64 // complex_matrix | 4 | * | 10 | * | |
|
767
|
65 // +---+---+----+----+ |
|
1
|
66 |
|
767
|
67 // -*- 1 -*- |
|
2086
|
68 octave_value |
|
1
|
69 xpow (double a, double b) |
|
|
70 { |
|
2800
|
71 if (a < 0.0 && static_cast<int> (b) != b) |
|
1
|
72 { |
|
4682
|
73 // XXX FIXME XXX -- avoid apparent GNU libm bug by converting |
|
|
74 // A and B to complex instead of just A. |
|
|
75 |
|
1
|
76 Complex atmp (a); |
|
4682
|
77 Complex btmp (b); |
|
|
78 |
|
|
79 return pow (atmp, btmp); |
|
1
|
80 } |
|
|
81 else |
|
1567
|
82 return pow (a, b); |
|
1
|
83 } |
|
|
84 |
|
767
|
85 // -*- 2 -*- |
|
2086
|
86 octave_value |
|
164
|
87 xpow (double a, const Matrix& b) |
|
1
|
88 { |
|
2086
|
89 octave_value retval; |
|
1
|
90 |
|
|
91 int nr = b.rows (); |
|
2365
|
92 int nc = b.cols (); |
|
1
|
93 |
|
|
94 if (nr == 0 || nc == 0 || nr != nc) |
|
|
95 error ("for x^A, A must be square"); |
|
|
96 else |
|
|
97 { |
|
|
98 EIG b_eig (b); |
|
|
99 ComplexColumnVector lambda (b_eig.eigenvalues ()); |
|
|
100 ComplexMatrix Q (b_eig.eigenvectors ()); |
|
|
101 |
|
|
102 for (int i = 0; i < nr; i++) |
|
|
103 { |
|
2305
|
104 Complex elt = lambda (i); |
|
1
|
105 if (imag (elt) == 0.0) |
|
2305
|
106 lambda (i) = pow (a, real (elt)); |
|
1
|
107 else |
|
2305
|
108 lambda (i) = pow (a, elt); |
|
1
|
109 } |
|
|
110 ComplexDiagMatrix D (lambda); |
|
|
111 |
|
1567
|
112 retval = ComplexMatrix (Q * D * Q.inverse ()); |
|
1
|
113 } |
|
|
114 |
|
|
115 return retval; |
|
|
116 } |
|
|
117 |
|
767
|
118 // -*- 3 -*- |
|
2086
|
119 octave_value |
|
164
|
120 xpow (double a, const Complex& b) |
|
1
|
121 { |
|
|
122 Complex result; |
|
|
123 Complex atmp (a); |
|
|
124 result = pow (atmp, b); |
|
1567
|
125 return result; |
|
1
|
126 } |
|
|
127 |
|
767
|
128 // -*- 4 -*- |
|
2086
|
129 octave_value |
|
164
|
130 xpow (double a, const ComplexMatrix& b) |
|
1
|
131 { |
|
2086
|
132 octave_value retval; |
|
1
|
133 |
|
|
134 int nr = b.rows (); |
|
2365
|
135 int nc = b.cols (); |
|
1
|
136 |
|
|
137 if (nr == 0 || nc == 0 || nr != nc) |
|
2365
|
138 error ("for x^A, A must be square"); |
|
1
|
139 else |
|
|
140 { |
|
|
141 EIG b_eig (b); |
|
|
142 ComplexColumnVector lambda (b_eig.eigenvalues ()); |
|
|
143 ComplexMatrix Q (b_eig.eigenvectors ()); |
|
|
144 |
|
|
145 for (int i = 0; i < nr; i++) |
|
|
146 { |
|
2305
|
147 Complex elt = lambda (i); |
|
1
|
148 if (imag (elt) == 0.0) |
|
2305
|
149 lambda (i) = pow (a, real (elt)); |
|
1
|
150 else |
|
2305
|
151 lambda (i) = pow (a, elt); |
|
1
|
152 } |
|
|
153 ComplexDiagMatrix D (lambda); |
|
|
154 |
|
1567
|
155 retval = ComplexMatrix (Q * D * Q.inverse ()); |
|
1
|
156 } |
|
|
157 |
|
|
158 return retval; |
|
|
159 } |
|
|
160 |
|
767
|
161 // -*- 5 -*- |
|
2086
|
162 octave_value |
|
164
|
163 xpow (const Matrix& a, double b) |
|
1
|
164 { |
|
2086
|
165 octave_value retval; |
|
1
|
166 |
|
|
167 int nr = a.rows (); |
|
2365
|
168 int nc = a.cols (); |
|
1
|
169 |
|
|
170 if (nr == 0 || nc == 0 || nr != nc) |
|
2365
|
171 error ("for A^b, A must be square"); |
|
1567
|
172 else |
|
1
|
173 { |
|
2800
|
174 if (static_cast<int> (b) == b) |
|
1
|
175 { |
|
2804
|
176 int btmp = static_cast<int> (b); |
|
1567
|
177 if (btmp == 0) |
|
|
178 { |
|
|
179 retval = DiagMatrix (nr, nr, 1.0); |
|
|
180 } |
|
|
181 else |
|
|
182 { |
|
|
183 // Too much copying? |
|
|
184 // XXX FIXME XXX -- we shouldn't do this if the exponent is |
|
|
185 // large... |
|
|
186 |
|
|
187 Matrix atmp; |
|
|
188 if (btmp < 0) |
|
|
189 { |
|
|
190 btmp = -btmp; |
|
1655
|
191 |
|
|
192 int info; |
|
|
193 double rcond = 0.0; |
|
|
194 |
|
1656
|
195 atmp = a.inverse (info, rcond, 1); |
|
1655
|
196 |
|
|
197 if (info == -1) |
|
|
198 warning ("inverse: matrix singular to machine\ |
|
|
199 precision, rcond = %g", rcond); |
|
1567
|
200 } |
|
|
201 else |
|
|
202 atmp = a; |
|
|
203 |
|
|
204 Matrix result (atmp); |
|
3178
|
205 |
|
|
206 btmp--; |
|
|
207 |
|
|
208 while (btmp > 0) |
|
|
209 { |
|
|
210 if (btmp & 1) |
|
|
211 result = result * atmp; |
|
|
212 |
|
|
213 btmp >>= 1; |
|
|
214 |
|
|
215 if (btmp > 0) |
|
|
216 atmp = atmp * atmp; |
|
|
217 } |
|
1567
|
218 |
|
|
219 retval = result; |
|
|
220 } |
|
1
|
221 } |
|
|
222 else |
|
|
223 { |
|
1567
|
224 EIG a_eig (a); |
|
|
225 ComplexColumnVector lambda (a_eig.eigenvalues ()); |
|
|
226 ComplexMatrix Q (a_eig.eigenvectors ()); |
|
|
227 |
|
|
228 for (int i = 0; i < nr; i++) |
|
2305
|
229 lambda (i) = pow (lambda (i), b); |
|
1567
|
230 |
|
|
231 ComplexDiagMatrix D (lambda); |
|
|
232 |
|
|
233 retval = ComplexMatrix (Q * D * Q.inverse ()); |
|
|
234 } |
|
|
235 } |
|
1358
|
236 |
|
1567
|
237 return retval; |
|
|
238 } |
|
1
|
239 |
|
1567
|
240 // -*- 6 -*- |
|
2086
|
241 octave_value |
|
1567
|
242 xpow (const Matrix& a, const Complex& b) |
|
|
243 { |
|
2086
|
244 octave_value retval; |
|
1
|
245 |
|
1567
|
246 int nr = a.rows (); |
|
2365
|
247 int nc = a.cols (); |
|
1567
|
248 |
|
|
249 if (nr == 0 || nc == 0 || nr != nc) |
|
2365
|
250 error ("for A^b, A must be square"); |
|
1
|
251 else |
|
|
252 { |
|
|
253 EIG a_eig (a); |
|
|
254 ComplexColumnVector lambda (a_eig.eigenvalues ()); |
|
|
255 ComplexMatrix Q (a_eig.eigenvectors ()); |
|
|
256 |
|
|
257 for (int i = 0; i < nr; i++) |
|
2305
|
258 lambda (i) = pow (lambda (i), b); |
|
1
|
259 |
|
|
260 ComplexDiagMatrix D (lambda); |
|
|
261 |
|
1567
|
262 retval = ComplexMatrix (Q * D * Q.inverse ()); |
|
1
|
263 } |
|
|
264 |
|
|
265 return retval; |
|
|
266 } |
|
|
267 |
|
767
|
268 // -*- 7 -*- |
|
2086
|
269 octave_value |
|
164
|
270 xpow (const Complex& a, double b) |
|
1
|
271 { |
|
|
272 Complex result; |
|
1567
|
273 |
|
|
274 if (xisint (b)) |
|
2800
|
275 result = pow (a, static_cast<int> (b)); |
|
1567
|
276 else |
|
|
277 result = pow (a, b); |
|
|
278 |
|
|
279 return result; |
|
1
|
280 } |
|
|
281 |
|
767
|
282 // -*- 8 -*- |
|
2086
|
283 octave_value |
|
164
|
284 xpow (const Complex& a, const Matrix& b) |
|
1
|
285 { |
|
2086
|
286 octave_value retval; |
|
1
|
287 |
|
|
288 int nr = b.rows (); |
|
2365
|
289 int nc = b.cols (); |
|
1
|
290 |
|
|
291 if (nr == 0 || nc == 0 || nr != nc) |
|
2365
|
292 error ("for x^A, A must be square"); |
|
1
|
293 else |
|
|
294 { |
|
|
295 EIG b_eig (b); |
|
|
296 ComplexColumnVector lambda (b_eig.eigenvalues ()); |
|
|
297 ComplexMatrix Q (b_eig.eigenvectors ()); |
|
|
298 |
|
|
299 for (int i = 0; i < nr; i++) |
|
|
300 { |
|
2305
|
301 Complex elt = lambda (i); |
|
1
|
302 if (imag (elt) == 0.0) |
|
2305
|
303 lambda (i) = pow (a, real (elt)); |
|
1
|
304 else |
|
2305
|
305 lambda (i) = pow (a, elt); |
|
1
|
306 } |
|
|
307 ComplexDiagMatrix D (lambda); |
|
|
308 |
|
1567
|
309 retval = ComplexMatrix (Q * D * Q.inverse ()); |
|
1
|
310 } |
|
|
311 |
|
|
312 return retval; |
|
|
313 } |
|
|
314 |
|
767
|
315 // -*- 9 -*- |
|
2086
|
316 octave_value |
|
164
|
317 xpow (const Complex& a, const Complex& b) |
|
1
|
318 { |
|
|
319 Complex result; |
|
|
320 result = pow (a, b); |
|
1567
|
321 return result; |
|
1
|
322 } |
|
|
323 |
|
767
|
324 // -*- 10 -*- |
|
2086
|
325 octave_value |
|
164
|
326 xpow (const Complex& a, const ComplexMatrix& b) |
|
1
|
327 { |
|
2086
|
328 octave_value retval; |
|
1
|
329 |
|
|
330 int nr = b.rows (); |
|
2365
|
331 int nc = b.cols (); |
|
1
|
332 |
|
|
333 if (nr == 0 || nc == 0 || nr != nc) |
|
2365
|
334 error ("for x^A, A must be square"); |
|
1
|
335 else |
|
|
336 { |
|
|
337 EIG b_eig (b); |
|
|
338 ComplexColumnVector lambda (b_eig.eigenvalues ()); |
|
|
339 ComplexMatrix Q (b_eig.eigenvectors ()); |
|
|
340 |
|
|
341 for (int i = 0; i < nr; i++) |
|
|
342 { |
|
2305
|
343 Complex elt = lambda (i); |
|
1
|
344 if (imag (elt) == 0.0) |
|
2305
|
345 lambda (i) = pow (a, real (elt)); |
|
1
|
346 else |
|
2305
|
347 lambda (i) = pow (a, elt); |
|
1
|
348 } |
|
|
349 ComplexDiagMatrix D (lambda); |
|
|
350 |
|
1567
|
351 retval = ComplexMatrix (Q * D * Q.inverse ()); |
|
1
|
352 } |
|
|
353 |
|
|
354 return retval; |
|
|
355 } |
|
|
356 |
|
767
|
357 // -*- 11 -*- |
|
2086
|
358 octave_value |
|
164
|
359 xpow (const ComplexMatrix& a, double b) |
|
1
|
360 { |
|
2086
|
361 octave_value retval; |
|
1
|
362 |
|
|
363 int nr = a.rows (); |
|
2365
|
364 int nc = a.cols (); |
|
1
|
365 |
|
|
366 if (nr == 0 || nc == 0 || nr != nc) |
|
2365
|
367 error ("for A^b, A must be square"); |
|
1567
|
368 else |
|
1
|
369 { |
|
2800
|
370 if (static_cast<int> (b) == b) |
|
1
|
371 { |
|
2804
|
372 int btmp = static_cast<int> (b); |
|
1567
|
373 if (btmp == 0) |
|
|
374 { |
|
|
375 retval = DiagMatrix (nr, nr, 1.0); |
|
|
376 } |
|
|
377 else |
|
|
378 { |
|
|
379 // Too much copying? |
|
|
380 // XXX FIXME XXX -- we shouldn't do this if the exponent is |
|
|
381 // large... |
|
|
382 |
|
|
383 ComplexMatrix atmp; |
|
|
384 if (btmp < 0) |
|
|
385 { |
|
|
386 btmp = -btmp; |
|
1655
|
387 |
|
|
388 int info; |
|
|
389 double rcond = 0.0; |
|
|
390 |
|
1656
|
391 atmp = a.inverse (info, rcond, 1); |
|
1655
|
392 |
|
|
393 if (info == -1) |
|
|
394 warning ("inverse: matrix singular to machine\ |
|
|
395 precision, rcond = %g", rcond); |
|
1567
|
396 } |
|
|
397 else |
|
|
398 atmp = a; |
|
|
399 |
|
|
400 ComplexMatrix result (atmp); |
|
3178
|
401 |
|
|
402 btmp--; |
|
|
403 |
|
|
404 while (btmp > 0) |
|
|
405 { |
|
|
406 if (btmp & 1) |
|
|
407 result = result * atmp; |
|
|
408 |
|
|
409 btmp >>= 1; |
|
|
410 |
|
|
411 if (btmp > 0) |
|
|
412 atmp = atmp * atmp; |
|
|
413 } |
|
1567
|
414 |
|
|
415 retval = result; |
|
|
416 } |
|
1
|
417 } |
|
|
418 else |
|
|
419 { |
|
1567
|
420 EIG a_eig (a); |
|
|
421 ComplexColumnVector lambda (a_eig.eigenvalues ()); |
|
|
422 ComplexMatrix Q (a_eig.eigenvectors ()); |
|
|
423 |
|
|
424 for (int i = 0; i < nr; i++) |
|
2305
|
425 lambda (i) = pow (lambda (i), b); |
|
1567
|
426 |
|
|
427 ComplexDiagMatrix D (lambda); |
|
|
428 |
|
|
429 retval = ComplexMatrix (Q * D * Q.inverse ()); |
|
|
430 } |
|
|
431 } |
|
1358
|
432 |
|
1567
|
433 return retval; |
|
|
434 } |
|
1
|
435 |
|
1567
|
436 // -*- 12 -*- |
|
2086
|
437 octave_value |
|
1567
|
438 xpow (const ComplexMatrix& a, const Complex& b) |
|
|
439 { |
|
2086
|
440 octave_value retval; |
|
1
|
441 |
|
1567
|
442 int nr = a.rows (); |
|
2365
|
443 int nc = a.cols (); |
|
1567
|
444 |
|
|
445 if (nr == 0 || nc == 0 || nr != nc) |
|
2365
|
446 error ("for A^b, A must be square"); |
|
1
|
447 else |
|
|
448 { |
|
|
449 EIG a_eig (a); |
|
|
450 ComplexColumnVector lambda (a_eig.eigenvalues ()); |
|
|
451 ComplexMatrix Q (a_eig.eigenvectors ()); |
|
|
452 |
|
|
453 for (int i = 0; i < nr; i++) |
|
2305
|
454 lambda (i) = pow (lambda (i), b); |
|
1
|
455 |
|
|
456 ComplexDiagMatrix D (lambda); |
|
|
457 |
|
1567
|
458 retval = ComplexMatrix (Q * D * Q.inverse ()); |
|
1
|
459 } |
|
|
460 |
|
|
461 return retval; |
|
|
462 } |
|
|
463 |
|
767
|
464 // Safer pow functions that work elementwise for matrices. |
|
|
465 // |
|
|
466 // op2 \ op1: s m cs cm |
|
|
467 // +-- +---+---+----+----+ |
|
|
468 // scalar | | * | 3 | * | 9 | |
|
|
469 // +---+---+----+----+ |
|
|
470 // matrix | 1 | 4 | 7 | 10 | |
|
|
471 // +---+---+----+----+ |
|
|
472 // complex_scalar | * | 5 | * | 11 | |
|
|
473 // +---+---+----+----+ |
|
|
474 // complex_matrix | 2 | 6 | 8 | 12 | |
|
|
475 // +---+---+----+----+ |
|
|
476 // |
|
|
477 // * -> not needed. |
|
1
|
478 |
|
3162
|
479 // XXX FIXME XXX -- these functions need to be fixed so that things |
|
|
480 // like |
|
|
481 // |
|
|
482 // a = -1; b = [ 0, 0.5, 1 ]; r = a .^ b |
|
|
483 // |
|
|
484 // and |
|
|
485 // |
|
|
486 // a = -1; b = [ 0, 0.5, 1 ]; for i = 1:3, r(i) = a .^ b(i), end |
|
|
487 // |
|
|
488 // produce identical results. Also, it would be nice if -1^0.5 |
|
|
489 // produced a pure imaginary result instead of a complex number with a |
|
|
490 // small real part. But perhaps that's really a problem with the math |
|
|
491 // library... |
|
|
492 |
|
767
|
493 // -*- 1 -*- |
|
2086
|
494 octave_value |
|
164
|
495 elem_xpow (double a, const Matrix& b) |
|
1
|
496 { |
|
2086
|
497 octave_value retval; |
|
1
|
498 |
|
|
499 int nr = b.rows (); |
|
2365
|
500 int nc = b.cols (); |
|
1
|
501 |
|
3162
|
502 double d1, d2; |
|
1358
|
503 |
|
3162
|
504 if (a < 0.0 && ! b.all_integers (d1, d2)) |
|
1
|
505 { |
|
|
506 Complex atmp (a); |
|
|
507 ComplexMatrix result (nr, nc); |
|
|
508 for (int j = 0; j < nc; j++) |
|
|
509 for (int i = 0; i < nr; i++) |
|
4153
|
510 { |
|
|
511 OCTAVE_QUIT; |
|
|
512 result (i, j) = pow (atmp, b (i, j)); |
|
|
513 } |
|
1
|
514 |
|
1567
|
515 retval = result; |
|
1
|
516 } |
|
|
517 else |
|
|
518 { |
|
|
519 Matrix result (nr, nc); |
|
|
520 for (int j = 0; j < nc; j++) |
|
|
521 for (int i = 0; i < nr; i++) |
|
4153
|
522 { |
|
|
523 OCTAVE_QUIT; |
|
|
524 result (i, j) = pow (a, b (i, j)); |
|
|
525 } |
|
1
|
526 |
|
1567
|
527 retval = result; |
|
1
|
528 } |
|
|
529 |
|
|
530 return retval; |
|
|
531 } |
|
|
532 |
|
767
|
533 // -*- 2 -*- |
|
2086
|
534 octave_value |
|
164
|
535 elem_xpow (double a, const ComplexMatrix& b) |
|
1
|
536 { |
|
|
537 int nr = b.rows (); |
|
2365
|
538 int nc = b.cols (); |
|
1
|
539 |
|
|
540 ComplexMatrix result (nr, nc); |
|
3125
|
541 Complex atmp (a); |
|
1
|
542 for (int j = 0; j < nc; j++) |
|
|
543 for (int i = 0; i < nr; i++) |
|
4153
|
544 { |
|
|
545 OCTAVE_QUIT; |
|
|
546 result (i, j) = pow (atmp, b (i, j)); |
|
|
547 } |
|
1
|
548 |
|
1567
|
549 return result; |
|
1
|
550 } |
|
|
551 |
|
767
|
552 // -*- 3 -*- |
|
2086
|
553 octave_value |
|
164
|
554 elem_xpow (const Matrix& a, double b) |
|
1
|
555 { |
|
2086
|
556 octave_value retval; |
|
1
|
557 |
|
|
558 int nr = a.rows (); |
|
2365
|
559 int nc = a.cols (); |
|
1
|
560 |
|
2800
|
561 if (static_cast<int> (b) != b && a.any_element_is_negative ()) |
|
1
|
562 { |
|
|
563 ComplexMatrix result (nr, nc); |
|
|
564 for (int j = 0; j < nc; j++) |
|
|
565 for (int i = 0; i < nr; i++) |
|
|
566 { |
|
4153
|
567 OCTAVE_QUIT; |
|
4790
|
568 |
|
|
569 // XXX FIXME XXX -- avoid apparent GNU libm bug by |
|
|
570 // converting A and B to complex instead of just A. |
|
|
571 |
|
2305
|
572 Complex atmp (a (i, j)); |
|
4790
|
573 Complex btmp (b); |
|
|
574 |
|
|
575 result (i, j) = pow (atmp, btmp); |
|
1
|
576 } |
|
|
577 |
|
1567
|
578 retval = result; |
|
1
|
579 } |
|
|
580 else |
|
|
581 { |
|
|
582 Matrix result (nr, nc); |
|
|
583 for (int j = 0; j < nc; j++) |
|
|
584 for (int i = 0; i < nr; i++) |
|
4153
|
585 { |
|
|
586 OCTAVE_QUIT; |
|
|
587 result (i, j) = pow (a (i, j), b); |
|
|
588 } |
|
1
|
589 |
|
1567
|
590 retval = result; |
|
1
|
591 } |
|
|
592 |
|
|
593 return retval; |
|
|
594 } |
|
|
595 |
|
767
|
596 // -*- 4 -*- |
|
2086
|
597 octave_value |
|
164
|
598 elem_xpow (const Matrix& a, const Matrix& b) |
|
1
|
599 { |
|
2086
|
600 octave_value retval; |
|
1567
|
601 |
|
1
|
602 int nr = a.rows (); |
|
2365
|
603 int nc = a.cols (); |
|
|
604 |
|
|
605 int b_nr = b.rows (); |
|
|
606 int b_nc = b.cols (); |
|
1
|
607 |
|
2365
|
608 if (nr != b_nr || nc != b_nc) |
|
|
609 { |
|
|
610 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
|
|
611 return octave_value (); |
|
|
612 } |
|
1
|
613 |
|
|
614 int convert_to_complex = 0; |
|
|
615 for (int j = 0; j < nc; j++) |
|
1321
|
616 for (int i = 0; i < nr; i++) |
|
1
|
617 { |
|
4153
|
618 OCTAVE_QUIT; |
|
2305
|
619 double atmp = a (i, j); |
|
|
620 double btmp = b (i, j); |
|
2800
|
621 if (atmp < 0.0 && static_cast<int> (btmp) != btmp) |
|
1
|
622 { |
|
|
623 convert_to_complex = 1; |
|
|
624 goto done; |
|
|
625 } |
|
|
626 } |
|
|
627 |
|
2365
|
628 done: |
|
1
|
629 |
|
|
630 if (convert_to_complex) |
|
|
631 { |
|
|
632 ComplexMatrix complex_result (nr, nc); |
|
|
633 |
|
1321
|
634 for (int j = 0; j < nc; j++) |
|
|
635 for (int i = 0; i < nr; i++) |
|
1
|
636 { |
|
4153
|
637 OCTAVE_QUIT; |
|
2305
|
638 Complex atmp (a (i, j)); |
|
|
639 Complex btmp (b (i, j)); |
|
|
640 complex_result (i, j) = pow (atmp, btmp); |
|
1
|
641 } |
|
1567
|
642 |
|
|
643 retval = complex_result; |
|
1
|
644 } |
|
|
645 else |
|
|
646 { |
|
|
647 Matrix result (nr, nc); |
|
|
648 |
|
1321
|
649 for (int j = 0; j < nc; j++) |
|
|
650 for (int i = 0; i < nr; i++) |
|
4153
|
651 { |
|
|
652 OCTAVE_QUIT; |
|
|
653 result (i, j) = pow (a (i, j), b (i, j)); |
|
|
654 } |
|
1
|
655 |
|
1567
|
656 retval = result; |
|
1
|
657 } |
|
1567
|
658 |
|
|
659 return retval; |
|
1
|
660 } |
|
|
661 |
|
767
|
662 // -*- 5 -*- |
|
2086
|
663 octave_value |
|
164
|
664 elem_xpow (const Matrix& a, const Complex& b) |
|
1
|
665 { |
|
|
666 int nr = a.rows (); |
|
2365
|
667 int nc = a.cols (); |
|
1
|
668 |
|
|
669 ComplexMatrix result (nr, nc); |
|
|
670 for (int j = 0; j < nc; j++) |
|
|
671 for (int i = 0; i < nr; i++) |
|
4153
|
672 { |
|
|
673 OCTAVE_QUIT; |
|
|
674 result (i, j) = pow (Complex (a (i, j)), b); |
|
|
675 } |
|
1
|
676 |
|
1567
|
677 return result; |
|
1
|
678 } |
|
|
679 |
|
767
|
680 // -*- 6 -*- |
|
2086
|
681 octave_value |
|
164
|
682 elem_xpow (const Matrix& a, const ComplexMatrix& b) |
|
1
|
683 { |
|
|
684 int nr = a.rows (); |
|
2365
|
685 int nc = a.cols (); |
|
|
686 |
|
|
687 int b_nr = b.rows (); |
|
|
688 int b_nc = b.cols (); |
|
1
|
689 |
|
2365
|
690 if (nr != b_nr || nc != b_nc) |
|
|
691 { |
|
|
692 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
|
|
693 return octave_value (); |
|
|
694 } |
|
1
|
695 |
|
|
696 ComplexMatrix result (nr, nc); |
|
|
697 for (int j = 0; j < nc; j++) |
|
|
698 for (int i = 0; i < nr; i++) |
|
4153
|
699 { |
|
|
700 OCTAVE_QUIT; |
|
|
701 result (i, j) = pow (Complex (a (i, j)), b (i, j)); |
|
|
702 } |
|
1
|
703 |
|
1567
|
704 return result; |
|
1
|
705 } |
|
|
706 |
|
767
|
707 // -*- 7 -*- |
|
2086
|
708 octave_value |
|
164
|
709 elem_xpow (const Complex& a, const Matrix& b) |
|
1
|
710 { |
|
|
711 int nr = b.rows (); |
|
2365
|
712 int nc = b.cols (); |
|
1
|
713 |
|
|
714 ComplexMatrix result (nr, nc); |
|
|
715 for (int j = 0; j < nc; j++) |
|
|
716 for (int i = 0; i < nr; i++) |
|
1567
|
717 { |
|
4153
|
718 OCTAVE_QUIT; |
|
2305
|
719 double btmp = b (i, j); |
|
1567
|
720 if (xisint (btmp)) |
|
2800
|
721 result (i, j) = pow (a, static_cast<int> (btmp)); |
|
1567
|
722 else |
|
2305
|
723 result (i, j) = pow (a, btmp); |
|
1567
|
724 } |
|
1
|
725 |
|
1567
|
726 return result; |
|
1
|
727 } |
|
|
728 |
|
767
|
729 // -*- 8 -*- |
|
2086
|
730 octave_value |
|
164
|
731 elem_xpow (const Complex& a, const ComplexMatrix& b) |
|
1
|
732 { |
|
|
733 int nr = b.rows (); |
|
2365
|
734 int nc = b.cols (); |
|
1
|
735 |
|
|
736 ComplexMatrix result (nr, nc); |
|
|
737 for (int j = 0; j < nc; j++) |
|
|
738 for (int i = 0; i < nr; i++) |
|
4153
|
739 { |
|
|
740 OCTAVE_QUIT; |
|
|
741 result (i, j) = pow (a, b (i, j)); |
|
|
742 } |
|
1
|
743 |
|
1567
|
744 return result; |
|
1
|
745 } |
|
|
746 |
|
767
|
747 // -*- 9 -*- |
|
2086
|
748 octave_value |
|
164
|
749 elem_xpow (const ComplexMatrix& a, double b) |
|
1
|
750 { |
|
|
751 int nr = a.rows (); |
|
2365
|
752 int nc = a.cols (); |
|
1
|
753 |
|
|
754 ComplexMatrix result (nr, nc); |
|
|
755 |
|
1567
|
756 if (xisint (b)) |
|
|
757 { |
|
|
758 for (int j = 0; j < nc; j++) |
|
|
759 for (int i = 0; i < nr; i++) |
|
4153
|
760 { |
|
|
761 OCTAVE_QUIT; |
|
|
762 result (i, j) = pow (a (i, j), static_cast<int> (b)); |
|
|
763 } |
|
1567
|
764 } |
|
|
765 else |
|
|
766 { |
|
|
767 for (int j = 0; j < nc; j++) |
|
|
768 for (int i = 0; i < nr; i++) |
|
4153
|
769 { |
|
|
770 OCTAVE_QUIT; |
|
|
771 result (i, j) = pow (a (i, j), b); |
|
|
772 } |
|
1567
|
773 } |
|
|
774 |
|
|
775 return result; |
|
1
|
776 } |
|
|
777 |
|
767
|
778 // -*- 10 -*- |
|
2086
|
779 octave_value |
|
164
|
780 elem_xpow (const ComplexMatrix& a, const Matrix& b) |
|
1
|
781 { |
|
|
782 int nr = a.rows (); |
|
2365
|
783 int nc = a.cols (); |
|
|
784 |
|
|
785 int b_nr = b.rows (); |
|
|
786 int b_nc = b.cols (); |
|
1
|
787 |
|
2365
|
788 if (nr != b_nr || nc != b_nc) |
|
|
789 { |
|
|
790 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
|
|
791 return octave_value (); |
|
|
792 } |
|
1
|
793 |
|
|
794 ComplexMatrix result (nr, nc); |
|
|
795 for (int j = 0; j < nc; j++) |
|
|
796 for (int i = 0; i < nr; i++) |
|
1567
|
797 { |
|
4153
|
798 OCTAVE_QUIT; |
|
2305
|
799 double btmp = b (i, j); |
|
1567
|
800 if (xisint (btmp)) |
|
2800
|
801 result (i, j) = pow (a (i, j), static_cast<int> (btmp)); |
|
1567
|
802 else |
|
2305
|
803 result (i, j) = pow (a (i, j), btmp); |
|
1567
|
804 } |
|
1
|
805 |
|
1567
|
806 return result; |
|
1
|
807 } |
|
|
808 |
|
767
|
809 // -*- 11 -*- |
|
2086
|
810 octave_value |
|
164
|
811 elem_xpow (const ComplexMatrix& a, const Complex& b) |
|
1
|
812 { |
|
|
813 int nr = a.rows (); |
|
2365
|
814 int nc = a.cols (); |
|
1
|
815 |
|
|
816 ComplexMatrix result (nr, nc); |
|
|
817 for (int j = 0; j < nc; j++) |
|
|
818 for (int i = 0; i < nr; i++) |
|
4153
|
819 { |
|
|
820 OCTAVE_QUIT; |
|
|
821 result (i, j) = pow (a (i, j), b); |
|
|
822 } |
|
1
|
823 |
|
1567
|
824 return result; |
|
1
|
825 } |
|
|
826 |
|
767
|
827 // -*- 12 -*- |
|
2086
|
828 octave_value |
|
164
|
829 elem_xpow (const ComplexMatrix& a, const ComplexMatrix& b) |
|
1
|
830 { |
|
|
831 int nr = a.rows (); |
|
2365
|
832 int nc = a.cols (); |
|
|
833 |
|
|
834 int b_nr = b.rows (); |
|
|
835 int b_nc = b.cols (); |
|
|
836 |
|
|
837 if (nr != b_nr || nc != b_nc) |
|
|
838 { |
|
|
839 gripe_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
|
|
840 return octave_value (); |
|
|
841 } |
|
1
|
842 |
|
|
843 ComplexMatrix result (nr, nc); |
|
|
844 for (int j = 0; j < nc; j++) |
|
|
845 for (int i = 0; i < nr; i++) |
|
4153
|
846 { |
|
|
847 OCTAVE_QUIT; |
|
|
848 result (i, j) = pow (a (i, j), b (i, j)); |
|
|
849 } |
|
1
|
850 |
|
1567
|
851 return result; |
|
1
|
852 } |
|
|
853 |
|
4543
|
854 // Safer pow functions that work elementwise for N-d arrays. |
|
|
855 // |
|
|
856 // op2 \ op1: s nd cs cnd |
|
|
857 // +-- +---+---+----+----+ |
|
|
858 // scalar | | * | 3 | * | 9 | |
|
|
859 // +---+---+----+----+ |
|
|
860 // N_d | 1 | 4 | 7 | 10 | |
|
|
861 // +---+---+----+----+ |
|
|
862 // complex_scalar | * | 5 | * | 11 | |
|
|
863 // +---+---+----+----+ |
|
|
864 // complex_N_d | 2 | 6 | 8 | 12 | |
|
|
865 // +---+---+----+----+ |
|
|
866 // |
|
|
867 // * -> not needed. |
|
|
868 |
|
|
869 // XXX FIXME XXX -- these functions need to be fixed so that things |
|
|
870 // like |
|
|
871 // |
|
|
872 // a = -1; b = [ 0, 0.5, 1 ]; r = a .^ b |
|
|
873 // |
|
|
874 // and |
|
|
875 // |
|
|
876 // a = -1; b = [ 0, 0.5, 1 ]; for i = 1:3, r(i) = a .^ b(i), end |
|
|
877 // |
|
|
878 // produce identical results. Also, it would be nice if -1^0.5 |
|
|
879 // produced a pure imaginary result instead of a complex number with a |
|
|
880 // small real part. But perhaps that's really a problem with the math |
|
|
881 // library... |
|
|
882 |
|
|
883 // -*- 1 -*- |
|
|
884 octave_value |
|
|
885 elem_xpow (double a, const NDArray& b) |
|
|
886 { |
|
|
887 octave_value retval; |
|
|
888 |
|
|
889 double d1, d2; |
|
|
890 |
|
|
891 if (a < 0.0 && ! b.all_integers (d1, d2)) |
|
|
892 { |
|
|
893 Complex atmp (a); |
|
|
894 ComplexNDArray result (b.dims ()); |
|
|
895 for (int i = 0; i < b.length (); i++) |
|
|
896 { |
|
|
897 OCTAVE_QUIT; |
|
|
898 result(i) = pow (atmp, b(i)); |
|
|
899 } |
|
|
900 |
|
|
901 retval = result; |
|
|
902 } |
|
|
903 else |
|
|
904 { |
|
|
905 NDArray result (b.dims ()); |
|
|
906 for (int i = 0; i < b.length (); i++) |
|
|
907 { |
|
|
908 OCTAVE_QUIT; |
|
|
909 result (i) = pow (a, b(i)); |
|
|
910 } |
|
|
911 |
|
|
912 retval = result; |
|
|
913 } |
|
|
914 |
|
|
915 return retval; |
|
|
916 } |
|
|
917 |
|
|
918 // -*- 2 -*- |
|
|
919 octave_value |
|
|
920 elem_xpow (double a, const ComplexNDArray& b) |
|
|
921 { |
|
|
922 ComplexNDArray result (b.dims ()); |
|
|
923 Complex atmp (a); |
|
|
924 for (int i = 0; i < b.length (); i++) |
|
|
925 { |
|
|
926 OCTAVE_QUIT; |
|
|
927 result(i) = pow (atmp, b(i)); |
|
|
928 } |
|
|
929 |
|
|
930 return result; |
|
|
931 } |
|
|
932 |
|
|
933 // -*- 3 -*- |
|
|
934 octave_value |
|
|
935 elem_xpow (const NDArray& a, double b) |
|
|
936 { |
|
|
937 octave_value retval; |
|
|
938 |
|
|
939 if (static_cast<int> (b) != b && a.any_element_is_negative ()) |
|
|
940 { |
|
|
941 ComplexNDArray result (a.dims ()); |
|
|
942 |
|
|
943 for (int i = 0; i < a.length (); i++) |
|
|
944 { |
|
|
945 OCTAVE_QUIT; |
|
4790
|
946 |
|
|
947 // XXX FIXME XXX -- avoid apparent GNU libm bug by |
|
|
948 // converting A and B to complex instead of just A. |
|
|
949 |
|
|
950 Complex atmp (a (i)); |
|
|
951 Complex btmp (b); |
|
|
952 |
|
|
953 result(i) = pow (atmp, btmp); |
|
4543
|
954 } |
|
|
955 |
|
|
956 retval = result; |
|
|
957 } |
|
|
958 else |
|
|
959 { |
|
|
960 NDArray result (a.dims ()); |
|
|
961 |
|
|
962 for (int i = 0; i < a.length (); i++) |
|
|
963 { |
|
|
964 OCTAVE_QUIT; |
|
|
965 result(i) = pow (a(i), b); |
|
|
966 } |
|
|
967 |
|
|
968 retval = result; |
|
|
969 } |
|
|
970 |
|
|
971 return retval; |
|
|
972 } |
|
|
973 |
|
|
974 // -*- 4 -*- |
|
|
975 octave_value |
|
|
976 elem_xpow (const NDArray& a, const NDArray& b) |
|
|
977 { |
|
|
978 octave_value retval; |
|
|
979 |
|
|
980 dim_vector a_dims = a.dims (); |
|
|
981 dim_vector b_dims = b.dims (); |
|
|
982 |
|
|
983 if (a_dims != b_dims) |
|
|
984 { |
|
|
985 gripe_nonconformant ("operator .^", a_dims, b_dims); |
|
|
986 return octave_value (); |
|
|
987 } |
|
|
988 |
|
|
989 int len = a.length (); |
|
|
990 |
|
|
991 bool convert_to_complex = false; |
|
|
992 |
|
|
993 for (int i = 0; i < len; i++) |
|
|
994 { |
|
|
995 OCTAVE_QUIT; |
|
|
996 double atmp = a(i); |
|
|
997 double btmp = b(i); |
|
|
998 if (atmp < 0.0 && static_cast<int> (btmp) != btmp) |
|
|
999 { |
|
|
1000 convert_to_complex = true; |
|
|
1001 goto done; |
|
|
1002 } |
|
|
1003 } |
|
|
1004 |
|
|
1005 done: |
|
|
1006 |
|
|
1007 if (convert_to_complex) |
|
|
1008 { |
|
|
1009 ComplexNDArray complex_result (a_dims); |
|
|
1010 |
|
|
1011 for (int i = 0; i < len; i++) |
|
|
1012 { |
|
|
1013 OCTAVE_QUIT; |
|
|
1014 Complex atmp (a(i)); |
|
|
1015 Complex btmp (b(i)); |
|
|
1016 complex_result(i) = pow (atmp, btmp); |
|
|
1017 } |
|
|
1018 |
|
|
1019 retval = complex_result; |
|
|
1020 } |
|
|
1021 else |
|
|
1022 { |
|
|
1023 NDArray result (a_dims); |
|
|
1024 |
|
|
1025 for (int i = 0; i < len; i++) |
|
|
1026 { |
|
|
1027 OCTAVE_QUIT; |
|
|
1028 result(i) = pow (a(i), b(i)); |
|
|
1029 } |
|
|
1030 |
|
|
1031 retval = result; |
|
|
1032 } |
|
|
1033 |
|
|
1034 return retval; |
|
|
1035 } |
|
|
1036 |
|
|
1037 // -*- 5 -*- |
|
|
1038 octave_value |
|
|
1039 elem_xpow (const NDArray& a, const Complex& b) |
|
|
1040 { |
|
|
1041 ComplexNDArray result (a.dims ()); |
|
|
1042 |
|
|
1043 for (int i = 0; i < a.length (); i++) |
|
|
1044 { |
|
|
1045 OCTAVE_QUIT; |
|
|
1046 result(i) = pow (Complex (a(i)), b); |
|
|
1047 } |
|
|
1048 |
|
|
1049 return result; |
|
|
1050 } |
|
|
1051 |
|
|
1052 // -*- 6 -*- |
|
|
1053 octave_value |
|
|
1054 elem_xpow (const NDArray& a, const ComplexNDArray& b) |
|
|
1055 { |
|
|
1056 dim_vector a_dims = a.dims (); |
|
|
1057 dim_vector b_dims = b.dims (); |
|
|
1058 |
|
|
1059 if (a_dims != b_dims) |
|
|
1060 { |
|
|
1061 gripe_nonconformant ("operator .^", a_dims, b_dims); |
|
|
1062 return octave_value (); |
|
|
1063 } |
|
|
1064 |
|
|
1065 ComplexNDArray result (a_dims); |
|
|
1066 for (int i = 0; i < a.length (); i++) |
|
|
1067 { |
|
|
1068 OCTAVE_QUIT; |
|
|
1069 result(i) = pow (Complex (a(i)), b(i)); |
|
|
1070 } |
|
|
1071 |
|
|
1072 return result; |
|
|
1073 } |
|
|
1074 |
|
|
1075 // -*- 7 -*- |
|
|
1076 octave_value |
|
|
1077 elem_xpow (const Complex& a, const NDArray& b) |
|
|
1078 { |
|
|
1079 ComplexNDArray result (b.dims ()); |
|
|
1080 for (int i = 0; i < b.length (); i++) |
|
|
1081 { |
|
|
1082 OCTAVE_QUIT; |
|
|
1083 double btmp = b(i); |
|
|
1084 if (xisint (btmp)) |
|
|
1085 result(i) = pow (a, static_cast<int> (btmp)); |
|
|
1086 else |
|
|
1087 result(i) = pow (a, btmp); |
|
|
1088 } |
|
|
1089 |
|
|
1090 return result; |
|
|
1091 } |
|
|
1092 |
|
|
1093 // -*- 8 -*- |
|
|
1094 octave_value |
|
|
1095 elem_xpow (const Complex& a, const ComplexNDArray& b) |
|
|
1096 { |
|
|
1097 ComplexNDArray result (b.dims ()); |
|
|
1098 for (int i = 0; i < b.length (); i++) |
|
|
1099 { |
|
|
1100 OCTAVE_QUIT; |
|
|
1101 result(i) = pow (a, b(i)); |
|
|
1102 } |
|
|
1103 |
|
|
1104 return result; |
|
|
1105 } |
|
|
1106 |
|
|
1107 // -*- 9 -*- |
|
|
1108 octave_value |
|
|
1109 elem_xpow (const ComplexNDArray& a, double b) |
|
|
1110 { |
|
|
1111 ComplexNDArray result (a.dims ()); |
|
|
1112 |
|
|
1113 if (xisint (b)) |
|
|
1114 { |
|
|
1115 for (int i = 0; i < a.length (); i++) |
|
|
1116 { |
|
|
1117 OCTAVE_QUIT; |
|
|
1118 result(i) = pow (a(i), static_cast<int> (b)); |
|
|
1119 } |
|
|
1120 } |
|
|
1121 else |
|
|
1122 { |
|
|
1123 for (int i = 0; i < a.length (); i++) |
|
|
1124 { |
|
|
1125 OCTAVE_QUIT; |
|
|
1126 result(i) = pow (a(i), b); |
|
|
1127 } |
|
|
1128 } |
|
|
1129 |
|
|
1130 return result; |
|
|
1131 } |
|
|
1132 |
|
|
1133 // -*- 10 -*- |
|
|
1134 octave_value |
|
|
1135 elem_xpow (const ComplexNDArray& a, const NDArray& b) |
|
|
1136 { |
|
|
1137 dim_vector a_dims = a.dims (); |
|
|
1138 dim_vector b_dims = b.dims (); |
|
|
1139 |
|
|
1140 if (a_dims != b_dims) |
|
|
1141 { |
|
|
1142 gripe_nonconformant ("operator .^", a_dims, b_dims); |
|
|
1143 return octave_value (); |
|
|
1144 } |
|
|
1145 |
|
|
1146 ComplexNDArray result (a_dims); |
|
|
1147 for (int i = 0; i < a.length (); i++) |
|
|
1148 { |
|
|
1149 OCTAVE_QUIT; |
|
|
1150 double btmp = b(i); |
|
|
1151 if (xisint (btmp)) |
|
|
1152 result(i) = pow (a(i), static_cast<int> (btmp)); |
|
|
1153 else |
|
|
1154 result(i) = pow (a(i), btmp); |
|
|
1155 } |
|
|
1156 |
|
|
1157 return result; |
|
|
1158 } |
|
|
1159 |
|
|
1160 // -*- 11 -*- |
|
|
1161 octave_value |
|
|
1162 elem_xpow (const ComplexNDArray& a, const Complex& b) |
|
|
1163 { |
|
|
1164 ComplexNDArray result (a.dims ()); |
|
|
1165 for (int i = 0; i < a.length (); i++) |
|
|
1166 { |
|
|
1167 OCTAVE_QUIT; |
|
|
1168 result(i) = pow (a(i), b); |
|
|
1169 } |
|
|
1170 |
|
|
1171 return result; |
|
|
1172 } |
|
|
1173 |
|
|
1174 // -*- 12 -*- |
|
|
1175 octave_value |
|
|
1176 elem_xpow (const ComplexNDArray& a, const ComplexNDArray& b) |
|
|
1177 { |
|
|
1178 dim_vector a_dims = a.dims (); |
|
|
1179 dim_vector b_dims = b.dims (); |
|
|
1180 |
|
|
1181 if (a_dims != b_dims) |
|
|
1182 { |
|
|
1183 gripe_nonconformant ("operator .^", a_dims, b_dims); |
|
|
1184 return octave_value (); |
|
|
1185 } |
|
|
1186 |
|
|
1187 ComplexNDArray result (a_dims); |
|
|
1188 for (int i = 0; i < a.length (); i++) |
|
|
1189 { |
|
|
1190 OCTAVE_QUIT; |
|
|
1191 result(i) = pow (a(i), b(i)); |
|
|
1192 } |
|
|
1193 |
|
|
1194 return result; |
|
|
1195 } |
|
|
1196 |
|
1
|
1197 /* |
|
|
1198 ;;; Local Variables: *** |
|
|
1199 ;;; mode: C++ *** |
|
|
1200 ;;; End: *** |
|
|
1201 */ |
