Mercurial > hg > octave-lyh

comparison libinterp/corefcn/ellipj.cc @ 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
comparison
equal deleted inserted replaced
16585:1a3bfb14b5da 16586:f423873d3275
1 /* 1 /*
2 Copyright (C) 2001 Leopoldo Cerbaro <redbliss@libero.it> 2
3 3 Copyright (C) 2001 Leopoldo Cerbaro <redbliss@libero.it>
4 This program is free software; you can redistribute it and/or modify 4
5 it under the terms of the GNU General Public License as published by 5 This file is part of Octave.
6 the Free Software Foundation; either version 3 of the License, or 6
7 (at your option) any later version. 7 Octave is free software; you can redistribute it and/or modify it
8 8 under the terms of the GNU General Public License as published by the
9 This program is distributed in the hope that it will be useful, 9 Free Software Foundation; either version 3 of the License, or (at your
10 but WITHOUT ANY WARRANTY; without even the implied warranty of 10 option) any later version.
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11
12 GNU General Public License for more details. 12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 You should have received a copy of the GNU General Public License 14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 along with this program; If not, see <http://www.gnu.org/licenses/>. 15 for more details.
16 16
17 Compute the Jacobi elliptic functions sn(u|m), cn(u|m) and dn(u|m) for 17 You should have received a copy of the GNU General Public License
18 argument u (real or complex) and parameter m. 18 along with Octave; see the file COPYING. If not, see
19 19 <http://www.gnu.org/licenses/>.
20 usage: [sn,cn,dn] = ellipj(u,m[,tol]) 20
21
22 u and can be complex.
23 m is restricted to 0 <= m <= 1.
24 They can be scalars, matrix and scalar, scalar and matrix,
25 column and row, conformant matrices.
26
27 modified so u can be complex. Leopoldo Cerbaro redbliss@libero.it
28
29 Ref: Abramowitz, Milton and Stegun, Irene A
30 Handbook of Mathematical Functions, Dover, 1965
31 Chapter 16 (Sections 16.4, 16.13 and 16.15)
32
33 Based upon ellipj.m made by David Billinghurst <David.Billinghurst@riotinto.com>
34 and besselj.cc
35
36 Author: Leopoldo Cerbaro <redbliss@libero.it>
37 Created: 15 December 2001
38 */ 21 */
39 22
40 #ifdef HAVE_CONFIG_H 23 #ifdef HAVE_CONFIG_H
41 #include <config.h> 24 #include <config.h>
42 #endif 25 #endif
44 #include "defun.h" 27 #include "defun.h"
45 #include "error.h" 28 #include "error.h"
46 #include "lo-ieee.h" 29 #include "lo-ieee.h"
47 30
48 static void 31 static void
49 gripe_ellipj_arg ( const char *arg) 32 gripe_ellipj_arg (const char *arg)
50 { 33 {
51 error ("ellipj: expecting scalar or matrix as %s argument", arg); 34 error ("ellipj: expecting scalar or matrix as %s argument", arg);
52 } 35 }
53 36
54 const double eps = 2.220446049e-16;
55 const int Nmax = 16;
56
57 static void 37 static void
58 sncndn ( double u, double m, double& sn, double& cn, double& dn, double& err) { 38 sncndn (double u, double m, double& sn, double& cn, double& dn, double& err)
59 /* real */ 39 {
60 double sqrt_eps, m1, t=0., si_u, co_u, se_u, ta_u, b, c[Nmax], a[Nmax], phi; 40 static const int Nmax = 16;
61 int n, Nn, ii; 41 double m1, t=0, si_u, co_u, se_u, ta_u, b, c[Nmax], a[Nmax], phi;
62 42 int n, Nn, ii;
63 if (m < 0. || m > 1.) { 43
64 warning ("ellipj: expecting 0. <= m <= 1."); /* -lc- */ 44 if (m < 0 || m > 1)
65 sn = cn = dn = lo_ieee_nan_value (); 45 {
66 return; 46 warning ("ellipj: expecting 0 <= m <= 1"); /* -lc- */
67 } 47 sn = cn = dn = lo_ieee_nan_value ();
68 sqrt_eps = sqrt(eps); 48 return;
69 if (m < sqrt_eps) { 49 }
70 /* # For small m, ( Abramowitz and Stegun, Section 16.13 ) */ 50
71 /*{{{*/ 51 double sqrt_eps = sqrt (std::numeric_limits<double>::epsilon ());
72 si_u = sin(u); 52 if (m < sqrt_eps)
73 co_u = cos(u); 53 {
74 t = 0.25*m*(u-si_u*co_u); 54 /* # For small m, ( Abramowitz and Stegun, Section 16.13 ) */
75 sn = si_u - t * co_u; 55 si_u = sin (u);
76 cn = co_u + t * si_u; 56 co_u = cos (u);
77 dn = 1.0 - 0.5*m*si_u*si_u; 57 t = 0.25*m*(u - si_u*co_u);
78 /*}}}*/ 58 sn = si_u - t * co_u;
79 } else if ( (1.0 - m) < sqrt_eps ) { 59 cn = co_u + t * si_u;
80 /* For m1 = (1-m) small ( Abramowitz and Stegun, Section 16.15 ) */ 60 dn = 1 - 0.5*m*si_u*si_u;
81 /*{{{*/ 61 }
82 m1 = 1.0-m; 62 else if ((1 - m) < sqrt_eps)
83 si_u = sinh(u); 63 {
84 co_u = cosh(u); 64 /* For m1 = (1-m) small ( Abramowitz and Stegun, Section 16.15 ) */
85 ta_u = tanh(u); 65 m1 = 1 - m;
86 se_u = 1.0/co_u; 66 si_u = sinh (u);
87 sn = ta_u + 0.25*m1*(si_u*co_u-u)*se_u*se_u; 67 co_u = cosh (u);
88 cn = se_u - 0.25*m1*(si_u*co_u-u)*ta_u*se_u; 68 ta_u = tanh (u);
89 dn = se_u + 0.25*m1*(si_u*co_u+u)*ta_u*se_u; 69 se_u = 1/co_u;
90 /*}}}*/ 70 sn = ta_u + 0.25*m1*(si_u*co_u - u)*se_u*se_u;
91 } else { 71 cn = se_u - 0.25*m1*(si_u*co_u - u)*ta_u*se_u;
92 /*{{{*/ 72 dn = se_u + 0.25*m1*(si_u*co_u + u)*ta_u*se_u;
93 /* 73 }
94 // Arithmetic-Geometric Mean (AGM) algorithm 74 else
95 // ( Abramowitz and Stegun, Section 16.4 ) 75 {
96 */ 76 /*
97 77 // Arithmetic-Geometric Mean (AGM) algorithm
98 a[0] = 1.0; 78 // ( Abramowitz and Stegun, Section 16.4 )
99 b = sqrt(1.0-m); 79 */
100 c[0] = sqrt(m); 80
101 for (n = 1; n<Nmax; ++n) { 81 a[0] = 1;
102 a[n] = (a[n-1]+b)/2; 82 b = sqrt (1 - m);
103 c[n] = (a[n-1]-b)/2; 83 c[0] = sqrt (m);
104 b = sqrt(a[n-1]*b); 84 for (n = 1; n < Nmax; ++n)
105 if ( c[n]/a[n] < eps) break; 85 {
86 a[n] = (a[n - 1] + b)/2;
87 c[n] = (a[n - 1] - b)/2;
88 b = sqrt (a[n - 1]*b);
89 if (c[n]/a[n] < std::numeric_limits<double>::epsilon ()) break;
106 } 90 }
107 if ( n >= Nmax-1) { 91 if (n >= Nmax - 1)
108 // fprintf(stderr, "Not enough workspace\n"); 92 {
109 err = 1.; 93 err = 1;
110 return; 94 return;
111 } 95 }
112 Nn = n; 96 Nn = n;
113 for ( ii = 1; n>0; ii = ii*2, --n) ; // pow(2, Nn) 97 for (ii = 1; n > 0; ii = ii*2, --n) ; // ii = pow(2,Nn)
114 phi = ii*a[Nn]*u; 98 phi = ii*a[Nn]*u;
115 for ( n = Nn; n > 0; --n) { 99 for (n = Nn; n > 0; --n)
100 {
116 t = phi; 101 t = phi;
117 phi = (asin((c[n]/a[n])* sin(phi))+phi)/2.; 102 phi = (asin ((c[n]/a[n])* sin (phi)) + phi)/2;
118 } 103 }
119 sn = sin(phi); 104 sn = sin (phi);
120 cn = cos(phi); 105 cn = cos (phi);
121 dn = cn/cos(t-phi); 106 dn = cn/cos (t - phi);
122 /*}}}*/ 107 }
123 }
124 return;
125 } 108 }
126 109
127 static void 110 static void
128 sncndn ( Complex& u, double m, 111 sncndn (Complex& u, double m, Complex& sn, Complex& cn, Complex& dn,
129 Complex& sn, Complex& cn, Complex& dn, double& err) { 112 double& err)
130 double m1 = 1.-m, ss1, cc1, dd1; 113 {
131 114 double m1 = 1 - m, ss1, cc1, dd1;
132 sncndn( imag(u), m1, ss1, cc1, dd1, err); 115
133 if ( real(u) == 0.) { 116 sncndn (imag (u), m1, ss1, cc1, dd1, err);
134 /* u is pure imag: Jacoby imag. transf. */ 117 if (real (u) == 0)
135 /*{{{*/ 118 {
136 sn = Complex (0. , ss1/cc1); 119 /* u is pure imag: Jacoby imag. transf. */
137 cn = 1/cc1; // cn.imag = 0.; 120 sn = Complex (0, ss1/cc1);
138 dn = dd1/cc1; // dn.imag = 0.; 121 cn = 1/cc1; // cn.imag = 0;
139 /*}}}*/ 122 dn = dd1/cc1; // dn.imag = 0;
140 } else { 123 }
141 /* u is generic complex */ 124 else
142 /*{{{*/ 125 {
143 double ss, cc, dd, ddd; 126 /* u is generic complex */
144 127 double ss, cc, dd, ddd;
145 sncndn( real(u), m, ss, cc, dd, err); 128
129 sncndn (real (u), m, ss, cc, dd, err);
146 ddd = cc1*cc1 + m*ss*ss*ss1*ss1; 130 ddd = cc1*cc1 + m*ss*ss*ss1*ss1;
147 sn = Complex (ss*dd1/ddd, cc*dd*ss1*cc1/ddd); 131 sn = Complex (ss*dd1/ddd, cc*dd*ss1*cc1/ddd);
148 cn = Complex (cc*cc1/ddd, -ss*dd*ss1*dd1/ddd); 132 cn = Complex (cc*cc1/ddd, -ss*dd*ss1*dd1/ddd);
149 dn = Complex (dd*cc1*dd1/ddd, -m*ss*cc*ss1/ddd); 133 dn = Complex (dd*cc1*dd1/ddd, -m*ss*cc*ss1/ddd);
150 /*}}}*/ 134 }
151 }
152 return;
153 } 135 }
154 136
155 DEFUN (ellipj, args, nargout, 137 DEFUN (ellipj, args, nargout,
156 "-*- texinfo -*-\n\ 138 "-*- texinfo -*-\n\
157 @deftypefn {Loadable Function} {[@var{sn}, @var{cn}, @var{dn}] =} \ 139 @deftypefn {Built-in Function} {[@var{sn}, @var{cn}, @var{dn}] =} \
158 ellipj (@var{u}, @var{m}, @var{err})\n\ 140 ellipj (@var{u}, @var{m}, @var{err})\n\
159 Compute the Jacobi elliptic functions sn, cn, dn of complex argument and real parameter.\n\ 141 Compute the Jacobi elliptic functions sn, cn, dn of complex argument and real parameter.\n\
160 \n\ 142 \n\
161 If @var{m} is a scalar, the results are the same size as @var{u}.\n\ 143 If @var{m} is a scalar, the results are the same size as @var{u}.\n\
162 If @var{u} is a scalar, the results are the same size as @var{m}.\n\ 144 If @var{u} is a scalar, the results are the same size as @var{m}.\n\
164 results are matrices with @code{length (@var{u})} rows and\n\ 146 results are matrices with @code{length (@var{u})} rows and\n\
165 @code{length (@var{m})} columns. Otherwise, @var{u} and\n\ 147 @code{length (@var{m})} columns. Otherwise, @var{u} and\n\
166 @var{m} must conform and the results will be the same size.\n\ 148 @var{m} must conform and the results will be the same size.\n\
167 \n\ 149 \n\
168 The value of @var{u} may be complex.\n\ 150 The value of @var{u} may be complex.\n\
169 The value of @var{m} must be 0 <= m <= 1. .\n\ 151 The value of @var{m} must be 0 <= m <= 1.\n\
170 \n\ 152 \n\
171 If requested, @var{err} contains the following status information\n\ 153 If requested, @var{err} contains the following status information\n\
172 and is the same size as the result.\n\ 154 and is the same size as the result.\n\
173 \n\ 155 \n\
174 @enumerate 0\n\ 156 @enumerate 0\n\
176 Normal return.\n\ 158 Normal return.\n\
177 @item\n\ 159 @item\n\
178 Error---no computation, algorithm termination condition not met,\n\ 160 Error---no computation, algorithm termination condition not met,\n\
179 return @code{NaN}.\n\ 161 return @code{NaN}.\n\
180 @end enumerate\n\ 162 @end enumerate\n\
163 Ref: Abramowitz, Milton and Stegun, Irene A\n\
164 Handbook of Mathematical Functions, Dover, 1965\n\
165 Chapter 16 (Sections 16.4, 16.13 and 16.15)\n\
181 @end deftypefn") 166 @end deftypefn")
182 { 167 {
183 octave_value_list retval; 168 octave_value_list retval;
184 169
185 int nargin = args.length (); 170 int nargin = args.length ();
186 171
187 if (nargin == 2 ) { 172 if (nargin != 2 )
188 octave_value u_arg = args(0); 173 {
189 octave_value m_arg = args(1); 174 print_usage ();
190 175 return retval;
191 if (m_arg.is_scalar_type ()) { // m is scalar 176 }
192 double m = args(1).double_value (); 177
193 178 octave_value u_arg = args(0);
194 if (! error_state) { 179 octave_value m_arg = args(1);
195 180
196 if (u_arg.is_scalar_type ()) { /* u scalar */ 181 if (m_arg.is_scalar_type ())
197 /*{{{*/ 182 {
198 if (u_arg.is_real_type ()) { // u real 183 double m = args(1).double_value ();
199 double u = args(0).double_value (); 184
200 185 if (error_state)
201 if (! error_state) { 186 {
202 double sn, cn, dn; 187 gripe_ellipj_arg ("second");
203 double err=0; 188 return retval;
204 octave_value result; 189 }
205 190
206 sncndn(u, m, sn, cn, dn, err); 191 if (u_arg.is_scalar_type ())
207 retval (0) = sn; 192 {
208 retval (1) = cn; 193 if (u_arg.is_real_type ())
209 retval (2) = dn; 194 { // u real
210 if (nargout > 3) 195 double u = args(0).double_value ();
211 retval(3) = err; 196
212 } else 197 if (error_state)
213 gripe_ellipj_arg ( "first"); 198 {
214 199 gripe_ellipj_arg ("first");
215 } else { // u complex 200 return retval;
201 }
202 double sn, cn, dn;
203 double err = 0;
204
205 sncndn (u, m, sn, cn, dn, err);
206 retval (0) = sn;
207 retval (1) = cn;
208 retval (2) = dn;
209 if (nargout > 3) retval(3) = err;
210 }
211 else
212 { // u complex
216 Complex u = u_arg.complex_value (); 213 Complex u = u_arg.complex_value ();
217 214
218 if (! error_state) { 215 if (error_state)
219 Complex sn, cn, dn; 216 {
220 double err; 217 gripe_ellipj_arg ("second");
221 octave_value result; 218 return retval;
222 219 }
223 sncndn( u, m, sn, cn, dn, err); 220
224 221 Complex sn, cn, dn;
225 retval (0) = sn; 222 double err;
226 retval (1) = cn; 223
227 retval (2) = dn; 224 sncndn (u, m, sn, cn, dn, err);
228 if (nargout > 3) retval(3) = err; 225
229 } else 226 retval (0) = sn;
230 gripe_ellipj_arg ( "second"); 227 retval (1) = cn;
228 retval (2) = dn;
229 if (nargout > 3) retval(3) = err;
231 } 230 }
232 /*}}}*/ 231 }
233 } else { /* u is matrix ( m is scalar ) */ 232 else
234 /*{{{*/ 233 { /* u is matrix ( m is scalar ) */
235 ComplexMatrix u = u_arg.complex_matrix_value (); 234 ComplexMatrix u = u_arg.complex_matrix_value ();
236 235
237 if (! error_state) { 236 if (error_state)
238 octave_value result; 237 {
239 int nr = u.rows (); 238 gripe_ellipj_arg ("first");
240 int nc = u.cols (); 239 return retval;
240 }
241
242 octave_idx_type nr = u.rows ();
243 octave_idx_type nc = u.cols ();
244
245 ComplexMatrix sn (nr, nc), cn (nr, nc), dn (nr, nc);
246 Matrix err (nr, nc);
247
248 for (octave_idx_type j = 0; j < nc; j++)
249 for (octave_idx_type i = 0; i < nr; i++)
250 sncndn (u(i,j), m, sn(i,j), cn(i,j), dn(i,j), err(i,j));
251
252 retval (0) = sn;
253 retval (1) = cn;
254 retval (2) = dn;
255 if (nargout > 3) retval(3) = err;
256 }
257 }
258 else
259 {
260 Matrix m = args(1).matrix_value ();
261
262 if (error_state)
263 {
264 gripe_ellipj_arg ("second");
265 return retval;
266 }
267
268 octave_idx_type mr = m.rows ();
269 octave_idx_type mc = m.cols ();
270
271 if (u_arg.is_scalar_type ())
272 { /* u is scalar */
273 octave_idx_type nr = m.rows ();
274 octave_idx_type nc = m.cols ();
275 Matrix err (nr, nc);
276
277 if (u_arg.is_real_type ())
278 {
279 double u = u_arg.double_value ();
280
281 if (error_state)
282 {
283 gripe_ellipj_arg ("first");
284 return retval;
285 }
286
287 Matrix sn (nr, nc), cn (nr, nc), dn (nr, nc);
288 for (octave_idx_type j = 0; j < nc; j++)
289 for (octave_idx_type i = 0; i < nr; i++)
290 sncndn (u, m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
291
292 retval (0) = sn;
293 retval (1) = cn;
294 retval (2) = dn;
295 if (nargout > 3) retval(3) = err;
296 }
297 else
298 {
299 Complex u = u_arg.complex_value ();
300 if (error_state)
301 {
302 gripe_ellipj_arg ("first");
303 return retval;
304 }
241 305
242 ComplexMatrix sn (nr, nc), cn (nr, nc), dn (nr, nc); 306 ComplexMatrix sn (nr, nc), cn (nr, nc), dn (nr, nc);
243 Matrix err (nr, nc); 307 for (octave_idx_type j = 0; j < nc; j++)
244 308 for (octave_idx_type i = 0; i < nr; i++)
245 for (int j = 0; j < nc; j++) 309 sncndn (u, m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
246 for (int i = 0; i < nr; i++) 310 retval (0) = sn;
247 sncndn (u(i,j), m, sn(i,j), cn(i,j), dn(i,j), err(i,j)); 311 retval (1) = cn;
248 312 retval (2) = dn;
249 retval (0) = sn; 313 if (nargout > 3) retval(3) = err;
250 retval (1) = cn; 314 }
251 retval (2) = dn; 315 }
252 if (nargout > 3) retval(3) = err; 316 else
253 } else 317 { // u is matrix (m is matrix)
254 gripe_ellipj_arg ( "first"); 318 if (u_arg.is_real_type ())
255 /*}}}*/ 319 { // u real matrix
256 }
257 } else
258 gripe_ellipj_arg ( "second");
259 } else { // m is matrix
260 Matrix m = args(1).matrix_value ();
261
262 if (! error_state) {
263 int mr = m.rows ();
264 int mc = m.cols ();
265
266 if (u_arg.is_scalar_type ()) { /* u is scalar */
267 /*{{{*/
268 octave_value result;
269 int nr = m.rows ();
270 int nc = m.cols ();
271 Matrix err (nr, nc);
272
273 if (u_arg.is_real_type ()) {
274 double u = u_arg.double_value ();
275 Matrix sn (nr, nc), cn (nr, nc), dn (nr, nc);
276 if (! error_state) {
277 for (int j = 0; j < nc; j++)
278 for (int i = 0; i < nr; i++)
279 sncndn (u, m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
280
281 retval (0) = sn;
282 retval (1) = cn;
283 retval (2) = dn;
284 if (nargout > 3) retval(3) = err;
285 } else
286 gripe_ellipj_arg ( "first");
287 } else {
288 Complex u = u_arg.complex_value ();
289 ComplexMatrix sn (nr, nc), cn (nr, nc), dn (nr, nc);
290 if (! error_state) {
291 for (int j = 0; j < nc; j++)
292 for (int i = 0; i < nr; i++)
293 sncndn (u, m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
294 retval (0) = sn;
295 retval (1) = cn;
296 retval (2) = dn;
297 if (nargout > 3) retval(3) = err;
298 } else
299 gripe_ellipj_arg ( "first");
300 }
301 /*}}}*/
302 } else { // u is matrix (m is matrix)
303 /*{{{*/
304 if (u_arg.is_real_type ()) { // u real matrix
305 320
306 Matrix u = u_arg.matrix_value (); 321 Matrix u = u_arg.matrix_value ();
307 if (! error_state) { 322 if (error_state)
308 int ur = u.rows (); 323 {
309 int uc = u.cols (); 324 gripe_ellipj_arg ("first ");
310 325 return retval;
311 if (mr == 1 && uc == 1) { // u column, m row 326 }
312 RowVector rm = m.row ((octave_idx_type)0); 327
313 ColumnVector cu = u.column ((octave_idx_type)0); 328 octave_idx_type ur = u.rows ();
314 329 octave_idx_type uc = u.cols ();
315 Matrix sn (ur, mc), cn (ur, mc), dn (ur, mc); 330
316 Matrix err(ur,mc); 331 if (mr == 1 && uc == 1)
317 // octave_value result; 332 { // u column, m row
318 333 RowVector rm = m.row (0);
319 for (int j = 0; j < mc; j++) 334 ColumnVector cu = u.column (0);
320 for (int i = 0; i < ur; i++) 335
321 sncndn (cu(i), rm(j), sn(i,j), cn(i,j), dn(i,j), err(i,j)); 336 Matrix sn (ur, mc), cn (ur, mc), dn (ur, mc);
322 337 Matrix err (ur,mc);
323 retval (0) = sn; 338
324 retval (1) = cn; 339 for (octave_idx_type j = 0; j < mc; j++)
325 retval (2) = dn; 340 for (octave_idx_type i = 0; i < ur; i++)
326 if (nargout > 3) retval(3) = err; 341 sncndn (cu(i), rm(j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
327 } else if (ur == mr && uc == mc) { 342
328 Matrix sn (ur, mc), cn (ur, mc), dn (ur, mc); 343 retval (0) = sn;
329 Matrix err(ur,mc); 344 retval (1) = cn;
330 // octave_value result; 345 retval (2) = dn;
331 346 if (nargout > 3) retval(3) = err;
332 for (int j = 0; j < uc; j++) 347 }
333 for (int i = 0; i < ur; i++) 348 else if (ur == mr && uc == mc)
334 sncndn (u(i,j), m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j)); 349 {
335 350 Matrix sn (ur, mc), cn (ur, mc), dn (ur, mc);
336 retval (0) = sn; 351 Matrix err (ur,mc);
337 retval (1) = cn; 352
338 retval (2) = dn; 353 for (octave_idx_type j = 0; j < uc; j++)
339 if (nargout > 3) retval(3) = err; 354 for (octave_idx_type i = 0; i < ur; i++)
340 } else 355 sncndn (u(i,j), m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
341 error("u m invalid"); 356
342 } else 357 retval (0) = sn;
343 gripe_ellipj_arg ( "first "); 358 retval (1) = cn;
344 } else { // u complex matrix 359 retval (2) = dn;
360 if (nargout > 3) retval(3) = err;
361 }
362 else
363 error ("u m invalid");
364 }
365 else
366 { // u complex matrix
345 ComplexMatrix u = u_arg.complex_matrix_value (); 367 ComplexMatrix u = u_arg.complex_matrix_value ();
346 if (! error_state) { 368 if (error_state)
347 int ur = u.rows (); 369 {
348 int uc = u.cols (); 370 gripe_ellipj_arg ("second");
349 371 return retval;
350 if (mr == 1 && uc == 1) { 372 }
351 RowVector rm = m.row ((octave_idx_type)0); 373
352 ComplexColumnVector cu = u.column ((octave_idx_type)0); 374 octave_idx_type ur = u.rows ();
353 375 octave_idx_type uc = u.cols ();
354 ComplexMatrix sn (ur, mc), cn (ur, mc), dn (ur, mc); 376
355 Matrix err(ur,mc); 377 if (mr == 1 && uc == 1)
356 // octave_value result; 378 {
357 379 RowVector rm = m.row (0);
358 for (int j = 0; j < mc; j++) 380 ComplexColumnVector cu = u.column (0);
359 for (int i = 0; i < ur; i++) 381
360 sncndn (cu(i), rm(j), sn(i,j), cn(i,j), dn(i,j), err(i,j)); 382 ComplexMatrix sn (ur, mc), cn (ur, mc), dn (ur, mc);
361 383 Matrix err (ur,mc);
362 retval (0) = sn; 384
363 retval (1) = cn; 385 for (octave_idx_type j = 0; j < mc; j++)
364 retval (2) = dn; 386 for (octave_idx_type i = 0; i < ur; i++)
365 if (nargout > 3) retval(3) = err; 387 sncndn (cu(i), rm(j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
366 } else if (ur == mr && uc == mc) { 388
367 389 retval (0) = sn;
368 ComplexMatrix sn (ur, mc), cn (ur, mc), dn (ur, mc); 390 retval (1) = cn;
369 Matrix err(ur,mc); 391 retval (2) = dn;
370 // octave_value result; 392 if (nargout > 3) retval(3) = err;
371 393 }
372 for (int j = 0; j < uc; j++) 394 else if (ur == mr && uc == mc)
373 for (int i = 0; i < ur; i++) 395 {
374 sncndn (u(i,j), m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j)); 396 ComplexMatrix sn (ur, mc), cn (ur, mc), dn (ur, mc);
375 397 Matrix err (ur,mc);
376 retval (0) = sn; 398
377 retval (1) = cn; 399 for (octave_idx_type j = 0; j < uc; j++)
378 retval (2) = dn; 400 for (octave_idx_type i = 0; i < ur; i++)
379 if (nargout > 3) retval(3) = err; 401 sncndn (u(i,j), m(i,j), sn(i,j), cn(i,j), dn(i,j), err(i,j));
380 } else 402
381 error("u m invalid"); 403 retval (0) = sn;
382 } else 404 retval (1) = cn;
383 gripe_ellipj_arg ( "second"); 405 retval (2) = dn;
406 if (nargout > 3) retval(3) = err;
407 }
408 else
409 error ("u m invalid");
384 } 410 }
385 /*}}}*/ 411 }
386 } 412 } // m matrix
387 } else 413
388 gripe_ellipj_arg ( "second"); 414 return retval;
389 } // m matrix
390 } else // wrong n. of argin
391 print_usage ();
392 return retval;
393 } 415 }
394 416
395 /* 417 /*
396 ## demos taken from inst/ellipj.m 418 ## demos taken from inst/ellipj.m
397 419
442 /* 464 /*
443 ## tests taken from inst/test_sncndn.m 465 ## tests taken from inst/test_sncndn.m
444 466
445 %!test 467 %!test
446 %! k = (tan(pi/8.))^2; m = k*k; 468 %! k = (tan(pi/8.))^2; m = k*k;
447 %! SN = [ 469 %! SN = [
448 %! -1. + I * 0. , -0.8392965923 + 0. * I 470 %! -1. + I * 0. , -0.8392965923 + 0. * I
449 %! -1. + I * 0.2 , -0.8559363407 + 0.108250955 * I 471 %! -1. + I * 0.2 , -0.8559363407 + 0.108250955 * I
450 %! -1. + I * 0.4 , -0.906529758 + 0.2204040232 * I 472 %! -1. + I * 0.4 , -0.906529758 + 0.2204040232 * I
451 %! -1. + I * 0.6 , -0.9931306727 + 0.3403783409 * I 473 %! -1. + I * 0.6 , -0.9931306727 + 0.3403783409 * I
452 %! -1. + I * 0.8 , -1.119268095 + 0.4720784944 * I 474 %! -1. + I * 0.8 , -1.119268095 + 0.4720784944 * I