changeset 17393:1a1d0cff1977

Merged with main repo
author Andrej Lojdl <andrej.lojdl@gmail.com>
date Thu, 05 Sep 2013 22:56:31 +0200
parents 609d93683f15 (current diff) 15e2ad6372f7 (diff)
children bdda4992ccf9
files configure.ac libinterp/Makefile.am libinterp/corefcn/txt-eng-ft.cc
diffstat 14 files changed, 327 insertions(+), 229 deletions(-) [+]
line wrap: on
line diff
--- a/configure.ac
+++ b/configure.ac
@@ -41,7 +41,11 @@
 AC_CONFIG_AUX_DIR([build-aux])
 AC_CONFIG_MACRO_DIR([m4])
 
-AM_INIT_AUTOMAKE([1.11 tar-ustar])
+AM_INIT_AUTOMAKE([1.11 tar-ustar subdir-objects])
+
+## Add the option to enable silent rules, available since Automake 1.11
+## and included by default starting with Automake 1.13.
+AM_SILENT_RULES
 
 OCTAVE_CANONICAL_HOST
 
--- a/libgui/Makefile.am
+++ b/libgui/Makefile.am
@@ -20,8 +20,6 @@
 
 include $(top_srcdir)/build-aux/common.mk
 
-AUTOMAKE_OPTIONS = subdir-objects
-
 MOC_CPPFLAGS =
 
 octlib_LTLIBRARIES = liboctgui.la
--- a/libgui/src/find-files-model.cc
+++ b/libgui/src/find-files-model.cc
@@ -111,13 +111,13 @@
 }
 
 int 
-find_files_model::rowCount (const QModelIndex & p) const
+find_files_model::rowCount (const QModelIndex &) const
 {
   return _files.size();
 }
 
 int 
-find_files_model::columnCount (const QModelIndex & p) const
+find_files_model::columnCount (const QModelIndex &) const
 {
   return _columnNames.size ();
 }
--- a/libinterp/Makefile.am
+++ b/libinterp/Makefile.am
@@ -20,8 +20,6 @@
 
 include $(top_srcdir)/build-aux/common.mk
 
-AUTOMAKE_OPTIONS = subdir-objects
-
 ## Search local directories before those specified by the user.
 AM_CPPFLAGS = \
   -I$(top_srcdir)/liboctave/cruft/misc \
@@ -258,7 +256,7 @@
 ## though we don't want it.  It would be super awesome if automake
 ## would allow users to choose the header file extension.
 .yy.cc:
-	$(am__skipyacc) $(SHELL) $(YLWRAP) $< y.tab.c $@ y.tab.h $*.h y.output $*.output -- $(YACCCOMPILE)
+	$(AM_V_YACC)$(am__skipyacc) $(SHELL) $(YLWRAP) $< y.tab.c $@ y.tab.h $*.h y.output $*.output -- $(YACCCOMPILE)
 
 ## Special rules:
 ## Mostly for sources which must be built before rest of compilation.
--- a/libinterp/corefcn/ellipj.cc
+++ b/libinterp/corefcn/ellipj.cc
@@ -1,6 +1,6 @@
 /*
 
-// Author: Leopoldo Cerbaro <redbliss@libero.it>
+Copyright (C) 2013 Leopoldo Cerbaro <redbliss@libero.it>
 
 This file is part of Octave.
 
@@ -35,7 +35,7 @@
 }
 
 static void
-sncndn (double u, double m, double& sn, double& cn, double& dn, double& err)
+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;
@@ -43,7 +43,7 @@
 
   if (m < 0 || m > 1)
     {
-      warning ("ellipj: expecting 0 <= m <= 1"); /* -lc- */
+      warning ("ellipj: expecting 0 <= M <= 1"); // -lc-
       sn = cn = dn = lo_ieee_nan_value ();
       return;
     }
@@ -51,7 +51,7 @@
   double sqrt_eps = sqrt (std::numeric_limits<double>::epsilon ());
   if (m < sqrt_eps)
     {
-      /*  # For small m, ( Abramowitz and Stegun, Section 16.13 ) */
+      // 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);
@@ -61,7 +61,7 @@
     }
   else if ((1 - m) < sqrt_eps)
     {
-      /*  For m1 = (1-m) small ( Abramowitz and Stegun, Section 16.15 ) */
+      // For m1 = (1-m) small ( Abramowitz and Stegun, Section 16.15 )
       m1 = 1 - m;
       si_u = sinh (u);
       co_u = cosh (u);
@@ -108,25 +108,25 @@
 }
 
 static void
-sncndn (Complex& u, double m, Complex& sn, Complex& cn, Complex& dn,
+do_ellipj (const Complex& u, const 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);
+  do_ellipj (imag (u), m1, ss1, cc1, dd1, err);
   if (real (u) == 0)
     {
-      /* u is pure imag: Jacoby imag. transf. */
+      // 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 */
+      // u is generic complex
       double ss, cc, dd, ddd;
 
-      sncndn (real (u), m, ss, cc, dd, err);
+      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);
@@ -146,13 +146,14 @@
 If @var{u} is a column vector and @var{m} is a row vector, the\n\
 results are matrices with @code{length (@var{u})} rows and\n\
 @code{length (@var{m})} columns.  Otherwise, @var{u} and\n\
-@var{m} must conform and the results will be the same size.\n\
+@var{m} must conform in size and the results will be the same size as the\n\
+inputs.\n\
 \n\
 The value of @var{u} may be complex.\n\
-The value of @var{m} must be 0 @leq{} m @leq{} 1.\n\
+The value of @var{m} must be 0 @leq{} @var{m} @leq{} 1.\n\
 \n\
-@var{tol} is currently ignored (@sc{matlab} uses this to allow faster,\n\
-less accurate approximation).\n\
+The optinoal input @var{tol} is currently ignored (@sc{matlab} uses this to\n\
+allow faster, less accurate approximation).\n\
 \n\
 If requested, @var{err} contains the following status information\n\
 and is the same size as the result.\n\
@@ -165,9 +166,11 @@
 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\
+\n\
+Reference: Milton Abramowitz and Irene A Stegun,\n\
+@cite{Handbook of Mathematical Functions}, Chapter 16 (Sections 16.4, 16.13,\n\
+and 16.15), Dover, 1965.\n\
+\n\
 @seealso{ellipke}\n\
 @end deftypefn")
 {
@@ -197,7 +200,7 @@
       if (u_arg.is_scalar_type ())
         {
           if (u_arg.is_real_type ())
-            {  // u real
+            {  // u real, m scalar
               double u = args(0).double_value ();
 
               if (error_state)
@@ -205,65 +208,75 @@
                   gripe_ellipj_arg ("first");
                   return retval;
                 }
+
               double sn, cn, dn;
               double err = 0;
 
-              sncndn (u, m, sn, cn, dn, err);
-              retval (0) = sn;
-              retval (1) = cn;
-              retval (2) = dn;
-              if (nargout > 3) retval(3) =  err;
+              do_ellipj (u, m, sn, cn, dn, err);
+
+              if (nargout > 3)
+                retval(3) = err;
+              retval(2) = dn;
+              retval(1) = cn;
+              retval(0) = sn;
             }
           else
-            {  // u complex
+            {  // u complex, m scalar
               Complex u = u_arg.complex_value ();
 
               if (error_state)
                 {
-                  gripe_ellipj_arg ("second");
+                  gripe_ellipj_arg ("first");
                   return retval;
                 }
 
               Complex sn, cn, dn;
-              double err;
+              double err = 0;
 
-              sncndn (u, m, sn, cn, dn, err);
+              do_ellipj (u, m, sn, cn, dn, err);
 
-              retval (0) = sn;
-              retval (1) = cn;
-              retval (2) = dn;
-              if (nargout > 3) retval(3) = err;
+              if (nargout > 3)
+                retval(3) = err;
+              retval(2) = dn;
+              retval(1) = cn;
+              retval(0) = sn;
             }
         }
       else
-        {  /* u is matrix ( m is scalar ) */
-          ComplexMatrix u = u_arg.complex_matrix_value ();
-
+        {  // u is matrix, m is scalar
+          ComplexNDArray u = u_arg.complex_array_value ();
+          
           if (error_state)
             {
               gripe_ellipj_arg ("first");
               return retval;
             }
 
-          octave_idx_type nr = u.rows ();
-          octave_idx_type nc = u.cols ();
+          dim_vector sz_u = u.dims ();
 
-          ComplexMatrix sn (nr, nc), cn (nr, nc), dn (nr, nc);
-          Matrix err (nr, nc);
+          ComplexNDArray sn (sz_u), cn (sz_u), dn (sz_u);
+          NDArray err (sz_u);
 
-          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));
+          const Complex *pu = u.data ();
+          Complex *psn = sn.fortran_vec ();
+          Complex *pcn = cn.fortran_vec ();
+          Complex *pdn = dn.fortran_vec ();
+          double *perr = err.fortran_vec ();
+          octave_idx_type nel = u.numel ();
 
-          retval (0) = sn;
-          retval (1) = cn;
-          retval (2) = dn;
-          if (nargout > 3) retval(3) = err;
+          for (octave_idx_type i = 0; i < nel; i++)
+            do_ellipj (pu[i], m, psn[i], pcn[i], pdn[i], perr[i]);
+
+          if (nargout > 3)
+            retval(3) = err;
+          retval(2) = dn;
+          retval(1) = cn;
+          retval(0) = sn;
         }
     }
   else
     {
-      Matrix m = args(1).matrix_value ();
+      NDArray m = args(1).array_value ();
 
       if (error_state)
         {
@@ -271,17 +284,12 @@
           return retval;
         }
 
-      octave_idx_type mr = m.rows ();
-      octave_idx_type mc = m.cols ();
+      dim_vector sz_m = m.dims ();
 
       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);
-
+        {  // u is scalar, m is array
           if (u_arg.is_real_type ())
-            {
+            {  // u is real scalar, m is array
               double u = u_arg.double_value ();
 
               if (error_state)
@@ -290,129 +298,176 @@
                   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));
+              NDArray sn (sz_m), cn (sz_m), dn (sz_m);
+              NDArray err (sz_m);
 
-              retval (0) = sn;
-              retval (1) = cn;
-              retval (2) = dn;
-              if (nargout > 3)  retval(3) = err;
+              const double *pm = m.data ();
+              double *psn = sn.fortran_vec ();
+              double *pcn = cn.fortran_vec ();
+              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]);
+
+              if (nargout > 3)
+                retval(3) = err;
+              retval(2) = dn;
+              retval(1) = cn;
+              retval(0) = sn;
             }
           else
-            {
+            {  // u is complex scalar, m is array
               Complex u = u_arg.complex_value ();
+
+              if (error_state)
+                {
+                  gripe_ellipj_arg ("first");
+                  return retval;
+                }
+
+              ComplexNDArray sn (sz_m), cn (sz_m), dn (sz_m);
+              NDArray err (sz_m);
+
+              const double *pm = m.data ();
+              Complex *psn = sn.fortran_vec ();
+              Complex *pcn = cn.fortran_vec ();
+              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]);
+
+              if (nargout > 3)
+                retval(3) = err;
+              retval(2) = dn;
+              retval(1) = cn;
+              retval(0) = sn;
+            }
+        }
+      else
+        {  // u is array, m is array
+          if (u_arg.is_real_type ())
+            {  // u is real array, m is array
+              NDArray u = u_arg.array_value ();
+
               if (error_state)
                 {
                   gripe_ellipj_arg ("first");
                   return retval;
                 }
 
-              ComplexMatrix 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 matrix)
-          if (u_arg.is_real_type ())
-            {  // u real matrix
+              dim_vector sz_u = u.dims ();
 
-              Matrix u = u_arg.matrix_value ();
-              if (error_state)
-                {
-                  gripe_ellipj_arg ("first ");
-                  return retval;
-                }
+              if (sz_u.length () == 2 && sz_m.length () == 2
+                  && sz_u(1) == 1 && sz_m(0) == 1)
+                {  // u is real column vector, m is row vector
+                  octave_idx_type ur = sz_u(0);
+                  octave_idx_type mc = sz_m(1);
+                  dim_vector sz_out (ur, mc);
 
-              octave_idx_type ur = u.rows ();
-              octave_idx_type uc = u.cols ();
+                  NDArray sn (sz_out), cn (sz_out), dn (sz_out);
+                  NDArray err (sz_out);
 
-              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);
+                  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++)
-                      sncndn (cu(i), rm(j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
+                      do_ellipj (pu[i], pm[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;
+                  if (nargout > 3)
+                    retval(3) = err;
+                  retval(2) = dn;
+                  retval(1) = cn;
+                  retval(0) = sn;
                 }
-              else if (ur == mr && uc == mc)
+              else if (sz_m == sz_u)
                 {
-                  Matrix sn (ur, mc), cn (ur, mc), dn (ur, mc);
-                  Matrix err (ur,mc);
+                  NDArray sn (sz_m), cn (sz_m), dn (sz_m);
+                  NDArray err (sz_m);
 
-                  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));
+                  const double *pu = u.data ();
+                  const double *pm = m.data ();
+                  double *psn = sn.fortran_vec ();
+                  double *pcn = cn.fortran_vec ();
+                  double *pdn = dn.fortran_vec ();
+                  double *perr = err.fortran_vec ();
+                  octave_idx_type nel = m.numel ();
 
-                  retval (0) = sn;
-                  retval (1) = cn;
-                  retval (2) = dn;
-                  if (nargout > 3)  retval(3) = err;
+                  for (octave_idx_type i = 0; i < nel; i++)
+                    do_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;
+                  retval(0) = sn;
                 }
               else
-                error ("u m invalid");
+                error ("ellipj: Invalid size combination for U and M");
             }
           else
-            {  // u complex matrix
-              ComplexMatrix u = u_arg.complex_matrix_value ();
+            {  // u is complex array, m is array
+              ComplexNDArray u = u_arg.complex_array_value ();
+
               if (error_state)
                 {
                   gripe_ellipj_arg ("second");
                   return retval;
                 }
 
-              octave_idx_type ur = u.rows ();
-              octave_idx_type uc = u.cols ();
+              dim_vector sz_u = u.dims ();
 
-              if (mr == 1 && uc == 1)
-                {
-                  RowVector rm = m.row (0);
-                  ComplexColumnVector cu = u.column (0);
+              if (sz_u.length () == 2 && sz_m.length () == 2
+                  && sz_u(1) == 1 && sz_m(0) == 1)
+                {  // u is complex column vector, m is row vector
+                  octave_idx_type ur = sz_u(0);
+                  octave_idx_type mc = sz_m(1);
+                  dim_vector sz_out (ur, mc);
 
-                  ComplexMatrix sn (ur, mc), cn (ur, mc), dn (ur, mc);
-                  Matrix err (ur,mc);
+                  ComplexNDArray sn (sz_out), cn (sz_out), dn (sz_out);
+                  NDArray err (sz_out);
+
+                  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++)
-                      sncndn (cu(i), rm(j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
+                      do_ellipj (pu[i], pm[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;
+                  if (nargout > 3)
+                    retval(3) = err;
+                  retval(2) = dn;
+                  retval(1) = cn;
+                  retval(0) = sn;
                 }
-              else if (ur == mr && uc == mc)
+              else if (sz_m == sz_u)
                 {
-                  ComplexMatrix sn (ur, mc), cn (ur, mc), dn (ur, mc);
-                  Matrix err (ur,mc);
+                  ComplexNDArray sn (sz_m), cn (sz_m), dn (sz_m);
+                  NDArray err (sz_m);
 
-                  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));
+                  const Complex *pu = u.data ();
+                  const double  *pm = m.data ();
+                  Complex *psn = sn.fortran_vec ();
+                  Complex *pcn = cn.fortran_vec ();
+                  Complex *pdn = dn.fortran_vec ();
+                  double *perr = err.fortran_vec ();
+                  octave_idx_type nel = m.numel ();
 
-                  retval (0) = sn;
-                  retval (1) = cn;
-                  retval (2) = dn;
-                  if (nargout > 3)  retval(3) = err;
+                  for (octave_idx_type i = 0; i < nel; i++)
+                    do_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;
+                  retval(0) = sn;
                 }
               else
-                error ("u m invalid");
+                error ("ellipj: Invalid size combination for U and M");
             }
         }
     }  // m matrix
@@ -434,7 +489,6 @@
 %! [sn, cn, dn] = ellipj (U,M);
 %!
 %! ## Plotting
-%! c = colormap (hot (64));
 %! data = {sn,cn,dn};
 %! dname = {"sn","cn","dn"};
 %! for i=1:3
@@ -444,7 +498,7 @@
 %!   image (m,u,32*data{i}+32);
 %!   title (dname{i});
 %! endfor
-%! colormap (c);
+%! colormap (hot (64));
 
 %!demo
 %! N = 200;
@@ -848,9 +902,9 @@
 %!     ui =  y * 0.2;
 %!     ii = 1 + y + x*11;
 %!     [sn, cn, dn] = ellipj (ur + I * ui, m);
-%!     assert (SN(ii, 2), sn, tol);
-%!     assert (CN(ii, 2), cn, tol);
-%!     assert (DN(ii, 2), dn, tol);
+%!     assert (sn, SN(ii, 2), tol);
+%!     assert (cn, CN(ii, 2), tol);
+%!     assert (dn, DN(ii, 2), tol);
 %!   endfor
 %! endfor
 
@@ -891,7 +945,7 @@
 %! u6 = 0.2 + 0.6i; m6 = tan(pi/8)^4;
 %! res6 = [ 0.2369100139 + 0.624633635i, ...
 %!          1.16200643   - 0.1273503824i, ...
-%!          1.004913944 - 0.004334880912i ];
+%!          1.004913944  - 0.004334880912i ];
 %! [sn,cn,dn] = ellipj (u6,m6);
 %! assert ([sn,cn,dn], res6, 1e-8);
 
@@ -904,15 +958,23 @@
 %! assert ([sn,cn,dn], res7, 1e-10);
 
 %!test
-%! u=[0,pi/6,pi/4,pi/2]; m=0;
+%! u = [0,pi/6,pi/4,pi/2]; m=0;
 %! res = [0,1/2,1/sqrt(2),1;1,cos(pi/6),1/sqrt(2),0;1,1,1,1];
 %! [sn,cn,dn] = ellipj (u,m);
-%! assert ([sn;cn;dn],res, 100*eps);
+%! assert ([sn;cn;dn], res, 100*eps);
 %! [sn,cn,dn] = ellipj (u',0);
-%! assert ([sn,cn,dn],res', 100*eps);
+%! assert ([sn,cn,dn], res', 100*eps);
+
+## FIXME: need to check [real,complex]x[scalar,rowvec,colvec,matrix]x[u,m]
 
-## XXX FIXME XXX
-## need to check [real,complex]x[scalar,rowvec,colvec,matrix]x[u,m]
+## One test for u column vector x m row vector
+%!test
+%! u = [0,pi/6,pi/4,pi/2]';  m = [0 0 0 0];
+%! res = [0,1/2,1/sqrt(2),1;1,cos(pi/6),1/sqrt(2),0;1,1,1,1]';
+%! [sn,cn,dn] = ellipj (u,m);
+%! assert (sn, repmat (res(:,1), [1,4]), 100*eps);
+%! assert (cn, repmat (res(:,2), [1,4]), 100*eps);
+%! assert (dn, repmat (res(:,3), [1,4]), 100*eps);
 
 %!test
 %! ## Test Jacobi elliptic functions
@@ -921,27 +983,46 @@
 %! ## 1 February 2001
 %! u = [ 0.25; 0.25; 0.20; 0.20; 0.672; 0.5];
 %! m = [ 0.0;  1.0;  0.19; 0.81; 0.36;  0.9999999999];
-%! S = [ sin(0.25); tanh(0.25);
+%! S = [ sin(0.25);
+%!       tanh(0.25);
 %!       0.19842311013970879516;
 %!       0.19762082367187648571;
 %!       0.6095196917919021945;
 %!       0.4621171572617320908 ];
-%! C = [ cos(0.25); sech(0.25);
+%! C = [ cos(0.25);
+%!       sech(0.25);
 %!       0.9801164570409401062;
 %!       0.9802785369736752032;
 %!       0.7927709286533560550;
 %!       0.8868188839691764094 ];
-%! D = [ 1.0;  sech(0.25);
+%! D = [ 1.0;
+%!       sech(0.25);
 %!       0.9962526643271134302;
 %!       0.9840560289645665155;
 %!       0.9307281387786906491;
 %!       0.8868188839812167635 ];
 %! [sn,cn,dn] = ellipj (u,m);
-%! assert (sn,S,8*eps);
-%! assert (cn,C,8*eps);
-%! assert (dn,D,8*eps);
+%! assert (sn, S, 8*eps);
+%! assert (cn, C, 8*eps);
+%! assert (dn, D, 8*eps);
 
 %!error ellipj ()
 %!error ellipj (1)
 %!error ellipj (1,2,3,4)
+%!warning <expecting 0 <= M <= 1> ellipj (1,2);
+## FIXME: errors commented out untill lasterr() truly returns the last error.
+%!#error <expecting scalar or matrix as second argument> ellipj (1, "1")
+%!#error <expecting scalar or matrix as first argument> ellipj ("1", 1)
+%!#error <expecting scalar or matrix as first argument> ellipj ({1}, 1)
+%!#error <expecting scalar or matrix as first argument> ellipj ({1, 2}, 1)
+%!#error <expecting scalar or matrix as second argument> ellipj (1, {1, 2})
+%!#error <expecting scalar or matrix as first argument> ellipj ("1", [1, 2])
+%!#error <expecting scalar or matrix as first argument> ellipj ({1}, [1, 2])
+%!#error <expecting scalar or matrix as first argument> ellipj ({1}, [1, 2])
+%!#error <expecting scalar or matrix as first argument> ellipj ("1,2", [1, 2])
+%!#error <expecting scalar or matrix as first argument> ellipj ({1, 2}, [1, 2])
+%!error <Invalid size combination for U and M> ellipj ([1:4], [1:3])
+%!error <Invalid size combination for U and M> ellipj (complex (1:4,1:4), [1:3])
+
 */
+
--- a/libinterp/corefcn/octave-link.cc
+++ b/libinterp/corefcn/octave-link.cc
@@ -255,7 +255,7 @@
 
               std::list<std::string>::iterator it = items_lst.begin ();
 
-              for (unsigned int idx = 0; idx < nel; idx++)
+              for (int idx = 0; idx < nel; idx++)
                 {
                   items.xelem (idx) = *it;
                   it++;
--- a/libinterp/corefcn/txt-eng-ft.cc
+++ b/libinterp/corefcn/txt-eng-ft.cc
@@ -42,8 +42,8 @@
 #include "pr-output.h"
 #include "txt-eng-ft.h"
 
-// FIXME -- maybe issue at most one warning per glyph/font/size/weight
-// combination.
+// FIXME: maybe issue at most one warning per glyph/font/size/weight
+//        combination.
 
 static void
 gripe_missing_glyph (FT_ULong c)
@@ -62,10 +62,10 @@
 }
 
 #ifdef _MSC_VER
-// This is just a trick to avoid multiply symbols definition.
+// This is just a trick to avoid multiple symbol definitions.
 // PermMatrix.h contains a dllexport'ed Array<octave_idx_type>
-// that will make MSVC not to generate new instantiation and
-// use the imported one.
+// that will cause MSVC not to generate a new instantiation and
+// use the imported one instead.
 #include "PermMatrix.h"
 #endif
 
@@ -154,8 +154,7 @@
         FT_Done_FreeType (library);
 
 #if defined (HAVE_FONTCONFIG)
-      // FIXME -- Skip the call to FcFini because it can trigger the
-      // assertion
+      // FIXME: Skip the call to FcFini because it can trigger the assertion
       //
       //   octave: fccache.c:507: FcCacheFini: Assertion 'fcCacheChains[i] == ((void *)0)' failed.
       //
@@ -225,7 +224,7 @@
               FcDefaultSubstitute (pat);
               match = FcFontMatch (0, pat, &res);
 
-              // FIXME -- originally, this test also required that
+              // FIXME: originally, this test also required that
               // res != FcResultNoMatch.  Is that really needed?
               if (match)
                 {
@@ -284,8 +283,7 @@
     {
       if (face->generic.data)
         {
-          ft_key* pkey =
-            reinterpret_cast<ft_key*> (face->generic.data);
+          ft_key* pkey = reinterpret_cast<ft_key*> (face->generic.data);
 
           cache.erase (*pkey);
           delete pkey;
@@ -635,7 +633,8 @@
 
       std::string str = e.string_value ();
       size_t n = str.length (), curr = 0;
-      mbstate_t ps = { 0 };
+      mbstate_t ps;
+      memset (&ps, 0, sizeof (ps));  // Initialize state to 0.
       wchar_t wc;
 
       while (n > 0)
--- a/libinterp/dldfcn/__magick_read__.cc
+++ b/libinterp/dldfcn/__magick_read__.cc
@@ -71,12 +71,17 @@
 // 1) A grayscale jpeg image can report being indexed even though the
 //    JPEG format has no support for indexed images. We can at least
 //    fix this one.
+// 2) A PNG file is only an indexed image if color type orig is 3 (value comes
+//    from libpng)
 static bool
 is_indexed (const Magick::Image& img)
 {
   bool retval = false;
-
-  if (img.classType () == Magick::PseudoClass && img.magick () != "JPEG")
+  const std::string format = img.magick ();
+  if (img.classType () == Magick::PseudoClass
+      && format != "JPEG"
+      && (format != "PNG"
+          || const_cast<Magick::Image&> (img).attribute ("PNG:IHDR.color-type-orig") == "3"))
     retval = true;
 
   return retval;
@@ -101,7 +106,7 @@
 get_depth (Magick::Image& img)
 {
   octave_idx_type depth = img.depth ();
-  if (depth != 1
+  if (depth == 8
       && img.channelDepth (Magick::RedChannel)     == 1
       && img.channelDepth (Magick::CyanChannel)    == 1
       && img.channelDepth (Magick::OpacityChannel) == 1
@@ -169,14 +174,14 @@
   // can't call colorMapSize on const Magick::Image
   const octave_idx_type mapsize = img.colorMapSize ();
   Matrix cmap                   = Matrix (mapsize, 3); // colormap
-  Matrix amap                   = Matrix (mapsize, 3); // alpha map
+  ColumnVector amap             = ColumnVector (mapsize); // alpha map
   for (octave_idx_type i = 0; i < mapsize; i++)
     {
       const Magick::ColorRGB c = img.colorMap (i);
       cmap(i,0) = c.red   ();
       cmap(i,1) = c.green ();
       cmap(i,2) = c.blue  ();
-      amap(i,0) = c.alpha ();
+      amap(i)   = c.alpha ();
     }
   octave_value_list maps;
   maps(0) = cmap;
--- a/libinterp/parse-tree/oct-parse.in.yy
+++ b/libinterp/parse-tree/oct-parse.in.yy
@@ -4436,7 +4436,7 @@
   return retval;
 }
 
-DEFUN (__parse_file__, args, nargout,
+DEFUN (__parse_file__, args, ,
   "-*- texinfo -*-\n\
 @deftypefn {Built-in Function} {} __parse_file__ (@var{file}, @var{verbose})\n\
 Undocumented internal function.\n\
--- a/liboctave/Makefile.am
+++ b/liboctave/Makefile.am
@@ -20,8 +20,6 @@
 
 include $(top_srcdir)/build-aux/common.mk
 
-AUTOMAKE_OPTIONS = subdir-objects
-
 ## Run cruft dir with stand-alone Makefile.
 ## Eventually this will use module.mk syntax.
 SUBDIRS = cruft
--- a/liboctave/cruft/Makefile.am
+++ b/liboctave/cruft/Makefile.am
@@ -20,8 +20,6 @@
 
 include $(top_srcdir)/build-aux/common.mk
 
-AUTOMAKE_OPTIONS = subdir-objects
-
 ## Search local directories before those specified by the user.
 AM_CPPFLAGS = \
   -I$(top_builddir)/libgnu -I$(top_srcdir)/libgnu
--- a/scripts/Makefile.am
+++ b/scripts/Makefile.am
@@ -20,8 +20,6 @@
 
 include $(top_srcdir)/build-aux/common.mk
 
-AUTOMAKE_OPTIONS = subdir-objects
-
 EXTRA_DIST =
 
 CLEANFILES =
--- a/scripts/image/imread.m
+++ b/scripts/image/imread.m
@@ -104,7 +104,7 @@
     filename{2} = varargin{2};
   endif
 
-  [img, varargout{2:nargout}] = imageIO (@__imread__, "read", filename, varargin{:});
+  [img, varargout{1:nargout-1}] = imageIO (@__imread__, "read", filename, varargin{:});
 endfunction
 
 
--- a/scripts/specfun/ellipke.m
+++ b/scripts/specfun/ellipke.m
@@ -19,19 +19,22 @@
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn  {Function File} {} ellipke (@var{m})
-## @deftypefnx {Function File} {} ellipke (@var{m}, @var{tol})
+## @deftypefn  {Function File} {@var{k} =} ellipke (@var{m})
+## @deftypefnx {Function File} {@var{k} =} ellipke (@var{m}, @var{tol})
 ## @deftypefnx {Function File} {[@var{k}, @var{e}] =} ellipke (@dots{})
-## Compute complete elliptic integral of the first K(@var{m}) and second
+## Compute complete elliptic integrals of the first K(@var{m}) and second
 ## E(@var{m}) kind.
 ##
-## @var{m} is either real array or scalar with 0 @leq{} m @leq{} 1.
+## @var{m} must be a scalar or real array with -Inf @leq{} @var{m} @leq{} 1.
+##
+## The optional input @var{tol} is currently ignored (@sc{matlab} uses this
+## to allow a faster, less accurate approximation).
 ##
-## @var{tol} is currently ignored (@sc{matlab} uses this to allow faster,
-## less accurate approximation).
+## Called with only one output, elliptic integrals of the first kind are
+## returned.
 ##
-## Ref: Abramowitz, Milton and Stegun, Irene A. Handbook of Mathematical
-## Functions, Dover, 1965, Chapter 17.
+## Reference: Milton Abramowitz and Irene A. Stegun,
+## @cite{Handbook of Mathematical Functions}, Chapter 17, Dover, 1965.
 ## @seealso{ellipj}
 ## @end deftypefn
 
@@ -45,55 +48,54 @@
     print_usage ();
   endif
 
-  k = e = zeros (size (m));
   m = m(:);
-  if (any (!isreal (m)))
-    error ("ellipke must have real m");
-  endif
-  if (any (m > 1))
-    error ("ellipke must have m <= 1");
+  if (! isreal (m))
+    error ("ellipke: M must be real");
+  elseif (any (m > 1))
+    error ("ellipke: M must be <= 1");
   endif
 
-  Nmax = 16;
-  idx = find (m == 1);
-  if (!isempty (idx))
-    k(idx) = Inf;
-    e(idx) = 1;
-  endif
+  k = e = zeros (size (m));
 
-  idx = find (m == -Inf);
-  if (!isempty (idx))
-    k(idx) = 0;
-    e(idx) = Inf;
-  endif
+  ## Handle extreme values
+  idx_1 = (m == 1);
+  k(idx_1) = Inf;
+  e(idx_1) = 1;
+
+  idx_neginf = (m == -Inf);
+  k(idx_neginf) = 0;
+  e(idx_neginf) = Inf;
 
   ## 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);
+  Nmax = 16;
+  idx = !idx_1 & !idx_neginf;
+  if (any (idx))
+    idx_neg = find (m < 0 & !idx_neginf);
     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);
+    m(idx_neg) = -m(idx_neg) ./ (1 - m(idx_neg));
+    a = ones (sum (idx), 1);
     b = sqrt (1 - m(idx));
     c = sqrt (m(idx));
     f = 0.5;
-    sum = f*c.*c;
-    for n = 2:Nmax
+    sum = f*c.^2;
+    n = 2;
+    do
       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
-    endfor
-    if (n >= Nmax) error ("ellipke: not enough workspace"); endif
-    k(idx) = 0.5*pi./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);
+      f *= 2;
+      sum += f*c.^2;
+    until (all (c./a < eps) || (++n > Nmax))
+    if (n >= Nmax)
+      error ("ellipke: algorithm did not converge in %d iterations", Nmax);
+    endif
+    k(idx) = 0.5*pi ./ 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
@@ -102,17 +104,16 @@
 ## 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 ];
+%! m = [0.0; 0.01; 0.1; 0.5; 0.9; 0.99; 1.0];
 %! [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 ];
+%!          Inf ];
 %! e_exp = [1.5707963267948966192;
 %!          1.5668619420216682912;
 %!          1.5307576368977632025;
@@ -120,7 +121,6 @@
 %!          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);
 
@@ -153,6 +153,25 @@
 %! assert (k, k_exp, 1e-15);
 %! assert (e, e_exp, 1e-8);
 
+## Test negative values against "exact" solution from Mathematica.
+%! m = [-0.01; -1; -5; -100; -1000; -Inf];
+%! [k,e] = ellipke (m);
+%!
+%! k_exp = [1.5668912730681963584;
+%!          1.3110287771460599052;
+%!          0.9555039270640439337;
+%!          0.3682192486091410329;
+%!          0.1530293349884987857;
+%!          0];
+%! e_exp = [1.5747159850169884130;
+%!          1.9100988945138560089;
+%!          2.8301982463458773125;
+%!          10.209260919814572009;
+%!          31.707204053711259719;
+%!          Inf ];
+%! assert (k, k_exp, 8*eps);
+%! assert (e, e_exp, 8*eps (e_exp));
+
 ## Test input validation
 %!error ellipke ()
 %!error ellipke (1,2,3)