[2563] | 1 | /*************************************************************************
|
---|
| 2 | Cephes Math Library Release 2.8: June, 2000
|
---|
| 3 | Copyright 1984, 1987, 2000 by Stephen L. Moshier
|
---|
| 4 |
|
---|
| 5 | Contributors:
|
---|
| 6 | * Sergey Bochkanov (ALGLIB project). Translation from C to
|
---|
| 7 | pseudocode.
|
---|
| 8 |
|
---|
| 9 | See subroutines comments for additional copyrights.
|
---|
| 10 |
|
---|
| 11 | >>> SOURCE LICENSE >>>
|
---|
| 12 | This program is free software; you can redistribute it and/or modify
|
---|
| 13 | it under the terms of the GNU General Public License as published by
|
---|
| 14 | the Free Software Foundation (www.fsf.org); either version 2 of the
|
---|
| 15 | License, or (at your option) any later version.
|
---|
| 16 |
|
---|
| 17 | This program is distributed in the hope that it will be useful,
|
---|
| 18 | but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 19 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 20 | GNU General Public License for more details.
|
---|
| 21 |
|
---|
| 22 | A copy of the GNU General Public License is available at
|
---|
| 23 | http://www.fsf.org/licensing/licenses
|
---|
| 24 |
|
---|
| 25 | >>> END OF LICENSE >>>
|
---|
| 26 | *************************************************************************/
|
---|
| 27 |
|
---|
| 28 | using System;
|
---|
| 29 |
|
---|
| 30 | namespace alglib
|
---|
| 31 | {
|
---|
| 32 | public class jacobianelliptic
|
---|
| 33 | {
|
---|
| 34 | /*************************************************************************
|
---|
| 35 | Jacobian Elliptic Functions
|
---|
| 36 |
|
---|
| 37 | Evaluates the Jacobian elliptic functions sn(u|m), cn(u|m),
|
---|
| 38 | and dn(u|m) of parameter m between 0 and 1, and real
|
---|
| 39 | argument u.
|
---|
| 40 |
|
---|
| 41 | These functions are periodic, with quarter-period on the
|
---|
| 42 | real axis equal to the complete elliptic integral
|
---|
| 43 | ellpk(1.0-m).
|
---|
| 44 |
|
---|
| 45 | Relation to incomplete elliptic integral:
|
---|
| 46 | If u = ellik(phi,m), then sn(u|m) = sin(phi),
|
---|
| 47 | and cn(u|m) = cos(phi). Phi is called the amplitude of u.
|
---|
| 48 |
|
---|
| 49 | Computation is by means of the arithmetic-geometric mean
|
---|
| 50 | algorithm, except when m is within 1e-9 of 0 or 1. In the
|
---|
| 51 | latter case with m close to 1, the approximation applies
|
---|
| 52 | only for phi < pi/2.
|
---|
| 53 |
|
---|
| 54 | ACCURACY:
|
---|
| 55 |
|
---|
| 56 | Tested at random points with u between 0 and 10, m between
|
---|
| 57 | 0 and 1.
|
---|
| 58 |
|
---|
| 59 | Absolute error (* = relative error):
|
---|
| 60 | arithmetic function # trials peak rms
|
---|
| 61 | IEEE phi 10000 9.2e-16* 1.4e-16*
|
---|
| 62 | IEEE sn 50000 4.1e-15 4.6e-16
|
---|
| 63 | IEEE cn 40000 3.6e-15 4.4e-16
|
---|
| 64 | IEEE dn 10000 1.3e-12 1.8e-14
|
---|
| 65 |
|
---|
| 66 | Peak error observed in consistency check using addition
|
---|
| 67 | theorem for sn(u+v) was 4e-16 (absolute). Also tested by
|
---|
| 68 | the above relation to the incomplete elliptic integral.
|
---|
| 69 | Accuracy deteriorates when u is large.
|
---|
| 70 |
|
---|
| 71 | Cephes Math Library Release 2.8: June, 2000
|
---|
| 72 | Copyright 1984, 1987, 2000 by Stephen L. Moshier
|
---|
| 73 | *************************************************************************/
|
---|
| 74 | public static void jacobianellipticfunctions(double u,
|
---|
| 75 | double m,
|
---|
| 76 | ref double sn,
|
---|
| 77 | ref double cn,
|
---|
| 78 | ref double dn,
|
---|
| 79 | ref double ph)
|
---|
| 80 | {
|
---|
| 81 | double ai = 0;
|
---|
| 82 | double b = 0;
|
---|
| 83 | double phi = 0;
|
---|
| 84 | double t = 0;
|
---|
| 85 | double twon = 0;
|
---|
| 86 | double[] a = new double[0];
|
---|
| 87 | double[] c = new double[0];
|
---|
| 88 | int i = 0;
|
---|
| 89 |
|
---|
| 90 | System.Diagnostics.Debug.Assert((double)(m)>=(double)(0) & (double)(m)<=(double)(1), "Domain error in JacobianEllipticFunctions: m<0 or m>1");
|
---|
| 91 | a = new double[8+1];
|
---|
| 92 | c = new double[8+1];
|
---|
| 93 | if( (double)(m)<(double)(1.0e-9) )
|
---|
| 94 | {
|
---|
| 95 | t = Math.Sin(u);
|
---|
| 96 | b = Math.Cos(u);
|
---|
| 97 | ai = 0.25*m*(u-t*b);
|
---|
| 98 | sn = t-ai*b;
|
---|
| 99 | cn = b+ai*t;
|
---|
| 100 | ph = u-ai;
|
---|
| 101 | dn = 1.0-0.5*m*t*t;
|
---|
| 102 | return;
|
---|
| 103 | }
|
---|
| 104 | if( (double)(m)>=(double)(0.9999999999) )
|
---|
| 105 | {
|
---|
| 106 | ai = 0.25*(1.0-m);
|
---|
| 107 | b = Math.Cosh(u);
|
---|
| 108 | t = Math.Tanh(u);
|
---|
| 109 | phi = 1.0/b;
|
---|
| 110 | twon = b*Math.Sinh(u);
|
---|
| 111 | sn = t+ai*(twon-u)/(b*b);
|
---|
| 112 | ph = 2.0*Math.Atan(Math.Exp(u))-1.57079632679489661923+ai*(twon-u)/b;
|
---|
| 113 | ai = ai*t*phi;
|
---|
| 114 | cn = phi-ai*(twon-u);
|
---|
| 115 | dn = phi+ai*(twon+u);
|
---|
| 116 | return;
|
---|
| 117 | }
|
---|
| 118 | a[0] = 1.0;
|
---|
| 119 | b = Math.Sqrt(1.0-m);
|
---|
| 120 | c[0] = Math.Sqrt(m);
|
---|
| 121 | twon = 1.0;
|
---|
| 122 | i = 0;
|
---|
| 123 | while( (double)(Math.Abs(c[i]/a[i]))>(double)(AP.Math.MachineEpsilon) )
|
---|
| 124 | {
|
---|
| 125 | if( i>7 )
|
---|
| 126 | {
|
---|
| 127 | System.Diagnostics.Debug.Assert(false, "Overflow in JacobianEllipticFunctions");
|
---|
| 128 | break;
|
---|
| 129 | }
|
---|
| 130 | ai = a[i];
|
---|
| 131 | i = i+1;
|
---|
| 132 | c[i] = 0.5*(ai-b);
|
---|
| 133 | t = Math.Sqrt(ai*b);
|
---|
| 134 | a[i] = 0.5*(ai+b);
|
---|
| 135 | b = t;
|
---|
| 136 | twon = twon*2.0;
|
---|
| 137 | }
|
---|
| 138 | phi = twon*a[i]*u;
|
---|
| 139 | do
|
---|
| 140 | {
|
---|
| 141 | t = c[i]*Math.Sin(phi)/a[i];
|
---|
| 142 | b = phi;
|
---|
| 143 | phi = (Math.Asin(t)+phi)/2.0;
|
---|
| 144 | i = i-1;
|
---|
| 145 | }
|
---|
| 146 | while( i!=0 );
|
---|
| 147 | sn = Math.Sin(phi);
|
---|
| 148 | t = Math.Cos(phi);
|
---|
| 149 | cn = t;
|
---|
| 150 | dn = t/Math.Cos(phi-b);
|
---|
| 151 | ph = phi;
|
---|
| 152 | }
|
---|
| 153 | }
|
---|
| 154 | }
|
---|