[2563] | 1 | /*************************************************************************
|
---|
| 2 | Cephes Math Library Release 2.8: June, 2000
|
---|
| 3 | Copyright 1984, 1987, 1995, 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 elliptic
|
---|
| 33 | {
|
---|
| 34 | /*************************************************************************
|
---|
| 35 | Complete elliptic integral of the first kind
|
---|
| 36 |
|
---|
| 37 | Approximates the integral
|
---|
| 38 |
|
---|
| 39 |
|
---|
| 40 |
|
---|
| 41 | pi/2
|
---|
| 42 | -
|
---|
| 43 | | |
|
---|
| 44 | | dt
|
---|
| 45 | K(m) = | ------------------
|
---|
| 46 | | 2
|
---|
| 47 | | | sqrt( 1 - m sin t )
|
---|
| 48 | -
|
---|
| 49 | 0
|
---|
| 50 |
|
---|
| 51 | using the approximation
|
---|
| 52 |
|
---|
| 53 | P(x) - log x Q(x).
|
---|
| 54 |
|
---|
| 55 | ACCURACY:
|
---|
| 56 |
|
---|
| 57 | Relative error:
|
---|
| 58 | arithmetic domain # trials peak rms
|
---|
| 59 | IEEE 0,1 30000 2.5e-16 6.8e-17
|
---|
| 60 |
|
---|
| 61 | Cephes Math Library, Release 2.8: June, 2000
|
---|
| 62 | Copyright 1984, 1987, 2000 by Stephen L. Moshier
|
---|
| 63 | *************************************************************************/
|
---|
| 64 | public static double ellipticintegralk(double m)
|
---|
| 65 | {
|
---|
| 66 | double result = 0;
|
---|
| 67 |
|
---|
| 68 | result = ellipticintegralkhighprecision(1.0-m);
|
---|
| 69 | return result;
|
---|
| 70 | }
|
---|
| 71 |
|
---|
| 72 |
|
---|
| 73 | /*************************************************************************
|
---|
| 74 | Complete elliptic integral of the first kind
|
---|
| 75 |
|
---|
| 76 | Approximates the integral
|
---|
| 77 |
|
---|
| 78 |
|
---|
| 79 |
|
---|
| 80 | pi/2
|
---|
| 81 | -
|
---|
| 82 | | |
|
---|
| 83 | | dt
|
---|
| 84 | K(m) = | ------------------
|
---|
| 85 | | 2
|
---|
| 86 | | | sqrt( 1 - m sin t )
|
---|
| 87 | -
|
---|
| 88 | 0
|
---|
| 89 |
|
---|
| 90 | where m = 1 - m1, using the approximation
|
---|
| 91 |
|
---|
| 92 | P(x) - log x Q(x).
|
---|
| 93 |
|
---|
| 94 | The argument m1 is used rather than m so that the logarithmic
|
---|
| 95 | singularity at m = 1 will be shifted to the origin; this
|
---|
| 96 | preserves maximum accuracy.
|
---|
| 97 |
|
---|
| 98 | K(0) = pi/2.
|
---|
| 99 |
|
---|
| 100 | ACCURACY:
|
---|
| 101 |
|
---|
| 102 | Relative error:
|
---|
| 103 | arithmetic domain # trials peak rms
|
---|
| 104 | IEEE 0,1 30000 2.5e-16 6.8e-17
|
---|
| 105 |
|
---|
| 106 | Àëãîðèòì âçÿò èç áèáëèîòåêè Cephes
|
---|
| 107 | *************************************************************************/
|
---|
| 108 | public static double ellipticintegralkhighprecision(double m1)
|
---|
| 109 | {
|
---|
| 110 | double result = 0;
|
---|
| 111 | double p = 0;
|
---|
| 112 | double q = 0;
|
---|
| 113 |
|
---|
| 114 | if( (double)(m1)<=(double)(AP.Math.MachineEpsilon) )
|
---|
| 115 | {
|
---|
| 116 | result = 1.3862943611198906188E0-0.5*Math.Log(m1);
|
---|
| 117 | }
|
---|
| 118 | else
|
---|
| 119 | {
|
---|
| 120 | p = 1.37982864606273237150E-4;
|
---|
| 121 | p = p*m1+2.28025724005875567385E-3;
|
---|
| 122 | p = p*m1+7.97404013220415179367E-3;
|
---|
| 123 | p = p*m1+9.85821379021226008714E-3;
|
---|
| 124 | p = p*m1+6.87489687449949877925E-3;
|
---|
| 125 | p = p*m1+6.18901033637687613229E-3;
|
---|
| 126 | p = p*m1+8.79078273952743772254E-3;
|
---|
| 127 | p = p*m1+1.49380448916805252718E-2;
|
---|
| 128 | p = p*m1+3.08851465246711995998E-2;
|
---|
| 129 | p = p*m1+9.65735902811690126535E-2;
|
---|
| 130 | p = p*m1+1.38629436111989062502E0;
|
---|
| 131 | q = 2.94078955048598507511E-5;
|
---|
| 132 | q = q*m1+9.14184723865917226571E-4;
|
---|
| 133 | q = q*m1+5.94058303753167793257E-3;
|
---|
| 134 | q = q*m1+1.54850516649762399335E-2;
|
---|
| 135 | q = q*m1+2.39089602715924892727E-2;
|
---|
| 136 | q = q*m1+3.01204715227604046988E-2;
|
---|
| 137 | q = q*m1+3.73774314173823228969E-2;
|
---|
| 138 | q = q*m1+4.88280347570998239232E-2;
|
---|
| 139 | q = q*m1+7.03124996963957469739E-2;
|
---|
| 140 | q = q*m1+1.24999999999870820058E-1;
|
---|
| 141 | q = q*m1+4.99999999999999999821E-1;
|
---|
| 142 | result = p-q*Math.Log(m1);
|
---|
| 143 | }
|
---|
| 144 | return result;
|
---|
| 145 | }
|
---|
| 146 |
|
---|
| 147 |
|
---|
| 148 | /*************************************************************************
|
---|
| 149 | Incomplete elliptic integral of the first kind F(phi|m)
|
---|
| 150 |
|
---|
| 151 | Approximates the integral
|
---|
| 152 |
|
---|
| 153 |
|
---|
| 154 |
|
---|
| 155 | phi
|
---|
| 156 | -
|
---|
| 157 | | |
|
---|
| 158 | | dt
|
---|
| 159 | F(phi_\m) = | ------------------
|
---|
| 160 | | 2
|
---|
| 161 | | | sqrt( 1 - m sin t )
|
---|
| 162 | -
|
---|
| 163 | 0
|
---|
| 164 |
|
---|
| 165 | of amplitude phi and modulus m, using the arithmetic -
|
---|
| 166 | geometric mean algorithm.
|
---|
| 167 |
|
---|
| 168 |
|
---|
| 169 |
|
---|
| 170 |
|
---|
| 171 | ACCURACY:
|
---|
| 172 |
|
---|
| 173 | Tested at random points with m in [0, 1] and phi as indicated.
|
---|
| 174 |
|
---|
| 175 | Relative error:
|
---|
| 176 | arithmetic domain # trials peak rms
|
---|
| 177 | IEEE -10,10 200000 7.4e-16 1.0e-16
|
---|
| 178 |
|
---|
| 179 | Cephes Math Library Release 2.8: June, 2000
|
---|
| 180 | Copyright 1984, 1987, 2000 by Stephen L. Moshier
|
---|
| 181 | *************************************************************************/
|
---|
| 182 | public static double incompleteellipticintegralk(double phi,
|
---|
| 183 | double m)
|
---|
| 184 | {
|
---|
| 185 | double result = 0;
|
---|
| 186 | double a = 0;
|
---|
| 187 | double b = 0;
|
---|
| 188 | double c = 0;
|
---|
| 189 | double e = 0;
|
---|
| 190 | double temp = 0;
|
---|
| 191 | double pio2 = 0;
|
---|
| 192 | double t = 0;
|
---|
| 193 | double k = 0;
|
---|
| 194 | int d = 0;
|
---|
| 195 | int md = 0;
|
---|
| 196 | int s = 0;
|
---|
| 197 | int npio2 = 0;
|
---|
| 198 |
|
---|
| 199 | pio2 = 1.57079632679489661923;
|
---|
| 200 | if( (double)(m)==(double)(0) )
|
---|
| 201 | {
|
---|
| 202 | result = phi;
|
---|
| 203 | return result;
|
---|
| 204 | }
|
---|
| 205 | a = 1-m;
|
---|
| 206 | if( (double)(a)==(double)(0) )
|
---|
| 207 | {
|
---|
| 208 | result = Math.Log(Math.Tan(0.5*(pio2+phi)));
|
---|
| 209 | return result;
|
---|
| 210 | }
|
---|
| 211 | npio2 = (int)Math.Floor(phi/pio2);
|
---|
| 212 | if( npio2%2!=0 )
|
---|
| 213 | {
|
---|
| 214 | npio2 = npio2+1;
|
---|
| 215 | }
|
---|
| 216 | if( npio2!=0 )
|
---|
| 217 | {
|
---|
| 218 | k = ellipticintegralk(1-a);
|
---|
| 219 | phi = phi-npio2*pio2;
|
---|
| 220 | }
|
---|
| 221 | else
|
---|
| 222 | {
|
---|
| 223 | k = 0;
|
---|
| 224 | }
|
---|
| 225 | if( (double)(phi)<(double)(0) )
|
---|
| 226 | {
|
---|
| 227 | phi = -phi;
|
---|
| 228 | s = -1;
|
---|
| 229 | }
|
---|
| 230 | else
|
---|
| 231 | {
|
---|
| 232 | s = 0;
|
---|
| 233 | }
|
---|
| 234 | b = Math.Sqrt(a);
|
---|
| 235 | t = Math.Tan(phi);
|
---|
| 236 | if( (double)(Math.Abs(t))>(double)(10) )
|
---|
| 237 | {
|
---|
| 238 | e = 1.0/(b*t);
|
---|
| 239 | if( (double)(Math.Abs(e))<(double)(10) )
|
---|
| 240 | {
|
---|
| 241 | e = Math.Atan(e);
|
---|
| 242 | if( npio2==0 )
|
---|
| 243 | {
|
---|
| 244 | k = ellipticintegralk(1-a);
|
---|
| 245 | }
|
---|
| 246 | temp = k-incompleteellipticintegralk(e, m);
|
---|
| 247 | if( s<0 )
|
---|
| 248 | {
|
---|
| 249 | temp = -temp;
|
---|
| 250 | }
|
---|
| 251 | result = temp+npio2*k;
|
---|
| 252 | return result;
|
---|
| 253 | }
|
---|
| 254 | }
|
---|
| 255 | a = 1.0;
|
---|
| 256 | c = Math.Sqrt(m);
|
---|
| 257 | d = 1;
|
---|
| 258 | md = 0;
|
---|
| 259 | while( (double)(Math.Abs(c/a))>(double)(AP.Math.MachineEpsilon) )
|
---|
| 260 | {
|
---|
| 261 | temp = b/a;
|
---|
| 262 | phi = phi+Math.Atan(t*temp)+md*Math.PI;
|
---|
| 263 | md = (int)((phi+pio2)/Math.PI);
|
---|
| 264 | t = t*(1.0+temp)/(1.0-temp*t*t);
|
---|
| 265 | c = 0.5*(a-b);
|
---|
| 266 | temp = Math.Sqrt(a*b);
|
---|
| 267 | a = 0.5*(a+b);
|
---|
| 268 | b = temp;
|
---|
| 269 | d = d+d;
|
---|
| 270 | }
|
---|
| 271 | temp = (Math.Atan(t)+md*Math.PI)/(d*a);
|
---|
| 272 | if( s<0 )
|
---|
| 273 | {
|
---|
| 274 | temp = -temp;
|
---|
| 275 | }
|
---|
| 276 | result = temp+npio2*k;
|
---|
| 277 | return result;
|
---|
| 278 | }
|
---|
| 279 |
|
---|
| 280 |
|
---|
| 281 | /*************************************************************************
|
---|
| 282 | Complete elliptic integral of the second kind
|
---|
| 283 |
|
---|
| 284 | Approximates the integral
|
---|
| 285 |
|
---|
| 286 |
|
---|
| 287 | pi/2
|
---|
| 288 | -
|
---|
| 289 | | | 2
|
---|
| 290 | E(m) = | sqrt( 1 - m sin t ) dt
|
---|
| 291 | | |
|
---|
| 292 | -
|
---|
| 293 | 0
|
---|
| 294 |
|
---|
| 295 | using the approximation
|
---|
| 296 |
|
---|
| 297 | P(x) - x log x Q(x).
|
---|
| 298 |
|
---|
| 299 | ACCURACY:
|
---|
| 300 |
|
---|
| 301 | Relative error:
|
---|
| 302 | arithmetic domain # trials peak rms
|
---|
| 303 | IEEE 0, 1 10000 2.1e-16 7.3e-17
|
---|
| 304 |
|
---|
| 305 | Cephes Math Library, Release 2.8: June, 2000
|
---|
| 306 | Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
|
---|
| 307 | *************************************************************************/
|
---|
| 308 | public static double ellipticintegrale(double m)
|
---|
| 309 | {
|
---|
| 310 | double result = 0;
|
---|
| 311 | double p = 0;
|
---|
| 312 | double q = 0;
|
---|
| 313 |
|
---|
| 314 | System.Diagnostics.Debug.Assert((double)(m)>=(double)(0) & (double)(m)<=(double)(1), "Domain error in EllipticIntegralE: m<0 or m>1");
|
---|
| 315 | m = 1-m;
|
---|
| 316 | if( (double)(m)==(double)(0) )
|
---|
| 317 | {
|
---|
| 318 | result = 1;
|
---|
| 319 | return result;
|
---|
| 320 | }
|
---|
| 321 | p = 1.53552577301013293365E-4;
|
---|
| 322 | p = p*m+2.50888492163602060990E-3;
|
---|
| 323 | p = p*m+8.68786816565889628429E-3;
|
---|
| 324 | p = p*m+1.07350949056076193403E-2;
|
---|
| 325 | p = p*m+7.77395492516787092951E-3;
|
---|
| 326 | p = p*m+7.58395289413514708519E-3;
|
---|
| 327 | p = p*m+1.15688436810574127319E-2;
|
---|
| 328 | p = p*m+2.18317996015557253103E-2;
|
---|
| 329 | p = p*m+5.68051945617860553470E-2;
|
---|
| 330 | p = p*m+4.43147180560990850618E-1;
|
---|
| 331 | p = p*m+1.00000000000000000299E0;
|
---|
| 332 | q = 3.27954898576485872656E-5;
|
---|
| 333 | q = q*m+1.00962792679356715133E-3;
|
---|
| 334 | q = q*m+6.50609489976927491433E-3;
|
---|
| 335 | q = q*m+1.68862163993311317300E-2;
|
---|
| 336 | q = q*m+2.61769742454493659583E-2;
|
---|
| 337 | q = q*m+3.34833904888224918614E-2;
|
---|
| 338 | q = q*m+4.27180926518931511717E-2;
|
---|
| 339 | q = q*m+5.85936634471101055642E-2;
|
---|
| 340 | q = q*m+9.37499997197644278445E-2;
|
---|
| 341 | q = q*m+2.49999999999888314361E-1;
|
---|
| 342 | result = p-q*m*Math.Log(m);
|
---|
| 343 | return result;
|
---|
| 344 | }
|
---|
| 345 |
|
---|
| 346 |
|
---|
| 347 | /*************************************************************************
|
---|
| 348 | Incomplete elliptic integral of the second kind
|
---|
| 349 |
|
---|
| 350 | Approximates the integral
|
---|
| 351 |
|
---|
| 352 |
|
---|
| 353 | phi
|
---|
| 354 | -
|
---|
| 355 | | |
|
---|
| 356 | | 2
|
---|
| 357 | E(phi_\m) = | sqrt( 1 - m sin t ) dt
|
---|
| 358 | |
|
---|
| 359 | | |
|
---|
| 360 | -
|
---|
| 361 | 0
|
---|
| 362 |
|
---|
| 363 | of amplitude phi and modulus m, using the arithmetic -
|
---|
| 364 | geometric mean algorithm.
|
---|
| 365 |
|
---|
| 366 | ACCURACY:
|
---|
| 367 |
|
---|
| 368 | Tested at random arguments with phi in [-10, 10] and m in
|
---|
| 369 | [0, 1].
|
---|
| 370 | Relative error:
|
---|
| 371 | arithmetic domain # trials peak rms
|
---|
| 372 | IEEE -10,10 150000 3.3e-15 1.4e-16
|
---|
| 373 |
|
---|
| 374 | Cephes Math Library Release 2.8: June, 2000
|
---|
| 375 | Copyright 1984, 1987, 1993, 2000 by Stephen L. Moshier
|
---|
| 376 | *************************************************************************/
|
---|
| 377 | public static double incompleteellipticintegrale(double phi,
|
---|
| 378 | double m)
|
---|
| 379 | {
|
---|
| 380 | double result = 0;
|
---|
| 381 | double pio2 = 0;
|
---|
| 382 | double a = 0;
|
---|
| 383 | double b = 0;
|
---|
| 384 | double c = 0;
|
---|
| 385 | double e = 0;
|
---|
| 386 | double temp = 0;
|
---|
| 387 | double lphi = 0;
|
---|
| 388 | double t = 0;
|
---|
| 389 | double ebig = 0;
|
---|
| 390 | int d = 0;
|
---|
| 391 | int md = 0;
|
---|
| 392 | int npio2 = 0;
|
---|
| 393 | int s = 0;
|
---|
| 394 |
|
---|
| 395 | pio2 = 1.57079632679489661923;
|
---|
| 396 | if( (double)(m)==(double)(0) )
|
---|
| 397 | {
|
---|
| 398 | result = phi;
|
---|
| 399 | return result;
|
---|
| 400 | }
|
---|
| 401 | lphi = phi;
|
---|
| 402 | npio2 = (int)Math.Floor(lphi/pio2);
|
---|
| 403 | if( npio2%2!=0 )
|
---|
| 404 | {
|
---|
| 405 | npio2 = npio2+1;
|
---|
| 406 | }
|
---|
| 407 | lphi = lphi-npio2*pio2;
|
---|
| 408 | if( (double)(lphi)<(double)(0) )
|
---|
| 409 | {
|
---|
| 410 | lphi = -lphi;
|
---|
| 411 | s = -1;
|
---|
| 412 | }
|
---|
| 413 | else
|
---|
| 414 | {
|
---|
| 415 | s = 1;
|
---|
| 416 | }
|
---|
| 417 | a = 1.0-m;
|
---|
| 418 | ebig = ellipticintegrale(m);
|
---|
| 419 | if( (double)(a)==(double)(0) )
|
---|
| 420 | {
|
---|
| 421 | temp = Math.Sin(lphi);
|
---|
| 422 | if( s<0 )
|
---|
| 423 | {
|
---|
| 424 | temp = -temp;
|
---|
| 425 | }
|
---|
| 426 | result = temp+npio2*ebig;
|
---|
| 427 | return result;
|
---|
| 428 | }
|
---|
| 429 | t = Math.Tan(lphi);
|
---|
| 430 | b = Math.Sqrt(a);
|
---|
| 431 |
|
---|
| 432 | //
|
---|
| 433 | // Thanks to Brian Fitzgerald <fitzgb@mml0.meche.rpi.edu>
|
---|
| 434 | // for pointing out an instability near odd multiples of pi/2
|
---|
| 435 | //
|
---|
| 436 | if( (double)(Math.Abs(t))>(double)(10) )
|
---|
| 437 | {
|
---|
| 438 |
|
---|
| 439 | //
|
---|
| 440 | // Transform the amplitude
|
---|
| 441 | //
|
---|
| 442 | e = 1.0/(b*t);
|
---|
| 443 |
|
---|
| 444 | //
|
---|
| 445 | // ... but avoid multiple recursions.
|
---|
| 446 | //
|
---|
| 447 | if( (double)(Math.Abs(e))<(double)(10) )
|
---|
| 448 | {
|
---|
| 449 | e = Math.Atan(e);
|
---|
| 450 | temp = ebig+m*Math.Sin(lphi)*Math.Sin(e)-incompleteellipticintegrale(e, m);
|
---|
| 451 | if( s<0 )
|
---|
| 452 | {
|
---|
| 453 | temp = -temp;
|
---|
| 454 | }
|
---|
| 455 | result = temp+npio2*ebig;
|
---|
| 456 | return result;
|
---|
| 457 | }
|
---|
| 458 | }
|
---|
| 459 | c = Math.Sqrt(m);
|
---|
| 460 | a = 1.0;
|
---|
| 461 | d = 1;
|
---|
| 462 | e = 0.0;
|
---|
| 463 | md = 0;
|
---|
| 464 | while( (double)(Math.Abs(c/a))>(double)(AP.Math.MachineEpsilon) )
|
---|
| 465 | {
|
---|
| 466 | temp = b/a;
|
---|
| 467 | lphi = lphi+Math.Atan(t*temp)+md*Math.PI;
|
---|
| 468 | md = (int)((lphi+pio2)/Math.PI);
|
---|
| 469 | t = t*(1.0+temp)/(1.0-temp*t*t);
|
---|
| 470 | c = 0.5*(a-b);
|
---|
| 471 | temp = Math.Sqrt(a*b);
|
---|
| 472 | a = 0.5*(a+b);
|
---|
| 473 | b = temp;
|
---|
| 474 | d = d+d;
|
---|
| 475 | e = e+c*Math.Sin(lphi);
|
---|
| 476 | }
|
---|
| 477 | temp = ebig/ellipticintegralk(m);
|
---|
| 478 | temp = temp*((Math.Atan(t)+md*Math.PI)/(d*a));
|
---|
| 479 | temp = temp+e;
|
---|
| 480 | if( s<0 )
|
---|
| 481 | {
|
---|
| 482 | temp = -temp;
|
---|
| 483 | }
|
---|
| 484 | result = temp+npio2*ebig;
|
---|
| 485 | return result;
|
---|
| 486 | return result;
|
---|
| 487 | }
|
---|
| 488 | }
|
---|
| 489 | }
|
---|