# HG changeset patch
# User John W. Eaton <jwe@octave.org>
# Date 1289282269 18000
# Node ID 94d9d412a2a07b902b1df7830edd893933dca167
# Parent  a44ba1cdfbb557c850ab54278e6d4dcbff559a28
improve Matlab compatibility of rem and mod

diff --git a/ChangeLog b/ChangeLog
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,7 @@
+2010-11-09  John W. Eaton  <jwe@octave.org>
+
+	* configure.ac: Don't check for trunc, copysign, or _copysign.
+
 2010-11-08  John W. Eaton  <jwe@octave.org>
 
 	* configure.ac (--without-opengl): New configure option.
diff --git a/configure.ac b/configure.ac
--- a/configure.ac
+++ b/configure.ac
@@ -1572,7 +1572,7 @@
   mkstemp pipe putenv \
   realpath resolvepath rindex roundl select setgrent setlocale \
   setpwent setvbuf siglongjmp \
-  strsignal tempnam tgammaf trunc umask \
+  strsignal tempnam tgammaf umask \
   uname utime waitpid \
   _chmod x_utime _utime32)
 
@@ -1744,8 +1744,8 @@
     AC_DEFINE(HAVE_ISNAN, 1, [Define if you have isnan().])
   ;;
   *)
-    AC_CHECK_FUNCS(finite isnan isinf copysign signbit)
-    AC_CHECK_FUNCS(_finite _isnan _copysign)
+    AC_CHECK_FUNCS(finite isnan isinf signbit)
+    AC_CHECK_FUNCS(_finite _isnan)
     AC_CHECK_DECLS(signbit, , , [#include <math.h>])
   ;;
 esac
diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,9 @@
+2010-11-09  John W. Eaton  <jwe@octave.org>
+
+	* lo-mappers.cc (xmod, xrem): New functions.
+	(mod): Delete.
+	* lo-mappers.h (xmod, xrem): Provide decls.
+
 2010-10-31  Michael Goffioul  <michael.goffioul@gmail.com>
 
 	* dim-vector.h (class dim_vector): tag with OCTAVE_API.
diff --git a/liboctave/lo-mappers.cc b/liboctave/lo-mappers.cc
--- a/liboctave/lo-mappers.cc
+++ b/liboctave/lo-mappers.cc
@@ -107,13 +107,85 @@
 }
 
 double
-mod (double x, double y)
+xmod (double x, double y)
 {
+  double retval;
+
   if (y == 0)
-    return x;
+    retval = x;
+  else
+    {
+      double q = x / y;
+
+      double n = floor (q);
+
+      if (D_NINT (y) != y)
+        {
+          if (D_NINT (q) == q)
+            n = q;
+          else
+            {
+              if (x >= -1 && x <= 1)
+                {
+                  if (std::abs (q - D_NINT (q)) < DBL_EPSILON)
+                    n = D_NINT (q);
+                }
+              else
+                {
+                  if (std::abs ((q - D_NINT (q))/ D_NINT (q)) < DBL_EPSILON)
+                    n = D_NINT (q);
+                }
+            }
+        }
+
+      retval = x - y * n;
+    }
+
+  if (x != y && y != 0)
+    retval = copysignf (retval, y);
+
+  return retval;
+}
 
-  double r = fmod (x, y);
-  return ((r < 0) != (y < 0)) ? y+r : r;
+double
+xrem (double x, double y)
+{
+  double retval;
+
+  if (y == 0)
+    retval = x;
+  else
+    {
+      double q = x / y;
+
+      double n = trunc (q);
+
+      if (D_NINT (y) != y)
+        {
+          if (D_NINT (q) == q)
+            n = q;
+          else
+            {
+              if (x >= -1 && x <= 1)
+                {
+                  if (std::abs (q - D_NINT (q)) < DBL_EPSILON)
+                    n = D_NINT (q);
+                }
+              else
+                {
+                  if (std::abs ((q - D_NINT (q))/ D_NINT (q)) < DBL_EPSILON)
+                    n = D_NINT (q);
+                }
+            }
+        }
+
+      retval = x - y * n;
+    }
+
+  if (x != y && y != 0)
+    retval = copysignf (retval, x);
+
+  return retval;
 }
 
 double
@@ -396,13 +468,85 @@
 }
 
 float
-mod (float x, float y)
+xmod (float x, float y)
 {
+  float retval;
+
   if (y == 0)
-    return x;
+    retval = x;
+  else
+    {
+      float q = x / y;
+
+      float n = floor (q);
+
+      if (F_NINT (y) != y)
+        {
+          if (F_NINT (q) == q)
+            n = q;
+          else
+            {
+              if (x >= -1 && x <= 1)
+                {
+                  if (std::abs (q - F_NINT (q)) < FLT_EPSILON)
+                    n = F_NINT (q);
+                }
+              else
+                {
+                  if (std::abs ((q - F_NINT (q))/ F_NINT (q)) < FLT_EPSILON)
+                    n = F_NINT (q);
+                }
+            }
+        }
+
+      retval = x - y * n;
+    }
+
+  if (x != y && y != 0)
+    retval = copysignf (retval, y);
+
+  return retval;
+}
 
-  float r = fmodf (x, y);
-  return ((r < 0) != (y < 0)) ? y+r : r;
+float
+xrem (float x, float y)
+{
+  float retval;
+
+  if (y == 0)
+    retval = x;
+  else
+    {
+      float q = x / y;
+
+      float n = truncf (q);
+
+      if (F_NINT (y) != y)
+        {
+          if (F_NINT (q) == q)
+            n = q;
+          else
+            {
+              if (x >= -1 && x <= 1)
+                {
+                  if (std::abs (q - F_NINT (q)) < FLT_EPSILON)
+                    n = F_NINT (q);
+                }
+              else
+                {
+                  if (std::abs ((q - F_NINT (q))/ F_NINT (q)) < FLT_EPSILON)
+                    n = F_NINT (q);
+                }
+            }
+        }
+
+      retval = x - y * n;
+    }
+
+  if (x != y && y != 0)
+    retval = copysignf (retval, x);
+
+  return retval;
 }
 
 float
diff --git a/liboctave/lo-mappers.h b/liboctave/lo-mappers.h
--- a/liboctave/lo-mappers.h
+++ b/liboctave/lo-mappers.h
@@ -38,7 +38,8 @@
 extern OCTAVE_API double xroundb (double x);
 extern OCTAVE_API double signum (double x);
 extern OCTAVE_API double xtrunc (double x);
-extern OCTAVE_API double mod (double x, double y);
+extern OCTAVE_API double xmod (double x, double y);
+extern OCTAVE_API double xrem (double x, double y);
 extern OCTAVE_API double xlog2 (double x); 
 extern OCTAVE_API Complex xlog2 (const Complex& x); 
 extern OCTAVE_API double xlog2 (double x, int& exp);
@@ -127,7 +128,8 @@
 extern OCTAVE_API float xroundb (float x);
 extern OCTAVE_API float signum (float x);
 extern OCTAVE_API float xtrunc (float x);
-extern OCTAVE_API float mod (float x, float y);
+extern OCTAVE_API float xmod (float x, float y);
+extern OCTAVE_API float xrem (float x, float y);
 extern OCTAVE_API float xlog2 (float x); 
 extern OCTAVE_API FloatComplex xlog2 (const FloatComplex& x); 
 extern OCTAVE_API float xlog2 (float x, int& exp);
diff --git a/src/ChangeLog b/src/ChangeLog
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,10 @@
+2010-11-09  John W. Eaton  <jwe@octave.org>
+
+	* gl-render.cc (make_marker_list): Call fmod instead of mod.
+
+	* data.cc (Frem): Use xrem instead of fmod and fmodf.
+	(Fmod): Use xmod instead of mod.
+
 2010-11-08  Kai Habel  <kai.habel@gmx.de>
 
 	* fltk-backend.cc (do_find_uimenu_children): Simplify code, Remove 
diff --git a/src/data.cc b/src/data.cc
--- a/src/data.cc
+++ b/src/data.cc
@@ -567,12 +567,12 @@
       else if (args(0).is_single_type () || args(1).is_single_type ())
         {
           if (args(0).is_scalar_type () && args(1).is_scalar_type ())
-            retval = fmodf (args(0).float_value (), args(1).float_value ());
+            retval = xrem (args(0).float_value (), args(1).float_value ());
           else
             {
               FloatNDArray a0 = args(0).float_array_value ();
               FloatNDArray a1 = args(1).float_array_value ();
-              retval = binmap<float> (a0, a1, fmodf, "rem");
+              retval = binmap<float> (a0, a1, xrem, "rem");
             }
         }
       else
@@ -580,19 +580,19 @@
           bool a0_scalar = args(0).is_scalar_type ();
           bool a1_scalar = args(1).is_scalar_type ();
           if (a0_scalar && a1_scalar)
-            retval = fmod (args(0).scalar_value (), args(1).scalar_value ());
+            retval = xrem (args(0).scalar_value (), args(1).scalar_value ());
           else if ((a0_scalar || args(0).is_sparse_type ()) 
                    && (a1_scalar || args(1).is_sparse_type ()))
             {
               SparseMatrix m0 = args(0).sparse_matrix_value ();
               SparseMatrix m1 = args(1).sparse_matrix_value ();
-              retval = binmap<double> (m0, m1, fmod, "rem");
+              retval = binmap<double> (m0, m1, xrem, "rem");
             }
           else
             {
               NDArray a0 = args(0).array_value ();
               NDArray a1 = args(1).array_value ();
-              retval = binmap<double> (a0, a1, fmod, "rem");
+              retval = binmap<double> (a0, a1, xrem, "rem");
             }
         }
     }
@@ -700,12 +700,12 @@
       else if (args(0).is_single_type () || args(1).is_single_type ())
         {
           if (args(0).is_scalar_type () && args(1).is_scalar_type ())
-            retval = mod (args(0).float_value (), args(1).float_value ());
+            retval = xmod (args(0).float_value (), args(1).float_value ());
           else
             {
               FloatNDArray a0 = args(0).float_array_value ();
               FloatNDArray a1 = args(1).float_array_value ();
-              retval = binmap<float> (a0, a1, mod, "mod");
+              retval = binmap<float> (a0, a1, xmod, "mod");
             }
         }
       else
@@ -713,19 +713,19 @@
           bool a0_scalar = args(0).is_scalar_type ();
           bool a1_scalar = args(1).is_scalar_type ();
           if (a0_scalar && a1_scalar)
-            retval = mod (args(0).scalar_value (), args(1).scalar_value ());
+            retval = xmod (args(0).scalar_value (), args(1).scalar_value ());
           else if ((a0_scalar || args(0).is_sparse_type ()) 
                    && (a1_scalar || args(1).is_sparse_type ()))
             {
               SparseMatrix m0 = args(0).sparse_matrix_value ();
               SparseMatrix m1 = args(1).sparse_matrix_value ();
-              retval = binmap<double> (m0, m1, mod, "mod");
+              retval = binmap<double> (m0, m1, xmod, "mod");
             }
           else
             {
               NDArray a0 = args(0).array_value ();
               NDArray a1 = args(1).array_value ();
-              retval = binmap<double> (a0, a1, mod, "mod");
+              retval = binmap<double> (a0, a1, xmod, "mod");
             }
         }
     }
diff --git a/src/gl-render.cc b/src/gl-render.cc
--- a/src/gl-render.cc
+++ b/src/gl-render.cc
@@ -3288,7 +3288,7 @@
         for (int i = 0; i < 2*5; i++)
           {
             ang = (-0.5 + double(i+1)/5) * M_PI;
-            r = 1.0 - (dr * mod(double(i+1), 2.0));
+            r = 1.0 - (dr * fmod(double(i+1), 2.0));
             glVertex2d (sz*r*cos(ang)/2, sz*r*sin(ang)/2);
           }
         glEnd ();
@@ -3304,7 +3304,7 @@
         for (int i = 0; i < 2*6; i++)
           {
             ang = (0.5 + double(i+1)/6.0) * M_PI;
-            r = 1.0 - (dr * mod(double(i+1), 2.0));
+            r = 1.0 - (dr * fmod(double(i+1), 2.0));
             glVertex2d (sz*r*cos(ang)/2, sz*r*sin(ang)/2);
           }
         glEnd ();