// Stijn "Adhemar" Vandamme // On the derivatives of powers of trigonometric and hyperbolic sine and cosine // // // To the extent possible under law, // the author Stijn "Adhemar" Vandamme // has dedicated all copyright and related and neighboring rights // to this software source code to the public domain worldwide. // This software source code is distributed without any warranty. // See . namespace Adhemar.Derivatives open System open System.Diagnostics.Contracts open System.Numerics open System.Text open Microsoft.FSharp.Core.Operators open Microsoft.FSharp.Reflection module Complex = let ofFloat64 f = Complex (f, 0.0) let ofInt32 n = ofFloat64 (float n) let ofUInt32 (n : uint32) = ofFloat64 (float n) let ofInt64 (n : int64) = ofFloat64 (float n) module IntegerMath = let isEven k = k % 2u = 0u let powerOfMinusOne p = if p % 2 = 0 then 1 else -1 let rec power n p = if p = 0u then 1 else n * power n (p - 1u) let rec factorial n = if n = 0u then 1u else n * factorial (n - 1u) let choose n m = if m > n then 0u else factorial n / (factorial m * factorial (n - m)) type ICoefficient = // InvalidCastExceptions may occur when adding or multiplying with ICoefficients of different types // ArgumentExceptions may occur when adding or multiplying with incompatible ICoefficients abstract member Zero : ICoefficient abstract member IsOne : bool abstract member Add : ICoefficient -> ICoefficient abstract member Scale64 : int64 -> ICoefficient abstract member Multiply : ICoefficient -> ICoefficient abstract member IsNegative : bool abstract member CoefficientString : bool * bool -> string abstract member ApplyFloat64 : Map -> float abstract member ApplyComplex : Map -> Complex type Coefficient64 (value : int64) = member __.Value = value static member Zero = Coefficient64 0L member __.IsOne = __.Value = 1L static member (+) (m : Coefficient64, n : Coefficient64) = Coefficient64 (m.Value + n.Value) static member (*) (scalingFactor : int64, m : Coefficient64) = Coefficient64 (scalingFactor * m.Value) static member (*) (m : Coefficient64, n : Coefficient64) = Coefficient64 (m.Value * n.Value) member __.IsNegative = __.Value < 0L static member (~-) (m : Coefficient64) = -1L * m override __.ToString () = sprintf "%d" __.Value member __.CoefficientString (first : bool, explicitOne : bool) = if explicitOne || (not __.IsOne && not (-__).IsOne) then if first then string __ elif __.IsNegative then sprintf " - %d" (-__).Value else sprintf " + %d" __.Value else if __.IsOne then if first then String.Empty else " + " else if first then "-" else " - " interface ICoefficient with member __.Zero = upcast Coefficient64.Zero member __.IsOne = __.IsOne member __.Add (t : ICoefficient) = upcast (__ + downcast t) member __.Scale64 (scalingFactor : int64) = upcast (scalingFactor * __) member __.Multiply (t : ICoefficient) = upcast (__ * downcast t) member __.IsNegative = __.IsNegative member __.CoefficientString (first : bool, explicitOne : bool) = __.CoefficientString (first, explicitOne) member __.ApplyFloat64 (_m : Map) = float __.Value member __.ApplyComplex (m : Map) = Complex.ofInt64 __.Value member __.IsZero = __ = Coefficient64.Zero static member (-) (m : Coefficient64, n : Coefficient64) = m + (-n) static member One = Coefficient64 1L static member MinusOne = -Coefficient64.One override __.GetHashCode () = __.Value.GetHashCode () member __.Equals (that : Coefficient64) = __.Value = that.Value interface IEquatable with member __.Equals (that : Coefficient64) = __.Equals that override __.Equals (that : obj) = match that with | null -> false | :? Coefficient64 as thatCoefficient64 -> __.Equals thatCoefficient64 | _ -> false static member op_Equality (r : Coefficient64, t : Coefficient64) = r.Equals t static member op_Inequality (r : Coefficient64, t : Coefficient64) = not (Coefficient64.op_Equality (r, t)) static member op_Implicit (value : int64) = Coefficient64 value type Polynomial<'T when 'T : equality and 'T :> ICoefficient> (x : string, ys : string list option, f : (string list option -> Map -> float) option, g : (string list option -> Map -> Complex) option, coefficients : 'T list) = static let isEmptyOrSingleton cs = match cs with | [] | [_] -> true | _ -> false static let isZero (t : 'T) = let x = t :> ICoefficient x = x.Zero static let rec exactlyOneNonZero (cs : 'T list) = match cs with | [] -> None | [c] -> Some c | h :: ts -> if isZero h then exactlyOneNonZero ts else None static let rec add (cs, ds) = match cs, ds with | _, [] -> cs | [], _ -> ds | g :: ts, h :: rs -> downcast ((g :> ICoefficient).Add h) :: add (ts, rs) static let rec shiftRank k (cs : 'T list) = if k <= 0 then List.skip (-k) cs elif List.isEmpty cs then List.empty else downcast (List.head cs).Zero :: (shiftRank (k - 1) cs) static let rec hasOnlyEvenIndexedItems (cs : 'T list) = match cs with | [] -> true | h :: ts -> hasOnlyOddIndexedItems ts and hasOnlyOddIndexedItems (cs : 'T list) = match cs with | [] -> true | h :: ts -> isZero h && hasOnlyEvenIndexedItems ts do let rec nonZeroHighestCoefficient cs = match cs with | [] -> true | [c] when isZero c -> false | _ :: ts -> nonZeroHighestCoefficient ts let isNoneOrAllAreNotEmpty zs = match zs with | None -> true | Some bs -> let isNeitherNullNorEmpty s = not (String.IsNullOrEmpty s) List.forall isNeitherNullNorEmpty bs let distinct v ws = let allDistinct mySequence = let folder state element = match state with | Some elementsSoFar when not (Set.contains element elementsSoFar) -> Some (Set.add element elementsSoFar) | _ -> None let initialState = Some Set.empty let scanning = Seq.scan folder initialState mySequence Seq.forall Option.isSome scanning match ws with | None -> true | Some zs -> allDistinct (v :: zs) if isNull x then nullArg "x" elif not (nonZeroHighestCoefficient coefficients) then invalidArg "coefficients" "Zero highest coefficient" elif not (isNoneOrAllAreNotEmpty ys) then invalidArg "ys" "At least one empty y in ys" elif isEmptyOrSingleton coefficients && x <> String.Empty then invalidArg "x" "Non-empty variable x for constant Polynomial" elif isEmptyOrSingleton coefficients && Option.isSome ys then invalidArg "ys" "Some arguments ys for constant Polynomial" elif isEmptyOrSingleton coefficients && Option.isSome f then invalidArg "ys" "Some function f for constant Polynomial" elif not (isEmptyOrSingleton coefficients) && x = String.Empty then invalidArg "x" "Empty x for non-constant Polynomial" elif not (distinct x ys) then invalidArg "ys" ("Some duplicate argument y" + " or argument y with same token as variable function x") elif Option.isSome ys && Option.isNone f then invalidArg "f" "Some argument y but no function f" static member Create (x : string, ys : string list option, f : (string list option -> Map -> float) option, g : (string list option -> Map -> Complex) option, coefficients : 'T list) = let rec deZeroHighestCoefficient cs = match cs with | [] -> List.empty | [c] when isZero c -> List.empty | h :: ts -> let deZeroTail = deZeroHighestCoefficient ts if isZero h && List.isEmpty (deZeroHighestCoefficient ts) then List.empty else h :: deZeroTail let deZeroCoefficients = deZeroHighestCoefficient coefficients if isEmptyOrSingleton deZeroCoefficients then Polynomial (String.Empty, None, None, None, deZeroCoefficients) else let emptyInsteadOfNull (s : string) = if isNull s then String.Empty else s Polynomial (emptyInsteadOfNull x, ys, f, g, coefficients) member __.X = x member __.Ys = ys member __.F = f member __.G = g member __.Coefficients = coefficients static member Zero = Polynomial (String.Empty, None, None, None, List.empty<'T>) member __.IsOne = match __.Coefficients with | [c] -> (c :> ICoefficient).IsOne | _ -> false member __.IsConstant = isEmptyOrSingleton __.Coefficients static member AreCompatible (p : Polynomial<'T>, q : Polynomial<'T>) = p.X = q.X && p.Ys = q.Ys || p.IsConstant || q.IsConstant static member Compatibleness (p : Polynomial<'T>, q : Polynomial<'T>) = if not (Polynomial.AreCompatible (p, q)) then invalidArg "q" "Polynomial p and q are incompatible" Contract.EndContractBlock () if p.IsConstant then q.X, q.Ys, q.F, q.G else p.X, p.Ys, p.F, p.G static member (+) (p : Polynomial<'T>, q : Polynomial<'T>) = if not (Polynomial.AreCompatible (p, q)) then invalidArg "q" "Polynomial p and q are incompatible" Contract.EndContractBlock () let x, ys, f, g = Polynomial.Compatibleness (p, q) Polynomial.Create (x, ys, f, g, add (p.Coefficients, q.Coefficients)) static member (*) (scalingFactor : int64, p : Polynomial<'T>) = let mapper f (t : 'T) = ((t :> ICoefficient).Scale64 f) :?> 'T Polynomial.Create (p.X, p.Ys, p.F, p.G, List.map (mapper scalingFactor) p.Coefficients) static member (*) (p : Polynomial<'T>, q : Polynomial<'T>) = if not (Polynomial.AreCompatible (p, q)) then invalidArg "q" "Polynomial p and q are incompatible" Contract.EndContractBlock () let rec multiply (cs : 'T list, ds : 'T list) = match cs with | [] -> List.empty | h :: ts -> let multiplyTerm g d = downcast ((g :> ICoefficient).Multiply d) add (List.map (multiplyTerm h) ds, shiftRank 1 (multiply (ts, ds))) let x, ys, f, g = Polynomial.Compatibleness (p, q) Polynomial.Create (x, ys, f, g, multiply (p.Coefficients, q.Coefficients)) member __.IsNegative = match exactlyOneNonZero __.Coefficients with | Some c -> (c :> ICoefficient).IsNegative | None -> false member __.IsZero = List.isEmpty __.Coefficients member __.Rank = List.length coefficients - 1 member __.IsMonomial = Option.isSome (exactlyOneNonZero __.Coefficients) static member (~-) (p : Polynomial<'T>) = -1L * p member __.YsString = match __.Ys with | None -> String.Empty | Some zs -> let separator = "," sprintf "(%s)" (String.Join (separator, zs)) override __.ToString () = if __.IsZero then "0" else let PolynomialTerms = seq { for k = __.Rank downto 0 do let d = List.item k __.Coefficients if not (isZero d) then yield (k = __.Rank), (k = 0), k, d } let stringBuilder = StringBuilder () for fi, eo, p, c in PolynomialTerms do let t = c.CoefficientString (fi, eo) ignore (stringBuilder.Append t) if p <> 0 then if not (String.IsNullOrEmpty t) && t <> "-" && not (t.EndsWith " ") then ignore (stringBuilder.Append " ") ignore (stringBuilder.Append __.X) if p <> 1 then ignore (stringBuilder.Append (sprintf "^%d" p)) ignore (stringBuilder.Append __.YsString) string stringBuilder member __.CoefficientString (first : bool, explicitOne : bool) = if explicitOne || (not __.IsOne && not (-__).IsOne) then if __.IsZero || __.IsMonomial then if __.IsNegative then if first then sprintf "-%O" -__ else sprintf " - %O" -__ else if first then sprintf "%O" __ else sprintf " + %O" __ else if first then if explicitOne then string __ else sprintf "(%O)" __ else sprintf " + (%O)" __ else if __.IsOne then if first then String.Empty else " + " else if first then "-" else " - " member __.ApplyFloat64 (m : Map) = if __.IsZero then 0.0 elif __.IsConstant then (List.exactlyOne __.Coefficients :> ICoefficient).ApplyFloat64 m else let mapper u a (i : int) c = ((c :> ICoefficient).ApplyFloat64 a) * Math.Pow (u, float i) let v = match Map.tryFind __.X m with | None -> (Option.get f) __.Ys m | Some w -> w Seq.sum (Seq.mapi (mapper v m) __.Coefficients) member __.ApplyComplex (m : Map) = if __.IsZero then Complex.Zero elif __.IsConstant then (List.exactlyOne __.Coefficients :> ICoefficient).ApplyComplex m else let seqAddComplex (zs : Complex seq) = Seq.fold (+) Complex.Zero zs let mapper u a (i : int) c = ((c :> ICoefficient).ApplyComplex a) * Complex.Pow (u, float i) let v = match Map.tryFind __.X m with | None -> (Option.get g) __.Ys m | Some w -> w seqAddComplex (Seq.mapi (mapper v m) __.Coefficients) interface ICoefficient with member __.Zero = upcast Polynomial<'T>.Zero member __.IsOne = __.IsOne member __.Add (t : ICoefficient) = upcast (__ + downcast t) member __.Scale64 (scalingFactor : int64) = upcast (scalingFactor * __) member __.Multiply (t : ICoefficient) = upcast (__ * downcast t) member __.IsNegative = __.IsNegative member __.CoefficientString (first : bool, explicitOne : bool) = __.CoefficientString (first, explicitOne) member __.ApplyFloat64 (m : Map) = __.ApplyFloat64 m member __.ApplyComplex (m : Map) = __.ApplyComplex m static member (-) (p : Polynomial<'T>, q : Polynomial<'T>) = p + (-q) static member (<<<) (p : Polynomial<'T>, k : int) = if p.IsConstant && not (p.IsZero) && k > 0 then invalidArg "k" ("Cannot shift rank of non-zero constant polynomial" + " (as constant polynomials have no variable)") Polynomial.Create (p.X, p.Ys, p.F, p.G, shiftRank k p.Coefficients) static member (>>>) (p : Polynomial<'T>, k : int) = p <<< -k member __.Derivative = if __.IsConstant then Polynomial<'T>.Zero else let mapper p (t : 'T) = (t :> ICoefficient).Scale64 (int64 p) :?> 'T let coefficients = List.mapi mapper __.Coefficients (Polynomial.Create (__.X, __.Ys, __.F, __.G, coefficients)) >>> 1 static member Constant (t : 'T) = Polynomial<'T>.Create (String.Empty, None, None, None, [t]) member __.HasOnlyEvenPowers = hasOnlyEvenIndexedItems __.Coefficients member __.HasOnlyOddPowers = hasOnlyOddIndexedItems __.Coefficients override __.GetHashCode () = __.Coefficients.GetHashCode () ^^^ (__.X.GetHashCode () <<< 1) ^^^ (__.Ys.GetHashCode () <<< 2) member __.Equals (that : Polynomial<'T>) = __.X = that.X && __.Ys = that.Ys && __.Coefficients = that.Coefficients interface IEquatable> with member __.Equals (that : Polynomial<'T>) = __.Equals that override __.Equals (that : obj) = match that with | null -> false | :? Polynomial<'T> as thatPolynomial -> __.Equals thatPolynomial | _ -> false static member op_Equality (p : Polynomial<'T>, q : Polynomial<'T>) = p.Equals q static member op_Inequality (p : Polynomial<'T>, q : Polynomial<'T>) = not (Polynomial.op_Equality (p, q)) module DiscriminatedUnions = let instance<'T> (caseInfo : UnionCaseInfo) = FSharpValue.MakeUnion (caseInfo, Array.empty) :?> 'T let cases<'T> = Seq.map instance<'T> (FSharpType.GetUnionCases typeof<'T>) module TrigHypes = let zero = Polynomial.Zero let one = Polynomial.Constant Coefficient64.One let minusOne = Polynomial.Constant Coefficient64.MinusOne let alwaysTrue _ = true type TrigHyp = | Sin | Cos | Sinh | Cosh with member __.Token = (sprintf "%A" __).ToLowerInvariant () member __.Float64Function = match __ with | TrigHyp.Sin -> Math.Sin | TrigHyp.Cos -> Math.Cos | TrigHyp.Sinh -> Math.Sinh | TrigHyp.Cosh -> Math.Cosh member __.ComplexFunction = match __ with | TrigHyp.Sin -> Complex.Sin | TrigHyp.Cos -> Complex.Cos | TrigHyp.Sinh -> Complex.Sinh | TrigHyp.Cosh -> Complex.Cosh member __.IsTrigonometric = List.contains __ [TrigHyp.Sin; TrigHyp.Cos] member __.IsHyperbolic = List.contains __ [TrigHyp.Sinh; TrigHyp.Cosh] member __.IsTrigonometricOrHyperbolicSine = List.contains __ [TrigHyp.Sin; TrigHyp.Sinh] member __.Other = match __ with | TrigHyp.Sin -> TrigHyp.Cos | TrigHyp.Cos -> TrigHyp.Sin | TrigHyp.Sinh -> TrigHyp.Cosh | TrigHyp.Cosh -> TrigHyp.Sinh member __.NegateInductionInOther = __ = TrigHyp.Cos member __.DerivativeFactorInOtherCoefficients = match __ with | TrigHyp.Cosh -> [TrigHypes.one; TrigHypes.zero; TrigHypes.one] | _ -> [TrigHypes.minusOne; TrigHypes.zero; TrigHypes.one] member __.SubstituteInSelfCoefficients = match __ with | TrigHyp.Sinh -> [TrigHypes.one; TrigHypes.zero; TrigHypes.one] | TrigHyp.Cosh -> [TrigHypes.minusOne; TrigHypes.zero; TrigHypes.one] | _ -> [TrigHypes.one; TrigHypes.zero; TrigHypes.minusOne] member __.NonNegativeConditionInOther = match __ with | TrigHyp.Cos -> IntegerMath.isEven | _ -> TrigHypes.alwaysTrue member __.NonNegativeConditionInSelf = let sinNonNegativeConditionInSelf k = IntegerMath.isEven (k / 2u) let cosNonNegativeConditionInSelf k = IntegerMath.isEven ((k + 1u) / 2u) match __ with | TrigHyp.Sin -> sinNonNegativeConditionInSelf | TrigHyp.Cos -> cosNonNegativeConditionInSelf | _ -> TrigHypes.alwaysTrue // The derivatives of the powers of trigonometric and hyperbolic sine and cosine // using the complex definitions and the binomial theorem module TrigHypComplexDefinitions = // For trigonometric sine and cosine, see also: // Feng Qi. Derivatives of tangent function and tangent numbers. // Applied Mathematics and Computation, volume 268, Thursday 1 October 2015, pages 844-858. // Specificaly pages 854-855. // // // let private seqAddComplex (zs : Complex seq) = Seq.fold (+) Complex.Zero zs // Evaluates the k'th derivative of trigHyp^n(x), evaluated at x (: float), // using the complex definition and the binomial theorem let derivativeOfPowerFloat64 trigHyp k n x = let term (t : TrigHyp) r m y q = let twoQMinusM = 2 * q - int m let factorCosine = int (IntegerMath.choose m (uint32 q)) * IntegerMath.power twoQMinusM r let factorInt32 = if t.IsTrigonometricOrHyperbolicSine then IntegerMath.powerOfMinusOne q * factorCosine else factorCosine let factor = float factorInt32 if t.IsTrigonometric then let argument = let halfPi = Math.PI / 2.0 if t = TrigHyp.Sin then float (r - m) * halfPi + float twoQMinusM * y else float r * halfPi + float twoQMinusM * y Complex (factor * Math.Cos argument, factor * Math.Sin argument) else Complex.ofFloat64 (factor * Math.Exp (float twoQMinusM * y)) let ns = int n let denominator = Complex.ofInt32 (IntegerMath.power 2 n) let unsigned = seqAddComplex (Seq.map (term trigHyp k n x) {0 .. ns}) / denominator if trigHyp.IsTrigonometricOrHyperbolicSine then let sign = Complex.ofInt32 (IntegerMath.powerOfMinusOne ns) sign * unsigned else unsigned // Evaluates the k'th derivative of trigHyp^n(x), evaluated at x (: Complex), // using the complex definition and the binomial theorem let derivativeOfPowerComplex trigHyp k n x = let term (t : TrigHyp) r m y q = let twoQMinusM = 2 * q - int m let factorCosine = int (IntegerMath.choose m (uint32 q)) * IntegerMath.power twoQMinusM r let factorInt32 = if t.IsTrigonometricOrHyperbolicSine then IntegerMath.powerOfMinusOne q * factorCosine else factorCosine let factor = Complex.ofInt32 factorInt32 if t.IsTrigonometric then let argument = let halfPi = Complex.ofFloat64 (Math.PI / 2.0) if t = TrigHyp.Sin then Complex.ofUInt32 (r - m) * halfPi + Complex.ofInt32 twoQMinusM * y else Complex.ofUInt32 r * halfPi + Complex.ofInt32 twoQMinusM * y factor * Complex.Cos argument + factor * Complex.Sin argument * Complex.ImaginaryOne else factor * Complex.Exp (Complex.ofInt32 twoQMinusM * y) let ns = int n let denominator = Complex.ofInt32 (IntegerMath.power 2 n) let unsigned = seqAddComplex (Seq.map (term trigHyp k n x) {0 .. ns}) / denominator if trigHyp.IsTrigonometricOrHyperbolicSine then let sign = Complex.ofInt32 (IntegerMath.powerOfMinusOne ns) sign * unsigned else unsigned // The derivatives of the powers of trigonometric and hyperbolic sine and cosine // using polynomials // both intermediate (polynomial in other) // and final (polynomial in self, possibly to be multiplied by a single additional factor other) module TrigHypPolynomials = let singleVariableFunction (h : float -> float) (ys : string list option) (m : Map) = h (Map.find (List.head (Option.get ys)) m) let singleVariableFunctionComplex (h : Complex -> Complex) (ys : string list option) (m : Map) = h (Map.find (List.head (Option.get ys)) m) // Intermediate polynomials: f_k(n)(u), g_k(n)(u), h_k(n)(u) let rec derivativePowerPolynomialInOther (trigHyp : TrigHyp) (k : uint32) = if k = 0u then Polynomial.Constant TrigHypes.one else let n = "n" let x = trigHyp.Other.Token let ys = Some ["x"] let f = Some (singleVariableFunction trigHyp.Other.Float64Function) let g = Some (singleVariableFunctionComplex trigHyp.Other.ComplexFunction) let zero = Polynomial.Zero let p (j : uint32) = let oneMinusJ = Coefficient64 (1L - int64 j) let nPlusOneMinusJ = Polynomial (n, None, None, None, [oneMinusJ; Coefficient64.One]) Polynomial (x, ys, f, g, [zero; nPlusOneMinusJ]) let q = Polynomial (x, ys, f, g, trigHyp.DerivativeFactorInOtherCoefficients) let previousPolynomial = derivativePowerPolynomialInOther trigHyp (k - 1u) let r = p k * previousPolynomial + q * previousPolynomial.Derivative if trigHyp.NegateInductionInOther then -r else r // Final polynomials: s_k(n)(v), c_k(n)(v), t_k(n)(v), d_k(n)(v) let derivativePowerPolynomialInSelf (trigHyp : TrigHyp) (k : uint32) = let derivativePowerPolynomialFromPolynomialInOther (g : TrigHyp) (p : Polynomial>) = if not p.HasOnlyEvenPowers then invalidArg "p" "Polynomial p has at least one odd power" elif not (p.X = g.Other.Token && Option.isSome p.Ys || p.IsConstant) then invalidArg "p" (sprintf "Polynomial p is not compatible with a polynomial in %s" g.Other.Token) Contract.EndContractBlock () let mapper (h : TrigHyp) i t = let rec powerOfSinSquaredMinusOne n = let one = Polynomial.Constant Coefficient64.One if n = 0u then Polynomial.Constant one else let x = h.Token let ys = p.Ys let f = Some (singleVariableFunction h.Float64Function) let g = Some (singleVariableFunctionComplex h.ComplexFunction) let r = Polynomial (x, ys, f, g, h.SubstituteInSelfCoefficients) if n = 1u then r else r * powerOfSinSquaredMinusOne (n - 1u) (Polynomial.Constant t) * powerOfSinSquaredMinusOne ((uint32 i) / 2u) Seq.sum (Seq.mapi (mapper g) p.Coefficients) let q = if IntegerMath.isEven k then derivativePowerPolynomialInOther trigHyp k else (derivativePowerPolynomialInOther trigHyp k) >>> 1 derivativePowerPolynomialFromPolynomialInOther trigHyp q // Prints the intermediate and the final expression for the k'th derivative of trigHyp^n(x) let printDerivative stream (trigHyp : TrigHyp) k = let otherOrEmpty = if IntegerMath.isEven k then String.Empty else sprintf " %s(x)" trigHyp.Other.Token let token = trigHyp.Token fprintfn stream "d^%d/dx^%d [%s^n(x)]" k k token let inOther = derivativePowerPolynomialInOther trigHyp k if trigHyp.NonNegativeConditionInOther k then fprintfn stream " = %s^(n-%d)(x) [%A]" token k inOther else fprintfn stream " = -%s^(n-%d)(x) [%A]" token k -inOther let inSelf = derivativePowerPolynomialInSelf trigHyp k if trigHyp.NonNegativeConditionInSelf k then fprintfn stream " = %s^(n-%d)(x)%s [%A]" token k otherOrEmpty inSelf else fprintfn stream " = -%s^(n-%d)(x)%s [%A]" token k otherOrEmpty -inSelf if k = 0u then fprintfn stream " = %s^n(x)" token let printDerivativesUpto stream (trigHyp : TrigHyp) (lastK : uint32) = fprintfn stream "Derivatives of %s^n(x)" trigHyp.Token fprintfn stream "" for k = 0 to int lastK do printDerivative stream trigHyp (uint32 k) // Evaluates the k'th derivative of trigHyp^n(x), evaluated at x (: float), using the intermediate step let derivativeOfPowerInOtherFloat64 (trigHyp : TrigHyp) k n x = let y = trigHyp.Float64Function x let o = trigHyp.Other.Float64Function x let p = derivativePowerPolynomialInOther trigHyp k let m = Map.ofList [trigHyp.Other.Token, o; "n", float n] Math.Pow (y, float (n - int k)) * p.ApplyFloat64 m // Evaluates the k'th derivative of trigHyp^n(x), evaluated at x (: Complex), using the intermediate step let derivativeOfPowerInOtherComplex (trigHyp : TrigHyp) k n z = let y = trigHyp.ComplexFunction z let o = trigHyp.Other.ComplexFunction z let p = derivativePowerPolynomialInOther trigHyp k let m = Map.ofList [trigHyp.Other.Token, o; "n", Complex.ofInt32 n] Complex.Pow (y, Complex.ofInt32 (n - int k)) * p.ApplyComplex m // Evaluates the k'th derivative of trigHyp^n(x), evaluated at x (: float), using the final step let derivativeOfPowerFloat64 (trigHyp : TrigHyp) k n x = let y = trigHyp.Float64Function x let p = derivativePowerPolynomialInSelf trigHyp k let m = Map.ofList [trigHyp.Token, y; "n", float n] let a = Math.Pow (y, float (n - int k)) * p.ApplyFloat64 m if IntegerMath.isEven k then a else trigHyp.Other.Float64Function x * a // Evaluates the k'th derivative of trigHyp^n(x), evaluated at x (: Complex), using the final step let derivativeOfPowerComplex (trigHyp : TrigHyp) k n z = let y = trigHyp.ComplexFunction z let p = derivativePowerPolynomialInSelf trigHyp k let m = Map.ofList [trigHyp.Token, y; "n", Complex.ofInt32 n] let a = Complex.Pow (y, Complex.ofInt32 (n - int k)) * p.ApplyComplex m if IntegerMath.isEven k then a else trigHyp.Other.ComplexFunction z * a module Derivatives = // Prints the k'th derivative of trigHyp^n(z) evaluated at x (: float), // calculated using the different expressions let printDerivativeOfPower stream (trigHyp : TrigHyp) k n x = let z = Complex.ofFloat64 x fprintfn stream "d^%d/dx^%d [%s^%d(x)] | x = %s" k k trigHyp.Token n (x.ToString "R") fprintfn stream " = %O" (TrigHypComplexDefinitions.derivativeOfPowerComplex trigHyp k (uint32 n) z) fprintfn stream " = %O" (TrigHypComplexDefinitions.derivativeOfPowerFloat64 trigHyp k (uint32 n) x) fprintfn stream " = %O" (TrigHypPolynomials.derivativeOfPowerInOtherComplex trigHyp k n z) fprintfn stream " = %s" ((TrigHypPolynomials.derivativeOfPowerInOtherFloat64 trigHyp k n x).ToString "R") fprintfn stream " = %O" (TrigHypPolynomials.derivativeOfPowerComplex trigHyp k n z) fprintfn stream " = %s" ((TrigHypPolynomials.derivativeOfPowerFloat64 trigHyp k n x).ToString "R") [] let main _args = let lastK = 6u let stream = Console.Out for th in DiscriminatedUnions.cases do TrigHypPolynomials.printDerivativesUpto stream th lastK fprintfn stream "" ignore (Console.ReadKey ()) 0