1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1996-2021 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25
26 #if defined (HAVE_CONFIG_H)
27 # include "config.h"
28 #endif
29
30 #include <istream>
31 #include <ostream>
32 #include <sstream>
33
34 #include "lo-ieee.h"
35 #include "lo-specfun.h"
36 #include "lo-mappers.h"
37
38 #include "mxarray.h"
39 #include "ovl.h"
40 #include "oct-hdf5.h"
41 #include "oct-stream.h"
42 #include "ops.h"
43 #include "ov-complex.h"
44 #include "ov-flt-complex.h"
45 #include "ov-base.h"
46 #include "ov-base-scalar.h"
47 #include "ov-base-scalar.cc"
48 #include "ov-cx-mat.h"
49 #include "ov-scalar.h"
50 #include "errwarn.h"
51 #include "pr-output.h"
52 #include "ops.h"
53
54 #include "ls-oct-text.h"
55 #include "ls-hdf5.h"
56
57 // Prevent implicit instantiations on some systems (Windows, others?)
58 // that can lead to duplicate definitions of static data members.
59
60 extern template class OCTINTERP_API octave_base_scalar<double>;
61 extern template class OCTINTERP_API octave_base_scalar<FloatComplex>;
62
63 template class octave_base_scalar<Complex>;
64
65 DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (octave_complex,
66 "complex scalar", "double");
67
68 namespace octave
69 {
70 // Complain if a complex value is used as a subscript.
71
72 class complex_index_exception : public index_exception
73 {
74 public:
75
complex_index_exception(const std::string & value)76 complex_index_exception (const std::string& value)
77 : index_exception (value)
78 {
79 // Virtual, but the one we want to call is defined in this class.
80 update_message ();
81 }
82
83 ~complex_index_exception (void) = default;
84
update_message(void)85 void update_message (void)
86 {
87 set_message (expression ()
88 + ": subscripts must be real (forgot to initialize i or j?)");
89 }
90
91 // ID of error to throw.
err_id(void) const92 const char * err_id (void) const
93 {
94 return "Octave:invalid-index";
95 }
96 };
97 }
98
99 static octave_base_value *
default_numeric_demotion_function(const octave_base_value & a)100 default_numeric_demotion_function (const octave_base_value& a)
101 {
102 const octave_complex& v = dynamic_cast<const octave_complex&> (a);
103
104 return new octave_float_complex (v.float_complex_value ());
105 }
106
107 octave_base_value::type_conv_info
numeric_demotion_function(void) const108 octave_complex::numeric_demotion_function (void) const
109 {
110 return
111 octave_base_value::type_conv_info (default_numeric_demotion_function,
112 octave_float_complex::static_type_id ());
113 }
114
115 octave_base_value *
try_narrowing_conversion(void)116 octave_complex::try_narrowing_conversion (void)
117 {
118 octave_base_value *retval = nullptr;
119
120 double im = scalar.imag ();
121
122 if (im == 0.0)
123 retval = new octave_scalar (scalar.real ());
124
125 return retval;
126 }
127
128 octave_value
do_index_op(const octave_value_list & idx,bool resize_ok)129 octave_complex::do_index_op (const octave_value_list& idx, bool resize_ok)
130 {
131 // FIXME: this doesn't solve the problem of
132 //
133 // a = i; a([1,1], [1,1], [1,1])
134 //
135 // and similar constructions. Hmm...
136
137 // FIXME: using this constructor avoids narrowing the
138 // 1x1 matrix back to a scalar value. Need a better solution
139 // to this problem.
140
141 octave_value tmp (new octave_complex_matrix (complex_matrix_value ()));
142
143 return tmp.do_index_op (idx, resize_ok);
144 }
145
146 // Can't make an index_vector from a complex number. Throw an error.
147 idx_vector
index_vector(bool) const148 octave_complex::index_vector (bool) const
149 {
150 std::ostringstream buf;
151 buf << scalar.real () << std::showpos << scalar.imag () << 'i';
152 octave::complex_index_exception e (buf.str ());
153
154 throw e;
155 }
156
157 double
double_value(bool force_conversion) const158 octave_complex::double_value (bool force_conversion) const
159 {
160 if (! force_conversion)
161 warn_implicit_conversion ("Octave:imag-to-real",
162 "complex scalar", "real scalar");
163
164 return scalar.real ();
165 }
166
167 float
float_value(bool force_conversion) const168 octave_complex::float_value (bool force_conversion) const
169 {
170 if (! force_conversion)
171 warn_implicit_conversion ("Octave:imag-to-real",
172 "complex scalar", "real scalar");
173
174 return scalar.real ();
175 }
176
177 Matrix
matrix_value(bool force_conversion) const178 octave_complex::matrix_value (bool force_conversion) const
179 {
180 Matrix retval;
181
182 if (! force_conversion)
183 warn_implicit_conversion ("Octave:imag-to-real",
184 "complex scalar", "real matrix");
185
186 retval = Matrix (1, 1, scalar.real ());
187
188 return retval;
189 }
190
191 FloatMatrix
float_matrix_value(bool force_conversion) const192 octave_complex::float_matrix_value (bool force_conversion) const
193 {
194 FloatMatrix retval;
195
196 if (! force_conversion)
197 warn_implicit_conversion ("Octave:imag-to-real",
198 "complex scalar", "real matrix");
199
200 retval = FloatMatrix (1, 1, scalar.real ());
201
202 return retval;
203 }
204
205 NDArray
array_value(bool force_conversion) const206 octave_complex::array_value (bool force_conversion) const
207 {
208 NDArray retval;
209
210 if (! force_conversion)
211 warn_implicit_conversion ("Octave:imag-to-real",
212 "complex scalar", "real matrix");
213
214 retval = NDArray (dim_vector (1, 1), scalar.real ());
215
216 return retval;
217 }
218
219 FloatNDArray
float_array_value(bool force_conversion) const220 octave_complex::float_array_value (bool force_conversion) const
221 {
222 FloatNDArray retval;
223
224 if (! force_conversion)
225 warn_implicit_conversion ("Octave:imag-to-real",
226 "complex scalar", "real matrix");
227
228 retval = FloatNDArray (dim_vector (1, 1), scalar.real ());
229
230 return retval;
231 }
232
233 Complex
complex_value(bool) const234 octave_complex::complex_value (bool) const
235 {
236 return scalar;
237 }
238
239 FloatComplex
float_complex_value(bool) const240 octave_complex::float_complex_value (bool) const
241 {
242 return static_cast<FloatComplex> (scalar);
243 }
244
245 ComplexMatrix
complex_matrix_value(bool) const246 octave_complex::complex_matrix_value (bool) const
247 {
248 return ComplexMatrix (1, 1, scalar);
249 }
250
251 FloatComplexMatrix
float_complex_matrix_value(bool) const252 octave_complex::float_complex_matrix_value (bool) const
253 {
254 return FloatComplexMatrix (1, 1, static_cast<FloatComplex> (scalar));
255 }
256
257 ComplexNDArray
complex_array_value(bool) const258 octave_complex::complex_array_value (bool /* force_conversion */) const
259 {
260 return ComplexNDArray (dim_vector (1, 1), scalar);
261 }
262
263 FloatComplexNDArray
float_complex_array_value(bool) const264 octave_complex::float_complex_array_value (bool /* force_conversion */) const
265 {
266 return FloatComplexNDArray (dim_vector (1, 1),
267 static_cast<FloatComplex> (scalar));
268 }
269
270 octave_value
resize(const dim_vector & dv,bool fill) const271 octave_complex::resize (const dim_vector& dv, bool fill) const
272 {
273 if (fill)
274 {
275 ComplexNDArray retval (dv, Complex (0));
276
277 if (dv.numel ())
278 retval(0) = scalar;
279
280 return retval;
281 }
282 else
283 {
284 ComplexNDArray retval (dv);
285
286 if (dv.numel ())
287 retval(0) = scalar;
288
289 return retval;
290 }
291 }
292
293 octave_value
as_double(void) const294 octave_complex::as_double (void) const
295 {
296 return scalar;
297 }
298
299 octave_value
as_single(void) const300 octave_complex::as_single (void) const
301 {
302 return FloatComplex (scalar);
303 }
304
305 octave_value
diag(octave_idx_type m,octave_idx_type n) const306 octave_complex::diag (octave_idx_type m, octave_idx_type n) const
307 {
308 return ComplexDiagMatrix (Array<Complex> (dim_vector (1, 1), scalar), m, n);
309 }
310
311 bool
save_ascii(std::ostream & os)312 octave_complex::save_ascii (std::ostream& os)
313 {
314 Complex c = complex_value ();
315
316 octave_write_complex (os, c);
317
318 os << "\n";
319
320 return true;
321 }
322
323 bool
load_ascii(std::istream & is)324 octave_complex::load_ascii (std::istream& is)
325 {
326 scalar = octave_read_value<Complex> (is);
327
328 if (! is)
329 error ("load: failed to load complex scalar constant");
330
331 return true;
332 }
333
334 bool
save_binary(std::ostream & os,bool)335 octave_complex::save_binary (std::ostream& os, bool /* save_as_floats */)
336 {
337 char tmp = static_cast<char> (LS_DOUBLE);
338 os.write (reinterpret_cast<char *> (&tmp), 1);
339 Complex ctmp = complex_value ();
340 os.write (reinterpret_cast<char *> (&ctmp), 16);
341
342 return true;
343 }
344
345 bool
load_binary(std::istream & is,bool swap,octave::mach_info::float_format fmt)346 octave_complex::load_binary (std::istream& is, bool swap,
347 octave::mach_info::float_format fmt)
348 {
349 char tmp;
350 if (! is.read (reinterpret_cast<char *> (&tmp), 1))
351 return false;
352
353 Complex ctmp;
354 read_doubles (is, reinterpret_cast<double *> (&ctmp),
355 static_cast<save_type> (tmp), 2, swap, fmt);
356
357 if (! is)
358 return false;
359
360 scalar = ctmp;
361 return true;
362 }
363
364 bool
save_hdf5(octave_hdf5_id loc_id,const char * name,bool)365 octave_complex::save_hdf5 (octave_hdf5_id loc_id, const char *name,
366 bool /* save_as_floats */)
367 {
368 bool retval = false;
369
370 #if defined (HAVE_HDF5)
371
372 hsize_t dimens[3];
373 hid_t space_hid, type_hid, data_hid;
374 space_hid = type_hid = data_hid = -1;
375
376 space_hid = H5Screate_simple (0, dimens, nullptr);
377 if (space_hid < 0)
378 return false;
379
380 type_hid = hdf5_make_complex_type (H5T_NATIVE_DOUBLE);
381 if (type_hid < 0)
382 {
383 H5Sclose (space_hid);
384 return false;
385 }
386 #if defined (HAVE_HDF5_18)
387 data_hid = H5Dcreate (loc_id, name, type_hid, space_hid,
388 octave_H5P_DEFAULT, octave_H5P_DEFAULT, octave_H5P_DEFAULT);
389 #else
390 data_hid = H5Dcreate (loc_id, name, type_hid, space_hid, octave_H5P_DEFAULT);
391 #endif
392 if (data_hid < 0)
393 {
394 H5Sclose (space_hid);
395 H5Tclose (type_hid);
396 return false;
397 }
398
399 Complex tmp = complex_value ();
400 retval = H5Dwrite (data_hid, type_hid, octave_H5S_ALL, octave_H5S_ALL,
401 octave_H5P_DEFAULT, &tmp) >= 0;
402
403 H5Dclose (data_hid);
404 H5Tclose (type_hid);
405 H5Sclose (space_hid);
406
407 #else
408 octave_unused_parameter (loc_id);
409 octave_unused_parameter (name);
410
411 warn_save ("hdf5");
412 #endif
413
414 return retval;
415 }
416
417 bool
load_hdf5(octave_hdf5_id loc_id,const char * name)418 octave_complex::load_hdf5 (octave_hdf5_id loc_id, const char *name)
419 {
420 bool retval = false;
421
422 #if defined (HAVE_HDF5)
423
424 #if defined (HAVE_HDF5_18)
425 hid_t data_hid = H5Dopen (loc_id, name, octave_H5P_DEFAULT);
426 #else
427 hid_t data_hid = H5Dopen (loc_id, name);
428 #endif
429 hid_t type_hid = H5Dget_type (data_hid);
430
431 hid_t complex_type = hdf5_make_complex_type (H5T_NATIVE_DOUBLE);
432
433 if (! hdf5_types_compatible (type_hid, complex_type))
434 {
435 H5Tclose (complex_type);
436 H5Dclose (data_hid);
437 return false;
438 }
439
440 hid_t space_id = H5Dget_space (data_hid);
441 hsize_t rank = H5Sget_simple_extent_ndims (space_id);
442
443 if (rank != 0)
444 {
445 H5Tclose (complex_type);
446 H5Sclose (space_id);
447 H5Dclose (data_hid);
448 return false;
449 }
450
451 // complex scalar:
452 Complex ctmp;
453 if (H5Dread (data_hid, complex_type, octave_H5S_ALL, octave_H5S_ALL,
454 octave_H5P_DEFAULT, &ctmp) >= 0)
455 {
456 retval = true;
457 scalar = ctmp;
458 }
459
460 H5Tclose (complex_type);
461 H5Sclose (space_id);
462 H5Dclose (data_hid);
463
464 #else
465 octave_unused_parameter (loc_id);
466 octave_unused_parameter (name);
467
468 warn_load ("hdf5");
469 #endif
470
471 return retval;
472 }
473
474 mxArray *
as_mxArray(void) const475 octave_complex::as_mxArray (void) const
476 {
477 mxArray *retval = new mxArray (mxDOUBLE_CLASS, 1, 1, mxCOMPLEX);
478
479 double *pr = static_cast<double *> (retval->get_data ());
480 double *pi = static_cast<double *> (retval->get_imag_data ());
481
482 pr[0] = scalar.real ();
483 pi[0] = scalar.imag ();
484
485 return retval;
486 }
487
488 octave_value
map(unary_mapper_t umap) const489 octave_complex::map (unary_mapper_t umap) const
490 {
491 switch (umap)
492 {
493 #define SCALAR_MAPPER(UMAP, FCN) \
494 case umap_ ## UMAP: \
495 return octave_value (FCN (scalar))
496
497 SCALAR_MAPPER (abs, std::abs);
498 SCALAR_MAPPER (acos, octave::math::acos);
499 SCALAR_MAPPER (acosh, octave::math::acosh);
500 SCALAR_MAPPER (angle, std::arg);
501 SCALAR_MAPPER (arg, std::arg);
502 SCALAR_MAPPER (asin, octave::math::asin);
503 SCALAR_MAPPER (asinh, octave::math::asinh);
504 SCALAR_MAPPER (atan, octave::math::atan);
505 SCALAR_MAPPER (atanh, octave::math::atanh);
506 SCALAR_MAPPER (erf, octave::math::erf);
507 SCALAR_MAPPER (erfc, octave::math::erfc);
508 SCALAR_MAPPER (erfcx, octave::math::erfcx);
509 SCALAR_MAPPER (erfi, octave::math::erfi);
510 SCALAR_MAPPER (dawson, octave::math::dawson);
511 SCALAR_MAPPER (ceil, octave::math::ceil);
512 SCALAR_MAPPER (conj, std::conj);
513 SCALAR_MAPPER (cos, std::cos);
514 SCALAR_MAPPER (cosh, std::cosh);
515 SCALAR_MAPPER (exp, std::exp);
516 SCALAR_MAPPER (expm1, octave::math::expm1);
517 SCALAR_MAPPER (fix, octave::math::fix);
518 SCALAR_MAPPER (floor, octave::math::floor);
519 SCALAR_MAPPER (imag, std::imag);
520 SCALAR_MAPPER (log, std::log);
521 SCALAR_MAPPER (log2, octave::math::log2);
522 SCALAR_MAPPER (log10, std::log10);
523 SCALAR_MAPPER (log1p, octave::math::log1p);
524 SCALAR_MAPPER (real, std::real);
525 SCALAR_MAPPER (round, octave::math::round);
526 SCALAR_MAPPER (roundb, octave::math::roundb);
527 SCALAR_MAPPER (signum, octave::math::signum);
528 SCALAR_MAPPER (sin, std::sin);
529 SCALAR_MAPPER (sinh, std::sinh);
530 SCALAR_MAPPER (sqrt, std::sqrt);
531 SCALAR_MAPPER (tan, std::tan);
532 SCALAR_MAPPER (tanh, std::tanh);
533 SCALAR_MAPPER (isfinite, octave::math::isfinite);
534 SCALAR_MAPPER (isinf, octave::math::isinf);
535 SCALAR_MAPPER (isna, octave::math::isna);
536 SCALAR_MAPPER (isnan, octave::math::isnan);
537
538 // Special cases for Matlab compatibility.
539 case umap_xtolower:
540 case umap_xtoupper:
541 return scalar;
542
543 default:
544 return octave_base_value::map (umap);
545 }
546 }
547