octave-nkf: libinterp/corefcn/ellipj.cc comparison

comparison libinterp/corefcn/ellipj.cc @ 17502:578805a293e5

ellipj: Move numerical code into liboctave * lo-specfun.cc, lo-specfun.h (ellipj): New functions, adapted from Fellipj. * ellipj.cc (Fellipj): Call ellipj. (do_ellipj): Delete.

author	Mike Miller <mtmiller@ieee.org>
date	Thu, 26 Sep 2013 21:34:26 -0400
parents	991c7c812e38
children	5b916efea542

comparison

equal deleted inserted replaced

-:55680de6a897
+:578805a293e5
 #include <config.h>
 #endif
 #include "defun.h"
 #include "error.h"
-#include "lo-ieee.h"
+#include "lo-specfun.h"
 static void
 gripe_ellipj_arg (const char *arg)
 {
 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\
 }
 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;
 retval(2) = dn;
 retval(1) = cn;
 }
 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;
 retval(2) = dn;
 retval(1) = cn;
 Complex *pdn = dn.fortran_vec ();
 double *perr = err.fortran_vec ();
 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;
 retval(2) = dn;
 retval(1) = cn;
 double *pdn = dn.fortran_vec ();
 double *perr = err.fortran_vec ();
 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;
 retval(2) = dn;
 retval(1) = cn;
 Complex *pdn = dn.fortran_vec ();
 double *perr = err.fortran_vec ();
 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;
 retval(2) = dn;
 retval(1) = cn;
 const double *pu = u.data ();
 const double *pm = m.data ();
 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;
 retval(2) = dn;
 retval(1) = cn;
 double *pdn = dn.fortran_vec ();
 double *perr = err.fortran_vec ();
 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;
 retval(2) = dn;
 retval(1) = cn;
 const Complex *pu = u.data ();
 const double  *pm = m.data ();
 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;
 retval(2) = dn;
 retval(1) = cn;
 Complex *pdn = dn.fortran_vec ();
 double *perr = err.fortran_vec ();
 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;
 retval(2) = dn;
 retval(1) = cn;

Mercurial > hg > octave-nkf

comparison libinterp/corefcn/ellipj.cc @ 17502:578805a293e5