--- a/libinterp/corefcn/ellipj.cc
+++ b/libinterp/corefcn/ellipj.cc
@@ -26,7 +26,7 @@
 
 #include "defun.h"
 #include "error.h"
-#include "lo-ieee.h"
+#include "lo-specfun.h"
 
 static void
 gripe_ellipj_arg (const char *arg)
@@ -34,106 +34,6 @@
   error ("ellipj: expecting scalar or matrix as %s argument", arg);
 }
 
-static void
-do_ellipj (const double u, const double m, double& sn, double& cn, double& dn, double& err)
-{
-  static const int Nmax = 16;
-  double m1, t=0, si_u, co_u, se_u, ta_u, b, c[Nmax], a[Nmax], phi;
-  int n, Nn, ii;
-
-  if (m < 0 || m > 1)
-    {
-      warning ("ellipj: expecting 0 <= M <= 1"); // -lc-
-      sn = cn = dn = lo_ieee_nan_value ();
-      return;
-    }
-
-  double sqrt_eps = sqrt (std::numeric_limits<double>::epsilon ());
-  if (m < sqrt_eps)
-    {
-      // For small m, ( Abramowitz and Stegun, Section 16.13 ) 
-      si_u = sin (u);
-      co_u = cos (u);
-      t = 0.25*m*(u - si_u*co_u);
-      sn = si_u - t * co_u;
-      cn = co_u + t * si_u;
-      dn = 1 - 0.5*m*si_u*si_u;
-    }
-  else if ((1 - m) < sqrt_eps)
-    {
-      // For m1 = (1-m) small ( Abramowitz and Stegun, Section 16.15 )
-      m1 = 1 - m;
-      si_u = sinh (u);
-      co_u = cosh (u);
-      ta_u = tanh (u);
-      se_u = 1/co_u;
-      sn = ta_u + 0.25*m1*(si_u*co_u - u)*se_u*se_u;
-      cn = se_u - 0.25*m1*(si_u*co_u - u)*ta_u*se_u;
-      dn = se_u + 0.25*m1*(si_u*co_u + u)*ta_u*se_u;
-    }
-  else
-    {
-      /*
-      //  Arithmetic-Geometric Mean (AGM) algorithm
-      //    ( Abramowitz and Stegun, Section 16.4 )
-      */
-
-      a[0] = 1;
-      b    = sqrt (1 - m);
-      c[0] = sqrt (m);
-      for (n = 1; n < Nmax; ++n)
-        {
-          a[n] = (a[n - 1] + b)/2;
-          c[n] = (a[n - 1] - b)/2;
-          b = sqrt (a[n - 1]*b);
-          if (c[n]/a[n] < std::numeric_limits<double>::epsilon ()) break;
-        }
-      if (n >= Nmax - 1)
-        {
-          err = 1;
-          return;
-        }
-      Nn = n;
-      for (ii = 1; n > 0; ii = ii*2, --n) ; // ii = pow(2,Nn)
-      phi = ii*a[Nn]*u;
-      for (n = Nn; n > 0; --n)
-        {
-          t = phi;
-          phi = (asin ((c[n]/a[n])* sin (phi)) + phi)/2;
-        }
-      sn = sin (phi);
-      cn = cos (phi);
-      dn = cn/cos (t - phi);
-    }
-}
-
-static void
-do_ellipj (const Complex& u, const double m, Complex& sn, Complex& cn, Complex& dn,
-        double& err)
-{
-  double m1 = 1 - m, ss1, cc1, dd1;
-
-  do_ellipj (imag (u), m1, ss1, cc1, dd1, err);
-  if (real (u) == 0)
-    {
-      // u is pure imag: Jacoby imag. transf.
-      sn = Complex (0, ss1/cc1);
-      cn = 1/cc1;         //    cn.imag = 0;
-      dn = dd1/cc1;       //    dn.imag = 0;
-    }
-  else
-    {
-      // u is generic complex
-      double ss, cc, dd, ddd;
-
-      do_ellipj (real (u), m, ss, cc, dd, err);
-      ddd = cc1*cc1 + m*ss*ss*ss1*ss1;
-      sn = Complex (ss*dd1/ddd, cc*dd*ss1*cc1/ddd);
-      cn = Complex (cc*cc1/ddd, -ss*dd*ss1*dd1/ddd);
-      dn = Complex (dd*cc1*dd1/ddd, -m*ss*cc*ss1/ddd);
-    }
-}
-
 DEFUN (ellipj, args, nargout,
   "-*- texinfo -*-\n\
 @deftypefn  {Built-in Function} {[@var{sn}, @var{cn}, @var{dn}, @var{err}] =} ellipj (@var{u}, @var{m})\n\
@@ -212,7 +112,7 @@
               double sn, cn, dn;
               double err = 0;
 
-              do_ellipj (u, m, sn, cn, dn, err);
+              ellipj (u, m, sn, cn, dn, err);
 
               if (nargout > 3)
                 retval(3) = err;
@@ -233,7 +133,7 @@
               Complex sn, cn, dn;
               double err = 0;
 
-              do_ellipj (u, m, sn, cn, dn, err);
+              ellipj (u, m, sn, cn, dn, err);
 
               if (nargout > 3)
                 retval(3) = err;
@@ -265,7 +165,7 @@
           octave_idx_type nel = u.numel ();
 
           for (octave_idx_type i = 0; i < nel; i++)
-            do_ellipj (pu[i], m, psn[i], pcn[i], pdn[i], perr[i]);
+            ellipj (pu[i], m, psn[i], pcn[i], pdn[i], perr[i]);
 
           if (nargout > 3)
             retval(3) = err;
@@ -309,7 +209,7 @@
               octave_idx_type nel = m.numel ();
 
               for (octave_idx_type i = 0; i < nel; i++)
-                do_ellipj (u, pm[i], psn[i], pcn[i], pdn[i], perr[i]);
+                ellipj (u, pm[i], psn[i], pcn[i], pdn[i], perr[i]);
 
               if (nargout > 3)
                 retval(3) = err;
@@ -338,7 +238,7 @@
               octave_idx_type nel = m.numel ();
 
               for (octave_idx_type i = 0; i < nel; i++)
-                do_ellipj (u, pm[i], psn[i], pcn[i], pdn[i], perr[i]);
+                ellipj (u, pm[i], psn[i], pcn[i], pdn[i], perr[i]);
 
               if (nargout > 3)
                 retval(3) = err;
@@ -376,7 +276,7 @@
 
                   for (octave_idx_type j = 0; j < mc; j++)
                     for (octave_idx_type i = 0; i < ur; i++)
-                      do_ellipj (pu[i], pm[j], sn(i,j), cn(i,j), dn(i,j), err(i,j));
+                      ellipj (pu[i], pm[j], sn(i,j), cn(i,j), dn(i,j), err(i,j));
 
                   if (nargout > 3)
                     retval(3) = err;
@@ -398,7 +298,7 @@
                   octave_idx_type nel = m.numel ();
 
                   for (octave_idx_type i = 0; i < nel; i++)
-                    do_ellipj (pu[i], pm[i], psn[i], pcn[i], pdn[i], perr[i]);
+                    ellipj (pu[i], pm[i], psn[i], pcn[i], pdn[i], perr[i]);
 
                   if (nargout > 3)
                     retval(3) = err;
@@ -436,7 +336,7 @@
 
                   for (octave_idx_type j = 0; j < mc; j++)
                     for (octave_idx_type i = 0; i < ur; i++)
-                      do_ellipj (pu[i], pm[j], sn(i,j), cn(i,j), dn(i,j), err(i,j));
+                      ellipj (pu[i], pm[j], sn(i,j), cn(i,j), dn(i,j), err(i,j));
 
                   if (nargout > 3)
                     retval(3) = err;
@@ -458,7 +358,7 @@
                   octave_idx_type nel = m.numel ();
 
                   for (octave_idx_type i = 0; i < nel; i++)
-                    do_ellipj (pu[i], pm[i], psn[i], pcn[i], pdn[i], perr[i]);
+                    ellipj (pu[i], pm[i], psn[i], pcn[i], pdn[i], perr[i]);
 
                   if (nargout > 3)
                     retval(3) = err;