GNU Octave  3.8.0
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-2013 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 #ifdef 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 
46  return (! (xisnan (rng_base) || xisnan (rng_inc))
47  && (NINTbig (rng_base) == rng_base || rng_nelem < 1)
48  && (NINTbig (rng_inc) == rng_inc || rng_nelem <= 1));
49 }
50 
51 Matrix
52 Range::matrix_value (void) const
53 {
54  if (rng_nelem > 0 && cache.nelem () == 0)
55  {
56  cache.resize (1, rng_nelem);
57  double b = rng_base;
58  double increment = rng_inc;
59  if (rng_nelem > 0)
60  {
61  // The first element must always be *exactly* the base.
62  // E.g, -0 would otherwise become +0 in the loop (-0 + 0*increment).
63  cache(0) = b;
64  for (octave_idx_type i = 1; i < rng_nelem; i++)
65  cache(i) = b + i * increment;
66  }
67 
68  // On some machines (x86 with extended precision floating point
69  // arithmetic, for example) it is possible that we can overshoot
70  // the limit by approximately the machine precision even though
71  // we were very careful in our calculation of the number of
72  // elements. The tests need equality (>= rng_limit or <= rng_limit)
73  // to have expressions like -5:1:-0 result in a -0 endpoint.
74 
75  if ((rng_inc > 0 && cache(rng_nelem-1) >= rng_limit)
76  || (rng_inc < 0 && cache(rng_nelem-1) <= rng_limit))
78  }
79 
80  return cache;
81 }
82 
83 double
85 {
86  if (i < 0 || i >= rng_nelem)
88 
89  if (i == 0)
90  return rng_base;
91  else if (i < rng_nelem - 1)
92  return rng_base + i * rng_inc;
93  else
94  {
95  double end = rng_base + i * rng_inc;
96  if ((rng_inc > 0 && end >= rng_limit)
97  || (rng_inc < 0 && end <= rng_limit))
98  return rng_limit;
99  else
100  return end;
101  }
102 }
103 
104 double
106 {
107 #if defined (BOUNDS_CHECKING)
108  return checkelem (i);
109 #else
110  if (i == 0)
111  return rng_base;
112  else if (i < rng_nelem - 1)
113  return rng_base + i * rng_inc;
114  else
115  {
116  double end = rng_base + i * rng_inc;
117  if ((rng_inc > 0 && end >= rng_limit)
118  || (rng_inc < 0 && end <= rng_limit))
119  return rng_limit;
120  else
121  return end;
122  }
123 #endif
124 }
125 
126 // Helper class used solely for idx_vector.loop () function call
128 {
129 public:
130  __rangeidx_helper (double *a, double b, double i, double l, octave_idx_type n)
131  : array (a), base (b), inc (i), limit (l), nmax (n-1) { }
132 
134  {
135  if (i == 0)
136  *array++ = base;
137  else if (i < nmax)
138  *array++ = base + i * inc;
139  else
140  {
141  double end = base + i * inc;
142  if ((inc > 0 && end >= limit) || (inc < 0 && end <= limit))
143  *array++ = limit;
144  else
145  *array++ = end;
146  }
147  }
148 
149 private:
150 
151  double *array, base, inc, limit;
153 
154 };
155 
157 Range::index (const idx_vector& i) const
158 {
159  Array<double> retval;
160 
162 
163  if (i.is_colon ())
164  {
165  retval = matrix_value ().reshape (dim_vector (rng_nelem, 1));
166  }
167  else
168  {
169  if (i.extent (n) != n)
170  gripe_index_out_of_range (1, 1, i.extent (n), n); // throws
171 
172  dim_vector rd = i.orig_dimensions ();
173  octave_idx_type il = i.length (n);
174 
175  // taken from Array.cc.
176  if (n != 1 && rd.is_vector ())
177  rd = dim_vector (1, il);
178 
179  retval.clear (rd);
180 
181  // idx_vector loop across all values in i,
182  // executing __rangeidx_helper (i) for each i
183  i.loop (n, __rangeidx_helper (retval.fortran_vec (),
185  }
186 
187  return retval;
188 }
189 
190 // NOTE: max and min only return useful values if nelem > 0.
191 // do_minmax_body() in max.cc avoids calling Range::min/max if nelem == 0.
192 
193 double
194 Range::min (void) const
195 {
196  double retval = 0.0;
197  if (rng_nelem > 0)
198  {
199  if (rng_inc > 0)
200  retval = rng_base;
201  else
202  {
203  retval = rng_base + (rng_nelem - 1) * rng_inc;
204 
205  // See the note in the matrix_value method above.
206  if (retval <= rng_limit)
207  retval = rng_limit;
208  }
209 
210  }
211  return retval;
212 }
213 
214 double
215 Range::max (void) const
216 {
217  double retval = 0.0;
218  if (rng_nelem > 0)
219  {
220  if (rng_inc > 0)
221  {
222  retval = rng_base + (rng_nelem - 1) * rng_inc;
223 
224  // See the note in the matrix_value method above.
225  if (retval >= rng_limit)
226  retval = rng_limit;
227  }
228  else
229  retval = rng_base;
230  }
231  return retval;
232 }
233 
234 void
235 Range::sort_internal (bool ascending)
236 {
237  if (ascending && rng_base > rng_limit && rng_inc < 0.0)
238  {
239  double tmp = rng_base;
240  rng_base = min ();
241  rng_limit = tmp;
242  rng_inc = -rng_inc;
243  clear_cache ();
244  }
245  else if (! ascending && rng_base < rng_limit && rng_inc > 0.0)
246  {
247  double tmp = rng_limit;
248  rng_limit = min ();
249  rng_base = tmp;
250  rng_inc = -rng_inc;
251  clear_cache ();
252  }
253 }
254 
255 void
257 {
258  octave_idx_type nel = nelem ();
259 
260  sidx.resize (dim_vector (1, nel));
261 
262  octave_idx_type *psidx = sidx.fortran_vec ();
263 
264  bool reverse = false;
265 
266  if (ascending && rng_base > rng_limit && rng_inc < 0.0)
267  {
268  double tmp = rng_base;
269  rng_base = min ();
270  rng_limit = tmp;
271  rng_inc = -rng_inc;
272  clear_cache ();
273  reverse = true;
274  }
275  else if (! ascending && rng_base < rng_limit && rng_inc > 0.0)
276  {
277  double tmp = rng_limit;
278  rng_limit = min ();
279  rng_base = tmp;
280  rng_inc = -rng_inc;
281  clear_cache ();
282  reverse = true;
283  }
284 
285  octave_idx_type tmp = reverse ? nel - 1 : 0;
286  octave_idx_type stp = reverse ? -1 : 1;
287 
288  for (octave_idx_type i = 0; i < nel; i++, tmp += stp)
289  psidx[i] = tmp;
290 
291 }
292 
293 Matrix
295 {
296  return matrix_value ().diag (k);
297 }
298 
299 Range
301 {
302  Range retval = *this;
303 
304  if (dim == 1)
305  {
306  if (mode == ASCENDING)
307  retval.sort_internal (true);
308  else if (mode == DESCENDING)
309  retval.sort_internal (false);
310  }
311  else if (dim != 0)
312  (*current_liboctave_error_handler) ("Range::sort: invalid dimension");
313 
314  return retval;
315 }
316 
317 Range
319  sortmode mode) const
320 {
321  Range retval = *this;
322 
323  if (dim == 1)
324  {
325  if (mode == ASCENDING)
326  retval.sort_internal (sidx, true);
327  else if (mode == DESCENDING)
328  retval.sort_internal (sidx, false);
329  }
330  else if (dim != 0)
331  (*current_liboctave_error_handler) ("Range::sort: invalid dimension");
332 
333  return retval;
334 }
335 
336 sortmode
338 {
339  if (rng_nelem > 1 && rng_inc > 0)
340  mode = (mode == DESCENDING) ? UNSORTED : ASCENDING;
341  else if (rng_nelem > 1 && rng_inc < 0)
342  mode = (mode == ASCENDING) ? UNSORTED : DESCENDING;
343  else
344  mode = mode ? mode : ASCENDING;
345 
346  return mode;
347 }
348 
349 std::ostream&
350 operator << (std::ostream& os, const Range& a)
351 {
352  double b = a.base ();
353  double increment = a.inc ();
354  octave_idx_type num_elem = a.nelem ();
355 
356  if (num_elem > 1)
357  {
358  // First element must be the base *exactly* (-0).
359  os << b << " ";
360  for (octave_idx_type i = 1; i < num_elem-1; i++)
361  os << b + i * increment << " ";
362  }
363 
364  // Prevent overshoot. See comment in the matrix_value method above.
365  os << (increment > 0 ? a.max () : a.min ()) << "\n";
366 
367  return os;
368 }
369 
370 std::istream&
371 operator >> (std::istream& is, Range& a)
372 {
373  is >> a.rng_base;
374  if (is)
375  {
376  is >> a.rng_limit;
377  if (is)
378  {
379  is >> a.rng_inc;
380  a.rng_nelem = a.nelem_internal ();
381  }
382  }
383 
384  return is;
385 }
386 
387 Range
389 {
390  return Range (-r.base (), -r.inc (), r.nelem ());
391 }
392 
393 Range operator + (double x, const Range& r)
394 {
395  Range result (x + r.base (), r.inc (), r.nelem ());
396  if (result.rng_nelem < 0)
397  result.cache = x + r.matrix_value ();
398 
399  return result;
400 }
401 
402 Range operator + (const Range& r, double x)
403 {
404  Range result (r.base () + x, r.inc (), r.nelem ());
405  if (result.rng_nelem < 0)
406  result.cache = r.matrix_value () + x;
407 
408  return result;
409 }
410 
411 Range operator - (double x, const Range& r)
412 {
413  Range result (x - r.base (), -r.inc (), r.nelem ());
414  if (result.rng_nelem < 0)
415  result.cache = x - r.matrix_value ();
416 
417  return result;
418 }
419 
420 Range operator - (const Range& r, double x)
421 {
422  Range result (r.base () - x, r.inc (), r.nelem ());
423  if (result.rng_nelem < 0)
424  result.cache = r.matrix_value () - x;
425 
426  return result;
427 }
428 
429 Range operator * (double x, const Range& r)
430 {
431  Range result (x * r.base (), x * r.inc (), r.nelem ());
432  if (result.rng_nelem < 0)
433  result.cache = x * r.matrix_value ();
434 
435  return result;
436 }
437 
438 Range operator * (const Range& r, double x)
439 {
440  Range result (r.base () * x, r.inc () * x, r.nelem ());
441  if (result.rng_nelem < 0)
442  result.cache = r.matrix_value () * x;
443 
444  return result;
445 }
446 
447 
448 // C See Knuth, Art Of Computer Programming, Vol. 1, Problem 1.2.4-5.
449 // C
450 // C===Tolerant FLOOR function.
451 // C
452 // C X - is given as a Double Precision argument to be operated on.
453 // C It is assumed that X is represented with M mantissa bits.
454 // C CT - is given as a Comparison Tolerance such that
455 // C 0.LT.CT.LE.3-SQRT(5)/2. If the relative difference between
456 // C X and A whole number is less than CT, then TFLOOR is
457 // C returned as this whole number. By treating the
458 // C floating-point numbers as a finite ordered set note that
459 // C the heuristic EPS=2.**(-(M-1)) and CT=3*EPS causes
460 // C arguments of TFLOOR/TCEIL to be treated as whole numbers
461 // C if they are exactly whole numbers or are immediately
462 // C adjacent to whole number representations. Since EPS, the
463 // C "distance" between floating-point numbers on the unit
464 // C interval, and M, the number of bits in X'S mantissa, exist
465 // C on every floating-point computer, TFLOOR/TCEIL are
466 // C consistently definable on every floating-point computer.
467 // C
468 // C For more information see the following references:
469 // C (1) P. E. Hagerty, "More On Fuzzy Floor And Ceiling," APL QUOTE
470 // C QUAD 8(4):20-24, June 1978. Note that TFLOOR=FL5.
471 // C (2) L. M. Breed, "Definitions For Fuzzy Floor And Ceiling", APL
472 // C QUOTE QUAD 8(3):16-23, March 1978. This paper cites FL1 through
473 // C FL5, the history of five years of evolutionary development of
474 // C FL5 - the seven lines of code below - by open collaboration
475 // C and corroboration of the mathematical-computing community.
476 // C
477 // C Penn State University Center for Academic Computing
478 // C H. D. Knoble - August, 1978.
479 
480 static inline double
481 tfloor (double x, double ct)
482 {
483 // C---------FLOOR(X) is the largest integer algebraically less than
484 // C or equal to X; that is, the unfuzzy FLOOR function.
485 
486 // DINT (X) = X - DMOD (X, 1.0);
487 // FLOOR (X) = DINT (X) - DMOD (2.0 + DSIGN (1.0, X), 3.0);
488 
489 // C---------Hagerty's FL5 function follows...
490 
491  double q = 1.0;
492 
493  if (x < 0.0)
494  q = 1.0 - ct;
495 
496  double rmax = q / (2.0 - ct);
497 
498  double t1 = 1.0 + gnulib::floor (x);
499  t1 = (ct / q) * (t1 < 0.0 ? -t1 : t1);
500  t1 = rmax < t1 ? rmax : t1;
501  t1 = ct > t1 ? ct : t1;
502  t1 = gnulib::floor (x + t1);
503 
504  if (x <= 0.0 || (t1 - x) < rmax)
505  return t1;
506  else
507  return t1 - 1.0;
508 }
509 
510 static inline double
511 tceil (double x, double ct)
512 {
513  return -tfloor (-x, ct);
514 }
515 
516 static inline bool
517 teq (double u, double v,
518  double ct = 3.0 * std::numeric_limits<double>::epsilon ())
519 {
520  double tu = fabs (u);
521  double tv = fabs (v);
522 
523  return fabs (u - v) < ((tu > tv ? tu : tv) * ct);
524 }
525 
528 {
529  octave_idx_type retval = -1;
530 
531  if (rng_inc == 0
532  || (rng_limit > rng_base && rng_inc < 0)
533  || (rng_limit < rng_base && rng_inc > 0))
534  {
535  retval = 0;
536  }
537  else
538  {
539  double ct = 3.0 * std::numeric_limits<double>::epsilon ();
540 
541  double tmp = tfloor ((rng_limit - rng_base + rng_inc) / rng_inc, ct);
542 
543  octave_idx_type n_elt = (tmp > 0.0 ? static_cast<octave_idx_type> (tmp)
544  : 0);
545 
546  // If the final element that we would compute for the range is
547  // equal to the limit of the range, or is an adjacent floating
548  // point number, accept it. Otherwise, try a range with one
549  // fewer element. If that fails, try again with one more
550  // element.
551  //
552  // I'm not sure this is very good, but it seems to work better than
553  // just using tfloor as above. For example, without it, the
554  // expression 1.8:0.05:1.9 fails to produce the expected result of
555  // [1.8, 1.85, 1.9].
556 
557  if (! teq (rng_base + (n_elt - 1) * rng_inc, rng_limit))
558  {
559  if (teq (rng_base + (n_elt - 2) * rng_inc, rng_limit))
560  n_elt--;
561  else if (teq (rng_base + n_elt * rng_inc, rng_limit))
562  n_elt++;
563  }
564 
566  ? n_elt : -1;
567  }
568 
569  return retval;
570 }