# HG changeset patch
# User Craig Hudson <c_hudson_phd@hotmail.com>
# Date 1377866707 -3600
# Node ID 44a293eafe37b1635298f8a98e620e26acf9e1fd
# Parent  975aff6fafe5adfe8ef4e4556f6287f50ec095ea
Improve complex mapper functions

diff -r 975aff6fafe5 -r 44a293eafe37 libgui/src/resource-manager.cc
--- a/libgui/src/resource-manager.cc	Tue Jul 02 14:17:33 2013 +0100
+++ b/libgui/src/resource-manager.cc	Fri Aug 30 13:45:07 2013 +0100
@@ -1546,6 +1546,7 @@
       "sqp "
       "sqrt "
       "sqrtm "
+      "square "
       "squeeze "
       "sscanf "
       "stairs "
diff -r 975aff6fafe5 -r 44a293eafe37 libinterp/corefcn/mappers.cc
--- a/libinterp/corefcn/mappers.cc	Tue Jul 02 14:17:33 2013 +0100
+++ b/libinterp/corefcn/mappers.cc	Fri Aug 30 13:45:07 2013 +0100
@@ -2017,6 +2017,25 @@
 %!error sqrt (1, 2)
 */
 
+DEFUN (square, args, ,
+    "-*- texinfo -*-\n\
+@deftypefn {Mapping Function} {} square (@var{x})\n\
+Compute the square of each element of @var{x}.  If @var{x} is complex,\n\
+overflow is avoided unless necessary, using the algorithm described by\n\
+W. Kahan in Branch Cuts for Complex Elementary Functions, or Much Ado\n\
+About Nothing's Sign Bit (from The State of the Art in Numerical Analysis,\n\
+eds. Iserles and Powell, Clarendon Press, Oxford, 1987)\n\
+@end deftypefn")
+{
+  octave_value retval;
+  if (args.length () == 1)
+    retval = args(0).square ();
+  else
+    print_usage ();
+
+  return retval;
+}
+
 DEFUN (tan, args, ,
     "-*- texinfo -*-\n\
 @deftypefn {Mapping Function} {} tan (@var{z})\n\
diff -r 975aff6fafe5 -r 44a293eafe37 libinterp/octave-value/ov-base.cc
--- a/libinterp/octave-value/ov-base.cc	Tue Jul 02 14:17:33 2013 +0100
+++ b/libinterp/octave-value/ov-base.cc	Fri Aug 30 13:45:07 2013 +0100
@@ -1221,6 +1221,7 @@
       "sin",
       "sinh",
       "sqrt",
+      "square",
       "tan",
       "tanh",
       "isalnum",
diff -r 975aff6fafe5 -r 44a293eafe37 libinterp/octave-value/ov-base.h
--- a/libinterp/octave-value/ov-base.h	Tue Jul 02 14:17:33 2013 +0100
+++ b/libinterp/octave-value/ov-base.h	Fri Aug 30 13:45:07 2013 +0100
@@ -721,6 +721,7 @@
       umap_sin,
       umap_sinh,
       umap_sqrt,
+      umap_square,
       umap_tan,
       umap_tanh,
       umap_xisalnum,
diff -r 975aff6fafe5 -r 44a293eafe37 libinterp/octave-value/ov-complex.cc
--- a/libinterp/octave-value/ov-complex.cc	Tue Jul 02 14:17:33 2013 +0100
+++ b/libinterp/octave-value/ov-complex.cc	Fri Aug 30 13:45:07 2013 +0100
@@ -75,10 +75,16 @@
 {
   octave_base_value *retval = 0;
 
+  double re = std::real (scalar);
   double im = std::imag (scalar);
 
-  if (im == 0.0)
-    retval = new octave_scalar (std::real (scalar));
+  // Do not cast complex expression to real if imaginary part is negative zero.
+  // Also, do not narrow if real and imaginary parts are both zero, as it
+  // prevents expressions such as 2 - 0*i from being evaluated as such, by
+  // casting 0*i to zero.
+  
+  if (im == 0.0 && signbit(im) == 0 && re != 0.0)
+    retval = new octave_scalar (re);
 
   return retval;
 }
@@ -422,11 +428,11 @@
     case umap_ ## UMAP: \
       return octave_value (FCN (scalar))
 
-      SCALAR_MAPPER (abs, std::abs);
+      SCALAR_MAPPER (abs, ::abs);
       SCALAR_MAPPER (acos, ::acos);
       SCALAR_MAPPER (acosh, ::acosh);
-      SCALAR_MAPPER (angle, std::arg);
-      SCALAR_MAPPER (arg, std::arg);
+      SCALAR_MAPPER (angle, ::arg);
+      SCALAR_MAPPER (arg, ::arg);
       SCALAR_MAPPER (asin, ::asin);
       SCALAR_MAPPER (asinh, ::asinh);
       SCALAR_MAPPER (atan, ::atan);
@@ -445,7 +451,7 @@
       SCALAR_MAPPER (fix, ::fix);
       SCALAR_MAPPER (floor, ::floor);
       SCALAR_MAPPER (imag, std::imag);
-      SCALAR_MAPPER (log, std::log);
+      SCALAR_MAPPER (log, ::log);
       SCALAR_MAPPER (log2, xlog2);
       SCALAR_MAPPER (log10, std::log10);
       SCALAR_MAPPER (log1p, ::log1p);
@@ -455,9 +461,10 @@
       SCALAR_MAPPER (signum, ::signum);
       SCALAR_MAPPER (sin, std::sin);
       SCALAR_MAPPER (sinh, std::sinh);
-      SCALAR_MAPPER (sqrt, std::sqrt);
-      SCALAR_MAPPER (tan, std::tan);
-      SCALAR_MAPPER (tanh, std::tanh);
+      SCALAR_MAPPER (sqrt, ::sqrt);
+      SCALAR_MAPPER (square, ::square);
+      SCALAR_MAPPER (tan, ::tan);
+      SCALAR_MAPPER (tanh, ::tanh);
       SCALAR_MAPPER (finite, xfinite);
       SCALAR_MAPPER (isinf, xisinf);
       SCALAR_MAPPER (isna, octave_is_NA);
diff -r 975aff6fafe5 -r 44a293eafe37 libinterp/octave-value/ov-cx-mat.cc
--- a/libinterp/octave-value/ov-cx-mat.cc	Tue Jul 02 14:17:33 2013 +0100
+++ b/libinterp/octave-value/ov-cx-mat.cc	Fri Aug 30 13:45:07 2013 +0100
@@ -86,8 +86,12 @@
   if (matrix.numel () == 1)
     {
       Complex c = matrix (0);
+      double im = std::imag (c);
 
-      if (std::imag (c) == 0.0)
+      // Do not cast complex expression to real if imaginary part is
+      // negative zero.
+  
+      if (im == 0.0 && signbit(im) == 0)
         retval = new octave_scalar (std::real (c));
       else
         retval = new octave_complex (c);
@@ -760,7 +764,6 @@
     case umap_ ## UMAP: \
       return octave_value (matrix.FCN ())
 
-      ARRAY_METHOD_MAPPER (abs, abs);
       ARRAY_METHOD_MAPPER (isnan, isnan);
       ARRAY_METHOD_MAPPER (isinf, isinf);
       ARRAY_METHOD_MAPPER (finite, isfinite);
@@ -769,10 +772,11 @@
     case umap_ ## UMAP: \
       return octave_value (matrix.map<TYPE> (FCN))
 
+      ARRAY_MAPPER (abs, double, ::abs);
       ARRAY_MAPPER (acos, Complex, ::acos);
       ARRAY_MAPPER (acosh, Complex, ::acosh);
-      ARRAY_MAPPER (angle, double, std::arg);
-      ARRAY_MAPPER (arg, double, std::arg);
+      ARRAY_MAPPER (angle, double, ::arg);
+      ARRAY_MAPPER (arg, double, ::arg);
       ARRAY_MAPPER (asin, Complex, ::asin);
       ARRAY_MAPPER (asinh, Complex, ::asinh);
       ARRAY_MAPPER (atan, Complex, ::atan);
@@ -789,7 +793,7 @@
       ARRAY_MAPPER (expm1, Complex, ::expm1);
       ARRAY_MAPPER (fix, Complex, ::fix);
       ARRAY_MAPPER (floor, Complex, ::floor);
-      ARRAY_MAPPER (log, Complex, std::log);
+      ARRAY_MAPPER (log, Complex, ::log);
       ARRAY_MAPPER (log2, Complex, xlog2);
       ARRAY_MAPPER (log10, Complex, std::log10);
       ARRAY_MAPPER (log1p, Complex, ::log1p);
@@ -798,9 +802,10 @@
       ARRAY_MAPPER (signum, Complex, ::signum);
       ARRAY_MAPPER (sin, Complex, std::sin);
       ARRAY_MAPPER (sinh, Complex, std::sinh);
-      ARRAY_MAPPER (sqrt, Complex, std::sqrt);
-      ARRAY_MAPPER (tan, Complex, std::tan);
-      ARRAY_MAPPER (tanh, Complex, std::tanh);
+      ARRAY_MAPPER (sqrt, Complex, ::sqrt);
+      ARRAY_MAPPER (square, Complex, ::square);
+      ARRAY_MAPPER (tan, Complex, ::tan);
+      ARRAY_MAPPER (tanh, Complex, ::tanh);
       ARRAY_MAPPER (isna, bool, octave_is_NA);
 
     default:
diff -r 975aff6fafe5 -r 44a293eafe37 libinterp/octave-value/ov-cx-sparse.cc
--- a/libinterp/octave-value/ov-cx-sparse.cc	Tue Jul 02 14:17:33 2013 +0100
+++ b/libinterp/octave-value/ov-cx-sparse.cc	Fri Aug 30 13:45:07 2013 +0100
@@ -72,8 +72,12 @@
           const SparseComplexMatrix tmp (matrix);
 
           Complex c = tmp (0, 0);
+          double im = std::imag (c);
 
-          if (std::imag (c) == 0.0)
+          // Do not cast complex expression to real if imaginary part is
+          // negative zero.
+  
+          if (im == 0.0 && signbit(im) == 0)
             retval = new octave_scalar (std::real (c));
           else
             retval = new octave_complex (c);
@@ -885,20 +889,15 @@
     case umap_imag:
       return ::imag (matrix);
 
-#define ARRAY_METHOD_MAPPER(UMAP, FCN) \
-    case umap_ ## UMAP: \
-      return octave_value (matrix.FCN ())
-
-      ARRAY_METHOD_MAPPER (abs, abs);
-
 #define ARRAY_MAPPER(UMAP, TYPE, FCN) \
     case umap_ ## UMAP: \
       return octave_value (matrix.map<TYPE> (FCN))
 
+      ARRAY_MAPPER (abs, double, ::abs);
       ARRAY_MAPPER (acos, Complex, ::acos);
       ARRAY_MAPPER (acosh, Complex, ::acosh);
-      ARRAY_MAPPER (angle, double, std::arg);
-      ARRAY_MAPPER (arg, double, std::arg);
+      ARRAY_MAPPER (angle, double, ::arg);
+      ARRAY_MAPPER (arg, double, ::arg);
       ARRAY_MAPPER (asin, Complex, ::asin);
       ARRAY_MAPPER (asinh, Complex, ::asinh);
       ARRAY_MAPPER (atan, Complex, ::atan);
@@ -916,7 +915,7 @@
       ARRAY_MAPPER (expm1, Complex, ::expm1);
       ARRAY_MAPPER (fix, Complex, ::fix);
       ARRAY_MAPPER (floor, Complex, ::floor);
-      ARRAY_MAPPER (log, Complex, std::log);
+      ARRAY_MAPPER (log, Complex, ::log);
       ARRAY_MAPPER (log2, Complex, xlog2);
       ARRAY_MAPPER (log10, Complex, std::log10);
       ARRAY_MAPPER (log1p, Complex, ::log1p);
@@ -925,9 +924,10 @@
       ARRAY_MAPPER (signum, Complex, ::signum);
       ARRAY_MAPPER (sin, Complex, std::sin);
       ARRAY_MAPPER (sinh, Complex, std::sinh);
-      ARRAY_MAPPER (sqrt, Complex, std::sqrt);
-      ARRAY_MAPPER (tan, Complex, std::tan);
-      ARRAY_MAPPER (tanh, Complex, std::tanh);
+      ARRAY_MAPPER (sqrt, Complex, ::sqrt);
+      ARRAY_MAPPER (square, Complex, ::square);
+      ARRAY_MAPPER (tan, Complex, ::tan);
+      ARRAY_MAPPER (tanh, Complex, ::tanh);
       ARRAY_MAPPER (isnan, bool, xisnan);
       ARRAY_MAPPER (isna, bool, octave_is_NA);
       ARRAY_MAPPER (isinf, bool, xisinf);
diff -r 975aff6fafe5 -r 44a293eafe37 libinterp/octave-value/ov-float.cc
--- a/libinterp/octave-value/ov-float.cc	Tue Jul 02 14:17:33 2013 +0100
+++ b/libinterp/octave-value/ov-float.cc	Fri Aug 30 13:45:07 2013 +0100
@@ -315,8 +315,8 @@
       SCALAR_MAPPER (roundb, xroundb);
       SCALAR_MAPPER (signum, ::signum);
       SCALAR_MAPPER (sin, ::sinf);
-      SCALAR_MAPPER (sinh, ::sinhf);
       SCALAR_MAPPER (sqrt, rc_sqrt);
+      SCALAR_MAPPER (square, ::square);
       SCALAR_MAPPER (tan, ::tanf);
       SCALAR_MAPPER (tanh, ::tanhf);
       SCALAR_MAPPER (finite, xfinite);
diff -r 975aff6fafe5 -r 44a293eafe37 libinterp/octave-value/ov-flt-complex.cc
--- a/libinterp/octave-value/ov-flt-complex.cc	Tue Jul 02 14:17:33 2013 +0100
+++ b/libinterp/octave-value/ov-flt-complex.cc	Fri Aug 30 13:45:07 2013 +0100
@@ -60,10 +60,16 @@
 {
   octave_base_value *retval = 0;
 
+  float re = std::real (scalar);
   float im = std::imag (scalar);
 
-  if (im == 0.0)
-    retval = new octave_float_scalar (std::real (scalar));
+  // Do not cast complex expression to real if imaginary part is negative zero.
+  // Also, do not narrow if real and imaginary parts are both zero, as it
+  // prevents expressions such as 2 - 0*i from being evaluated as such, by
+  // casting 0*i to zero.
+  
+  if (im == 0.0 && signbit(im) == 0 && re != 0.0)
+    retval = new octave_float_scalar (re);
 
   return retval;
 }
@@ -407,11 +413,11 @@
     case umap_ ## UMAP: \
       return octave_value (FCN (scalar))
 
-      SCALAR_MAPPER (abs, std::abs);
+      SCALAR_MAPPER (abs, ::abs);
       SCALAR_MAPPER (acos, ::acos);
       SCALAR_MAPPER (acosh, ::acosh);
-      SCALAR_MAPPER (angle, std::arg);
-      SCALAR_MAPPER (arg, std::arg);
+      SCALAR_MAPPER (angle, ::arg);
+      SCALAR_MAPPER (arg, ::arg);
       SCALAR_MAPPER (asin, ::asin);
       SCALAR_MAPPER (asinh, ::asinh);
       SCALAR_MAPPER (atan, ::atan);
@@ -430,7 +436,7 @@
       SCALAR_MAPPER (fix, ::fix);
       SCALAR_MAPPER (floor, ::floor);
       SCALAR_MAPPER (imag, std::imag);
-      SCALAR_MAPPER (log, std::log);
+      SCALAR_MAPPER (log, ::log);
       SCALAR_MAPPER (log2, xlog2);
       SCALAR_MAPPER (log10, std::log10);
       SCALAR_MAPPER (log1p, ::log1p);
@@ -440,9 +446,10 @@
       SCALAR_MAPPER (signum, ::signum);
       SCALAR_MAPPER (sin, std::sin);
       SCALAR_MAPPER (sinh, std::sinh);
-      SCALAR_MAPPER (sqrt, std::sqrt);
-      SCALAR_MAPPER (tan, std::tan);
-      SCALAR_MAPPER (tanh, std::tanh);
+      SCALAR_MAPPER (sqrt, ::sqrt);
+      SCALAR_MAPPER (square, ::square);
+      SCALAR_MAPPER (tan, ::tan);
+      SCALAR_MAPPER (tanh, ::tanh);
       SCALAR_MAPPER (finite, xfinite);
       SCALAR_MAPPER (isinf, xisinf);
       SCALAR_MAPPER (isna, octave_is_NA);
diff -r 975aff6fafe5 -r 44a293eafe37 libinterp/octave-value/ov-flt-cx-mat.cc
--- a/libinterp/octave-value/ov-flt-cx-mat.cc	Tue Jul 02 14:17:33 2013 +0100
+++ b/libinterp/octave-value/ov-flt-cx-mat.cc	Fri Aug 30 13:45:07 2013 +0100
@@ -75,8 +75,12 @@
   if (matrix.numel () == 1)
     {
       FloatComplex c = matrix (0);
+      float im = std::imag (c);
 
-      if (std::imag (c) == 0.0)
+      // Do not cast complex expression to real if imaginary part is
+      // negative zero.
+  
+      if (im == 0.0 && signbit(im) == 0)
         retval = new octave_float_scalar (std::real (c));
       else
         retval = new octave_float_complex (c);
@@ -727,7 +731,6 @@
     case umap_ ## UMAP: \
       return octave_value (matrix.FCN ())
 
-      ARRAY_METHOD_MAPPER (abs, abs);
       ARRAY_METHOD_MAPPER (isnan, isnan);
       ARRAY_METHOD_MAPPER (isinf, isinf);
       ARRAY_METHOD_MAPPER (finite, isfinite);
@@ -736,10 +739,11 @@
     case umap_ ## UMAP: \
       return octave_value (matrix.map<TYPE> (FCN))
 
+      ARRAY_MAPPER (abs, float, ::abs);
       ARRAY_MAPPER (acos, FloatComplex, ::acos);
       ARRAY_MAPPER (acosh, FloatComplex, ::acosh);
-      ARRAY_MAPPER (angle, float, std::arg);
-      ARRAY_MAPPER (arg, float, std::arg);
+      ARRAY_MAPPER (angle, float, ::arg);
+      ARRAY_MAPPER (arg, float, ::arg);
       ARRAY_MAPPER (asin, FloatComplex, ::asin);
       ARRAY_MAPPER (asinh, FloatComplex, ::asinh);
       ARRAY_MAPPER (atan, FloatComplex, ::atan);
@@ -756,7 +760,7 @@
       ARRAY_MAPPER (expm1, FloatComplex, ::expm1);
       ARRAY_MAPPER (fix, FloatComplex, ::fix);
       ARRAY_MAPPER (floor, FloatComplex, ::floor);
-      ARRAY_MAPPER (log, FloatComplex, std::log);
+      ARRAY_MAPPER (log, FloatComplex, ::log);
       ARRAY_MAPPER (log2, FloatComplex, xlog2);
       ARRAY_MAPPER (log10, FloatComplex, std::log10);
       ARRAY_MAPPER (log1p, FloatComplex, ::log1p);
@@ -765,9 +769,10 @@
       ARRAY_MAPPER (signum, FloatComplex, ::signum);
       ARRAY_MAPPER (sin, FloatComplex, std::sin);
       ARRAY_MAPPER (sinh, FloatComplex, std::sinh);
-      ARRAY_MAPPER (sqrt, FloatComplex, std::sqrt);
-      ARRAY_MAPPER (tan, FloatComplex, std::tan);
-      ARRAY_MAPPER (tanh, FloatComplex, std::tanh);
+      ARRAY_MAPPER (sqrt, FloatComplex, ::sqrt);
+      ARRAY_MAPPER (square, FloatComplex, ::square);
+      ARRAY_MAPPER (tan, FloatComplex, ::tan);
+      ARRAY_MAPPER (tanh, FloatComplex, ::tanh);
       ARRAY_MAPPER (isna, bool, octave_is_NA);
 
     default:
diff -r 975aff6fafe5 -r 44a293eafe37 libinterp/octave-value/ov-flt-re-mat.cc
--- a/libinterp/octave-value/ov-flt-re-mat.cc	Tue Jul 02 14:17:33 2013 +0100
+++ b/libinterp/octave-value/ov-flt-re-mat.cc	Fri Aug 30 13:45:07 2013 +0100
@@ -803,8 +803,8 @@
       ARRAY_MAPPER (roundb, float, xroundb);
       ARRAY_MAPPER (signum, float, ::signum);
       ARRAY_MAPPER (sin, float, ::sinf);
-      ARRAY_MAPPER (sinh, float, ::sinhf);
       RC_ARRAY_MAPPER (sqrt, FloatComplex, rc_sqrt);
+      ARRAY_MAPPER (square, float, ::square);
       ARRAY_MAPPER (tan, float, ::tanf);
       ARRAY_MAPPER (tanh, float, ::tanhf);
       ARRAY_MAPPER (isna, bool, octave_is_NA);
diff -r 975aff6fafe5 -r 44a293eafe37 libinterp/octave-value/ov-re-mat.cc
--- a/libinterp/octave-value/ov-re-mat.cc	Tue Jul 02 14:17:33 2013 +0100
+++ b/libinterp/octave-value/ov-re-mat.cc	Fri Aug 30 13:45:07 2013 +0100
@@ -928,6 +928,7 @@
       ARRAY_MAPPER (sin, double, ::sin);
       ARRAY_MAPPER (sinh, double, ::sinh);
       RC_ARRAY_MAPPER (sqrt, Complex, rc_sqrt);
+      ARRAY_MAPPER (square, double, ::square);
       ARRAY_MAPPER (tan, double, ::tan);
       ARRAY_MAPPER (tanh, double, ::tanh);
       ARRAY_MAPPER (isna, bool, octave_is_NA);
diff -r 975aff6fafe5 -r 44a293eafe37 libinterp/octave-value/ov-re-sparse.cc
--- a/libinterp/octave-value/ov-re-sparse.cc	Tue Jul 02 14:17:33 2013 +0100
+++ b/libinterp/octave-value/ov-re-sparse.cc	Fri Aug 30 13:45:07 2013 +0100
@@ -935,6 +935,7 @@
       ARRAY_MAPPER (sin, double, ::sin);
       ARRAY_MAPPER (sinh, double, ::sinh);
       ARRAY_MAPPER (sqrt, Complex, rc_sqrt);
+      ARRAY_MAPPER (square, double, ::square);
       ARRAY_MAPPER (tan, double, ::tan);
       ARRAY_MAPPER (tanh, double, ::tanh);
       ARRAY_MAPPER (isnan, bool, xisnan);
diff -r 975aff6fafe5 -r 44a293eafe37 libinterp/octave-value/ov-scalar.cc
--- a/libinterp/octave-value/ov-scalar.cc	Tue Jul 02 14:17:33 2013 +0100
+++ b/libinterp/octave-value/ov-scalar.cc	Fri Aug 30 13:45:07 2013 +0100
@@ -333,6 +333,7 @@
       SCALAR_MAPPER (sin, ::sin);
       SCALAR_MAPPER (sinh, ::sinh);
       SCALAR_MAPPER (sqrt, rc_sqrt);
+      SCALAR_MAPPER (square, ::square);
       SCALAR_MAPPER (tan, ::tan);
       SCALAR_MAPPER (tanh, ::tanh);
       SCALAR_MAPPER (finite, xfinite);
diff -r 975aff6fafe5 -r 44a293eafe37 libinterp/octave-value/ov.h
--- a/libinterp/octave-value/ov.h	Tue Jul 02 14:17:33 2013 +0100
+++ b/libinterp/octave-value/ov.h	Fri Aug 30 13:45:07 2013 +0100
@@ -1169,6 +1169,7 @@
   MAPPER_FORWARD (sin)
   MAPPER_FORWARD (sinh)
   MAPPER_FORWARD (sqrt)
+  MAPPER_FORWARD (square)
   MAPPER_FORWARD (tan)
   MAPPER_FORWARD (tanh)
 
diff -r 975aff6fafe5 -r 44a293eafe37 liboctave/numeric/lo-mappers.cc
--- a/liboctave/numeric/lo-mappers.cc	Tue Jul 02 14:17:33 2013 +0100
+++ b/liboctave/numeric/lo-mappers.cc	Fri Aug 30 13:45:07 2013 +0100
@@ -26,6 +26,7 @@
 #endif
 
 #include <cfloat>
+#include <algorithm>
 
 #include "lo-error.h"
 #include "lo-ieee.h"
@@ -189,48 +190,487 @@
 
 // complex -> complex mappers.
 
-Complex
-acos (const Complex& x)
+// Complex mappers using code in 'Branch Cuts for Complex Elementary Functions
+// (or: Much Ado About Nothing's Sign Bit)' by W Kahan (from 'The State of the
+// Art in Numerical Analysis,' eds. Iserles and Powell, 1987).
+// Comments in these functions have been taken from the paper, except where
+// context implies otherwise.
+
+double
+abs (const Complex& z)
 {
-  static Complex i (0, 1);
+// The obvious formula can produce errors bigger than one ulp, and could
+// over/underflow spuriously. Not so for what follows.
 
-  return -i * (log (x + i * (sqrt (1.0 - x*x))));
+// These constants must be correctly rounded to working precision;
+// consequently r2pl + t2pl = 1 + sqrt(2) to double that precision.
+  static const double r2 = sqrt(2.0);
+  static const double r2p1 = r2 + 1.0;
+  static const double t2p1 = 1.0 + sqrt(2.0) - r2p1;
+  
+  double x = fabs (z.real ());
+  double y = fabs (z.imag ());
+  double s = 0.0;
+
+  if (x < y)
+    {
+      double temp = x;
+      x = y;
+      y = temp;
+    }
+    
+  if (isinf (y))
+    x = y;
+  
+  double t = x - y;
+  
+  if (!(isinf (x)) && t != x)
+    {
+      if (t > y)
+        {
+          s = x/y;
+          s += sqrt (1.0 + s*s);
+        }
+      else
+        {
+          s = t/y;
+          t = (2.0 + s) * s;
+          s = ((t2p1 + t / (r2 + sqrt (2.0 + t))) + s) + r2p1;
+        }
+      s = y/s;
+    }
+    
+  return x + s;
 }
 
 Complex
-acosh (const Complex& x)
+acos (const Complex& z)
 {
-  return log (x + sqrt (x*x - 1.0));
+  static const double pi = 3.14159265358979323846;
+
+  double x, y;
+  
+  Complex sqrt_1mz = sqrt (1.0 - z);
+  double re_sqrt_1pz = real (sqrt (1.0 + z));
+  
+  if (re_sqrt_1pz == 0.0)
+    x = pi;
+  else
+    x = 2.0 * atan (real (sqrt_1mz) / re_sqrt_1pz);
+  
+  y = asinh (imag (sqrt (conj (1.0 + z)) * sqrt_1mz));
+
+  return Complex (x, y);
 }
 
 Complex
-asin (const Complex& x)
+acosh (const Complex& z)
 {
-  static Complex i (0, 1);
+  static const double pi = 3.14159265358979323846;
 
-  return -i * log (i*x + sqrt (1.0 - x*x));
+  double x, y;
+  
+  Complex sqrt_1pz = sqrt (1.0 + z);
+  double re_sqrt_1pz = real (sqrt_1pz);
+  
+  x = asinh (real (sqrt (conj (z - 1.0)) * sqrt_1pz));
+
+  if (re_sqrt_1pz == 0.0)
+    y = pi;
+  else
+    y = 2.0 * atan (imag (sqrt (z - 1.0)) / re_sqrt_1pz);
+  
+  return Complex (x, y);
+}
+
+static Complex
+cbox (const Complex& z)
+{
+// Defined only for zero and infinite arguments.
+  
+  double x = z.real ();
+  double y = z.imag ();
+  
+  if (x == 0.0 && y == 0.0)
+    return Complex (copysign (1.0, x), y);
+  else
+    {
+      assert (isinf (x) || isinf (y));
+      if (isinf (x))
+        {
+          if (isinf (y))
+            return Complex (copysign (1.0, x), copysign (1.0, y));
+          else
+            // This line of code is different from Kahan's, as the latter
+            // gives the wrong sign for the imaginary part
+            return Complex (copysign (1.0, x), copysign (0.0, y));
+        }
+      else
+        // This line of code is different from Kahan's, as the latter
+        // gives the wrong sign for the real part
+        return Complex (copysign (0.0, x), copysign (1.0, y));
+    }
+}
+
+double
+arg (const Complex& z)
+{
+  static const double pi = 3.14159265358979323846;
+  static const double pi_2 = 1.57079632679489661923;
+
+  double x = z.real ();
+  double y = z.imag ();
+
+  if (x == 0.0 && y == 0.0)
+    x = copysign (1.0, x);
+  
+  if (isinf (x) || isinf (y))
+    {
+      Complex z1 = cbox (z); // leaves signs unchanged
+      x = z1.real ();
+      y = z1.imag ();
+    }
+  
+  if (fabs (y) > fabs (x))
+    return copysign (pi_2, y) - atan (x/y);
+  else if (x < 0.0)
+    return copysign (pi, y) + atan (y/x);
+  else
+    return atan (y/x);
+}
+
+static Complex
+xasin (const Complex& z, bool is_hyperbolic)
+{
+// Code for inverse circular and hyperbolic sines. The latter is calculated
+// using the identity arsinh(z) = -i.arcsin(iz). Note that iz etc. are not
+// coded using z * Complex(0,1) etc. as they produce spurious NaN's when
+// z = Inf+0i and variations thereof.
+
+  static const double pi_2 = 1.57079632679489661923;
+
+  Complex z1;
+
+  if (is_hyperbolic)
+    z1 = Complex (-z.imag (), z.real ());
+  else
+    z1 = z;
+  
+  Complex sqrt_1pz = sqrt (1.0 + z1);
+  double re_sqrt_1mz2 = real (sqrt (1.0 - z1) * sqrt_1pz);
+  
+  double x, y;
+  
+  if (re_sqrt_1mz2 == 0.0)
+    x = pi_2;
+  else
+    x = atan (real (z1) / re_sqrt_1mz2);
+  
+  y = asinh (imag (sqrt (conj (1.0 - z1)) * sqrt_1pz));
+
+  if (is_hyperbolic)
+    return Complex (y, -x);
+  else
+    return Complex (x, y);
 }
 
 Complex
-asinh (const Complex& x)
+asin (const Complex& z)
 {
-  return log (x + sqrt (x*x + 1.0));
+  return xasin (z, false);
 }
 
 Complex
-atan (const Complex& x)
+asinh (const Complex& z)
 {
-  static Complex i (0, 1);
+  return xasin (z, true);
+}
 
-  return i * log ((i + x) / (i - x)) / 2.0;
+static Complex
+xatan (const Complex& z, bool is_hyperbolic)
+{
+// Code for inverse circular and hyperbolic tangents. The former is calculated
+// using the identity arctan(z) = -i.artanh(iz). Note that iz etc. are not
+// coded using z * Complex(0,1) etc. as they produce spurious NaN's when
+// z = Inf+0i and variations thereof.
+
+  static const double cutoff = sqrt (DBL_MAX) / 4.0;
+  static const double pi_2 = 1.57079632679489661923;
+  static const double pi_4 = 0.78539816339744830962;
+
+  double x, y;
+  
+  if (is_hyperbolic)
+    {
+      x = z.real ();
+      y = z.imag ();
+    }
+  else
+    {
+      x = -z.imag ();
+      y = z.real ();
+    }
+
+  double abs_x = fabs (x);
+  double abs_y = fabs (y);
+
+  double x1, y1;
+  
+  if (abs_x > cutoff || abs_y > cutoff) // to avoid overflow
+    {
+      x1 = real (1.0 / z);
+      y1 = copysign (pi_2, y);
+    }
+  else if (abs_x == 1.0)
+    {
+      if (abs_y == 0.0)
+        {
+          // This is an addition to Kahan's code, as the latter gave a
+          // non-infinite value for the real part.
+          x1 = copysign (octave_Inf, x);
+          y1 = copysign (pi_4, y);
+        }
+      else
+        {
+          // Kahan's calculation adds a minuscule amount (=1/t) to the
+          // denominator of the division in the next line of code, and in
+          // several other places. This is unnecessary, and also gives
+          // inaccurate results, especially in the region of +/- 1 +/- 0i.
+          x1 = copysign (log (sqrt (sqrt (4.0 + y*y)) / sqrt (abs_y)), x);
+          y1 = copysign (pi_2 + atan ((abs_y) / 2.0), y) / 2.0;
+        }
+    }
+  else // Normal case. Using log1p(u) = log(1+u) accurately even if u is tiny.
+    {
+      double u = 1.0 - x;
+      double y2 = y*y;
+      
+      x1 = log1p (4.0 * x / (u*u + y2)) / 4.0;
+      y1 = arg (Complex (u * (1.0 + x) - y2, 2.0 * y)) / 2.0;
+    }
+    
+  if (is_hyperbolic)
+    return Complex (x1, y1);
+  else
+    return Complex (y1, -x1);
 }
 
 Complex
-atanh (const Complex& x)
+atan (const Complex& z)
 {
-  return log ((1.0 + x) / (1.0 - x)) / 2.0;
+  return xatan (z, false);
 }
 
+Complex
+atanh (const Complex& z)
+{
+  return xatan (z, true);
+}
+
+static double
+ssqs (const double x, const double y, int& k)
+{
+// To compute r = |(x + iy) / 2^k|^2 scaled to avoid over/underflow
+
+  k = 0;
+  double r;
+
+  if (isinf (x) || isinf (y))
+    r = octave_Inf;
+  else
+    {
+      r = x*x + y*y;
+      if (isinf (r) || (r < DBL_MIN/DBL_EPSILON))
+        {
+          k = ilogb (std::max (fabs (x), fabs (y)));
+          double rx = scalbln (x, -k);
+          double ry = scalbln (y, -k);
+          r = rx*rx + ry*ry;
+        }
+    }
+  return r;
+}
+
+static Complex
+logs (const Complex& z, int j)
+{
+// For use with j != O only when |x+iy| is huge. This program is particularly
+// helpful for inverse trigonometric and hyperbolic functions that behave like
+// ln(2z) for huge |z|. This program uses a subprogram log1p(x) = ln(1+x)
+// presumed to be available with full relative accuracy for all tiny real x.
+// The accuracy of log1p influences the choice of thresholds tO, t1 and t2.
+
+  static const double t0 = 0.70710678118654752440, t1 = 1.25, t2 = 3.0;
+  static const double log_2 = 0.693147180559945309417232121458176568075;
+  
+  double x = z.real ();
+  double y = z.imag ();
+
+  int k = 0;
+  double r = ssqs (x, y, k); // Sum of squares scaled: See above
+  
+  double b = std::max (fabs (x), fabs (y));
+  double theta = std::min (fabs (x), fabs (y));
+  
+  if (k == 0 && t0 < b && (b <= t1 || r < t2))
+    r = log1p ((b - 1.0) * (b + 1.0) + theta * theta) / 2.0;
+  else
+    r = log (r) / 2.0 + (k + j) * log_2;
+  theta = arg (z);
+
+  return Complex (r, theta);
+}
+
+Complex
+log (const Complex& z)
+{
+  return logs (z, 0);
+}
+
+Complex
+square (const Complex& z)
+{
+  double x = z.real ();
+  double y = z.imag ();
+  
+  double x1 = (x - y) * (x + y); // Not x*x - y*y
+  double y1 = x*y;
+  y1 += y1; // One multiply, one add
+
+  // If a spurious NaN is created by overflow, it gets removed thus:
+
+  if (isnan (x1))
+    {
+      if (isinf (y1))
+        x1 = copysign (0.0, x);
+      else if (isinf (y))
+        x1 = -octave_Inf;
+      else if (isinf (x))
+        x1 = octave_Inf;
+    }
+  else if (isnan (y1) && isinf (x1))
+    {
+      // Kahan just makes y1 equal to copysign (0.0, y1), but for some reason
+      // signed 0 * signed Inf doesn't give the right sign to the NaN result.
+      y1 = copysign (0.0, x);
+      y1 *= copysign (1.0, y);
+    }
+
+  return Complex (x1, y1);
+}
+
+Complex
+sqrt (const Complex& z)
+{
+  double x = z.real ();
+  double y = z.imag ();
+
+  int k = 0;
+  double r = ssqs (x, y, k); // Sum of squares scaled: See above
+  
+  if (!(isnan (x)))
+    r = scalbln (fabs (x), -k) + sqrt (r);
+  if (k % 2)
+    k = (k - 1) / 2;
+  else
+    {
+      k = k/2 - 1;
+      r += r;
+    }
+  r = scalbln (sqrt (r), k);
+  // r is now sqrt (( |x + iy| + |x| ) / 2) without over/underflow
+  
+  double x1 = r;
+  double y1 = y;
+  
+  if (r != 0.0)
+    {
+      if (!(isinf (y1)))
+        y1 = (y1/r) / 2.0;
+      if (x < 0.0)
+        {
+          x1 = fabs (y1);
+          y1 = copysign (r, y);
+        }
+    }
+    
+  return Complex (x1, y1);
+}
+
+static Complex
+xtan (const Complex &z, bool is_hyperbolic)
+{
+// Code for circular and hyperbolic tangents. The former is calculated
+// using the identity tan(z) = -i.tanh(iz). Note that iz etc. are not
+// coded using z * Complex(0,1) etc. as they produce spurious NaN's when
+// z = Inf+0i and variations thereof.
+
+  static const double pi = 3.14159265358979323846;
+  static const double cutoff = asinh (DBL_MAX) / 4.0;
+
+  double x, y;
+  
+  if (is_hyperbolic)
+    {
+      x = z.real ();
+      y = z.imag ();
+    }
+  else
+    {
+      x = -z.imag ();
+      y = z.real ();
+    }
+  
+  Complex z1;
+
+  if (fabs (x) > cutoff) // Avoid overflow
+    z1 = Complex (copysign (1.0, x), copysign (0.0, y));
+  else
+    {
+      double s, r;
+  
+      s = sinh (x);
+      r = sqrt (1.0 + s*s); // = cosh(x)
+
+      double n; // NULL rather than &n in the next line should work, but doesn't
+  
+      if (modf (fabs (y) / pi, &n) == 0.5)
+        {
+          // Addition to Kahan's code, as a call to tan ((2n+1)*pi/2)
+          // doesn't give infinity as expected.
+          if (s == 0.0)
+            z1 = Complex (copysign (0.0, x), copysign (octave_Inf, y));
+          else
+            z1 = Complex (r / s, copysign (0.0, y));
+        }
+      else
+        {
+          double t = tan (y);
+          double b = 1.0 + t*t; // = 1/cos^2(y)
+          double d = 1.0 + b*s*s;
+          z1 = Complex (b * r * s / d, t / d);
+        }
+    }
+
+  if (is_hyperbolic)
+    return z1;
+  else
+    return Complex (z1.imag (), -z1.real ());
+}
+
+Complex
+tan (const Complex &z)
+{
+  return xtan (z, false);
+}
+
+Complex
+tanh (const Complex &z)
+{
+  return xtan (z, true);
+}
+
+
 // complex -> bool mappers.
 
 bool
@@ -418,48 +858,489 @@
 
 // complex -> complex mappers.
 
-FloatComplex
-acos (const FloatComplex& x)
+// Complex mappers using code in 'Branch Cuts for Complex Elementary Functions
+// (or: Much Ado About Nothing's Sign Bit)' by W Kahan (from 'The State of the
+// Art in Numerical Analysis,' eds. Iserles and Powell, 1987).
+// Comments in these functions have been taken from the paper, except where
+// context implies otherwise.
+
+float
+abs (const FloatComplex& z)
 {
-  static FloatComplex i (0, 1);
+// The obvious formula can produce errors bigger than one ulp, and could
+// over/underflow spuriously. Not so for what follows.
 
-  return -i * (log (x + i * (sqrt (static_cast<float>(1.0) - x*x))));
+// These constants must be correctly rounded to working precision;
+// consequently r2pl + t2pl = 1 + sqrt(2) to double that precision.
+  static const float r2 = sqrtf(2.0f);
+  static const float r2p1 = r2 + 1.0f;
+  static const float t2p1 = 1.0f + sqrtf(2.0f) - r2p1;
+  
+  float x = fabsf (z.real ());
+  float y = fabsf (z.imag ());
+  float s = 0.0f;
+
+  if (x < y)
+    {
+      float temp = x;
+      x = y;
+      y = temp;
+    }
+    
+  if (isinf (y))
+    x = y;
+  
+  float t = x - y;
+  
+  if (!(isinf (x)) && t != x)
+    {
+      if (t > y)
+        {
+          s = x/y;
+          s += sqrtf (1.0f + s*s);
+        }
+      else
+        {
+          s = t/y;
+          t = (2.0f + s) * s;
+          s = ((t2p1 + t / (r2 + sqrtf (2.0f + t))) + s) + r2p1;
+        }
+      s = y/s;
+    }
+    
+  return x + s;
 }
 
 FloatComplex
-acosh (const FloatComplex& x)
+acos (const FloatComplex& z)
 {
-  return log (x + sqrt (x*x - static_cast<float>(1.0)));
+  static const float pi = 3.14159265358979323846f;
+
+  float x, y;
+  
+  FloatComplex sqrt_1mz = sqrt (1.0f - z);
+  float re_sqrt_1pz = real (sqrt (1.0f + z));
+  
+  if (re_sqrt_1pz == 0.0f)
+    x = pi;
+  else
+    x = 2.0f * atanf (real (sqrt_1mz) / re_sqrt_1pz);
+  
+  y = asinhf (imag (sqrt (conj (1.0f + z)) * sqrt_1mz));
+
+  return FloatComplex (x, y);
 }
 
 FloatComplex
-asin (const FloatComplex& x)
+acosh (const FloatComplex& z)
 {
-  static FloatComplex i (0, 1);
+  static const float pi = 3.14159265358979323846f;
 
-  return -i * log (i*x + sqrt (static_cast<float>(1.0) - x*x));
+  float x, y;
+  
+  FloatComplex sqrt_1pz = sqrt (1.0f + z);
+  float re_sqrt_1pz = real (sqrt_1pz);
+  
+  x = asinhf (real (sqrt (conj (z - 1.0f)) * sqrt_1pz));
+
+  if (re_sqrt_1pz == 0.0f)
+    y = pi;
+  else
+    y = 2.0f * atanf (imag (sqrt (z - 1.0f)) / re_sqrt_1pz);
+  
+  return FloatComplex (x, y);
+}
+
+static FloatComplex
+cbox (const FloatComplex& z)
+{
+// Defined only for zero and infinite arguments.
+  
+  float x = z.real ();
+  float y = z.imag ();
+  
+  if (x == 0.0f && y == 0.0f)
+    return FloatComplex (copysignf (1.0f, x), y);
+  else
+    {
+      assert (isinf (x) || isinf (y));
+      if (isinf (x))
+        {
+          if (isinf (y))
+            return FloatComplex (copysignf (1.0f, x), copysignf (1.0f, y));
+          else
+            // This line of code is different from Kahan's, as the latter
+            // gives the wrong sign for the imaginary part
+            return FloatComplex (copysignf (1.0f, x), copysignf (0.0f, y));
+        }
+      else
+        // This line of code is different from Kahan's, as the latter
+        // gives the wrong sign for the real part
+        return FloatComplex (copysignf (0.0f, x), copysignf (1.0f, y));
+    }
+}
+
+float
+arg (const FloatComplex& z)
+{
+  static const float pi = 3.14159265358979323846f;
+  static const float pi_2 = 1.57079632679489661923f;
+
+  float x = z.real ();
+  float y = z.imag ();
+
+  if (x == 0.0f && y == 0.0f)
+    x = copysignf (1.0f, x);
+  
+  if (isinf (x) || isinf (y))
+    {
+      FloatComplex z1 = cbox (z); // leaves signs unchanged
+      x = z1.real ();
+      y = z1.imag ();
+    }
+  
+  if (fabsf (y) > fabsf (x))
+    return copysignf (pi_2, y) - atanf (x/y);
+  else if (x < 0.0f)
+    return copysignf (pi, y) + atanf (y/x);
+  else
+    return atanf (y/x);
+}
+
+static FloatComplex
+xasin (const FloatComplex& z, bool is_hyperbolic)
+{
+// Code for inverse circular and hyperbolic sines. The latter is calculated
+// using the identity arsinh(z) = -i.arcsin(iz). Note that iz etc. are not
+// coded using z * Complex(0,1) etc. as they produce spurious NaN's when
+// z = Inf+0i and variations thereof.
+
+  static const float pi_2 = 1.57079632679489661923f;
+
+  FloatComplex z1;
+
+  if (is_hyperbolic)
+    z1 = FloatComplex (-z.imag (), z.real ());
+  else
+    z1 = z;
+  
+  FloatComplex sqrt_1pz = sqrt (1.0f + z1);
+  float re_sqrt_1mz2 = real (sqrt (1.0f - z1) * sqrt_1pz);
+  
+  float x, y;
+  
+  if (re_sqrt_1mz2 == 0.0f)
+    x = pi_2;
+  else
+    x = atanf (real (z1) / re_sqrt_1mz2);
+  
+  y = asinhf (imag (sqrt (conj (1.0f - z1)) * sqrt_1pz));
+
+  if (is_hyperbolic)
+    return FloatComplex (y, -x);
+  else
+    return FloatComplex (x, y);
 }
 
 FloatComplex
-asinh (const FloatComplex& x)
+asin (const FloatComplex& z)
 {
-  return log (x + sqrt (x*x + static_cast<float>(1.0)));
+  return xasin (z, false);
 }
 
 FloatComplex
-atan (const FloatComplex& x)
+asinh (const FloatComplex& z)
 {
-  static FloatComplex i (0, 1);
+  return xasin (z, true);
+}
 
-  return i * log ((i + x) / (i - x)) / static_cast<float>(2.0);
+static FloatComplex
+xatan (const FloatComplex& z, bool is_hyperbolic)
+{
+// Code for inverse circular and hyperbolic tangents. The former is calculated
+// using the identity arctan(z) = -i.artanh(iz). Note that iz etc. are not
+// coded using z * Complex(0,1) etc. as they produce spurious NaN's when
+// z = Inf+0i and variations thereof.
+
+  static const float cutoff = sqrtf (FLT_MAX) / 4.0f;
+  static const float pi_2 = 1.57079632679489661923f;
+  static const float pi_4 = 0.78539816339744830962f;
+
+  float x, y;
+  
+  if (is_hyperbolic)
+    {
+      x = z.real ();
+      y = z.imag ();
+    }
+  else
+    {
+      x = -z.imag ();
+      y = z.real ();
+    }
+
+  float abs_x = fabsf (x);
+  float abs_y = fabsf (y);
+
+  float x1, y1;
+  
+  if (abs_x > cutoff || abs_y > cutoff) // to avoid overflow
+    {
+      x1 = real (1.0f / z);
+      y1 = copysignf (pi_2, y);
+    }
+  else if (abs_x == 1.0f)
+    {
+      if (abs_y == 0.0f)
+        {
+          // This is an addition to Kahan's code, as the latter gave a
+          // non-infinite value for the real part.
+          x1 = copysign (octave_Float_Inf, x);
+          y1 = copysign (pi_4, y);
+        }
+      else
+        {
+          // Kahan's calculation adds a minuscule amount (=1/t) to the
+          // denominator of the division in the next line of code, and in
+          // several other places. This is unnecessary, and also gives
+          // inaccurate results, especially in the region of +/- 1 +/- 0i.
+          x1 = copysignf (logf (sqrtf (sqrtf (4.0f + y*y)) / sqrtf (abs_y)),
+                          x);
+          y1 = copysignf (pi_2 + atanf ((abs_y) / 2.0f), y) / 2.0f;
+        }
+    }
+  else // Normal case. Using log1p(u) = log(1+u) accurately even if u is tiny.
+    {
+      float u = 1.0f - x;
+      float y2 = y*y;
+      
+      x1 = log1p (4.0f * x / (u*u + y2)) / 4.0f;
+      y1 = arg (FloatComplex (u * (1.0f + x) - y2, 2.0f * y)) / 2.0f;
+    }
+    
+  if (is_hyperbolic)
+    return FloatComplex (x1, y1);
+  else
+    return FloatComplex (y1, -x1);
 }
 
 FloatComplex
-atanh (const FloatComplex& x)
+atan (const FloatComplex& z)
 {
-  return log ((static_cast<float>(1.0) + x) / (static_cast<float>(1.0) - x)) / static_cast<float>(2.0);
+  return xatan (z, false);
 }
 
+FloatComplex
+atanh (const FloatComplex& z)
+{
+  return xatan (z, true);
+}
+
+static float
+ssqs (const float x, const float y, int& k)
+{
+// To compute r = |(x + iy) / 2^k|^2 scaled to avoid over/underflow
+
+  k = 0;
+  float r;
+
+  if (isinf (x) || isinf (y))
+    r = octave_Float_Inf;
+  else
+    {
+      r = x*x + y*y;
+      if (isinf (r) || (r < FLT_MIN/FLT_EPSILON))
+        {
+          k = ilogbf (std::max (fabsf (x), fabsf (y)));
+          double rx = scalblnf (x, -k);
+          double ry = scalblnf (y, -k);
+          r = rx*rx + ry*ry;
+        }
+    }
+  return r;
+}
+
+static FloatComplex
+logs (const FloatComplex& z, int j)
+{
+// For use with j != O only when |x+iy| is huge. This program is particularly
+// helpful for inverse trigonometric and hyperbolic functions that behave like
+// ln(2z) for huge |z|. This program uses a subprogram log1p(x) = ln(1+x)
+// presumed to be available with full relative accuracy for all tiny real x.
+// The accuracy of log1p influences the choice of thresholds tO, t1 and t2.
+
+  static const float t0 = 0.70710678118654752440f, t1 = 1.25f, t2 = 3.0f;
+  static const float log_2 = 0.693147180559945309417232121458176568075f;
+  
+  float x = z.real ();
+  float y = z.imag ();
+
+  int k = 0;
+  float r = ssqs (x, y, k); // Sum of squares scaled: See above
+  
+  float b = std::max (fabsf (x), fabsf (y));
+  float theta = std::min (fabsf (x), fabsf (y));
+  
+  if (k == 0 && t0 < b && (b <= t1 || r < t2))
+    r = log1pf ((b - 1.0f) * (b + 1.0f) + theta * theta) / 2.0f;
+  else
+    r = logf (r) / 2.0f + (k + j) * log_2;
+  theta = arg (z);
+
+  return FloatComplex (r, theta);
+}
+
+FloatComplex
+log (const FloatComplex& z)
+{
+  return logs (z, 0);
+}
+
+FloatComplex
+square (const FloatComplex& z)
+{
+  float x = z.real ();
+  float y = z.imag ();
+  
+  float x1 = (x - y) * (x + y); // Not x*x - y*y
+  float y1 = x*y;
+  y1 += y1; // One multiply, one add
+
+  // If a spurious NaN is created by overflow, it gets removed thus:
+
+  if (isnan (x1))
+    {
+      if (isinf (y1))
+        x1 = copysignf (0.0f, x);
+      else if (isinf (y))
+        x1 = -octave_Float_Inf;
+      else if (isinf (x))
+        x1 = octave_Float_Inf;
+    }
+  else if (isnan (y1) && isinf (x1))
+    {
+      // Kahan just makes y1 equal to copysign (0.0, y1), but for some reason
+      // signed 0 * signed Inf doesn't give the right sign to the NaN result.
+      y1 = copysignf (0.0f, x);
+      y1 *= copysignf (1.0f, y);
+    }
+
+  return FloatComplex (x1, y1);
+}
+
+FloatComplex
+sqrt (const FloatComplex& z)
+{
+  float x = z.real ();
+  float y = z.imag ();
+
+  int k = 0;
+  float r = ssqs (x, y, k); // Sum of squares scaled: See above
+  
+  if (!(isnan (x)))
+    r = scalblnf (fabsf (x), -k) + sqrtf (r);
+  if (k % 2)
+    k = (k - 1) / 2;
+  else
+    {
+      k = k/2 - 1;
+      r += r;
+    }
+  r = scalblnf (sqrtf (r), k);
+  // r is now sqrt (( |x + iy| + |x| ) / 2) without over/underflow
+  
+  float x1 = r;
+  float y1 = y;
+  
+  if (r != 0.0f)
+    {
+      if (!(isinf (y1)))
+        y1 = (y1/r) / 2.0f;
+      if (x < 0.0f)
+        {
+          x1 = fabsf (y1);
+          y1 = copysignf (r, y);
+        }
+    }
+    
+  return FloatComplex (x1, y1);
+}
+
+static FloatComplex
+xtan (const FloatComplex &z, bool is_hyperbolic)
+{
+// Code for circular and hyperbolic tangents. The former is calculated
+// using the identity tan(z) = -i.tanh(iz). Note that iz etc. are not
+// coded using z * Complex(0,1) etc. as they produce spurious NaN's when
+// z = Inf+0i and variations thereof.
+
+  static const float pi = 3.14159265358979323846f;
+  static const float cutoff = asinhf (FLT_MAX) / 4.0f;
+
+  float x, y;
+  
+  if (is_hyperbolic)
+    {
+      x = z.real ();
+      y = z.imag ();
+    }
+  else
+    {
+      x = -z.imag ();
+      y = z.real ();
+    }
+  
+  FloatComplex z1;
+
+  if (fabsf (x) > cutoff) // Avoid overflow
+    z1 = FloatComplex (copysignf (1.0f, x), copysign (0.0f, y));
+  else
+    {
+      float s, r;
+  
+      s = sinhf (x);
+      r = sqrtf (1.0f + s*s); // = cosh(x)
+
+      float n; // NULL rather than &n in the next line should work, but doesn't
+  
+      if (modff (fabsf (y) / pi, &n) == 0.5f)
+        {
+          // Addition to Kahan's code, as a call to tan ((2n+1)*pi/2)
+          // doesn't give infinity as expected.
+          if (s == 0.0f)
+            z1 = FloatComplex (copysignf (0.0f, x),
+                               copysignf (octave_Float_Inf, y));
+          else
+            z1 = FloatComplex (r / s, copysignf (0.0f, y));
+        }
+      else
+        {
+          float t = tanf (y);
+          float b = 1.0f + t*t; // = 1/cos^2(y)
+          float d = 1.0f + b*s*s;
+          z1 = FloatComplex (b * r * s / d, t / d);
+        }
+    }
+
+  if (is_hyperbolic)
+    return z1;
+  else
+    return FloatComplex (z1.imag (), -z1.real ());
+}
+
+FloatComplex
+tan (const FloatComplex &z)
+{
+  return xtan (z, false);
+}
+
+FloatComplex
+tanh (const FloatComplex &z)
+{
+  return xtan (z, true);
+}
+
+
 // complex -> bool mappers.
 
 bool
@@ -583,13 +1464,13 @@
 Complex
 rc_sqrt (double x)
 {
-  return x < 0.0 ? Complex (0.0, sqrt (-x)) : Complex (sqrt (x));
+  return x < 0.0 || x == -0.0 ? Complex (0.0, sqrt (-x)) : Complex (sqrt (x));
 }
 
 FloatComplex
 rc_sqrt (float x)
 {
-  return x < 0.0f ? FloatComplex (0.0f, sqrtf (-x)) : FloatComplex (sqrtf (x));
+  return x < 0.0f || x == -0.0f ? FloatComplex (0.0f, sqrtf (-x)) : FloatComplex (sqrtf (x));
 }
 
 bool
diff -r 975aff6fafe5 -r 44a293eafe37 liboctave/numeric/lo-mappers.h
--- a/liboctave/numeric/lo-mappers.h	Tue Jul 02 14:17:33 2013 +0100
+++ b/liboctave/numeric/lo-mappers.h	Fri Aug 30 13:45:07 2013 +0100
@@ -39,6 +39,7 @@
 inline double fix (double x) { return xtrunc (x); }
 inline double imag (double) { return 0.0; }
 inline double real (double x) { return x; }
+inline double square (double x) { return x*x; }
 extern OCTAVE_API double xround (double x);
 extern OCTAVE_API double xroundb (double x);
 extern OCTAVE_API double signum (double x);
@@ -102,12 +103,19 @@
   return xisnan (y) ? x : (x >= y ? x : y);
 }
 
-extern OCTAVE_API Complex acos (const Complex& x);
-extern OCTAVE_API Complex acosh (const Complex& x);
-extern OCTAVE_API Complex asin (const Complex& x);
-extern OCTAVE_API Complex asinh (const Complex& x);
-extern OCTAVE_API Complex atan (const Complex& x);
-extern OCTAVE_API Complex atanh (const Complex& x);
+extern OCTAVE_API double abs (const Complex& z);
+extern OCTAVE_API Complex acos (const Complex& z);
+extern OCTAVE_API Complex acosh (const Complex& z);
+extern OCTAVE_API double arg (const Complex& z);
+extern OCTAVE_API Complex asin (const Complex& z);
+extern OCTAVE_API Complex asinh (const Complex& z);
+extern OCTAVE_API Complex atan (const Complex& z);
+extern OCTAVE_API Complex atanh (const Complex& z);
+extern OCTAVE_API Complex log (const Complex& z);
+extern OCTAVE_API Complex square (const Complex& z);
+extern OCTAVE_API Complex sqrt (const Complex& z);
+extern OCTAVE_API Complex tan (const Complex& z);
+extern OCTAVE_API Complex tanh (const Complex& z);
 
 extern OCTAVE_API bool octave_is_NA (const Complex& x);
 extern OCTAVE_API bool octave_is_NaN_or_NA (const Complex& x);
@@ -125,6 +133,7 @@
 inline float fix (float x) { return xtrunc (x); }
 inline float imag (float) { return 0.0f; }
 inline float real (float x) { return x; }
+inline float square (float x) { return x*x; }
 extern OCTAVE_API float xround (float x);
 extern OCTAVE_API float xroundb (float x);
 extern OCTAVE_API float signum (float x);
@@ -168,12 +177,19 @@
   return xisnan (y) ? x : (x >= y ? x : y);
 }
 
-extern OCTAVE_API FloatComplex acos (const FloatComplex& x);
-extern OCTAVE_API FloatComplex acosh (const FloatComplex& x);
-extern OCTAVE_API FloatComplex asin (const FloatComplex& x);
-extern OCTAVE_API FloatComplex asinh (const FloatComplex& x);
-extern OCTAVE_API FloatComplex atan (const FloatComplex& x);
-extern OCTAVE_API FloatComplex atanh (const FloatComplex& x);
+extern OCTAVE_API float abs (const FloatComplex& z);
+extern OCTAVE_API FloatComplex acos (const FloatComplex& z);
+extern OCTAVE_API FloatComplex acosh (const FloatComplex& z);
+extern OCTAVE_API float arg (const FloatComplex& z);
+extern OCTAVE_API FloatComplex asin (const FloatComplex& z);
+extern OCTAVE_API FloatComplex asinh (const FloatComplex& z);
+extern OCTAVE_API FloatComplex atan (const FloatComplex& z);
+extern OCTAVE_API FloatComplex atanh (const FloatComplex& z);
+extern OCTAVE_API FloatComplex log (const FloatComplex& z);
+extern OCTAVE_API FloatComplex square (const FloatComplex& z);
+extern OCTAVE_API FloatComplex sqrt (const FloatComplex& z);
+extern OCTAVE_API FloatComplex tan (const FloatComplex& z);
+extern OCTAVE_API FloatComplex tanh (const FloatComplex& z);
 
 extern OCTAVE_API bool octave_is_NA (const FloatComplex& x);
 extern OCTAVE_API bool octave_is_NaN_or_NA (const FloatComplex& x);
diff -r 975aff6fafe5 -r 44a293eafe37 liboctave/operators/mx-inlines.cc
--- a/liboctave/operators/mx-inlines.cc	Tue Jul 02 14:17:33 2013 +0100
+++ b/liboctave/operators/mx-inlines.cc	Fri Aug 30 13:45:07 2013 +0100
@@ -234,7 +234,11 @@
 {
   for (size_t i = 0; i < n; i++)
     {
-      if (x[i].imag () != 0)
+      T im = x[i].imag ();
+
+      // Do not return true if imaginary part of element is negative zero.
+  
+      if (im != 0.0 || signbit(im) != 0)
         return false;
     }
 
