--- a/liboctave/ChangeLog	Wed May 21 19:25:08 2008 +0200
+++ b/liboctave/ChangeLog	Mon May 19 10:22:38 2008 +0200
@@ -1,3 +1,16 @@
+2008-05-21  David Bateman  <dbateman@free.fr>
+
+	* Quad-opts.in: Handle single precision tolerances.
+	* Quad.cc (float_user_fcn): New static variable.
+	(quad_float_fcn_ptr): New typedef.
+	(qagp, quagi): New QUADPACK decls.
+	(float_user_function): New function.
+	(DefQuad::do_integrate, IndefQuad::do_integrate): Float versions.
+	(FloatDefQuad::do_integrate, FloatIndefQuad::do_integrate):
+	New functions.
+	* Quad.h (class Quad): Handle float type.
+	(class FloatDefQuad, class FloatIndefQuad): New classes.
+
 2008-05-21  Jaroslav Hajek <highegg@gmail.com>
 
 	* fCMatrix.h (xgemm): Provide decl.

--- a/liboctave/Quad-opts.in	Wed May 21 19:25:08 2008 +0200
+++ b/liboctave/Quad-opts.in	Mon May 19 10:22:38 2008 +0200
@@ -39,3 +39,26 @@
   INIT_VALUE = "::sqrt (DBL_EPSILON)"
   SET_EXPR = "val"
 END_OPTION
+
+OPTION
+  NAME = "single precision absolute tolerance"
+  DOC_ITEM
+Absolute tolerance for single precision; may be zero for pure relative 
+error test.
+  END_DOC_ITEM
+  TYPE = "float"
+  INIT_VALUE = "::sqrt (FLT_EPSILON)"
+  SET_EXPR = "val"
+END_OPTION
+
+OPTION
+  NAME = "single precision relative tolerance"
+  DOC_ITEM
+Nonnegative relative tolerance for single precision.  If the absolute
+tolerance is zero, the relative tolerance must be greater than or equal to 
+@code{max (50*eps, 0.5e-28)}.
+  END_DOC_ITEM
+  TYPE = "float"
+  INIT_VALUE = "::sqrt (FLT_EPSILON)"
+  SET_EXPR = "val"
+END_OPTION

--- a/liboctave/Quad.cc	Wed May 21 19:25:08 2008 +0200
+++ b/liboctave/Quad.cc	Mon May 19 10:22:38 2008 +0200
@@ -32,6 +32,7 @@
 #include "sun-utils.h"
 
 static integrand_fcn user_fcn;
+static float_integrand_fcn float_user_fcn;
 
 // FIXME -- would be nice to not have to have this global
 // variable.
@@ -40,6 +41,7 @@
 int quad_integration_error = 0;
 
 typedef octave_idx_type (*quad_fcn_ptr) (double*, int&, double*);
+typedef octave_idx_type (*quad_float_fcn_ptr) (float*, int&, float*);
 			      
 extern "C"
 {
@@ -55,6 +57,19 @@
 			   const double&, const double&, double&,
 			   double&, octave_idx_type&, octave_idx_type&, const octave_idx_type&,
 			   const octave_idx_type&, octave_idx_type&, octave_idx_type*, double*); 
+
+  F77_RET_T
+  F77_FUNC (qagp, QAGP) (quad_float_fcn_ptr, const float&, const float&,
+			 const octave_idx_type&, const float*, const float&,
+			 const float&, float&, float&, octave_idx_type&,
+			 octave_idx_type&, const octave_idx_type&, const octave_idx_type&, octave_idx_type&, octave_idx_type*,
+			 float*);
+
+  F77_RET_T
+  F77_FUNC (qagi, QAGI) (quad_float_fcn_ptr, const float&, const octave_idx_type&,
+			 const float&, const float&, float&,
+			 float&, octave_idx_type&, octave_idx_type&, const octave_idx_type&,
+			 const octave_idx_type&, octave_idx_type&, octave_idx_type*, float*); 
 }
 
 static octave_idx_type
@@ -86,6 +101,35 @@
   return 0;
 }
 
+static octave_idx_type
+float_user_function (float *x, int& ierr, float *result)
+{
+  BEGIN_INTERRUPT_WITH_EXCEPTIONS;
+
+#if defined (sun) && defined (__GNUC__)
+  float xx = access_float (x);
+#else
+  float xx = *x;
+#endif
+
+  quad_integration_error = 0;
+
+  float xresult = (*float_user_fcn) (xx);
+
+#if defined (sun) && defined (__GNUC__)
+  assign_float (result, xresult);
+#else
+  *result = xresult;
+#endif
+
+  if (quad_integration_error)
+    ierr = -1;
+
+  END_INTERRUPT_WITH_EXCEPTIONS;
+
+  return 0;
+}
+
 double
 DefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr)
 {
@@ -115,6 +159,13 @@
   return result;
 }
 
+float
+DefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr)
+{
+  (*current_liboctave_error_handler) ("incorrect integration function called");
+  return 0.0;
+}
+
 double
 IndefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr)
 {
@@ -161,6 +212,102 @@
   return result;
 }
 
+float
+IndefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr)
+{
+  (*current_liboctave_error_handler) ("incorrect integration function called");
+  return 0.0;
+}
+
+double
+FloatDefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr)
+{
+  (*current_liboctave_error_handler) ("incorrect integration function called");
+  return 0.0;
+}
+
+float
+FloatDefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr)
+{
+  octave_idx_type npts = singularities.capacity () + 2;
+  float *points = singularities.fortran_vec ();
+  float result = 0.0;
+
+  octave_idx_type leniw = 183*npts - 122;
+  Array<octave_idx_type> iwork (leniw);
+  octave_idx_type *piwork = iwork.fortran_vec ();
+
+  octave_idx_type lenw = 2*leniw - npts;
+  Array<float> work (lenw);
+  float *pwork = work.fortran_vec ();
+
+  float_user_fcn = ff;
+  octave_idx_type last;
+
+  float abs_tol = single_precision_absolute_tolerance ();
+  float rel_tol = single_precision_relative_tolerance ();
+
+  F77_XFCN (qagp, QAGP, (float_user_function, lower_limit, upper_limit,
+			 npts, points, abs_tol, rel_tol, result,
+			 abserr, neval, ier, leniw, lenw, last,
+			 piwork, pwork));
+
+  return result;
+}
+
+double
+FloatIndefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr)
+{
+  (*current_liboctave_error_handler) ("incorrect integration function called");
+  return 0.0;
+}
+
+float
+FloatIndefQuad::do_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr)
+{
+  float result = 0.0;
+
+  octave_idx_type leniw = 128;
+  Array<octave_idx_type> iwork (leniw);
+  octave_idx_type *piwork = iwork.fortran_vec ();
+
+  octave_idx_type lenw = 8*leniw;
+  Array<float> work (lenw);
+  float *pwork = work.fortran_vec ();
+
+  float_user_fcn = ff;
+  octave_idx_type last;
+
+  octave_idx_type inf;
+  switch (type)
+    {
+    case bound_to_inf:
+      inf = 1;
+      break;
+
+    case neg_inf_to_bound:
+      inf = -1;
+      break;
+
+    case doubly_infinite:
+      inf = 2;
+      break;
+
+    default:
+      assert (0);
+      break;
+    }
+
+  float abs_tol = single_precision_absolute_tolerance ();
+  float rel_tol = single_precision_relative_tolerance ();
+
+  F77_XFCN (qagi, QAGI, (float_user_function, bound, inf, abs_tol, rel_tol,
+			 result, abserr, neval, ier, leniw, lenw,
+			 last, piwork, pwork));
+
+  return result;
+}
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***

--- a/liboctave/Quad.h	Wed May 21 19:25:08 2008 +0200
+++ b/liboctave/Quad.h	Mon May 19 10:22:38 2008 +0200
@@ -27,12 +27,14 @@
 #include <cfloat>
 
 #include "dColVector.h"
+#include "fColVector.h"
 #include "lo-math.h"
 
 #if !defined (octave_Quad_typedefs)
 #define octave_Quad_typedefs 1
 
 typedef double (*integrand_fcn) (double x);
+typedef float (*float_integrand_fcn) (float x);
 
 #endif
 
@@ -53,6 +55,9 @@
   Quad (integrand_fcn fcn)
     : Quad_options (), f (fcn) { }
 
+  Quad (float_integrand_fcn fcn)
+    : Quad_options (), ff (fcn) { }
+
   virtual ~Quad (void) { }
 
   virtual double integrate (void)
@@ -62,6 +67,13 @@
       return do_integrate (ier, neval, abserr);
     }
 
+  virtual float float_integrate (void)
+    {
+      octave_idx_type ier, neval;
+      float abserr;
+      return do_integrate (ier, neval, abserr);
+    }
+
   virtual double integrate (octave_idx_type& ier)
     {
       octave_idx_type neval;
@@ -69,22 +81,43 @@
       return do_integrate (ier, neval, abserr);
     }
 
+  virtual float float_integrate (octave_idx_type& ier)
+    {
+      octave_idx_type neval;
+      float abserr;
+      return do_integrate (ier, neval, abserr);
+    }
+
   virtual double integrate (octave_idx_type& ier, octave_idx_type& neval)
     {
       double abserr;
       return do_integrate (ier, neval, abserr);
     }
 
+  virtual float float_integrate (octave_idx_type& ier, octave_idx_type& neval)
+    {
+      float abserr;
+      return do_integrate (ier, neval, abserr);
+    }
+
   virtual double integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr)
     {
       return do_integrate (ier, neval, abserr);
     }
 
+  virtual float float_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr)
+    {
+      return do_integrate (ier, neval, abserr);
+    }
+
   virtual double do_integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr) = 0;
 
+  virtual float do_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr) = 0;
+
  protected:
 
   integrand_fcn f;
+  float_integrand_fcn ff;
 };
 
 class
@@ -112,6 +145,8 @@
 
   double do_integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr);
 
+  float do_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr);
+
  private:
 
   double lower_limit;
@@ -138,6 +173,8 @@
 
   double do_integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr);
 
+  float do_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr);
+
  private:
 
   double bound;
@@ -145,6 +182,68 @@
   int integration_error;
 };
 
+class
+OCTAVE_API
+FloatDefQuad : public Quad
+{
+ public:
+
+  FloatDefQuad (float_integrand_fcn fcn)
+    : Quad (fcn), lower_limit (0.0), upper_limit (1.0), singularities () { }
+
+  FloatDefQuad (float_integrand_fcn fcn, float ll, float ul)
+    : Quad (fcn), lower_limit (ll), upper_limit (ul), singularities () { }
+
+  FloatDefQuad (float_integrand_fcn fcn, float ll, float ul,
+	   const FloatColumnVector& sing)
+    : Quad (fcn), lower_limit (ll), upper_limit (ul),
+      singularities (sing) { }
+
+  FloatDefQuad (float_integrand_fcn fcn, const FloatColumnVector& sing)
+    : Quad (fcn), lower_limit (0.0), upper_limit (1.0),
+      singularities (sing) { }
+
+  ~FloatDefQuad (void) { }
+
+  double do_integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr);
+
+  float do_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr);
+
+ private:
+
+  float lower_limit;
+  float upper_limit;
+
+  FloatColumnVector singularities;
+};
+
+class
+OCTAVE_API
+FloatIndefQuad : public Quad
+{
+ public:
+
+  enum IntegralType { bound_to_inf, neg_inf_to_bound, doubly_infinite };
+
+  FloatIndefQuad (float_integrand_fcn fcn)
+    : Quad (fcn), bound (0.0), type (bound_to_inf) { }
+
+  FloatIndefQuad (float_integrand_fcn fcn, double b, IntegralType t)
+    : Quad (fcn), bound (b), type (t) { }
+
+  ~FloatIndefQuad (void) { }
+
+  double do_integrate (octave_idx_type& ier, octave_idx_type& neval, double& abserr);
+
+  float do_integrate (octave_idx_type& ier, octave_idx_type& neval, float& abserr);
+
+ private:
+
+  float bound;
+  IntegralType type;
+  int integration_error;
+};
+
 #endif
 
 /*

--- a/src/ChangeLog	Wed May 21 19:25:08 2008 +0200
+++ b/src/ChangeLog	Mon May 19 10:22:38 2008 +0200
@@ -1,3 +1,9 @@
+2008-05-21  David Bateman  <dbateman@free.fr>
+
+	* DLD-FUNCTIONS/quad.cc (quad_float_user_function): New function.
+	(Fquad): Handle float type.
+	New tests.
+
 2008-05-21  Jaroslav Hajek <highegg@gmail.com>
 
 	* OPERATORS/op-fcm-fcm.cc (trans_mul, mul_trans, herm_mul, mul_herm):

--- a/src/DLD-FUNCTIONS/quad.cc	Wed May 21 19:25:08 2008 +0200
+++ b/src/DLD-FUNCTIONS/quad.cc	Mon May 19 10:22:38 2008 +0200
@@ -103,6 +103,51 @@
   return retval;
 }
 
+float
+quad_float_user_function (float x)
+{
+  float retval = 0.0;
+
+  octave_value_list args;
+  args(0) = x;
+
+  if (quad_fcn)
+    {
+      octave_value_list tmp = quad_fcn->do_multi_index_op (1, args);
+
+      if (error_state)
+	{
+	  quad_integration_error = 1;  // FIXME
+	  gripe_user_supplied_eval ("quad");
+	  return retval;
+	}
+
+      if (tmp.length () && tmp(0).is_defined ())
+	{
+	  if (! warned_imaginary && tmp(0).is_complex_type ())
+	    {
+	      warning ("quad: ignoring imaginary part returned from user-supplied function");
+	      warned_imaginary = true;
+	    }
+
+	  retval = tmp(0).float_value ();
+
+	  if (error_state)
+	    {
+	      quad_integration_error = 1;  // FIXME
+	      gripe_user_supplied_eval ("quad");
+	    }
+	}
+      else
+	{
+	  quad_integration_error = 1;  // FIXME
+	  gripe_user_supplied_eval ("quad");
+	}
+    }
+
+  return retval;
+}
+
 #define QUAD_ABORT() \
   do \
     { \
@@ -203,110 +248,222 @@
       if (! quad_fcn)
 	QUAD_ABORT ();
 
-      double a = args(1).double_value ();
-
-      if (error_state)
-	QUAD_ABORT1 ("expecting second argument to be a scalar");
-
-      double b = args(2).double_value ();
-
-      if (error_state)
-	QUAD_ABORT1 ("expecting third argument to be a scalar");
-
-      int indefinite = 0;
-      IndefQuad::IntegralType indef_type = IndefQuad::doubly_infinite;
-      double bound = 0.0;
-      if (xisinf (a) && xisinf (b))
-	{
-	  indefinite = 1;
-	  indef_type = IndefQuad::doubly_infinite;
-	}
-      else if (xisinf (a))
+      if (args(1).is_single_type () || args(2).is_single_type ())
 	{
-	  indefinite = 1;
-	  bound = b;
-	  indef_type = IndefQuad::neg_inf_to_bound;
-	}
-      else if (xisinf (b))
-	{
-	  indefinite = 1;
-	  bound = a;
-	  indef_type = IndefQuad::bound_to_inf;
-	}
+	  float a = args(1).float_value ();
 
-      octave_idx_type ier = 0;
-      octave_idx_type nfun = 0;
-      double abserr = 0.0;
-      double val = 0.0;
-      bool have_sing = false;
-      ColumnVector sing;
-      ColumnVector tol;
+	  if (error_state)
+	    QUAD_ABORT1 ("expecting second argument to be a scalar");
 
-      switch (nargin)
-	{
-	case 5:
-	  if (indefinite)
-	    QUAD_ABORT1 ("singularities not allowed on infinite intervals");
-
-	  have_sing = true;
-
-	  sing = ColumnVector (args(4).vector_value ());
+	  float b = args(2).float_value ();
 
 	  if (error_state)
-	    QUAD_ABORT1 ("expecting vector of singularities as fourth argument");
+	    QUAD_ABORT1 ("expecting third argument to be a scalar");
 
-	case 4:
-	  tol = ColumnVector (args(3).vector_value ());
+	  int indefinite = 0;
+	  FloatIndefQuad::IntegralType indef_type = FloatIndefQuad::doubly_infinite;
+	  float bound = 0.0;
+	  if (xisinf (a) && xisinf (b))
+	    {
+	      indefinite = 1;
+	      indef_type = FloatIndefQuad::doubly_infinite;
+	    }
+	  else if (xisinf (a))
+	    {
+	      indefinite = 1;
+	      bound = b;
+	      indef_type = FloatIndefQuad::neg_inf_to_bound;
+	    }
+	  else if (xisinf (b))
+	    {
+	      indefinite = 1;
+	      bound = a;
+	      indef_type = FloatIndefQuad::bound_to_inf;
+	    }
 
-	  if (error_state)
-	    QUAD_ABORT1 ("expecting vector of tolerances as fifth argument");
+	  octave_idx_type ier = 0;
+	  octave_idx_type nfun = 0;
+	  float abserr = 0.0;
+	  float val = 0.0;
+	  bool have_sing = false;
+	  FloatColumnVector sing;
+	  FloatColumnVector tol;
+
+	  switch (nargin)
+	    {
+	    case 5:
+	      if (indefinite)
+		QUAD_ABORT1 ("singularities not allowed on infinite intervals");
+
+	      have_sing = true;
+
+	      sing = FloatColumnVector (args(4).float_vector_value ());
 
-	  switch (tol.capacity ())
-	    {
-	    case 2:
-	      quad_opts.set_relative_tolerance (tol (1));
+	      if (error_state)
+		QUAD_ABORT1 ("expecting vector of singularities as fourth argument");
+
+	    case 4:
+	      tol = FloatColumnVector (args(3).float_vector_value ());
+
+	      if (error_state)
+		QUAD_ABORT1 ("expecting vector of tolerances as fifth argument");
+
+	      switch (tol.capacity ())
+		{
+		case 2:
+		  quad_opts.set_single_precision_relative_tolerance (tol (1));
+
+		case 1:
+		  quad_opts.set_single_precision_absolute_tolerance (tol (0));
+		  break;
+
+		default:
+		  QUAD_ABORT1 ("expecting tol to contain no more than two values");
+		}
 
-	    case 1:
-	      quad_opts.set_absolute_tolerance (tol (0));
+	    case 3:
+	      if (indefinite)
+		{
+		  FloatIndefQuad iq (quad_float_user_function, bound, 
+				     indef_type);
+		  iq.set_options (quad_opts);
+		  val = iq.float_integrate (ier, nfun, abserr);
+		}
+	      else
+		{
+		  if (have_sing)
+		    {
+		      FloatDefQuad dq (quad_float_user_function, a, b, sing);
+		      dq.set_options (quad_opts);
+		      val = dq.float_integrate (ier, nfun, abserr);
+		    }
+		  else
+		    {
+		      FloatDefQuad dq (quad_float_user_function, a, b);
+		      dq.set_options (quad_opts);
+		      val = dq.float_integrate (ier, nfun, abserr);
+		    }
+		}
 	      break;
 
 	    default:
-	      QUAD_ABORT1 ("expecting tol to contain no more than two values");
+	      panic_impossible ();
+	      break;
+	    }
+
+	  retval(3) = abserr;
+	  retval(2) = static_cast<double> (nfun);
+	  retval(1) = static_cast<double> (ier);
+	  retval(0) = val;
+
+	}
+      else
+	{
+	  double a = args(1).double_value ();
+
+	  if (error_state)
+	    QUAD_ABORT1 ("expecting second argument to be a scalar");
+
+	  double b = args(2).double_value ();
+
+	  if (error_state)
+	    QUAD_ABORT1 ("expecting third argument to be a scalar");
+
+	  int indefinite = 0;
+	  IndefQuad::IntegralType indef_type = IndefQuad::doubly_infinite;
+	  double bound = 0.0;
+	  if (xisinf (a) && xisinf (b))
+	    {
+	      indefinite = 1;
+	      indef_type = IndefQuad::doubly_infinite;
+	    }
+	  else if (xisinf (a))
+	    {
+	      indefinite = 1;
+	      bound = b;
+	      indef_type = IndefQuad::neg_inf_to_bound;
+	    }
+	  else if (xisinf (b))
+	    {
+	      indefinite = 1;
+	      bound = a;
+	      indef_type = IndefQuad::bound_to_inf;
 	    }
 
-	case 3:
-	  if (indefinite)
+	  octave_idx_type ier = 0;
+	  octave_idx_type nfun = 0;
+	  double abserr = 0.0;
+	  double val = 0.0;
+	  bool have_sing = false;
+	  ColumnVector sing;
+	  ColumnVector tol;
+
+	  switch (nargin)
 	    {
-	      IndefQuad iq (quad_user_function, bound, indef_type);
-	      iq.set_options (quad_opts);
-	      val = iq.integrate (ier, nfun, abserr);
-	    }
-	  else
-	    {
-	      if (have_sing)
+	    case 5:
+	      if (indefinite)
+		QUAD_ABORT1 ("singularities not allowed on infinite intervals");
+
+	      have_sing = true;
+
+	      sing = ColumnVector (args(4).vector_value ());
+
+	      if (error_state)
+		QUAD_ABORT1 ("expecting vector of singularities as fourth argument");
+
+	    case 4:
+	      tol = ColumnVector (args(3).vector_value ());
+
+	      if (error_state)
+		QUAD_ABORT1 ("expecting vector of tolerances as fifth argument");
+
+	      switch (tol.capacity ())
 		{
-		  DefQuad dq (quad_user_function, a, b, sing);
-		  dq.set_options (quad_opts);
-		  val = dq.integrate (ier, nfun, abserr);
+		case 2:
+		  quad_opts.set_relative_tolerance (tol (1));
+
+		case 1:
+		  quad_opts.set_absolute_tolerance (tol (0));
+		  break;
+
+		default:
+		  QUAD_ABORT1 ("expecting tol to contain no more than two values");
+		}
+
+	    case 3:
+	      if (indefinite)
+		{
+		  IndefQuad iq (quad_user_function, bound, indef_type);
+		  iq.set_options (quad_opts);
+		  val = iq.integrate (ier, nfun, abserr);
 		}
 	      else
 		{
-		  DefQuad dq (quad_user_function, a, b);
-		  dq.set_options (quad_opts);
-		  val = dq.integrate (ier, nfun, abserr);
+		  if (have_sing)
+		    {
+		      DefQuad dq (quad_user_function, a, b, sing);
+		      dq.set_options (quad_opts);
+		      val = dq.integrate (ier, nfun, abserr);
+		    }
+		  else
+		    {
+		      DefQuad dq (quad_user_function, a, b);
+		      dq.set_options (quad_opts);
+		      val = dq.integrate (ier, nfun, abserr);
+		    }
 		}
-	    }
-	  break;
+	      break;
 
-	default:
-	  panic_impossible ();
-	  break;
-	}
+	    default:
+	      panic_impossible ();
+	      break;
+	    }
 
-      retval(3) = abserr;
-      retval(2) = static_cast<double> (nfun);
-      retval(1) = static_cast<double> (ier);
-      retval(0) = val;
+	  retval(3) = abserr;
+	  retval(2) = static_cast<double> (nfun);
+	  retval(1) = static_cast<double> (ier);
+	  retval(0) = val;
+	}
 
       if (fcn_name.length())
 	clear_function (fcn_name);
@@ -327,12 +484,19 @@
 %! [v, ier, nfun, err] = quad ("f", 0, 5);
 %! assert(ier == 0 && abs (v - 17.5) < sqrt (eps) && nfun > 0 && 
 %!        err < sqrt (eps))
+%!test
+%! [v, ier, nfun, err] = quad ("f", single(0), single(5));
+%! assert(ier == 0 && abs (v - 17.5) < sqrt (eps ("single")) && nfun > 0 && 
+%!        err < sqrt (eps ("single")))
 
 %!function y = f (x)
 %!  y = x .* sin (1 ./ x) .* sqrt (abs (1 - x));
 %!test
 %!  [v, ier, nfun, err] = quad ("f", 0.001, 3);
 %! assert((ier == 0 || ier == 1) && abs (v - 1.98194120273598) < sqrt (eps) && nfun > 0);
+%!test
+%!  [v, ier, nfun, err] = quad ("f", single(0.001), single(3));
+%! assert((ier == 0 || ier == 1) && abs (v - 1.98194120273598) < sqrt (eps ("single")) && nfun > 0);
 
 %!error <Invalid call to quad.*> quad ();