# HG changeset patch
# User jwe
# Date 813037166 0
# Node ID a272c4056bab95c261dd4f9ced45f4e8bad86e92
# Parent  8bdfa6fe5d6978d09806a897260f8d24f695a2de
[project @ 1995-10-07 03:38:40 by jwe]

diff --git a/liboctave/DASSL.cc b/liboctave/DASSL.cc
--- a/liboctave/DASSL.cc
+++ b/liboctave/DASSL.cc
@@ -102,7 +102,7 @@
     info [i] = 0;
 }
 
-DAE::DAE (const Vector& state, double time, DAEFunc& f)
+DAE::DAE (const ColumnVector& state, double time, DAEFunc& f)
 {
   n = state.capacity ();
   t = time;
@@ -129,7 +129,8 @@
     info [i] = 0;
 }
 
-DAE::DAE (const Vector& state, const Vector& deriv, double time, DAEFunc& f)
+DAE::DAE (const ColumnVector& state, const ColumnVector& deriv,
+	  double time, DAEFunc& f)
 {
   if (deriv.capacity () != state.capacity ())
     {
@@ -169,14 +170,14 @@
   delete iwork;
 }
 
-Vector
+ColumnVector
 DAE::deriv (void)
 {
   return xdot;
 }
 
 void
-DAE::initialize (const Vector& state, double time)
+DAE::initialize (const ColumnVector& state, double time)
 {
   integration_error = 0;
   restart = 1;
@@ -187,7 +188,8 @@
 }
 
 void
-DAE::initialize (const Vector& state, const Vector& deriv, double time)
+DAE::initialize (const ColumnVector& state,
+		 const ColumnVector& deriv, double time)
 {
   integration_error = 0;
   restart = 1;
@@ -200,9 +202,9 @@
 ddassl_f (const double& time, double *state, double *deriv,
 	  double *delta, int& ires, double *, int *)
 {
-  Vector tmp_deriv (nn);
-  Vector tmp_state (nn);
-  Vector tmp_delta (nn);
+  ColumnVector tmp_deriv (nn);
+  ColumnVector tmp_state (nn);
+  ColumnVector tmp_delta (nn);
 
   for (int i = 0; i < nn; i++)
     {
@@ -227,8 +229,8 @@
 ddassl_j (const double& time, double *, double *, double *pd, const
 	  double& cj, double *, int *)
 {
-  Vector tmp_state (nn);
-  Vector tmp_deriv (nn);
+  ColumnVector tmp_state (nn);
+  ColumnVector tmp_deriv (nn);
 
   // XXX FIXME XXX
 
@@ -252,7 +254,7 @@
   return 0;
 }
 
-Vector
+ColumnVector
 DAE::integrate (double tout)
 {
   integration_error = 0;
@@ -361,7 +363,7 @@
 }
 
 Matrix
-DAE::integrate (const Vector& tout, Matrix& xdot_out)
+DAE::integrate (const ColumnVector& tout, Matrix& xdot_out)
 {
   Matrix retval;
   int n_out = tout.capacity ();
@@ -396,7 +398,8 @@
 }
 
 Matrix
-DAE::integrate (const Vector& tout, Matrix& xdot_out, const Vector& tcrit) 
+DAE::integrate (const ColumnVector& tout, Matrix& xdot_out,
+		const ColumnVector& tcrit) 
 {
   Matrix retval;
   int n_out = tout.capacity ();
diff --git a/liboctave/LinConst.cc b/liboctave/LinConst.cc
--- a/liboctave/LinConst.cc
+++ b/liboctave/LinConst.cc
@@ -42,8 +42,8 @@
   (*current_liboctave_error_handler) ("fatal LinConst error: %s", msg);
 }
 
-LinConst::LinConst (const Matrix& a_eq, const Vector& b_eq,
-		    const Matrix& a_ineq, const Vector& b_ineq)
+LinConst::LinConst (const Matrix& a_eq, const ColumnVector& b_eq,
+		    const Matrix& a_ineq, const ColumnVector& b_ineq)
 {
   // Need some checks here.
 
@@ -123,10 +123,10 @@
   return retval;
 }
 
-Vector
+ColumnVector
 LinConst::eq_constraint_vector (void) const
 {
-  Vector retval (nb);
+  ColumnVector retval (nb);
   int count = 0;
   for (int i = 0; i < nb; i++)
     {
@@ -140,10 +140,10 @@
   return retval;
 }
 
-Vector
+ColumnVector
 LinConst::ineq_constraint_vector (void) const
 {
-  Vector retval (2*nb);
+  ColumnVector retval (2*nb);
   int count = 0;
   for (int i = 0; i < nb; i++)
     {
diff --git a/liboctave/Quad.cc b/liboctave/Quad.cc
--- a/liboctave/Quad.cc
+++ b/liboctave/Quad.cc
@@ -29,9 +29,6 @@
 #include <config.h>
 #endif
 
-#include <cfloat>
-#include <cmath>
-
 #include "Quad.h"
 #include "f77-uscore.h"
 #include "sun-utils.h"
diff --git a/liboctave/Range.cc b/liboctave/Range.cc
--- a/liboctave/Range.cc
+++ b/liboctave/Range.cc
@@ -29,7 +29,9 @@
 #include <config.h>
 #endif
 
+#include <cfloat>
 #include <climits>
+#include <cmath>
 
 #include <iostream.h>
 
@@ -136,6 +138,81 @@
   return is;
 }
 
+// C  See Knuth, Art Of Computer Programming, Vol. 1, Problem 1.2.4-5.
+// C
+// C===Tolerant FLOOR function.
+// C
+// C    X  -  is given as a Double Precision argument to be operated on.
+// C          It is assumed that X is represented with M mantissa bits.
+// C    CT -  is   given   as   a   Comparison   Tolerance   such   that
+// C          0.LT.CT.LE.3-SQRT(5)/2. If the relative difference between
+// C          X and A whole number is  less  than  CT,  then  TFLOOR  is
+// C          returned   as   this   whole   number.   By  treating  the
+// C          floating-point numbers as a finite ordered set  note  that
+// C          the  heuristic  EPS=2.**(-(M-1))   and   CT=3*EPS   causes
+// C          arguments  of  TFLOOR/TCEIL to be treated as whole numbers
+// C          if they are  exactly  whole  numbers  or  are  immediately
+// C          adjacent to whole number representations.  Since EPS,  the
+// C          "distance"  between  floating-point  numbers  on  the unit
+// C          interval, and M, the number of bits in X'S mantissa, exist
+// C          on  every  floating-point   computer,   TFLOOR/TCEIL   are
+// C          consistently definable on every floating-point computer.
+// C
+// C          For more information see the following references:
+// C    (1) P. E. Hagerty, "More On Fuzzy Floor And Ceiling," APL  QUOTE
+// C        QUAD 8(4):20-24, June 1978. Note that TFLOOR=FL5.
+// C    (2) L. M. Breed, "Definitions For Fuzzy Floor And Ceiling",  APL
+// C        QUOTE QUAD 8(3):16-23, March 1978. This paper cites FL1 through
+// C        FL5, the history of five years of evolutionary development of
+// C        FL5 - the seven lines of code below - by open collaboration
+// C        and corroboration of the mathematical-computing community.
+// C
+// C  Penn State University Center for Academic Computing
+// C  H. D. Knoble - August, 1978.
+
+static inline double
+tfloor (double x, double ct)
+{
+// C---------FLOOR(X) is the largest integer algebraically less than
+// C         or equal to X; that is, the unfuzzy FLOOR function.
+
+//  DINT (X) = X - DMOD (X, 1.0);
+//  FLOOR (X) = DINT (X) - DMOD (2.0 + DSIGN (1.0, X), 3.0);
+
+// C---------Hagerty's FL5 function follows...
+
+  double q = 1.0;
+
+  if (x < 0.0)
+    q = 1.0 - ct;
+
+  double rmax = q / (2.0 - ct);
+
+  double t1 = 1.0 + floor (x);
+  t1 = (ct / q) * (t1 < 0.0 ? -t1 : t1);
+  t1 = rmax < t1 ? rmax : t1;
+  t1 = ct > t1 ? ct : t1;
+  t1 = floor (x + t1);
+
+  if (x <= 0.0 || t1 - x < rmax)
+    return t1;
+  else
+    return t1 - 1.0;
+}
+
+static inline double
+tceil (double x, double ct)
+{
+  return -tfloor (-x, ct);
+}
+
+static inline double
+round (double x, double ct)
+{
+  return tfloor (x+0.5, ct);
+}
+
+
 // Find an approximate number of intervals, then do the best we can to
 // find the number of intervals that we would get if we had done
 // something like
@@ -152,65 +229,13 @@
 int
 Range::nelem_internal (void) const
 {
-  // We can't have more than INT_MAX elements in the range.
-
-  double d_n_intervals = (rng_limit - rng_base) / rng_inc;
-  int max_intervals = INT_MAX - 1;
-  double d_max_val = (double) max_intervals;
-
-  if (d_n_intervals > d_max_val)
-    return -1;
+  double ct = 3.0 * DBL_EPSILON;
 
-  int n_intervals = (d_n_intervals > 0)
-    ? ((int) (d_n_intervals + 0.5))
-    : ((int) (d_n_intervals - 0.5)); 
-
-  if (rng_limit > rng_base && rng_inc > 0)
-    {
-      // Our approximation may have been too big.
-
-      while (rng_base + n_intervals * rng_inc > rng_limit && n_intervals > 0)
-	n_intervals--;
-
-      // Now that we are close, get the actual number.  Try to avoid
-      // problems with extended precision registers.
+  double tmp = round ((rng_limit - rng_base + rng_inc) / rng_inc, ct);
 
-      for (;;)
-	{
-	  volatile double tmp_inc = (n_intervals + 1) * rng_inc;
-	  volatile double tmp_val = rng_base + tmp_inc;
-	  if (tmp_val <= rng_limit && n_intervals < max_intervals)
-	    n_intervals++;
-	  else
-	    break;
-	}
-    }
-  else if (rng_limit < rng_base && rng_inc < 0)
-    {
-      // Our approximation may have been too big.
-
-      while (rng_base + n_intervals * rng_inc < rng_limit && n_intervals > 0)
-	n_intervals--;
+  int n_intervals = (int) (tmp > 0.0 ? tmp : 0);
 
-      // Now that we are close, get the actual number.  Try to avoid
-      // problems with extended precision registers.
-
-      for (;;)
-	{
-	  volatile double tmp_inc = (n_intervals + 1) * rng_inc;
-	  volatile double tmp_val = rng_base + tmp_inc;
-	  if (tmp_val >= rng_limit && n_intervals < max_intervals)
-	    n_intervals++;
-	  else
-	    break;
-	}
-    }
-  else if (rng_limit == rng_base)
-    n_intervals = 0;
-  else
-    n_intervals = -1;
-
-  return (n_intervals >= max_intervals) ? -1 : n_intervals + 1;
+  return (n_intervals >= INT_MAX - 1) ? -1 : n_intervals;
 }
 
 /*