GNU Octave  4.2.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Range.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1993-2017 John W. Eaton
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 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.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #if defined (HAVE_CONFIG_H)
24 # include "config.h"
25 #endif
26 
27 #include <cfloat>
28 
29 #include <iostream>
30 #include <limits>
31 
32 #include "Range.h"
33 #include "lo-error.h"
34 #include "lo-mappers.h"
35 #include "lo-math.h"
36 #include "lo-utils.h"
37 #include "Array-util.h"
38 
39 bool
41 {
42  // If the base and increment are ints, the final value in the range
43  // will also be an integer, even if the limit is not. If there is one
44  // or fewer elements only the base needs to be an integer
45 
49 }
50 
51 Matrix
52 Range::matrix_value (void) const
53 {
54  if (rng_numel > 0 && cache.is_empty ())
55  {
56  cache.resize (1, rng_numel);
57 
58  // The first element must always be *exactly* the base.
59  // E.g, -0 would otherwise become +0 in the loop (-0 + 0*increment).
60  cache(0) = rng_base;
61 
62  double b = rng_base;
63  double increment = rng_inc;
64  for (octave_idx_type i = 1; i < rng_numel - 1; i++)
65  cache(i) = b + i * increment;
66 
67  cache(rng_numel - 1) = rng_limit;
68  }
69 
70  return cache;
71 }
72 
73 double
75 {
78 
79  if (i == 0)
80  return rng_base;
81  else if (i < rng_numel - 1)
82  return rng_base + i * rng_inc;
83  else
84  return rng_limit;
85 }
86 
87 double
89 {
90 #if defined (OCTAVE_ENABLE_BOUNDS_CHECK)
91  return checkelem (i);
92 #else
93  if (i == 0)
94  return rng_base;
95  else if (i < rng_numel - 1)
96  return rng_base + i * rng_inc;
97  else
98  return rng_limit;
99 #endif
100 }
101 
102 // Helper class used solely for idx_vector.loop () function call
104 {
105 public:
106  __rangeidx_helper (double *a, double b, double i, double l, octave_idx_type n)
107  : array (a), base (b), inc (i), limit (l), nmax (n-1) { }
108 
110  {
111  if (i == 0)
112  *array++ = base;
113  else if (i < nmax)
114  *array++ = base + i * inc;
115  else
116  *array++ = limit;
117  }
118 
119 private:
120 
121  double *array, base, inc, limit;
123 
124 };
125 
127 Range::index (const idx_vector& i) const
128 {
130 
132 
133  if (i.is_colon ())
134  {
135  retval = matrix_value ().reshape (dim_vector (rng_numel, 1));
136  }
137  else
138  {
139  if (i.extent (n) != n)
140  octave::err_index_out_of_range (1, 1, i.extent (n), n); // throws
141 
142  dim_vector rd = i.orig_dimensions ();
143  octave_idx_type il = i.length (n);
144 
145  // taken from Array.cc.
146  if (n != 1 && rd.is_vector ())
147  rd = dim_vector (1, il);
148 
149  retval.clear (rd);
150 
151  // idx_vector loop across all values in i,
152  // executing __rangeidx_helper (i) for each i
153  i.loop (n, __rangeidx_helper (retval.fortran_vec (),
155  }
156 
157  return retval;
158 }
159 
160 // NOTE: max and min only return useful values if numel > 0.
161 // do_minmax_body() in max.cc avoids calling Range::min/max if numel == 0.
162 
163 double
164 Range::min (void) const
165 {
166  double retval = 0.0;
167  if (rng_numel > 0)
168  {
169  if (rng_inc > 0)
170  retval = rng_base;
171  else
172  {
173  retval = rng_base + (rng_numel - 1) * rng_inc;
174 
175  // Require '<=' test. See note in max ().
176  if (retval <= rng_limit)
177  retval = rng_limit;
178  }
179 
180  }
181  return retval;
182 }
183 
184 double
185 Range::max (void) const
186 {
187  double retval = 0.0;
188  if (rng_numel > 0)
189  {
190  if (rng_inc > 0)
191  {
192  retval = rng_base + (rng_numel - 1) * rng_inc;
193 
194  // On some machines (x86 with extended precision floating point
195  // arithmetic, for example) it is possible that we can overshoot the
196  // limit by approximately the machine precision even though we were
197  // very careful in our calculation of the number of elements.
198  // Therefore, we clip the result to the limit if it overshoots.
199  // The test also includes equality (>= rng_limit) to have expressions
200  // such as -5:1:-0 result in a -0 endpoint.
201  if (retval >= rng_limit)
202  retval = rng_limit;
203  }
204  else
205  retval = rng_base;
206  }
207  return retval;
208 }
209 
210 void
211 Range::sort_internal (bool ascending)
212 {
213  if ((ascending && rng_base > rng_limit && rng_inc < 0.0)
214  || (! ascending && rng_base < rng_limit && rng_inc > 0.0))
215  {
217  rng_inc = -rng_inc;
218  clear_cache ();
219  }
220 }
221 
222 void
224 {
225  octave_idx_type nel = numel ();
226 
227  sidx.resize (dim_vector (1, nel));
228 
229  octave_idx_type *psidx = sidx.fortran_vec ();
230 
231  bool reverse = false;
232 
233  if ((ascending && rng_base > rng_limit && rng_inc < 0.0)
234  || (! ascending && rng_base < rng_limit && rng_inc > 0.0))
235  {
237  rng_inc = -rng_inc;
238  clear_cache ();
239  reverse = true;
240  }
241 
242  octave_idx_type tmp = reverse ? nel - 1 : 0;
243  octave_idx_type stp = reverse ? -1 : 1;
244 
245  for (octave_idx_type i = 0; i < nel; i++, tmp += stp)
246  psidx[i] = tmp;
247 
248 }
249 
250 Matrix
252 {
253  return matrix_value ().diag (k);
254 }
255 
256 Range
258 {
259  Range retval = *this;
260 
261  if (dim == 1)
262  {
263  if (mode == ASCENDING)
264  retval.sort_internal (true);
265  else if (mode == DESCENDING)
266  retval.sort_internal (false);
267  }
268  else if (dim != 0)
269  (*current_liboctave_error_handler) ("Range::sort: invalid dimension");
270 
271  return retval;
272 }
273 
274 Range
276  sortmode mode) const
277 {
278  Range retval = *this;
279 
280  if (dim == 1)
281  {
282  if (mode == ASCENDING)
283  retval.sort_internal (sidx, true);
284  else if (mode == DESCENDING)
285  retval.sort_internal (sidx, false);
286  }
287  else if (dim != 0)
288  (*current_liboctave_error_handler) ("Range::sort: invalid dimension");
289 
290  return retval;
291 }
292 
293 sortmode
295 {
296  if (rng_numel > 1 && rng_inc > 0)
297  mode = (mode == DESCENDING) ? UNSORTED : ASCENDING;
298  else if (rng_numel > 1 && rng_inc < 0)
299  mode = (mode == ASCENDING) ? UNSORTED : DESCENDING;
300  else
301  mode = mode ? mode : ASCENDING;
302 
303  return mode;
304 }
305 
306 void
308 {
309  if (rng_base != b)
310  {
311  rng_base = b;
312 
313  init ();
314  }
315 }
316 
317 void
319 {
320  if (rng_limit != l)
321  {
322  rng_limit = l;
323 
324  init ();
325  }
326 }
327 
328 void
330 {
331  if (rng_inc != i)
332  {
333  rng_inc = i;
334 
335  init ();
336  }
337 }
338 
339 std::ostream&
340 operator << (std::ostream& os, const Range& a)
341 {
342  double b = a.base ();
343  double increment = a.inc ();
344  octave_idx_type num_elem = a.numel ();
345 
346  if (num_elem > 1)
347  {
348  // First element must be the base *exactly* (-0).
349  os << b << " ";
350  for (octave_idx_type i = 1; i < num_elem-1; i++)
351  os << b + i * increment << " ";
352  }
353 
354  // Print out exactly the last element, rather than a calculated last element.
355  os << a.rng_limit << "\n";
356 
357  return os;
358 }
359 
360 std::istream&
361 operator >> (std::istream& is, Range& a)
362 {
363  is >> a.rng_base;
364  if (is)
365  {
366  double tmp_rng_limit;
367  is >> tmp_rng_limit;
368 
369  if (is)
370  is >> a.rng_inc;
371 
372  // Clip the rng_limit to the true limit and rebuild numel, clear cache
373  a.set_limit (tmp_rng_limit);
374  }
375 
376  return is;
377 }
378 
379 Range
381 {
382  return Range (-r.base (), -r.limit (), -r.inc (), r.numel ());
383 }
384 
385 Range operator + (double x, const Range& r)
386 {
387  Range result (x + r.base (), x + r.limit (), r.inc (), r.numel ());
388  if (result.rng_numel < 0)
389  result.cache = x + r.matrix_value ();
390 
391  return result;
392 }
393 
394 Range operator + (const Range& r, double x)
395 {
396  Range result (r.base () + x, r.limit () + x, r.inc (), r.numel ());
397  if (result.rng_numel < 0)
398  result.cache = r.matrix_value () + x;
399 
400  return result;
401 }
402 
403 Range operator - (double x, const Range& r)
404 {
405  Range result (x - r.base (), x - r.limit (), -r.inc (), r.numel ());
406  if (result.rng_numel < 0)
407  result.cache = x - r.matrix_value ();
408 
409  return result;
410 }
411 
412 Range operator - (const Range& r, double x)
413 {
414  Range result (r.base () - x, r.limit () - x, r.inc (), r.numel ());
415  if (result.rng_numel < 0)
416  result.cache = r.matrix_value () - x;
417 
418  return result;
419 }
420 
421 Range operator * (double x, const Range& r)
422 {
423  Range result (x * r.base (), x * r.limit (), x * r.inc (), r.numel ());
424  if (result.rng_numel < 0)
425  result.cache = x * r.matrix_value ();
426 
427  return result;
428 }
429 
430 Range operator * (const Range& r, double x)
431 {
432  Range result (r.base () * x, r.limit () * x, r.inc () * x, r.numel ());
433  if (result.rng_numel < 0)
434  result.cache = r.matrix_value () * x;
435 
436  return result;
437 }
438 
439 // C See Knuth, Art Of Computer Programming, Vol. 1, Problem 1.2.4-5.
440 // C
441 // C===Tolerant FLOOR function.
442 // C
443 // C X - is given as a Double Precision argument to be operated on.
444 // C It is assumed that X is represented with M mantissa bits.
445 // C CT - is given as a Comparison Tolerance such that
446 // C 0.LT.CT.LE.3-SQRT(5)/2. If the relative difference between
447 // C X and A whole number is less than CT, then TFLOOR is
448 // C returned as this whole number. By treating the
449 // C floating-point numbers as a finite ordered set note that
450 // C the heuristic EPS=2.**(-(M-1)) and CT=3*EPS causes
451 // C arguments of TFLOOR/TCEIL to be treated as whole numbers
452 // C if they are exactly whole numbers or are immediately
453 // C adjacent to whole number representations. Since EPS, the
454 // C "distance" between floating-point numbers on the unit
455 // C interval, and M, the number of bits in X'S mantissa, exist
456 // C on every floating-point computer, TFLOOR/TCEIL are
457 // C consistently definable on every floating-point computer.
458 // C
459 // C For more information see the following references:
460 // C (1) P. E. Hagerty, "More On Fuzzy Floor And Ceiling," APL QUOTE
461 // C QUAD 8(4):20-24, June 1978. Note that TFLOOR=FL5.
462 // C (2) L. M. Breed, "Definitions For Fuzzy Floor And Ceiling", APL
463 // C QUOTE QUAD 8(3):16-23, March 1978. This paper cites FL1 through
464 // C FL5, the history of five years of evolutionary development of
465 // C FL5 - the seven lines of code below - by open collaboration
466 // C and corroboration of the mathematical-computing community.
467 // C
468 // C Penn State University Center for Academic Computing
469 // C H. D. Knoble - August, 1978.
470 
471 static inline double
472 tfloor (double x, double ct)
473 {
474 // C---------FLOOR(X) is the largest integer algebraically less than
475 // C or equal to X; that is, the unfuzzy FLOOR function.
476 
477 // DINT (X) = X - DMOD (X, 1.0);
478 // FLOOR (X) = DINT (X) - DMOD (2.0 + DSIGN (1.0, X), 3.0);
479 
480 // C---------Hagerty's FL5 function follows...
481 
482  double q = 1.0;
483 
484  if (x < 0.0)
485  q = 1.0 - ct;
486 
487  double rmax = q / (2.0 - ct);
488 
489  double t1 = 1.0 + std::floor (x);
490  t1 = (ct / q) * (t1 < 0.0 ? -t1 : t1);
491  t1 = rmax < t1 ? rmax : t1;
492  t1 = ct > t1 ? ct : t1;
493  t1 = std::floor (x + t1);
494 
495  if (x <= 0.0 || (t1 - x) < rmax)
496  return t1;
497  else
498  return t1 - 1.0;
499 }
500 
501 static inline double
502 tceil (double x, double ct)
503 {
504  return -tfloor (-x, ct);
505 }
506 
507 static inline bool
508 teq (double u, double v,
509  double ct = 3.0 * std::numeric_limits<double>::epsilon ())
510 {
511  double tu = fabs (u);
512  double tv = fabs (v);
513 
514  return fabs (u - v) < ((tu > tv ? tu : tv) * ct);
515 }
516 
519 {
520  octave_idx_type retval = -1;
521 
522  if (rng_inc == 0
523  || (rng_limit > rng_base && rng_inc < 0)
524  || (rng_limit < rng_base && rng_inc > 0))
525  {
526  retval = 0;
527  }
528  else
529  {
530  double ct = 3.0 * std::numeric_limits<double>::epsilon ();
531 
532  double tmp = tfloor ((rng_limit - rng_base + rng_inc) / rng_inc, ct);
533 
534  octave_idx_type n_elt = (tmp > 0.0 ? static_cast<octave_idx_type> (tmp)
535  : 0);
536 
537  // If the final element that we would compute for the range is
538  // equal to the limit of the range, or is an adjacent floating
539  // point number, accept it. Otherwise, try a range with one
540  // fewer element. If that fails, try again with one more
541  // element.
542  //
543  // I'm not sure this is very good, but it seems to work better than
544  // just using tfloor as above. For example, without it, the
545  // expression 1.8:0.05:1.9 fails to produce the expected result of
546  // [1.8, 1.85, 1.9].
547 
548  if (! teq (rng_base + (n_elt - 1) * rng_inc, rng_limit))
549  {
550  if (teq (rng_base + (n_elt - 2) * rng_inc, rng_limit))
551  n_elt--;
552  else if (teq (rng_base + n_elt * rng_inc, rng_limit))
553  n_elt++;
554  }
555 
557  ? n_elt : -1;
558  }
559 
560  return retval;
561 }
562 
563 double
565 {
566  double tmp_limit = rng_limit;
567 
568  if (rng_inc > 0)
569  tmp_limit = max ();
570  else
571  tmp_limit = min ();
572 
573  return (tmp_limit != rng_limit) ? tmp_limit : rng_limit;
574 }
575 
576 void
578 {
581 
582  clear_cache ();
583 }
Matrix diag(octave_idx_type k=0) const
Definition: dMatrix.cc:2471
static bool teq(double u, double v, double ct=3.0 *std::numeric_limits< double >::epsilon())
Definition: Range.cc:508
bool is_empty(void) const
Definition: Array.h:575
Matrix cache
Definition: Range.h:142
std::ostream & operator<<(std::ostream &os, const Range &a)
Definition: Range.cc:340
octave_idx_type length(octave_idx_type n=0) const
Definition: idx-vector.h:541
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition: dMatrix.h:145
bool is_vector(void) const
Definition: dim-vector.h:458
double rng_limit
Definition: Range.h:137
Range operator-(const Range &r)
Definition: Range.cc:380
sortmode
Definition: oct-sort.h:105
double limit
Definition: Range.cc:121
Return the CPU time used by your Octave session The first output is the total time spent executing your process and is equal to the sum of second and third which are the number of CPU seconds spent executing in user mode and the number of CPU seconds spent executing in system mode
Definition: data.cc:6386
bool isnan(double x)
Definition: lo-mappers.cc:347
for large enough k
Definition: lu.cc:606
Range operator+(double x, const Range &r)
Definition: Range.cc:385
Definition: Range.h:33
void err_index_out_of_range(int nd, int dim, octave_idx_type idx, octave_idx_type ext)
double limit_internal(void) const
Definition: Range.cc:564
u
Definition: lu.cc:138
double rng_base
Definition: Range.h:136
void init(void)
Definition: Range.cc:577
void operator()(octave_idx_type i)
Definition: Range.cc:109
double min(void) const
Definition: Range.cc:164
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:398
bool swap
Definition: load-save.cc:725
Range sort(octave_idx_type dim=0, sortmode mode=ASCENDING) const
Definition: Range.cc:257
void loop(octave_idx_type n, Functor body) const
Definition: idx-vector.h:827
dim_vector orig_dimensions(void) const
Definition: idx-vector.h:583
void sort_internal(bool ascending=true)
Definition: Range.cc:211
double max(void) const
Definition: Range.cc:185
Range operator*(double x, const Range &r)
Definition: Range.cc:421
double limit(void) const
Definition: Range.h:79
void set_base(double b)
Definition: Range.cc:307
double inc(void) const
Definition: Range.h:80
octave_idx_type numel(void) const
Definition: Range.h:85
Matrix diag(octave_idx_type k=0) const
Definition: Range.cc:251
void resize(const dim_vector &dv, const T &rfv)
Definition: Array.cc:1028
double tmp
Definition: data.cc:6300
double elem(octave_idx_type i) const
Definition: Range.cc:88
octave_value retval
Definition: data.cc:6294
std::istream & operator>>(std::istream &is, Range &a)
Definition: Range.cc:361
sortmode is_sorted(sortmode mode=ASCENDING) const
Definition: Range.cc:294
Matrix matrix_value(void) const
Definition: Range.cc:52
void set_limit(double l)
Definition: Range.cc:318
Definition: dMatrix.h:37
octave_idx_type rng_numel
Definition: Range.h:140
With real return the complex result
Definition: data.cc:3375
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:228
Array< double > index(const idx_vector &i) const
Definition: Range.cc:127
void clear(void)
Definition: Array.cc:95
octave_idx_type extent(octave_idx_type n) const
Definition: idx-vector.h:544
void set_inc(double i)
Definition: Range.cc:329
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
octave_idx_type nint_big(double x)
Definition: lo-mappers.cc:409
octave_idx_type nmax
Definition: Range.cc:122
static double tfloor(double x, double ct)
Definition: Range.cc:472
bool all_elements_are_ints(void) const
Definition: Range.cc:40
octave_idx_type numel_internal(void) const
Definition: Range.cc:518
b
Definition: cellfun.cc:398
double base(void) const
Definition: Range.h:78
double floor(double x)
Definition: lo-mappers.cc:330
double rng_inc
Definition: Range.h:138
double * array
Definition: Range.cc:121
const T * fortran_vec(void) const
Definition: Array.h:584
MArray< T > reshape(const dim_vector &new_dims) const
Definition: MArray.h:91
double checkelem(octave_idx_type i) const
Definition: Range.cc:74
write the output to stdout if nargout is
Definition: load-save.cc:1576
static double tceil(double x, double ct)
Definition: Range.cc:502
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
__rangeidx_helper(double *a, double b, double i, double l, octave_idx_type n)
Definition: Range.cc:106
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE * x
bool is_colon(void) const
Definition: idx-vector.h:565
void clear_cache(void) const
Definition: Range.h:150