--- liboctave/ChangeLog.cvs	Tue Sep  9 11:03:22 2003
+++ liboctave/ChangeLog	Tue Sep  9 11:02:06 2003
@@ -1,3 +1,8 @@
+2003-09-09  David Bateman <address@hidden>
+
+	* lo-specfun.cc (zbesj, zbesy, zbesi, zbesk, zbesh1, zbesh2, 
+	airy, biry): Used scaled versions for overflow/underflow
+
 2003-09-04  John W. Eaton  <address@hidden>
 
 	* lo-specfun.cc (xlgamma): Require nonnegative argument.
--- liboctave/lo-specfun.cc.cvs	Tue Sep  9 10:57:37 2003
+++ liboctave/lo-specfun.cc	Tue Sep  9 10:46:28 2003
@@ -225,6 +225,16 @@
 
       F77_FUNC (zbesj, ZBESJ) (zr, zi, alpha, kode, 1, &yr, &yi, nz, ierr);
 
+      if (kode == 1 && ierr)
+	{
+	  // Overflow/Underflow. Try scaling
+	  F77_FUNC (zbesj, ZBESJ) (zr, zi, alpha, 2, 1, &yr, &yi, nz, ierr);
+
+	  double expz = exp(std::abs(zi)); 
+	  yr *= expz;
+	  yi *= expz;
+	}
+
       if (zi == 0.0 && zr >= 0.0)
 	yi = 0.0;
 
@@ -278,6 +288,17 @@
 	  F77_FUNC (zbesy, ZBESY) (zr, zi, alpha, kode, 1, &yr, &yi, nz,
 				  &wr, &wi, ierr);
 
+	  if (kode == 1 && ierr)
+	    {
+	      // Overflow/Underflow. Try scaling
+	      F77_FUNC (zbesy, ZBESY) (zr, zi, alpha, 2, 1, &yr, &yi, nz,
+				   &wr, &wi, ierr);
+	      
+	      double expz = exp(std::abs(zi));
+	      yr *= expz;
+	      yi *= expz;
+	    }
+
 	  if (zi == 0.0 && zr >= 0.0)
 	    yi = 0.0;
 	}
@@ -320,6 +341,16 @@
 
       F77_FUNC (zbesi, ZBESI) (zr, zi, alpha, kode, 1, &yr, &yi, nz, ierr);
 
+      if (kode == 1 && ierr)
+	{
+	  // Overflow/Underflow. Try scaling
+	  F77_FUNC (zbesi, ZBESI) (zr, zi, alpha, 2, 1, &yr, &yi, nz, ierr);
+	  
+	  double expz = exp(std::abs(zr));
+	  yr *= expz;
+	  yi *= expz;
+	}
+
       if (zi == 0.0 && zr >= 0.0)
 	yi = 0.0;
 
@@ -371,6 +402,19 @@
 	{
 	  F77_FUNC (zbesk, ZBESK) (zr, zi, alpha, kode, 1, &yr, &yi, nz, ierr);
 
+	  if (kode == 1 && ierr)
+	    {
+	      // Overflow/Underflow. Try scaling
+	      F77_FUNC (zbesk, ZBESK) (zr, zi, alpha, 2, 1, &yr, &yi, 
+				       nz, ierr);
+	      Complex expz = exp(-z);
+	      double rexpz = real(expz);
+	      double iexpz = imag(expz);
+	      double tmp = yr*rexpz - yi*iexpz;
+	      yi = yr*iexpz + yi*rexpz;
+	      yr = tmp;
+	    }
+
 	  if (zi == 0.0 && zr >= 0.0)
 	    yi = 0.0;
 	}
@@ -404,6 +448,19 @@
 
       F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, kode, 1, 1, &yr, &yi, nz, ierr);
 
+      if (kode == 1 && ierr)
+	{
+	  // Overflow/Underflow. Try scaling
+	  F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, 2, 1, 1, &yr, &yi, nz, ierr);
+
+	  Complex expz = exp( Complex ( 0., 1.) * z);
+	  double rexpz = real(expz);
+	  double iexpz = imag(expz);
+	  double tmp = yr*rexpz - yi*iexpz;
+	  yi = yr*iexpz + yi*rexpz;
+	  yr = tmp;
+	}
+
       retval = bessel_return_value (Complex (yr, yi), ierr);
     }
   else
@@ -437,6 +494,20 @@
 
       F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, kode, 2, 1, &yr, &yi, nz, ierr);
 
+      if (kode == 1 && ierr)
+	{
+	  // Overflow/Underflow. Try scaling
+	  F77_FUNC (zbesh, ZBESH) (zr, zi, alpha, 2, 2, 1, &yr, &yi, nz, 
+				   ierr);
+
+	  Complex expz = exp( - Complex ( 0., 1.) * z);
+	  double rexpz = real(expz);
+	  double iexpz = imag(expz);
+	  double tmp = yr*rexpz - yi*iexpz;
+	  yi = yr*iexpz + yi*rexpz;
+	  yr = tmp;
+	}
+
       retval = bessel_return_value (Complex (yr, yi), ierr);
     }
   else
@@ -622,6 +693,19 @@
 
   F77_FUNC (zairy, ZAIRY) (zr, zi, id, kode, ar, ai, nz, ierr);
 
+  if (kode == 1 && ierr)
+    {
+      // Overflow/Underflow. Try scaling
+      F77_FUNC (zairy, ZAIRY) (zr, zi, id, 2, ar, ai, nz, ierr);
+
+      Complex expz = exp( - 2. / 3. * z * sqrt(z));
+      double rexpz = real(expz);
+      double iexpz = imag(expz);
+      double tmp = ar*rexpz - ai*iexpz;
+      ai = ar*iexpz + ai*rexpz;
+      ar = tmp;
+    }
+
   if (zi == 0.0 && (! scaled || zr >= 0.0))
     ai = 0.0;
 
@@ -642,6 +726,19 @@
   int kode = scaled ? 2 : 1;
 
   F77_FUNC (zbiry, ZBIRY) (zr, zi, id, kode, ar, ai, ierr);
+
+  if (kode == 1 && ierr)
+    {
+      // Overflow/Underflow. Try scaling
+      F77_FUNC (zbiry, ZBIRY) (zr, zi, id, 2, ar, ai, ierr);
+
+      Complex expz = exp( std::abs( real( 2. / 3. * z * sqrt(z))));
+      double rexpz = real(expz);
+      double iexpz = imag(expz);
+      double tmp = ar*rexpz - ai*iexpz;
+      ai = ar*iexpz + ai*rexpz;
+      ar = tmp;
+    }
 
   if (zi == 0.0 && (! scaled || zr >= 0.0))
     ai = 0.0;