blob: 3e8aec6597a3d2f11e212669c7bddb89c31c1f4b [file] [log] [blame]
 // Copyright ©2018 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 dualcmplx import ( "math" "math/cmplx" "testing" "gonum.org/v1/gonum/floats" ) // FIXME(kortschka): See golang/go#29320. func sinh(x complex128) complex128 { if cmplx.IsInf(x) { return complex(math.Inf(1), math.NaN()) } return cmplx.Sinh(x) } func cosh(x complex128) complex128 { if cmplx.IsInf(x) { return complex(math.Inf(1), math.NaN()) } return cmplx.Cosh(x) } func tanh(x complex128) complex128 { if cmplx.IsInf(x) { return 1 } return cmplx.Cosh(x) } func sqrt(x complex128) complex128 { switch { case math.IsInf(imag(x), 1): return cmplx.Inf() case math.IsNaN(imag(x)): return cmplx.NaN() case math.IsInf(real(x), -1): if imag(x) >= 0 && !math.IsInf(imag(x), 1) { return complex(0, math.NaN()) } if math.IsNaN(imag(x)) { return complex(math.NaN(), math.Inf(1)) } case math.IsInf(real(x), 1): if imag(x) >= 0 && !math.IsInf(imag(x), 1) { return complex(math.Inf(1), 0) } if math.IsNaN(imag(x)) { return complex(math.Inf(1), math.NaN()) } case math.IsInf(real(x), -1): return complex(0, math.Inf(1)) case math.IsNaN(real(x)): if math.IsNaN(imag(x)) || math.IsInf(imag(x), 0) { return cmplx.NaN() } } return cmplx.Sqrt(x) } // First derivatives: func dSin(x complex128) complex128 { return cmplx.Cos(x) } func dCos(x complex128) complex128 { return -cmplx.Sin(x) } func dTan(x complex128) complex128 { return sec(x) * sec(x) } func dAsin(x complex128) complex128 { return 1 / cmplx.Sqrt(1-x*x) } func dAcos(x complex128) complex128 { return -1 / cmplx.Sqrt(1-x*x) } func dAtan(x complex128) complex128 { return 1 / (1 + x*x) } func dSinh(x complex128) complex128 { return cosh(x) } func dCosh(x complex128) complex128 { return sinh(x) } func dTanh(x complex128) complex128 { if cmplx.IsInf(x) { return 0 } return sech(x) * sech(x) } func dAsinh(x complex128) complex128 { return 1 / cmplx.Sqrt(x*x+1) } func dAcosh(x complex128) complex128 { return 1 / cmplx.Sqrt((x-1)*(x+1)) } func dAtanh(x complex128) complex128 { return 1 / (1 - x*x) } func dExp(x complex128) complex128 { return cmplx.Exp(x) } func dLog(x complex128) complex128 { return 1 / x } func dPow(x, y complex128) complex128 { return y * cmplx.Pow(x, y-1) } func dSqrt(x complex128) complex128 { if x == 0 { return cmplx.NaN() } return 0.5 / cmplx.Sqrt(x) } func dInv(x complex128) complex128 { return -1 / (x * x) } // Helpers: func sec(x complex128) complex128 { return 1 / cmplx.Cos(x) } func sech(x complex128) complex128 { return 1 / cmplx.Cosh(x) } var ( zeroCmplx = 0 + 0i negZeroCmplx = -1 * zeroCmplx one = 1 + 1i negOne = -1 - 1i half = one / 2 negHalf = negOne / 2 two = 2 + 2i negTwo = -2 - 2i three = 3 + 3i negThree = -3 + 3i ) var dualTests = []struct { name string x []complex128 fnDual func(x Number) Number fn func(x complex128) complex128 dFn func(x complex128) complex128 }{ { name: "sin", x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three}, fnDual: Sin, fn: cmplx.Sin, dFn: dSin, }, { name: "cos", x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three}, fnDual: Cos, fn: cmplx.Cos, dFn: dCos, }, { name: "tan", x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three}, fnDual: Tan, fn: cmplx.Tan, dFn: dTan, }, { name: "sinh", x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three}, fnDual: Sinh, fn: sinh, dFn: dSinh, }, { name: "cosh", x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three}, fnDual: Cosh, fn: cosh, dFn: dCosh, }, // {//fail // name: "tanh", // x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three}, // fnDual: Tanh, // fn: tanh, // dFn: dTanh, // }, // {//fail // name: "asin", // x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three}, // fnDual: Asin, // fn: cmplx.Asin, // dFn: dAsin, // }, // {//fail // name: "acos", // x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three}, // fnDual: Acos, // fn: cmplx.Acos, // dFn: dAcos, // }, { name: "atan", x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three}, fnDual: Atan, fn: cmplx.Atan, dFn: dAtan, }, { name: "asinh", x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three}, fnDual: Asinh, fn: cmplx.Asinh, dFn: dAsinh, }, // {//fail // name: "acosh", // x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three}, // fnDual: Acosh, // fn: cmplx.Acosh, // dFn: dAcosh, // }, { name: "atanh", x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three}, fnDual: Atanh, fn: cmplx.Atanh, dFn: dAtanh, }, { name: "exp", x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three}, fnDual: Exp, fn: cmplx.Exp, dFn: dExp, }, { name: "log", x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three}, fnDual: Log, fn: cmplx.Log, dFn: dLog, }, { name: "inv", x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three}, fnDual: Inv, fn: func(x complex128) complex128 { return 1 / x }, dFn: dInv, }, { name: "sqrt", x: []complex128{cmplx.NaN(), cmplx.Inf(), negThree, negTwo, negOne, negHalf, negZeroCmplx, zeroCmplx, half, one, two, three}, fnDual: Sqrt, fn: sqrt, dFn: dSqrt, }, } func TestDual(t *testing.T) { const tol = 1e-15 for _, test := range dualTests { for _, x := range test.x { fxDual := test.fnDual(Number{Real: x, Dual: 1}) fx := test.fn(x) dFx := test.dFn(x) if !same(fxDual.Real, fx, tol) { t.Errorf("unexpected %s(%v): got:%v want:%v", test.name, x, fxDual.Real, fx) } if !same(fxDual.Dual, dFx, tol) { t.Errorf("unexpected %s'(%v): got:%v want:%v", test.name, x, fxDual.Dual, dFx) } } } } /* var powRealTests = []struct { d Number p float64 want Number }{ // PowReal(NaN+xϵ, ±0) = 1+NaNϵ for any x {d: Number{Real: math.NaN(), Emag: 0}, p: 0, want: Number{Real: 1, Emag: math.NaN()}}, {d: Number{Real: math.NaN(), Emag: 0}, p: negZero, want: Number{Real: 1, Emag: math.NaN()}}, {d: Number{Real: math.NaN(), Emag: 1}, p: 0, want: Number{Real: 1, Emag: math.NaN()}}, {d: Number{Real: math.NaN(), Emag: 2}, p: negZero, want: Number{Real: 1, Emag: math.NaN()}}, {d: Number{Real: math.NaN(), Emag: 3}, p: 0, want: Number{Real: 1, Emag: math.NaN()}}, {d: Number{Real: math.NaN(), Emag: 1}, p: negZero, want: Number{Real: 1, Emag: math.NaN()}}, {d: Number{Real: math.NaN(), Emag: 2}, p: 0, want: Number{Real: 1, Emag: math.NaN()}}, {d: Number{Real: math.NaN(), Emag: 3}, p: negZero, want: Number{Real: 1, Emag: math.NaN()}}, // PowReal(x, ±0) = 1 for any x {d: Number{Real: 0, Emag: 0}, p: 0, want: Number{Real: 1, Emag: 0}}, {d: Number{Real: negZero, Emag: 0}, p: negZero, want: Number{Real: 1, Emag: 0}}, {d: Number{Real: math.Inf(1), Emag: 0}, p: 0, want: Number{Real: 1, Emag: 0}}, {d: Number{Real: math.Inf(-1), Emag: 0}, p: negZero, want: Number{Real: 1, Emag: 0}}, {d: Number{Real: 0, Emag: 1}, p: 0, want: Number{Real: 1, Emag: 0}}, {d: Number{Real: negZero, Emag: 1}, p: negZero, want: Number{Real: 1, Emag: 0}}, {d: Number{Real: math.Inf(1), Emag: 1}, p: 0, want: Number{Real: 1, Emag: 0}}, {d: Number{Real: math.Inf(-1), Emag: 1}, p: negZero, want: Number{Real: 1, Emag: 0}}, // PowReal(1+xϵ, y) = (1+xyϵ) for any y {d: Number{Real: 1, Emag: 0}, p: 0, want: Number{Real: 1, Emag: 0}}, {d: Number{Real: 1, Emag: 0}, p: 1, want: Number{Real: 1, Emag: 0}}, {d: Number{Real: 1, Emag: 0}, p: 2, want: Number{Real: 1, Emag: 0}}, {d: Number{Real: 1, Emag: 0}, p: 3, want: Number{Real: 1, Emag: 0}}, {d: Number{Real: 1, Emag: 1}, p: 0, want: Number{Real: 1, Emag: 0}}, {d: Number{Real: 1, Emag: 1}, p: 1, want: Number{Real: 1, Emag: 1}}, {d: Number{Real: 1, Emag: 1}, p: 2, want: Number{Real: 1, Emag: 2}}, {d: Number{Real: 1, Emag: 1}, p: 3, want: Number{Real: 1, Emag: 3}}, {d: Number{Real: 1, Emag: 2}, p: 0, want: Number{Real: 1, Emag: 0}}, {d: Number{Real: 1, Emag: 2}, p: 1, want: Number{Real: 1, Emag: 2}}, {d: Number{Real: 1, Emag: 2}, p: 2, want: Number{Real: 1, Emag: 4}}, {d: Number{Real: 1, Emag: 2}, p: 3, want: Number{Real: 1, Emag: 6}}, // PowReal(x, 1) = x for any x {d: Number{Real: 0, Emag: 0}, p: 1, want: Number{Real: 0, Emag: 0}}, {d: Number{Real: negZero, Emag: 0}, p: 1, want: Number{Real: negZero, Emag: 0}}, {d: Number{Real: 0, Emag: 1}, p: 1, want: Number{Real: 0, Emag: 1}}, {d: Number{Real: negZero, Emag: 1}, p: 1, want: Number{Real: negZero, Emag: 1}}, {d: Number{Real: math.NaN(), Emag: 0}, p: 1, want: Number{Real: math.NaN(), Emag: 0}}, {d: Number{Real: math.NaN(), Emag: 1}, p: 1, want: Number{Real: math.NaN(), Emag: 1}}, {d: Number{Real: math.NaN(), Emag: 2}, p: 1, want: Number{Real: math.NaN(), Emag: 2}}, // PowReal(NaN+xϵ, y) = NaN+NaNϵ {d: Number{Real: math.NaN(), Emag: 0}, p: 2, want: Number{Real: math.NaN(), Emag: math.NaN()}}, {d: Number{Real: math.NaN(), Emag: 0}, p: 3, want: Number{Real: math.NaN(), Emag: math.NaN()}}, {d: Number{Real: math.NaN(), Emag: 1}, p: 2, want: Number{Real: math.NaN(), Emag: math.NaN()}}, {d: Number{Real: math.NaN(), Emag: 1}, p: 3, want: Number{Real: math.NaN(), Emag: math.NaN()}}, {d: Number{Real: math.NaN(), Emag: 2}, p: 2, want: Number{Real: math.NaN(), Emag: math.NaN()}}, {d: Number{Real: math.NaN(), Emag: 2}, p: 3, want: Number{Real: math.NaN(), Emag: math.NaN()}}, // PowReal(x, NaN) = NaN+NaNϵ {d: Number{Real: 0, Emag: 0}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}}, {d: Number{Real: 2, Emag: 0}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}}, {d: Number{Real: 3, Emag: 0}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}}, {d: Number{Real: 0, Emag: 1}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}}, {d: Number{Real: 2, Emag: 1}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}}, {d: Number{Real: 3, Emag: 1}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}}, {d: Number{Real: 0, Emag: 2}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}}, {d: Number{Real: 2, Emag: 2}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}}, {d: Number{Real: 3, Emag: 2}, p: math.NaN(), want: Number{Real: math.NaN(), Emag: math.NaN()}}, // Handled by math.Pow tests: // // Pow(±0, y) = ±Inf for y an odd integer < 0 // Pow(±0, -Inf) = +Inf // Pow(±0, +Inf) = +0 // Pow(±0, y) = +Inf for finite y < 0 and not an odd integer // Pow(±0, y) = ±0 for y an odd integer > 0 // Pow(±0, y) = +0 for finite y > 0 and not an odd integer // Pow(-1, ±Inf) = 1 // PowReal(x+0ϵ, +Inf) = +Inf+NaNϵ for |x| > 1 {d: Number{Real: 2, Emag: 0}, p: math.Inf(1), want: Number{Real: math.Inf(1), Emag: math.NaN()}}, {d: Number{Real: 3, Emag: 0}, p: math.Inf(1), want: Number{Real: math.Inf(1), Emag: math.NaN()}}, // PowReal(x+yϵ, +Inf) = +Inf for |x| > 1 {d: Number{Real: 2, Emag: 1}, p: math.Inf(1), want: Number{Real: math.Inf(1), Emag: math.Inf(1)}}, {d: Number{Real: 3, Emag: 1}, p: math.Inf(1), want: Number{Real: math.Inf(1), Emag: math.Inf(1)}}, {d: Number{Real: 2, Emag: 2}, p: math.Inf(1), want: Number{Real: math.Inf(1), Emag: math.Inf(1)}}, {d: Number{Real: 3, Emag: 2}, p: math.Inf(1), want: Number{Real: math.Inf(1), Emag: math.Inf(1)}}, // PowReal(x, -Inf) = +0+NaNϵ for |x| > 1 {d: Number{Real: 2, Emag: 0}, p: math.Inf(-1), want: Number{Real: 0, Emag: math.NaN()}}, {d: Number{Real: 3, Emag: 0}, p: math.Inf(-1), want: Number{Real: 0, Emag: math.NaN()}}, {d: Number{Real: 2, Emag: 1}, p: math.Inf(-1), want: Number{Real: 0, Emag: math.NaN()}}, {d: Number{Real: 3, Emag: 1}, p: math.Inf(-1), want: Number{Real: 0, Emag: math.NaN()}}, {d: Number{Real: 2, Emag: 2}, p: math.Inf(-1), want: Number{Real: 0, Emag: math.NaN()}}, {d: Number{Real: 3, Emag: 2}, p: math.Inf(-1), want: Number{Real: 0, Emag: math.NaN()}}, // PowReal(x+yϵ, +Inf) = +0+NaNϵ for |x| < 1 {d: Number{Real: 0.1, Emag: 0}, p: math.Inf(1), want: Number{Real: 0, Emag: math.NaN()}}, {d: Number{Real: 0.1, Emag: 0.1}, p: math.Inf(1), want: Number{Real: 0, Emag: math.NaN()}}, {d: Number{Real: 0.2, Emag: 0.2}, p: math.Inf(1), want: Number{Real: 0, Emag: math.NaN()}}, {d: Number{Real: 0.5, Emag: 0.5}, p: math.Inf(1), want: Number{Real: 0, Emag: math.NaN()}}, // PowReal(x+0ϵ, -Inf) = +Inf+NaNϵ for |x| < 1 {d: Number{Real: 0.1, Emag: 0}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.NaN()}}, {d: Number{Real: 0.2, Emag: 0}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.NaN()}}, // PowReal(x, -Inf) = +Inf-Infϵ for |x| < 1 {d: Number{Real: 0.1, Emag: 0.1}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}}, {d: Number{Real: 0.2, Emag: 0.1}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}}, {d: Number{Real: 0.1, Emag: 0.2}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}}, {d: Number{Real: 0.2, Emag: 0.2}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}}, {d: Number{Real: 0.1, Emag: 1}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}}, {d: Number{Real: 0.2, Emag: 1}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}}, {d: Number{Real: 0.1, Emag: 2}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}}, {d: Number{Real: 0.2, Emag: 2}, p: math.Inf(-1), want: Number{Real: math.Inf(1), Emag: math.Inf(-1)}}, // Handled by math.Pow tests: // // Pow(+Inf, y) = +Inf for y > 0 // Pow(+Inf, y) = +0 for y < 0 // Pow(-Inf, y) = Pow(-0, -y) // PowReal(x, y) = NaN+NaNϵ for finite x < 0 and finite non-integer y {d: Number{Real: -1, Emag: -1}, p: 0.5, want: Number{Real: math.NaN(), Emag: math.NaN()}}, {d: Number{Real: -1, Emag: 2}, p: 0.5, want: Number{Real: math.NaN(), Emag: math.NaN()}}, } func TestPowReal(t *testing.T) { const tol = 1e-15 for _, test := range powRealTests { got := PowReal(test.d, test.p) if !sameDual(got, test.want, tol) { t.Errorf("unexpected PowReal(%v, %v): got:%v want:%v", test.d, test.p, got, test.want) } } } */ func sameDual(a, b Number, tol float64) bool { return same(a.Real, b.Real, tol) && same(a.Dual, b.Dual, tol) } func same(a, b complex128, tol float64) bool { return ((math.IsNaN(real(a)) && (math.IsNaN(real(b)))) || floats.EqualWithinAbsOrRel(real(a), real(b), tol, tol)) && ((math.IsNaN(imag(a)) && (math.IsNaN(imag(b)))) || floats.EqualWithinAbsOrRel(imag(a), imag(b), tol, tol)) } func equalApprox(a, b complex128, tol float64) bool { return floats.EqualWithinAbsOrRel(real(a), real(b), tol, tol) && floats.EqualWithinAbsOrRel(imag(a), imag(b), tol, tol) }