mirror of
https://github.com/gonum/gonum.git
synced 2025-10-26 00:30:27 +08:00
802755713a
* Implemented Carlson's and Legendre's elliptic integrals. Added Carlson's symmetric elliptic integrals RF and RD. Added Legendre's elliptic integrals of the 1st and 2nd kinds. * Combined const declarations. Combined const decls. Used 1 as a seed for the RNG. * Renamed CarsonRF, CarlsonRD to EllipticRF, EllipticRD. Renamed CarsonRF, CarlsonRD to EllipticRF, EllipticRD. Updated the docs as per the reviewer's suggestions. Used bit shifts to define 2^-1022. * Improved the docs. Improved the docs. Added an ArXiv preprint link. * Improved doc formatting. Used two spaces instead of tabs. Removed extra blank lines. * Improved the docs. Elaborated on the origins of the RF and RD integrals. * Multiple fixes as suggested by the reviewer. Multiple fixes in the docs and code. Added spot checks for several precomputed values of EllipticF and EllipticE.
104 lines
3.5 KiB
Go
104 lines
3.5 KiB
Go
// Copyright ©2017 The gonum Authors. All rights reserved.
|
|
// Use of this source code is governed by a BSD-style
|
|
// license that can be found in the LICENSE file.
|
|
|
|
package mathext
|
|
|
|
import (
|
|
"math"
|
|
"math/rand"
|
|
"testing"
|
|
)
|
|
|
|
// Testing EllipticF (and EllipticRF) using the addition theorems from http://dlmf.nist.gov/19.11.i
|
|
func TestEllipticF(t *testing.T) {
|
|
const tol = 1.0e-14
|
|
rnd := rand.New(rand.NewSource(1))
|
|
|
|
// The following EllipticF(pi/3,m), m=0.1(0.1)0.9 was computed in Maxima 5.38.0 using Bigfloat arithmetic.
|
|
vF := [...]float64{
|
|
1.0631390181954904767742338285104637431858016483079,
|
|
1.0803778062523490005579242592072579594037132891908,
|
|
1.0991352230920430074586978843452269008747645822123,
|
|
1.1196949183404746257742176145632376703505764745654,
|
|
1.1424290580457772555013955266260457822322036529624,
|
|
1.1678400583161860445148860686430780757517286094732,
|
|
1.1966306515644649360767197589467723191317720122309,
|
|
1.2298294422249382706933871574135731278765534034979,
|
|
1.2690359140762658660446752406901433173504503955036,
|
|
}
|
|
phi := math.Pi / 3
|
|
for m := 1; m <= 9; m++ {
|
|
mf := float64(m) / 10
|
|
delta := math.Abs(EllipticF(phi, mf) - vF[m-1])
|
|
if delta > tol {
|
|
t.Fatalf("EllipticF(pi/3,m) test fail for m=%v", mf)
|
|
}
|
|
}
|
|
|
|
for test := 0; test < 100; test++ {
|
|
alpha := rnd.Float64() * math.Pi / 4
|
|
beta := rnd.Float64() * math.Pi / 4
|
|
for mi := 0; mi < 9999; mi++ {
|
|
m := float64(mi) / 10000
|
|
Fa := EllipticF(alpha, m)
|
|
Fb := EllipticF(beta, m)
|
|
sina, cosa := math.Sincos(alpha)
|
|
sinb, cosb := math.Sincos(beta)
|
|
tan := (sina*math.Sqrt(1-m*sinb*sinb) + sinb*math.Sqrt(1-m*sina*sina)) / (cosa + cosb)
|
|
gamma := 2 * math.Atan(tan)
|
|
Fg := EllipticF(gamma, m)
|
|
delta := math.Abs(Fa + Fb - Fg)
|
|
if delta > tol {
|
|
t.Fatalf("EllipticF test fail for m=%v, alpha=%v, beta=%v", m, alpha, beta)
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Testing EllipticE (and EllipticRF, EllipticRD) using the addition theorems from http://dlmf.nist.gov/19.11.i
|
|
func TestEllipticE(t *testing.T) {
|
|
const tol = 1.0e-14
|
|
rnd := rand.New(rand.NewSource(1))
|
|
|
|
// The following EllipticE(pi/3,m), m=0.1(0.1)0.9 was computed in Maxima 5.38.0 using Bigfloat arithmetic.
|
|
vE := [...]float64{
|
|
1.0316510822817691068014397636905610074934300946730,
|
|
1.0156973658341766636288643556414001451527597364432,
|
|
9.9929636467826398814855428365155224243586391115108e-1,
|
|
9.8240033979859736941287149003648737502960015189033e-1,
|
|
9.6495145764299257550956863602992167490195750321518e-1,
|
|
9.4687829659158090935158610908054896203271861698355e-1,
|
|
9.2809053417715769009517654522979827392794124845027e-1,
|
|
9.0847044378047233264777277954768245721857017157916e-1,
|
|
8.8785835036531301307661603341327881634688308777383e-1,
|
|
}
|
|
phi := math.Pi / 3
|
|
for m := 1; m <= 9; m++ {
|
|
mf := float64(m) / 10
|
|
delta := math.Abs(EllipticE(phi, mf) - vE[m-1])
|
|
if delta > tol {
|
|
t.Fatalf("EllipticE(pi/3,m) test fail for m=%v", mf)
|
|
}
|
|
}
|
|
|
|
for test := 0; test < 100; test++ {
|
|
alpha := rnd.Float64() * math.Pi / 4
|
|
beta := rnd.Float64() * math.Pi / 4
|
|
for mi := 0; mi < 9999; mi++ {
|
|
m := float64(mi) / 10000
|
|
Ea := EllipticE(alpha, m)
|
|
Eb := EllipticE(beta, m)
|
|
sina, cosa := math.Sincos(alpha)
|
|
sinb, cosb := math.Sincos(beta)
|
|
tan := (sina*math.Sqrt(1-m*sinb*sinb) + sinb*math.Sqrt(1-m*sina*sina)) / (cosa + cosb)
|
|
gamma := 2 * math.Atan(tan)
|
|
Eg := EllipticE(gamma, m)
|
|
delta := math.Abs(Ea + Eb - Eg - m*sina*sinb*math.Sin(gamma))
|
|
if delta > tol {
|
|
t.Fatalf("EllipticE test fail for m=%v, alpha=%v, beta=%v", m, alpha, beta)
|
|
}
|
|
}
|
|
}
|
|
}
|