Mercurial > octave-nkf

changeset 16586:f423873d3275

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Style fixes for ellipj.cc, ellipke.m, and expint.m * ellipj.cc, ellipke.m, expint.m: Style fixes.
author Mike Miller <mtmiller@ieee.org>
date Sun, 28 Apr 2013 18:50:51 -0400
parents 1a3bfb14b5da
children a3fdd6041e64
files libinterp/corefcn/ellipj.cc scripts/specfun/ellipke.m scripts/specfun/expint.m
diffstat 3 files changed, 444 insertions(+), 400 deletions(-) [+]
line wrap: on
line diff
--- a/libinterp/corefcn/ellipj.cc	Sun Apr 28 17:08:11 2013 -0400
+++ b/libinterp/corefcn/ellipj.cc	Sun Apr 28 18:50:51 2013 -0400
@@ -1,40 +1,23 @@
 /*
- Copyright (C) 2001 Leopoldo Cerbaro <redbliss@libero.it>
-
- This program is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 3 of the License, or
- (at your option) any later version.
 
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- GNU General Public License for more details.
+Copyright (C) 2001 Leopoldo Cerbaro <redbliss@libero.it>
 
- You should have received a copy of the GNU General Public License
- along with this program; If not, see <http://www.gnu.org/licenses/>.
-
- Compute the Jacobi elliptic functions sn(u|m), cn(u|m) and dn(u|m) for 
- argument u (real or complex) and parameter m. 
+This file is part of Octave.
 
- usage: [sn,cn,dn] = ellipj(u,m[,tol])
- 
- u and can be complex.
- m is restricted to 0 <= m <= 1.
- They can be scalars, matrix and scalar, scalar and matrix,
- column and row, conformant matrices.
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
 
- modified so u can be complex.   Leopoldo Cerbaro redbliss@libero.it
- 
- Ref: Abramowitz, Milton and Stegun, Irene A
-      Handbook of Mathematical Functions, Dover, 1965
-      Chapter 16 (Sections 16.4, 16.13 and 16.15)
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
 
- Based upon ellipj.m  made by David Billinghurst <David.Billinghurst@riotinto.com>
- and besselj.cc
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, see
+<http://www.gnu.org/licenses/>.
 
- Author: Leopoldo Cerbaro <redbliss@libero.it>
- Created: 15 December 2001
 */
 
 #ifdef HAVE_CONFIG_H
@@ -46,115 +29,114 @@
 #include "lo-ieee.h"
 
 static void
-gripe_ellipj_arg ( const char *arg)
+gripe_ellipj_arg (const char *arg)
 {
   error ("ellipj: expecting scalar or matrix as %s argument", arg);
 }
 
-const double  eps = 2.220446049e-16;
-const int  Nmax = 16;
+static void
+sncndn (double u, 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;
 
-static void
-sncndn ( double u, double m, double& sn, double& cn, double& dn, double& err) {
-/* real */
-double sqrt_eps, 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;
+    }
 
-  if (m < 0. || m > 1.) {
-     warning ("ellipj: expecting 0. <= m <= 1."); /* -lc- */
-     sn = cn = dn = lo_ieee_nan_value ();
-     return;
-  }
-  sqrt_eps = sqrt(eps);
-  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 - 0.5*m*si_u*si_u;
-/*}}}*/
-  } else if ( (1.0 - m) < sqrt_eps ) {
-    /*  For m1 = (1-m) small ( Abramowitz and Stegun, Section 16.15 ) */
-    /*{{{*/
-        m1 = 1.0-m;
-        si_u = sinh(u);
-        co_u = cosh(u);
-        ta_u = tanh(u);
-        se_u = 1.0/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.0;
-        b    = sqrt(1.0-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] < eps) break; 
+  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) {
-           // fprintf(stderr, "Not enough workspace\n");
-           err = 1.;
-           return;
+      if (n >= Nmax - 1)
+        {
+          err = 1;
+          return;
         }
-        Nn = n;
-        for ( ii = 1;  n>0; ii = ii*2, --n) ; // pow(2, Nn)
-        phi = ii*a[Nn]*u;
-        for ( n = Nn; n > 0; --n) {
+      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.;
+          phi = (asin ((c[n]/a[n])* sin (phi)) + phi)/2;
         }
-        sn = sin(phi);
-        cn = cos(phi);
-        dn = cn/cos(t-phi);
-/*}}}*/
-  }
- return;
+      sn = sin (phi);
+      cn = cos (phi);
+      dn = cn/cos (t - phi);
+    }
 }
 
 static void
-sncndn ( Complex& u, double m, 
-         Complex& sn, Complex& cn, Complex& dn, double& err) {
-double m1 = 1.-m, ss1, cc1, dd1;
+sncndn (Complex& u, double m, Complex& sn, Complex& cn, Complex& dn,
+        double& err)
+{
+  double m1 = 1 - m, ss1, cc1, dd1;
 
-  sncndn( 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;
+  sncndn (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;
 
-    sncndn( real(u), m, ss, cc, dd, err);
+      sncndn (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); 
+      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);
-    /*}}}*/
-  }
- return;
+    }
 }
 
 DEFUN (ellipj, args, nargout,
   "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {[@var{sn}, @var{cn}, @var{dn}] =} \
+@deftypefn {Built-in Function} {[@var{sn}, @var{cn}, @var{dn}] =} \
 ellipj (@var{u}, @var{m}, @var{err})\n\
 Compute the Jacobi elliptic functions sn, cn, dn of complex argument and real parameter.\n\
 \n\
@@ -166,7 +148,7 @@
 @var{m} must conform and the results will be the same size.\n\
 \n\
 The value of @var{u} may be complex.\n\
-The value of @var{m} must be 0 <= m <= 1. .\n\
+The value of @var{m} must be 0 <= m <= 1.\n\
 \n\
 If requested, @var{err} contains the following status information\n\
 and is the same size as the result.\n\
@@ -178,218 +160,258 @@
 Error---no computation, algorithm termination condition not met,\n\
 return @code{NaN}.\n\
 @end enumerate\n\
+ Ref: Abramowitz, Milton and Stegun, Irene A\n\
+      Handbook of Mathematical Functions, Dover, 1965\n\
+      Chapter 16 (Sections 16.4, 16.13 and 16.15)\n\
 @end deftypefn")
 {
   octave_value_list retval;
 
   int nargin = args.length ();
 
-  if (nargin == 2 ) {
-      octave_value u_arg = args(0);
-      octave_value m_arg = args(1);
+  if (nargin != 2 )
+    {
+      print_usage ();
+      return retval;
+    }
 
-      if (m_arg.is_scalar_type ()) {  // m is scalar
-        double  m = args(1).double_value ();
+  octave_value u_arg = args(0);
+  octave_value m_arg = args(1);
 
-        if (! error_state) {
+  if (m_arg.is_scalar_type ())
+    {
+      double m = args(1).double_value ();
 
-          if (u_arg.is_scalar_type ()) {   /*  u scalar */
-            /*{{{*/
-            if (u_arg.is_real_type ()) {  // u real
-              double  u = args(0).double_value ();
+      if (error_state)
+        {
+          gripe_ellipj_arg ("second");
+          return retval;
+        }
 
-              if (! error_state) {
-                double sn, cn, dn;
-                double err=0;
-                octave_value result;
+      if (u_arg.is_scalar_type ())
+        {
+          if (u_arg.is_real_type ())
+            {  // u real
+              double u = args(0).double_value ();
 
-                sncndn(u, m, sn, cn, dn, err);
-                retval (0) = sn;
-                retval (1) = cn;
-                retval (2) = dn;
-                if (nargout > 3)
-                  retval(3) =  err;
-            } else 
-                gripe_ellipj_arg ( "first");
+              if (error_state)
+                {
+                  gripe_ellipj_arg ("first");
+                  return retval;
+                }
+              double sn, cn, dn;
+              double err = 0;
 
-            } else {  // u complex
+              sncndn (u, m, sn, cn, dn, err);
+              retval (0) = sn;
+              retval (1) = cn;
+              retval (2) = dn;
+              if (nargout > 3) retval(3) =  err;
+            }
+          else
+            {  // u complex
               Complex u = u_arg.complex_value ();
 
-              if (! error_state) {
-                Complex sn, cn, dn;
-                double err;
-                octave_value result;
+              if (error_state)
+                {
+                  gripe_ellipj_arg ("second");
+                  return retval;
+                }
+
+              Complex sn, cn, dn;
+              double err;
+
+              sncndn (u, m, sn, cn, dn, err);
 
-                sncndn( u, m, sn, cn, dn, err);
+              retval (0) = sn;
+              retval (1) = cn;
+              retval (2) = dn;
+              if (nargout > 3) retval(3) = err;
+            }
+        }
+      else
+        {  /* u is matrix ( m is scalar ) */
+          ComplexMatrix u = u_arg.complex_matrix_value ();
+
+          if (error_state)
+            {
+              gripe_ellipj_arg ("first");
+              return retval;
+            }
+
+          octave_idx_type nr = u.rows ();
+          octave_idx_type nc = u.cols ();
+
+          ComplexMatrix sn (nr, nc), cn (nr, nc), dn (nr, nc);
+          Matrix err (nr, nc);
+
+          for (octave_idx_type j = 0; j < nc; j++)
+            for (octave_idx_type i = 0; i < nr; i++)
+              sncndn (u(i,j), m, sn(i,j), cn(i,j), dn(i,j), err(i,j));
 
-                retval (0) = sn;
-                retval (1) = cn;
-                retval (2) = dn;
-                if (nargout > 3)  retval(3) = err;
-              } else
-                gripe_ellipj_arg ( "second");
+          retval (0) = sn;
+          retval (1) = cn;
+          retval (2) = dn;
+          if (nargout > 3) retval(3) = err;
+        }
+    }
+  else
+    {
+      Matrix m = args(1).matrix_value ();
+
+      if (error_state)
+        {
+          gripe_ellipj_arg ("second");
+          return retval;
+        }
+
+      octave_idx_type mr = m.rows ();
+      octave_idx_type mc = m.cols ();
+
+      if (u_arg.is_scalar_type ())
+        {    /* u is scalar */
+          octave_idx_type nr = m.rows ();
+          octave_idx_type nc = m.cols ();
+          Matrix err (nr, nc);
+
+          if (u_arg.is_real_type ())
+            {
+              double u = u_arg.double_value ();
+
+              if (error_state)
+                {
+                  gripe_ellipj_arg ("first");
+                  return retval;
+                }
+
+              Matrix sn (nr, nc), cn (nr, nc), dn (nr, nc);
+              for (octave_idx_type j = 0; j < nc; j++)
+                for (octave_idx_type i = 0; i < nr; i++)
+                  sncndn (u, m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
+
+              retval (0) = sn;
+              retval (1) = cn;
+              retval (2) = dn;
+              if (nargout > 3)  retval(3) = err;
             }
-            /*}}}*/
-          } else {  /* u is matrix ( m is scalar ) */
-            /*{{{*/
-            ComplexMatrix u = u_arg.complex_matrix_value ();
-
-            if (! error_state) {
-              octave_value result;
-              int nr = u.rows ();
-              int nc = u.cols ();
+          else
+            {
+              Complex u = u_arg.complex_value ();
+              if (error_state)
+                {
+                  gripe_ellipj_arg ("first");
+                  return retval;
+                }
 
               ComplexMatrix sn (nr, nc), cn (nr, nc), dn (nr, nc);
-              Matrix err (nr, nc);
-
-              for (int j = 0; j < nc; j++)
-                for (int i = 0; i < nr; i++)
-                  sncndn (u(i,j), m, sn(i,j), cn(i,j), dn(i,j), err(i,j));
-
-                retval (0) = sn;
-                retval (1) = cn;
-                retval (2) = dn;
-                if (nargout > 3)  retval(3) = err;
-            } else
-                gripe_ellipj_arg ( "first");
-            /*}}}*/
-          }
-        } else
-            gripe_ellipj_arg ( "second");
-     } else { // m is matrix
-       Matrix m = args(1).matrix_value ();
-
-       if (! error_state) {
-         int mr = m.rows ();
-         int mc = m.cols ();
-
-         if (u_arg.is_scalar_type ()) {    /* u is scalar */
-           /*{{{*/
-           octave_value result;
-           int nr = m.rows ();
-           int nc = m.cols ();
-           Matrix err (nr, nc);
-
-           if (u_arg.is_real_type ()) {
-             double  u = u_arg.double_value ();
-             Matrix sn (nr, nc), cn (nr, nc), dn (nr, nc);
-             if (! error_state) {
-               for (int j = 0; j < nc; j++)
-                 for (int i = 0; i < nr; i++)
-                   sncndn (u, m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
-
-               retval (0) = sn;
-               retval (1) = cn;
-               retval (2) = dn;
-               if (nargout > 3)  retval(3) = err;
-             } else
-               gripe_ellipj_arg ( "first");
-           } else {
-             Complex u = u_arg.complex_value ();
-             ComplexMatrix sn (nr, nc), cn (nr, nc), dn (nr, nc);
-             if (! error_state) {
-               for (int j = 0; j < nc; j++)
-                 for (int i = 0; i < nr; i++)
-                   sncndn (u, m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
-               retval (0) = sn;
-               retval (1) = cn;
-               retval (2) = dn;
-               if (nargout > 3)  retval(3) = err;
-             } else
-               gripe_ellipj_arg ( "first");
-           }
-           /*}}}*/
-         } else {    // u is matrix  (m is matrix)
-           /*{{{*/
-           if (u_arg.is_real_type ()) {  // u real matrix
+              for (octave_idx_type j = 0; j < nc; j++)
+                for (octave_idx_type i = 0; i < nr; i++)
+                  sncndn (u, m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
+              retval (0) = sn;
+              retval (1) = cn;
+              retval (2) = dn;
+              if (nargout > 3)  retval(3) = err;
+            }
+        }
+      else
+        {    // u is matrix  (m is matrix)
+          if (u_arg.is_real_type ())
+            {  // u real matrix
 
               Matrix u = u_arg.matrix_value ();
-              if (! error_state) {
-                int ur = u.rows ();
-                int uc = u.cols ();
-
-              if (mr == 1 && uc == 1) {  // u column, m row
-                RowVector rm = m.row ((octave_idx_type)0);
-                ColumnVector cu = u.column ((octave_idx_type)0);
+              if (error_state)
+                {
+                  gripe_ellipj_arg ("first ");
+                  return retval;
+                }
 
-                Matrix sn (ur, mc), cn (ur, mc), dn (ur, mc);
-                Matrix err(ur,mc);
-//               octave_value result;
+              octave_idx_type ur = u.rows ();
+              octave_idx_type uc = u.cols ();
 
-                for (int j = 0; j < mc; j++)
-                  for (int i = 0; i < ur; i++)
-                    sncndn (cu(i), rm(j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
+              if (mr == 1 && uc == 1)
+                {  // u column, m row
+                  RowVector rm = m.row (0);
+                  ColumnVector cu = u.column (0);
+
+                  Matrix sn (ur, mc), cn (ur, mc), dn (ur, mc);
+                  Matrix err (ur,mc);
 
-                retval (0) = sn;
-                retval (1) = cn;
-                retval (2) = dn;
-                if (nargout > 3)  retval(3) = err;
-              } else if (ur == mr && uc == mc) {
-                Matrix sn (ur, mc), cn (ur, mc), dn (ur, mc);
-                Matrix err(ur,mc);
-//               octave_value result;
-
-                for (int j = 0; j < uc; j++)
-                 for (int i = 0; i < ur; i++)
-                  sncndn (u(i,j), m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
+                  for (octave_idx_type j = 0; j < mc; j++)
+                    for (octave_idx_type i = 0; i < ur; i++)
+                      sncndn (cu(i), rm(j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
 
-                retval (0) = sn;
-                retval (1) = cn;
-                retval (2) = dn;
-                if (nargout > 3)  retval(3) = err;
-              } else
-                 error("u m invalid");
-              } else
-                gripe_ellipj_arg ( "first ");
-            } else {  // u complex matrix
-              ComplexMatrix u = u_arg.complex_matrix_value ();
-              if (! error_state) {
-                int ur = u.rows ();
-                int uc = u.cols ();
+                  retval (0) = sn;
+                  retval (1) = cn;
+                  retval (2) = dn;
+                  if (nargout > 3)  retval(3) = err;
+                }
+              else if (ur == mr && uc == mc)
+                {
+                  Matrix sn (ur, mc), cn (ur, mc), dn (ur, mc);
+                  Matrix err (ur,mc);
+
+                  for (octave_idx_type j = 0; j < uc; j++)
+                    for (octave_idx_type i = 0; i < ur; i++)
+                      sncndn (u(i,j), m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
 
-              if (mr == 1 && uc == 1) {
-                RowVector rm = m.row ((octave_idx_type)0);
-                ComplexColumnVector cu = u.column ((octave_idx_type)0);
-
-                ComplexMatrix sn (ur, mc), cn (ur, mc), dn (ur, mc);
-                Matrix err(ur,mc);
-//               octave_value result;
+                  retval (0) = sn;
+                  retval (1) = cn;
+                  retval (2) = dn;
+                  if (nargout > 3)  retval(3) = err;
+                }
+              else
+                error ("u m invalid");
+            }
+          else
+            {  // u complex matrix
+              ComplexMatrix u = u_arg.complex_matrix_value ();
+              if (error_state)
+                {
+                  gripe_ellipj_arg ("second");
+                  return retval;
+                }
 
-                for (int j = 0; j < mc; j++)
-                  for (int i = 0; i < ur; i++)
-                    sncndn (cu(i), rm(j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
+              octave_idx_type ur = u.rows ();
+              octave_idx_type uc = u.cols ();
 
-                retval (0) = sn;
-                retval (1) = cn;
-                retval (2) = dn;
-                if (nargout > 3)  retval(3) = err;
-              } else if (ur == mr && uc == mc) {
+              if (mr == 1 && uc == 1)
+                {
+                  RowVector rm = m.row (0);
+                  ComplexColumnVector cu = u.column (0);
 
-                ComplexMatrix sn (ur, mc), cn (ur, mc), dn (ur, mc);
-                Matrix err(ur,mc);
-//               octave_value result;
+                  ComplexMatrix sn (ur, mc), cn (ur, mc), dn (ur, mc);
+                  Matrix err (ur,mc);
 
-                for (int j = 0; j < uc; j++)
-                 for (int i = 0; i < ur; i++)
-                  sncndn (u(i,j), m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
+                  for (octave_idx_type j = 0; j < mc; j++)
+                    for (octave_idx_type i = 0; i < ur; i++)
+                      sncndn (cu(i), rm(j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
+
+                  retval (0) = sn;
+                  retval (1) = cn;
+                  retval (2) = dn;
+                  if (nargout > 3)  retval(3) = err;
+                }
+              else if (ur == mr && uc == mc)
+                {
+                  ComplexMatrix sn (ur, mc), cn (ur, mc), dn (ur, mc);
+                  Matrix err (ur,mc);
 
-                retval (0) = sn;
-                retval (1) = cn;
-                retval (2) = dn;
-                if (nargout > 3)  retval(3) = err;
-              } else
-                 error("u m invalid");
-              } else
-                gripe_ellipj_arg ( "second");
+                  for (octave_idx_type j = 0; j < uc; j++)
+                    for (octave_idx_type i = 0; i < ur; i++)
+                      sncndn (u(i,j), m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
+
+                  retval (0) = sn;
+                  retval (1) = cn;
+                  retval (2) = dn;
+                  if (nargout > 3)  retval(3) = err;
+                }
+              else
+                error ("u m invalid");
             }
-           /*}}}*/
-         }
-       } else
-          gripe_ellipj_arg ( "second");
-     }  // m matrix
-   } else  // wrong n. of argin
-       print_usage ();
-   return retval;
+        }
+    }  // m matrix
+
+  return retval;
 }
 
 /*
@@ -444,7 +466,7 @@
 
 %!test
 %! k = (tan(pi/8.))^2; m = k*k;
-%! SN = [ 
+%! SN = [
 %! -1. + I * 0. ,  -0.8392965923 + 0. * I
 %! -1. + I * 0.2 ,  -0.8559363407 + 0.108250955 * I
 %! -1. + I * 0.4 ,  -0.906529758 + 0.2204040232 * I
--- a/scripts/specfun/ellipke.m	Sun Apr 28 17:08:11 2013 -0400
+++ b/scripts/specfun/ellipke.m	Sun Apr 28 18:50:51 2013 -0400
@@ -2,25 +2,30 @@
 ## Copyright (C) 2001 Paul Kienzle <pkienzle@users.sf.net>
 ## Copyright (C) 2003 Jaakko Ruohio
 ##
-## This program is free software; you can redistribute it and/or modify it under
-## the terms of the GNU General Public License as published by the Free Software
-## Foundation; either version 3 of the License, or (at your option) any later
-## version.
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
 ##
-## This program is distributed in the hope that it will be useful, but WITHOUT
-## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-## details.
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
 ##
-## You should have received a copy of the GNU General Public License along with
-## this program; if not, see <http://www.gnu.org/licenses/>.
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {[@var{k}, @var{e}] =} ellipke (@var{m}[,@var{tol}])
+## @deftypefn  {Function File} {} ellipke (@var{m})
+## @deftypefnx {Function File} {} ellipke (@var{m}, @var{tol})
+## @deftypefnx {Function File} {[@var{k}, @var{e}] =} ellipke (@dots{})
 ## Compute complete elliptic integral of first K(@var{m}) and second E(@var{m}).
 ##
 ## @var{m} is either real array or scalar with 0 <= m <= 1
-## 
+##
 ## @var{tol} will be ignored (@sc{Matlab} uses this to allow faster, less
 ## accurate approximation)
 ##
@@ -29,91 +34,89 @@
 ## @seealso{ellipj}
 ## @end deftypefn
 
-function [k,e] = ellipke( m )
+function [k, e] = ellipke (m)
 
   if (nargin < 1 || nargin > 2)
-    print_usage;
+    print_usage ();
   endif
 
-  k = e = zeros(size(m));
+  k = e = zeros (size (m));
   m = m(:);
-  if any(~isreal(m))
-    error("ellipke must have real m"); 
+  if any (!isreal (m))
+    error ("ellipke must have real m");
   endif
-  if any(m>1)
-    error("ellipke must have m <= 1");
+  if any (m > 1)
+    error ("ellipke must have m <= 1");
   endif
 
   Nmax = 16;
-  idx = find(m == 1);
-  if (!isempty(idx))
+  idx = find (m == 1);
+  if (!isempty (idx))
     k(idx) = Inf;
-    e(idx) = 1.0;
+    e(idx) = 1;
   endif
-      
-  idx = find(m == -Inf);
-  if (!isempty(idx))
-    k(idx) = 0.0;
+
+  idx = find (m == -Inf);
+  if (!isempty (idx))
+    k(idx) = 0;
     e(idx) = Inf;
   endif
 
   ## Arithmetic-Geometric Mean (AGM) algorithm
   ## ( Abramowitz and Stegun, Section 17.6 )
-  idx = find(m != 1 & m != -Inf);
-  if (!isempty(idx))
-    idx_neg = find(m < 0 & m != -Inf);
-    mult_k = 1./sqrt(1-m(idx_neg));
-    mult_e = sqrt(1-m(idx_neg));
-    m(idx_neg) = -m(idx_neg)./(1-m(idx_neg));
-    a = ones(length(idx),1);
-    b = sqrt(1.0-m(idx));
-    c = sqrt(m(idx));
+  idx = find (m != 1 & m != -Inf);
+  if (!isempty (idx))
+    idx_neg = find (m < 0 & m != -Inf);
+    mult_k = 1./sqrt (1 - m(idx_neg));
+    mult_e = sqrt (1 - m(idx_neg));
+    m(idx_neg) = -m(idx_neg)./(1 - m(idx_neg));
+    a = ones (length (idx), 1);
+    b = sqrt (1 - m(idx));
+    c = sqrt (m(idx));
     f = 0.5;
     sum = f*c.*c;
     for n = 2:Nmax
-      t = (a+b)/2;
-      c = (a-b)/2;
-      b = sqrt(a.*b);
+      t = (a + b)/2;
+      c = (a - b)/2;
+      b = sqrt (a.*b);
       a = t;
       f = f * 2;
       sum = sum + f*c.*c;
-      if all(c./a < eps), break; endif
+      if (all (c./a < eps)) break; endif
     endfor
-    if n >= Nmax, error("ellipke: not enough workspace"); endif
+    if (n >= Nmax) error ("ellipke: not enough workspace"); endif
     k(idx) = 0.5*pi./a;
-    e(idx) = 0.5*pi.*(1.0-sum)./a;
+    e(idx) = 0.5*pi.*(1 - sum)./a;
     k(idx_neg) = mult_k.*k(idx_neg);
     e(idx_neg) = mult_e.*e(idx_neg);
   endif
 
 endfunction
 
-## Test complete elliptic functions of first and second kind
-## against "exact" solution from Mathematica 3.0
+%% Test complete elliptic functions of first and second kind
+%% against "exact" solution from Mathematica 3.0
 %!test
 %! m = [0.0; 0.01; 0.1; 0.5; 0.9; 0.99; 1.0 ];
-%! [k,e] = ellipke(m);
+%! [k,e] = ellipke (m);
 %!
 %! # K(1.0) is really infinity - see below
-%! k_exp = [ 
-%!  1.5707963267948966192;
-%!  1.5747455615173559527;
-%!  1.6124413487202193982;
-%!  1.8540746773013719184;
-%!  2.5780921133481731882;
-%!  3.6956373629898746778;
-%!  0.0 ];
-%! e_exp = [
-%!  1.5707963267948966192;
-%!  1.5668619420216682912;
-%!  1.5307576368977632025;
-%!  1.3506438810476755025;
-%!  1.1047747327040733261;
-%!  1.0159935450252239356;
-%!  1.0 ];
-%! if k(7)==Inf, k(7)=0.0; endif;
-%! assert(k,k_exp,k,8*eps);
-%! assert(e,e_exp,8*eps);
+%! k_exp = [1.5707963267948966192;
+%!          1.5747455615173559527;
+%!          1.6124413487202193982;
+%!          1.8540746773013719184;
+%!          2.5780921133481731882;
+%!          3.6956373629898746778;
+%!          0.0 ];
+%! e_exp = [1.5707963267948966192;
+%!          1.5668619420216682912;
+%!          1.5307576368977632025;
+%!          1.3506438810476755025;
+%!          1.1047747327040733261;
+%!          1.0159935450252239356;
+%!          1.0 ];
+%! if k(7)==Inf, k(7)=0; endif;
+%! assert (k, k_exp, 8*eps);
+%! assert (e, e_exp, 8*eps);
 
 %% Test against A&S Table 17.1
 %!test
--- a/scripts/specfun/expint.m	Sun Apr 28 17:08:11 2013 -0400
+++ b/scripts/specfun/expint.m	Sun Apr 28 18:50:51 2013 -0400
@@ -1,20 +1,23 @@
 ## Copyright (C) 2006 Sylvain Pelissier <sylvain.pelissier@gmail.com>
 ##
-## This program is free software; you can redistribute it and/or modify it under
-## the terms of the GNU General Public License as published by the Free Software
-## Foundation; either version 3 of the License, or (at your option) any later
-## version.
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
 ##
-## This program is distributed in the hope that it will be useful, but WITHOUT
-## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-## details.
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
 ##
-## You should have received a copy of the GNU General Public License along with
-## this program; if not, see <http://www.gnu.org/licenses/>.
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {@var{y} =} expint (@var{x})
+## @deftypefn {Function File} {} expint (@var{x})
 ## Compute the exponential integral,
 ## @verbatim
 ##                    infinity
@@ -23,14 +26,16 @@
 ##                   /
 ##                  x
 ## @end verbatim
-## @seealso{expint_E1, expint_Ei}
 ## @end deftypefn
 
-function y = expint(x)
+function y = expint (x)
+
   if (nargin != 1)
-    print_usage;
+    print_usage ();
   endif
-  y = expint_E1(x);
+
+  y = expint_E1 (x);
+
 endfunction
 
 ## -*- texinfo -*-
@@ -43,18 +48,28 @@
 ##                   /
 ##                  x
 ## @end verbatim
-## @seealso{expint, expint_Ei}
 ## @end deftypefn
 
-function y = expint_E1(x)
+function y = expint_E1 (x)
+
   if (nargin != 1)
-    print_usage;
+    print_usage ();
   endif
+
   y = x;
-  y(imag(x) > 0 & imag(x) != 0) = -expint_Ei(-y(imag(x) > 0 & imag(x) != 0)) -i.*pi;
-  y(imag(x) < 0 & imag(x) != 0) = -expint_Ei(-y(imag(x) < 0 & imag(x) != 0)) +i.*pi;
-  y(real(x) >= 0 & imag(x)==0) = -expint_Ei(-y(real(x) >= 0 & imag(x)==0));
-  y(real(x) < 0 & imag(x)==0) = -expint_Ei(-y(real(x) < 0 & imag(x)==0)) -i.*pi;
+
+  idx = (imag (x) > 0 & imag (x) != 0);
+  y(idx) = -expint_Ei (-y(idx)) - i.*pi;
+
+  idx = (imag (x) < 0 & imag (x) != 0);
+  y(idx) = -expint_Ei (-y(idx)) + i.*pi;
+
+  idx = (real (x) >= 0 & imag (x) == 0);
+  y(idx) = -expint_Ei (-y(idx));
+
+  idx = (real (x) < 0 & imag (x) == 0);
+  y(idx) = -expint_Ei (-y(idx)) - i.*pi;
+
 endfunction
 
 ## -*- texinfo -*-
@@ -67,39 +82,43 @@
 ##                     /
 ##                     -x
 ## @end verbatim
-## @seealso{expint, expint_E1}
 ## @end deftypefn
 
-function y = expint_Ei(x)
+function y = expint_Ei (x)
+
   if (nargin != 1)
-    print_usage;
+    print_usage ();
   endif
-  y = zeros(size(x));
-  F = @(x) exp(-x)./x;
-  s = prod(size(x));
+
+  y = zeros (size (x));
+  F = @(x) exp (-x)./x;
+  s = prod (size (x));
+
   for t = 1:s;
-    if(x(t)<0 && imag(x(t)) == 0)
-      y(t) = -quad(F,-x(t),Inf);
+    if (x(t) < 0 && imag (x(t)) == 0)
+      y(t) = -quad (F, -x(t), Inf);
     else
-      if(abs(x(t)) > 2 && imag(x(t)) == 0)
-        y(t) = expint_Ei(2) - quad(F,-x(t),-2);
+      if (abs (x(t)) > 2 && imag (x(t)) == 0)
+        y(t) = expint_Ei (2) - quad (F, -x(t), -2);
       else
-        if(abs(x(t)) >= 10)
-          if(imag(x(t)) <= 0)
+        if (abs (x(t)) >= 10)
+          if (imag (x(t)) <= 0)
             a1 = 4.03640;
             a2 = 1.15198;
             b1 = 5.03637;
             b2 = 4.19160;
-            y(t) = -(x(t).^2 - a1.*x(t) + a2)./((x(t).^2-b1.*x(t)+b2).*(-x(t)).*exp(-x(t)))-i.*pi;
+            y(t) = -(x(t).^2 - a1.*x(t) + a2) ...
+                   ./ ((x(t).^2 - b1.*x(t) + b2) .* (-x(t)) .* exp (-x(t))) ...
+                   - i.*pi;
           else
-            y(t) = conj(expint_Ei(conj(x(t))));
+            y(t) = conj (expint_Ei (conj (x(t))));
           endif;
         ## Serie Expansion
         else
           for k = 1:100;
-            y(t) = y(t) + x(t).^k./(k.*factorial(k));
+            y(t) = y(t) + x(t).^k ./ (k.*factorial (k));
           endfor
-          y(t) = 0.577215664901532860606512090082402431 + log(x(t)) + y(t);
+          y(t) = 0.577215664901532860606512090082402431 + log (x(t)) + y(t);
         endif
       endif
     endif