--- a/src/DLD-FUNCTIONS/quad.cc
+++ b/src/DLD-FUNCTIONS/quad.cc
@@ -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 ();