blob: 75b99fd19624304cca1a004a410796cb86fcbfd6 [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 ( "fmt" "math" "math/cmplx" "strings" ) // Number is a float64 precision anti-commutative dual complex number. type Number struct { Real, Dual complex128 } // Format implements fmt.Formatter. func (d Number) Format(fs fmt.State, c rune) { prec, pOk := fs.Precision() if !pOk { prec = -1 } width, wOk := fs.Width() if !wOk { width = -1 } switch c { case 'v': if fs.Flag('#') { fmt.Fprintf(fs, "%T{Real:%#v, Dual:%#v}", d, d.Real, d.Dual) return } if fs.Flag('+') { fmt.Fprintf(fs, "{Real:%+v, Dual:%+v}", d.Real, d.Dual) return } c = 'g' prec = -1 fallthrough case 'e', 'E', 'f', 'F', 'g', 'G': fre := fmtString(fs, c, prec, width, false) fim := fmtString(fs, c, prec, width, true) fmt.Fprintf(fs, fmt.Sprintf("(%s+%[2]sϵ)", fre, fim), d.Real, d.Dual) default: fmt.Fprintf(fs, "%%!%c(%T=%[2]v)", c, d) return } } // This is horrible, but it's what we have. func fmtString(fs fmt.State, c rune, prec, width int, wantPlus bool) string { var b strings.Builder b.WriteByte('%') for _, f := range "0+- " { if fs.Flag(int(f)) || (f == '+' && wantPlus) { b.WriteByte(byte(f)) } } if width >= 0 { fmt.Fprint(&b, width) } if prec >= 0 { b.WriteByte('.') if prec > 0 { fmt.Fprint(&b, prec) } } b.WriteRune(c) return b.String() } // Add returns the sum of x and y. func Add(x, y Number) Number { return Number{ Real: x.Real + y.Real, Dual: x.Dual + y.Dual, } } // Sub returns the difference of x and y, x-y. func Sub(x, y Number) Number { return Number{ Real: x.Real - y.Real, Dual: x.Dual - y.Dual, } } // Mul returns the dual product of x and y, x×y. func Mul(x, y Number) Number { return Number{ Real: x.Real * y.Real, Dual: x.Real*y.Dual + x.Dual*cmplx.Conj(y.Real), } } // Inv returns the dual inverse of d. func Inv(d Number) Number { return Number{ Real: 1 / d.Real, Dual: -d.Dual / (d.Real * cmplx.Conj(d.Real)), } } // Conj returns the conjugate of d₁+d₂ϵ, d̅₁+d₂ϵ. func Conj(d Number) Number { return Number{ Real: cmplx.Conj(d.Real), Dual: d.Dual, } } // Scale returns d scaled by f. func Scale(f float64, d Number) Number { return Number{Real: complex(f, 0) * d.Real, Dual: complex(f, 0) * d.Dual} } // Abs returns the absolute value of d. func Abs(d Number) float64 { return cmplx.Abs(d.Real) } // PowReal returns d**p, the base-d exponential of p. // // Special cases are (in order): // PowReal(NaN+xϵ, ±0) = 1+NaNϵ for any x // Pow(0+xϵ, y) = 0+Infϵ for all y < 1. // Pow(0+xϵ, y) = 0 for all y > 1. // PowReal(x, ±0) = 1 for any x // PowReal(1+xϵ, y) = 1+xyϵ for any y // Pow(Inf, y) = +Inf+NaNϵ for y > 0 // Pow(Inf, y) = +0+NaNϵ for y < 0 // PowReal(x, 1) = x for any x // PowReal(NaN+xϵ, y) = NaN+NaNϵ // PowReal(x, NaN) = NaN+NaNϵ // PowReal(-1, ±Inf) = 1 // PowReal(x+0ϵ, +Inf) = +Inf+NaNϵ for |x| > 1 // PowReal(x+yϵ, +Inf) = +Inf for |x| > 1 // PowReal(x, -Inf) = +0+NaNϵ for |x| > 1 // PowReal(x, +Inf) = +0+NaNϵ for |x| < 1 // PowReal(x+0ϵ, -Inf) = +Inf+NaNϵ for |x| < 1 // PowReal(x, -Inf) = +Inf-Infϵ for |x| < 1 // PowReal(+Inf, y) = +Inf for y > 0 // PowReal(+Inf, y) = +0 for y < 0 // PowReal(-Inf, y) = Pow(-0, -y) func PowReal(d Number, p float64) Number { switch { case p == 0: switch { case cmplx.IsNaN(d.Real): return Number{Real: 1, Dual: cmplx.NaN()} case d.Real == 0, cmplx.IsInf(d.Real): return Number{Real: 1} } case p == 1: if cmplx.IsInf(d.Real) { d.Dual = cmplx.NaN() } return d case math.IsInf(p, 1): if d.Real == -1 { return Number{Real: 1, Dual: cmplx.NaN()} } if Abs(d) > 1 { if d.Dual == 0 { return Number{Real: cmplx.Inf(), Dual: cmplx.NaN()} } return Number{Real: cmplx.Inf(), Dual: cmplx.Inf()} } return Number{Real: 0, Dual: cmplx.NaN()} case math.IsInf(p, -1): if d.Real == -1 { return Number{Real: 1, Dual: cmplx.NaN()} } if Abs(d) > 1 { return Number{Real: 0, Dual: cmplx.NaN()} } if d.Dual == 0 { return Number{Real: cmplx.Inf(), Dual: cmplx.NaN()} } return Number{Real: cmplx.Inf(), Dual: cmplx.Inf()} case math.IsNaN(p): return Number{Real: cmplx.NaN(), Dual: cmplx.NaN()} case d.Real == 0: if p < 1 { return Number{Real: d.Real, Dual: cmplx.Inf()} } return Number{Real: d.Real} case cmplx.IsInf(d.Real): if p < 0 { return Number{Real: 0, Dual: cmplx.NaN()} } return Number{Real: cmplx.Inf(), Dual: cmplx.NaN()} } return Pow(d, Number{Real: complex(p, 0)}) } // Pow returns d**p, the base-d exponential of p. func Pow(d, p Number) Number { return Exp(Mul(p, Log(d))) } // Sqrt returns the square root of d. // // Special cases are: // Sqrt(+Inf) = +Inf // Sqrt(±0) = (±0+Infϵ) // Sqrt(x < 0) = NaN // Sqrt(NaN) = NaN func Sqrt(d Number) Number { return PowReal(d, 0.5) } // Exp returns e**q, the base-e exponential of d. // // Special cases are: // Exp(+Inf) = +Inf // Exp(NaN) = NaN // Very large values overflow to 0 or +Inf. // Very small values underflow to 1. func Exp(d Number) Number { fn := cmplx.Exp(d.Real) if imag(d.Real) == 0 { return Number{Real: fn, Dual: fn * d.Dual} } conj := cmplx.Conj(d.Real) return Number{ Real: fn, Dual: ((fn - cmplx.Exp(conj)) / (d.Real - conj)) * d.Dual, } } // Log returns the natural logarithm of d. // // Special cases are: // Log(+Inf) = (+Inf+0ϵ) // Log(0) = (-Inf±Infϵ) // Log(x < 0) = NaN // Log(NaN) = NaN func Log(d Number) Number { fn := cmplx.Log(d.Real) switch { case d.Real == 0: return Number{ Real: fn, Dual: complex(math.Copysign(math.Inf(1), real(d.Real)), math.NaN()), } case imag(d.Real) == 0: return Number{ Real: fn, Dual: d.Dual / d.Real, } case cmplx.IsInf(d.Real): return Number{ Real: fn, Dual: 0, } } conj := cmplx.Conj(d.Real) return Number{ Real: fn, Dual: ((fn - cmplx.Log(conj)) / (d.Real - conj)) * d.Dual, } }