Back to index/Generated with Clerk from src/emmy/series.cljc@f807468
(ns emmy.series
"This namespace contains an implementation of two data types:
- [[Series]], which represents a generic infinite series of arbitrary values, and
- [[PowerSeries]], a series that represents a power series in a single
variable; in other words, a series where the nth entry is interpreted as
the coefficient of $x^n$:
```
$$[a b c d ...] == $a + bx + cx^2 + dx^3 + ...$$
```
Many of the functions below draw on the [[emmy.series.impl]] namespace,
which implements many of these operations on bare Clojure sequences.
The implementation follows Doug McIlroy's beautiful paper, [\"Power Series,
Power
Serious\"](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.333.3156&rep=rep1&type=pdf).
Doug also has a 10-line version in Haskell on [his
website](https://www.cs.dartmouth.edu/~doug/powser.html)."
(:refer-clojure :exclude [identity])
(:require [emmy.differential :as d]
[emmy.expression] ;; for the effect of the defmethod of zero? for Literals
[emmy.function :as f]
[emmy.generic :as g]
[emmy.series.impl :as i]
[emmy.util :as u]
[emmy.value :as v])
#?(:clj
(:import (clojure.lang AFn IFn IObj Seqable Sequential))))
(declare fmap s-zero s-one s-identity series-value)
(0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 10+ more elided)
(deftype Series [xs m]
f/IArity
(arity [_] (f/arity (first xs)))
d/IPerturbed
(perturbed? [_] false)
(replace-tag [s old new] (fmap #(d/replace-tag % old new) s))
(extract-tangent [s tag] (fmap #(d/extract-tangent % tag) s))
v/IKind
(kind [_] ::series)
Object
(toString [S] (str (g/freeze S)))
#?@
(:clj
[IObj
(meta [_] m)
(withMeta [_ meta] (Series. xs meta))
Sequential
Seqable
(seq [_] xs)
IFn
;; Invoking a series uses `series-value` to generate a new series.
(invoke [_]
(Series. (series-value xs []) nil))
(invoke [_ a]
(Series. (series-value xs [a]) nil))
(invoke [_ a b]
(Series. (series-value xs [a b]) nil))
(invoke [_ a b c]
(Series. (series-value xs [a b c]) nil))
(invoke [_ a b c d]
(Series. (series-value xs [a b c d]) nil))
(invoke [_ a b c d e]
(Series. (series-value xs [a b c d e]) nil))
(invoke [_ a b c d e f]
(Series. (series-value xs [a b c d e f]) nil))
(invoke [_ a b c d e f g]
(Series. (series-value xs [a b c d e f g]) nil))
(invoke [_ a b c d e f g h]
(Series. (series-value xs [a b c d e f g h]) nil))
(invoke [_ a b c d e f g h i]
(Series. (series-value xs [a b c d e f g h i]) nil))
(invoke [_ a b c d e f g h i j]
(Series. (series-value xs [a b c d e f g h i j]) nil))
(invoke [_ a b c d e f g h i j k]
(Series. (series-value xs [a b c d e f g h i j k]) nil))
(invoke [_ a b c d e f g h i j k l]
(Series. (series-value xs [a b c d e f g h i j k l]) nil))
(invoke [_ a b c d e f g h i j k l m]
(Series. (series-value xs [a b c d e f g h i j k l m]) nil))
(invoke [_ a b c d e f g h i j k l m n]
(Series. (series-value xs [a b c d e f g h i j k l m n]) nil))
(invoke [_ a b c d e f g h i j k l m n o]
(Series. (series-value xs [a b c d e f g h i j k l m n o]) nil))
(invoke [_ a b c d e f g h i j k l m n o p]
(Series. (series-value xs [a b c d e f g h i j k l m n o p]) nil))
(invoke [_ a b c d e f g h i j k l m n o p q]
(Series. (series-value xs [a b c d e f g h i j k l m n o p q]) nil))
(invoke [_ a b c d e f g h i j k l m n o p q r]
(Series. (series-value xs [a b c d e f g h i j k l m n o p q r]) nil))
(invoke [_ a b c d e f g h i j k l m n o p q r s]
(Series. (series-value xs [a b c d e f g h i j k l m n o p q r s]) nil))
(invoke [_ a b c d e f g h i j k l m n o p q r s t]
(Series. (series-value xs [a b c d e f g h i j k l m n o p q r s t]) nil))
(invoke [_ a b c d e f g h i j k l m n o p q r s t rest]
(Series. (series-value xs (concat [a b c d e f g h i j k l m n o p q r s t] rest)) nil))
(applyTo [s xs] (AFn/applyToHelper s xs))]
:cljs
[IMeta
(-meta [_] m)
IWithMeta
(-with-meta [_ meta] (Series. xs meta))
ISequential
ISeqable
(-seq [_] xs)
IPrintWithWriter
(-pr-writer [x writer _]
(write-all writer
"#object[emmy.series.Series \""
(.toString x)
"\"]"))
IFn
(-invoke [_]
(Series. (series-value xs []) nil))
(-invoke [_ a]
(Series. (series-value xs [a]) nil))
(-invoke [_ a b]
(Series. (series-value xs [a b]) nil))
(-invoke [_ a b c]
(Series. (series-value xs [a b c]) nil))
(-invoke [_ a b c d]
(Series. (series-value xs [a b c d]) nil))
(-invoke [_ a b c d e]
(Series. (series-value xs [a b c d e]) nil))
(-invoke [_ a b c d e f]
(Series. (series-value xs [a b c d e f]) nil))
(-invoke [_ a b c d e f g]
(Series. (series-value xs [a b c d e f g]) nil))
(-invoke [_ a b c d e f g h]
(Series. (series-value xs [a b c d e f g h]) nil))
(-invoke [_ a b c d e f g h i]
(Series. (series-value xs [a b c d e f g h i]) nil))
(-invoke [_ a b c d e f g h i j]
(Series. (series-value xs [a b c d e f g h i j]) nil))
(-invoke [_ a b c d e f g h i j k]
(Series. (series-value xs [a b c d e f g h i j k]) nil))
(-invoke [_ a b c d e f g h i j k l]
(Series. (series-value xs [a b c d e f g h i j k l]) nil))
(-invoke [_ a b c d e f g h i j k l m]
(Series. (series-value xs [a b c d e f g h i j k l m]) nil))
(-invoke [_ a b c d e f g h i j k l m n]
(Series. (series-value xs [a b c d e f g h i j k l m n]) nil))
(-invoke [_ a b c d e f g h i j k l m n o]
(Series. (series-value xs [a b c d e f g h i j k l m n o]) nil))
(-invoke [_ a b c d e f g h i j k l m n o p]
(Series. (series-value xs [a b c d e f g h i j k l m n o p]) nil))
(-invoke [_ a b c d e f g h i j k l m n o p q]
(Series. (series-value xs [a b c d e f g h i j k l m n o p q]) nil))
(-invoke [_ a b c d e f g h i j k l m n o p q r]
(Series. (series-value xs [a b c d e f g h i j k l m n o p q r]) nil))
(-invoke [_ a b c d e f g h i j k l m n o p q r s]
(Series. (series-value xs [a b c d e f g h i j k l m n o p q r s]) nil))
(-invoke [_ a b c d e f g h i j k l m n o p q r s t]
(Series. (series-value xs [a b c d e f g h i j k l m n o p q r s t]) nil))
(-invoke [_ a b c d e f g h i j k l m n o p q r s t rest]
(Series. (series-value xs (concat [a b c d e f g h i j k l m n o p q r s t] rest)) nil))]))
#object[emmy.series$eval88223$__GT_Series__88229 0x120cd6d7 "
emmy.series$eval88223$__GT_Series__88229@120cd6d7"
]

Unmap the auto-generated constructor and replace it with a better one.

#?(:clj
(defmethod print-method Series [^Series s ^java.io.Writer w]
(.write w (str "#object[emmy.series.Series \""
(.toString s)
"\"]"))))
#object[clojure.lang.MultiFn 0x5877451f "
clojure.lang.MultiFn@5877451f"
]

Power Series

The primary difference from Series is the IFn implementation; application of a power series multiples each coefficient by a successively larger power of its (single, for now) argument.

TODO Modify this description once we implement multivariable power series!

(declare zero one identity power-series-value)
(1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 10+ more elided)
(deftype PowerSeries [xs m]
f/IArity
(arity [_] [:exactly 1])
d/IPerturbed
(perturbed? [_] false)
(replace-tag [s old new] (fmap #(d/replace-tag % old new) s))
(extract-tangent [s tag] (fmap #(d/extract-tangent % tag) s))
v/IKind
(kind [_] ::power-series)
Object
(toString [S] (str (g/freeze S)))
#?@
(:clj
[IObj
(meta [_] m)
(withMeta [_ meta] (PowerSeries. xs meta))
Sequential
Seqable
(seq [_] xs)
IFn
(invoke [_ a] (Series. (power-series-value xs a) nil))
(applyTo [s xs] (AFn/applyToHelper s xs))]
:cljs
[IMeta
(-meta [_] m)
IWithMeta
(-with-meta [_ meta] (PowerSeries. xs meta))
ISequential
ISeqable
(-seq [_] xs)
IFn
(-invoke [_ a] (Series. (power-series-value xs a) nil))
IPrintWithWriter
(-pr-writer [this writer _]
(write-all writer
"#object[emmy.series.PowerSeries \""
(.toString this)
"\"]"))]))
#object[emmy.series$eval88236$__GT_PowerSeries__88242 0x671bd847 "
emmy.series$eval88236$__GT_PowerSeries__88242@671bd847"
]
#?(:clj
(defmethod print-method PowerSeries [^PowerSeries s ^java.io.Writer w]
(.write w (str "#object[emmy.series.PowerSeries \""
(.toString s)
"\"]"))))
#object[clojure.lang.MultiFn 0x5877451f "
clojure.lang.MultiFn@5877451f"
]

Constructors

(defn series?
"Returns true if `s` is either a [[Series]] or a [[PowerSeries]], false
otherwise."
[s]
(or (instance? Series s)
(instance? PowerSeries s)))
#object[emmy.series$series_QMARK_ 0x5031a736 "
emmy.series$series_QMARK_@5031a736"
]
(defn power-series?
"Returns true if `s` is specifically a [[PowerSeries]], false otherwise."
[s]
(instance? PowerSeries s))
#object[emmy.series$power_series_QMARK_ 0x6e398b9a "
emmy.series$power_series_QMARK_@6e398b9a"
]
(defn- -make
"Takes a [[series?]]-true object and returns the appropriate, more specific
constructor."
[s]
(if (power-series? s)
->PowerSeries
->Series))
#object[emmy.series$_make 0x1af752e4 "
emmy.series$_make@1af752e4"
]
(defn- kind->make
"Takes a keyword - either `::series` or `::power-series` - and returns the
appropriate series constructor. Throws if neither of these are supplied."
[kind]
(case kind
::series ->Series
::power-series ->PowerSeries
(u/illegal (str "Unsupported kind: " kind))))
#object[emmy.series$kind__GT_make 0x2397681d "
emmy.series$kind__GT_make@2397681d"
]
(defn series*
"Given a sequence, returns a new [[Series]] object that wraps that
sequence (potentially padding its tail with zeros if it's finite)."
[prefix]
(->Series (i/->series prefix) nil))
#object[emmy.series$series_STAR_ 0x7e42634e "
emmy.series$series_STAR_@7e42634e"
]
(defn series
"Return a [[Series]] starting with the supplied values. The remainder of the
series will be filled with the zero-value corresponding to the first of the
given values.
If you have a sequence already, prefer [[series*]]."
[& prefix]
(series* prefix))
#object[emmy.series$series 0x5c855247 "
emmy.series$series@5c855247"
]
(defn power-series*
"Given a sequence, returns a new [[PowerSeries]] object that wraps that
sequence (potentially padding its tail with zeros if it's finite)."
[prefix]
(->PowerSeries (i/->series prefix) nil))
#object[emmy.series$power_series_STAR_ 0x7c455b41 "
emmy.series$power_series_STAR_@7c455b41"
]
(defn power-series
"Return a [[PowerSeries]] starting with the supplied values. The remainder of
the series will be filled with the zero-value corresponding to the first of
the given values.
If you have a sequence already, prefer [[power-series*]]."
[& prefix]
(power-series* prefix))
#object[emmy.series$power_series 0x648a0354 "
emmy.series$power_series@648a0354"
]
(def ^:private s-zero (series* [0]))
(0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 10+ more elided)
(def ^:private s-one (series* [1]))
(1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 10+ more elided)
(def ^:private s-identity (series* [0 1]))
(0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 10+ more elided)

These exposed objects are PowerSeries instances, because the concepts of zero, one and identity don't make sense unless you interpret these as coefficients on a power series.

(def zero
"[[PowerSeries]] instance representing the constant 0."
(power-series* [0]))
(0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 10+ more elided)
(def one
"[[PowerSeries]] instance representing the constant 1."
(power-series* [1]))
(1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 10+ more elided)
(def identity
"[[PowerSeries]] instance representing the identity function."
(power-series* [0 1]))
(0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 10+ more elided)
(defn constant
"Returns a [[PowerSeries]] representing the supplied constant term.
Optionally, pass `kind` of either `::series` or `::power-series` to specify
the type of series returned."
([c] (power-series* [c]))
([c kind] ((kind->make kind) (i/->series [c]) nil)))
#object[emmy.series$constant 0x9bf4498 "
emmy.series$constant@9bf4498"
]
(defn xpow
"Returns a [[PowerSeries]] instance representing $x^n$."
[n]
{:pre [(>= n 0)]}
(power-series* (concat (repeat n 0) [1])))
#object[emmy.series$xpow 0x7f5e3882 "
emmy.series$xpow@7f5e3882"
]
(defn generate
"Returns a [[PowerSeries]] generated by `(f i)` for `i` in `0, 1, ...`
Optionally, pass `kind` of either `::series` or `::power-series` to specify
the type of series returned."
([f] (->PowerSeries (map f (range)) nil))
([f kind]
((kind->make kind) (map f (range)) nil)))
#object[emmy.series$generate 0x2997eb7b "
emmy.series$generate@2997eb7b"
]
(defn ->function
"Accepts a [[Series]] or [[PowerSeries]] and coerces the input to
a [[PowerSeries]] without any application. Returns the coerced [[PowerSeries]]
instance.
Supplying a non-series will throw."
[s]
(cond (power-series? s) s
(series? s) (->PowerSeries (seq s) (meta s))
:else (u/illegal "non-series provided to ->function.")))
#object[emmy.series$__GT_function 0x1d33cdbe "
emmy.series$__GT_function@1d33cdbe"
]

To go the other way we need Taylor's theorem to give us a power series:

(defn function->
"Returns a [[PowerSeries]] representing the [Taylor
series](https://en.wikipedia.org/wiki/Taylor_series) expansion of `f` at the
point specified by `xs`. Multiple arguments are allowed. If no arguments `xs`
are supplied, the expansion point defaults to 0.
The expansion at 0 is also called a 'Maclaurin series'.
NOTE: this function takes derivatives internally, so if you pass a function
make sure you require [[emmy.calculus.derivative]] to install the
derivative implementation for functions. If you pass some other callable,
differentiable function-like thing, like a polynomial, this is not necessary.
NOTE: The typical definition of a Taylor series of `f` expanded around some
point `x` is
$$T(p) = f(x) + \\frac{f'(x)}{1!}(p-x) + \\frac{f''(x)}{2!} (p-x)^2 + \\ldots,$$
where `p` is the evaluation point. When `(= p x)`, all derivatives of the
Taylor series expansion of `f` will exactly match the derivatives of `f`
itself.
The Taylor series returned here (call it $T'$) is actually a function of `dx`,
where
$$T'(dx) = T(x+dx) = f(x) + \\frac{f'(x)}{1!}(dx) + \\frac{f''(x)}{2!} (dx)^2 + \\ldots.$$"
([f] (function-> f 0))
([f & xs]
(letfn [(gen [i f fact-n]
(lazy-seq
(cons (g// (apply f xs) fact-n)
(gen (inc i)
(g/partial-derivative f [])
(* fact-n i)))))]
(->PowerSeries (gen 1 f 1) nil))))
#object[emmy.series$function__GT_ 0x6e5bc104 "
emmy.series$function__GT_@6e5bc104"
]

Application

Given some power series Loading..., we can "apply" the series to a value Loading... by multiplying each entry Loading... by Loading...:

(defn- power-series-value
"Evaluates the power series, and converts it back down to a normal series."
[f x]
(let [one (g/one-like x)
powers (iterate #(g/* x %) one)]
(map g/* f powers)))
#object[emmy.series$power_series_value 0x60cf7e65 "
emmy.series$power_series_value@60cf7e65"
]

Once a PowerSeries has been applied, it becomes a Series.

What does it mean to apply a Series? The concept only makes sense if the series contains "applicables", or objects that can act as functions themselves.

If it does, then application of a series to some argument list xs means applying each series element to xs.

One further wrinkle occurs if the applicable in some position returns a series. value should combine all of these resulting series, with each series shifted by its initial position in the first series. Concretely, suppose that Loading... has the form:

Loading...

Then, this series applied to x should yield the series of values (A1, (+ A2 B1), (+ A3 B2 C1), ...)

Here's the implementation:

(defn- series-value [f xs]
(letfn [(collect [f]
(let [result (apply (first f) xs)]
(if (series? result)
(lazy-seq
(let [[r & r-tail] result]
(cons r (i/seq:+ r-tail (collect (rest f))))))
;; note that we have already realized first-result,
;; so it does not need to be behind lazy-seq.
(cons result (lazy-seq (collect (rest f)))))))]
(collect (seq f))))
#object[emmy.series$series_value 0x76158d34 "
emmy.series$series_value@76158d34"
]
(defn value
"Returns the value of the supplied [[Series]] or [[PowerSeries]] applied to `xs`.
If a [[PowerSeries]] is supplied, `xs` (despite its name) must be a single
value. Returns a [[Series]] generated by multiplying each `i`th term in `s` by
$x^i$, where $x$ is the `xs` argument.
If a [[Series]] `s` is supplied:
Assumes that `s` is a series of applicables of arity equal to the count of
`xs`. If, in fact, `s` is a series of series-valued applicables, then the
result will be a sort of layered sum of the values.
Concretely, suppose that `s` has the form:
```
[x => [A1 A2 A3...], x => [B1 B2 B3...], x => [C1 C2 C3...], ...]
```
Then, this series applied to x will yield the new series:
```
[A1 (+ A2 B1) (+ A3 B2 C1) ...]
```
The way to think about this is, that if a power series has some other series
as the coefficient of the $x^n$ term, the series must shift by $n$ positions
before being added into the final total."
[s xs]
(cond (power-series? s) (power-series-value s xs)
(series? s) (series-value s xs)
:else
(u/illegal (str "value only works on `Series` or `PowerSeries`; received " s))))
#object[emmy.series$value 0x7548e075 "
emmy.series$value@7548e075"
]
(defn fmap
"Returns a new series generated by applying the supplied `f` to each element in
the input series `s`. The returned series will be the same type as the input
series, either [[Series]] or [[PowerSeries]].
NOTE scmutils calls this `series:elementwise`."
[f s]
((-make s) (map f s) (meta s)))
#object[emmy.series$fmap 0x45ab5c39 "
emmy.series$fmap@45ab5c39"
]
(defn inflate
"Accepts an input series `s` and an exponent `n`, and expands the series in the
`n`th power of its argument. Every term `i` maps to position `i*n`, with zeros
padded in the new missing slots.
For example:
```clojure
(inflate identity 3)
;; => (series 0 0 0 1)
(take 6 (inflate (generate inc) 3))
;; => (1 0 2 0 3 0)
```
NOTE this operation makes sense as described for a [[PowerSeries]], where each
entry represents the coefficient of some power of `x`; functionally it still
works with [[Series]] objects."
[s n]
(if (<= n 1)
s
(let [zero (g/zero-like (first s))
zeros (repeat (dec n) zero)]
((-make s)
(->> (map cons s (repeat zeros))
(apply concat))
(meta s)))))
#object[emmy.series$inflate 0x5e8c2ebf "
emmy.series$inflate@5e8c2ebf"
]
(defn partial-sums
"Returns a series (of the same type as the input) of partial sums of the terms
in the supplied series `s`."
[s]
((-make s) (reductions g/+ s) (meta s)))
#object[emmy.series$partial_sums 0x42a6477 "
emmy.series$partial_sums@42a6477"
]
(defn sum
"Returns the sum of all elements in the input series `s` up to order
`n` (inclusive). For example:
```clojure
(sum (series 1 1 1 1 1 1 1) 3)
;; => 4
```
NOTE that [[sum]] sums the first `n + 1` terms, since a series starts with an
order 0 term."
[s n]
(transduce (take (inc n)) g/+ s))
#object[emmy.series$sum 0x365367cf "
emmy.series$sum@365367cf"
]

Power Series Specific Functions

(defn compose
"Returns a new [[PowerSeries]] $U$ that represents the composition of the two
input power series $S$ and $T$, where $U$ evaluates like:
```
$$U(x) = S(T(x))$$
```"
[s t]
{:pre [(power-series? s)
(power-series? t)]}
(->PowerSeries (i/compose (seq s) (seq t))
nil))
#object[emmy.series$compose 0x2986bfa8 "
emmy.series$compose@2986bfa8"
]
(defn revert
"Returns a new [[PowerSeries]] $U$ that represents the compositional inverse (the
'reversion') of the input power series $S$, satisfying:
```
$$S(U(x)) = x$$
```"
[s]
{:pre [(power-series? s)]}
(->PowerSeries (i/revert (seq s))
(meta s)))
#object[emmy.series$revert 0x16c085f1 "
emmy.series$revert@16c085f1"
]
(defn integral
"Returns a [[PowerSeries]] $U$ that represents the definite integral of the
input power series $S$ with constant term $c$:
```
$$U = c + \\int_0^{\\infty} S$$
```"
([s] (integral s 0))
([s constant]
{:pre [(power-series? s)]}
(->PowerSeries (i/integral (seq s) constant)
(meta s))))
#object[emmy.series$integral 0x419e5282 "
emmy.series$integral@419e5282"
]
(defn arg-scale
"Given a univariate [[PowerSeries]] and a singleton sequence of `factors`,
returns a new [[PowerSeries]] that scales its argument by `(first factor)` on
application.
Given a [[Series]], recursively applies [[arg-scale]] to each element, making
this ONLY appropriate in its current form for a [[Series]] of [[PowerSeries]]
instances."
[s factors]
(if (power-series? s)
(do (assert (= (count factors) 1) "Only univariate [[PowerSeries]] are allowed.")
(compose s (power-series* [0 (first factors)])))
(fmap #(arg-scale % factors) s)))
#object[emmy.series$arg_scale 0x6c5e2861 "
emmy.series$arg_scale@6c5e2861"
]
(defn arg-shift
"Given a univariate [[PowerSeries]] and a singleton sequence of `shifts`,
returns a function that, when applied, returns a value equivalent to calling
the original `s` with its argument shifted by `(first shifts)`.
NOTE: [[arg-shift]] can't return a [[PowerSeries]] instance because the
implementation of [[compose]] does not currently allow a constant element in
the right-hand series.
Given a [[Series]], recursively applies [[arg-shift]] to each element, making
this ONLY appropriate in its current form for a [[Series]] of [[PowerSeries]]
instances. Returns a [[Series]] of functions."
[s shifts]
(if (power-series? s)
(do (assert (= (count shifts) 1) "Only univariate [[PowerSeries]] are allowed.")
(apply f/arg-shift s shifts))
(fmap #(arg-shift % shifts) s)))
#object[emmy.series$arg_shift 0x70b9b1bd "
emmy.series$arg_shift@70b9b1bd"
]

Built In Series

The following section defines a number of built in series that come up often enough to be included. There are, of course, far more! Please feel free to open a PR if you have a series you think should be included.

(def exp-series (->PowerSeries i/expx nil))
(1 1 1/2 1/6 1/24 1/120 1/720 1/5040 1/40320 1/362880 1/3628800 1/39916800 1/479001600 1/6227020800 1/87178291200 1/1307674368000 1/20922789888000 1/355687428096000 1/6402373705728000 1/121645100408832000 10+ more elided)
(def sin-series (->PowerSeries i/sinx nil))
(0 1 0 -1/6 0 1/120 0 -1/5040 0 1/362880 0 -1/39916800 0 1/6227020800 0 -1/1307674368000 0 1/355687428096000 0 -1/121645100408832000 10+ more elided)
(def cos-series (->PowerSeries i/cosx nil))
(1 0 -1/2 0 1/24 0 -1/720 0 1/40320 0 -1/3628800 0 1/479001600 0 -1/87178291200 0 1/20922789888000 0 -1/6402373705728000 0 10+ more elided)
(def tan-series (->PowerSeries i/tanx nil))
(0 1 0 1/3 0 2/15 0 17/315 0 62/2835 0 1382/155925 0 21844/6081075 0 929569/638512875 0 6404582/10854718875 0 443861162/1856156927625 10+ more elided)
(def sec-series (->PowerSeries i/secx nil))
(1 0 1/2 0 5/24 0 61/720 0 277/8064 0 50521/3628800 0 540553/95800320 0 199360981/87178291200 0 3878302429/4184557977600 0 2404879675441/6402373705728000 0 10+ more elided)
(def asin-series (->PowerSeries i/asinx nil))
(0 1 0 1/6 0 3/40 0 5/112 0 35/1152 0 63/2816 0 231/13312 0 143/10240 0 6435/557056 0 12155/1245184 10+ more elided)
(def acos-series (->PowerSeries i/acosx nil))
(1.5707963267948966 -1 0 -1/6 0 -3/40 0 -5/112 0 -35/1152 0 -63/2816 0 -231/13312 0 -143/10240 0 -6435/557056 0 -12155/1245184 10+ more elided)
(def atan-series (->PowerSeries i/atanx nil))
(0 1 0 -1/3 0 1/5 0 -1/7 0 1/9 0 -1/11 0 1/13 0 -1/15 0 1/17 0 -1/19 10+ more elided)
(def acot-series (->PowerSeries i/acotx nil))
(1.5707963267948966 -1 0 1/3 0 -1/5 0 1/7 0 -1/9 0 1/11 0 -1/13 0 1/15 0 -1/17 0 1/19 10+ more elided)
(def sinh-series (->PowerSeries i/sinhx nil))
(0 1 0 1/6 0 1/120 0 1/5040 0 1/362880 0 1/39916800 0 1/6227020800 0 1/1307674368000 0 1/355687428096000 0 1/121645100408832000 10+ more elided)
(def cosh-series (->PowerSeries i/coshx nil))
(1 0 1/2 0 1/24 0 1/720 0 1/40320 0 1/3628800 0 1/479001600 0 1/87178291200 0 1/20922789888000 0 1/6402373705728000 0 10+ more elided)
(def tanh-series (->PowerSeries i/tanhx nil))
(0 1 0 -1/3 0 2/15 0 -17/315 0 62/2835 0 -1382/155925 0 21844/6081075 0 -929569/638512875 0 6404582/10854718875 0 -443861162/1856156927625 10+ more elided)
(def asinh-series (->PowerSeries i/asinhx nil))
(0 1 0 -1/6 0 3/40 0 -5/112 0 35/1152 0 -63/2816 0 231/13312 0 -143/10240 0 6435/557056 0 -12155/1245184 10+ more elided)
(def atanh-series (->PowerSeries i/atanhx nil))
(0 1 0 1/3 0 1/5 0 1/7 0 1/9 0 1/11 0 1/13 0 1/15 0 1/17 0 1/19 10+ more elided)
(def log1+x-series (->PowerSeries i/log1+x nil))
(0 1 -1/2 1/3 -1/4 1/5 -1/6 1/7 -1/8 1/9 -1/10 1/11 -1/12 1/13 -1/14 1/15 -1/16 1/17 -1/18 1/19 10+ more elided)
(def log1-x-series (->PowerSeries i/log1-x nil))
(0 -1 -1/2 -1/3 -1/4 -1/5 -1/6 -1/7 -1/8 -1/9 -1/10 -1/11 -1/12 -1/13 -1/14 -1/15 -1/16 -1/17 -1/18 -1/19 10+ more elided)
(defn binomial-series
"Returns a [[PowerSeries]] instance representing a
[Binomial series](https://en.wikipedia.org/wiki/Binomial_series), i.e., the
taylor series of the function $f$ given by
```
$$f(x) = (1 + x)^\\alpha$$
```"
[alpha]
(->PowerSeries (i/binomial alpha) nil))
#object[emmy.series$binomial_series 0x62445d56 "
emmy.series$binomial_series@62445d56"
]

Series (vs PowerSeries)

These are interesting sequences, not taylor series, but nice to have in a library like Emmy.

(def fib-series (->Series i/fib nil))
(0 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 10+ more elided)
(def catalan-series (->Series i/catalan nil))
(1 1 2 5 14 42 132 429 1430 4862 16796 58786 208012 742900 2674440 9694845 35357670 129644790 477638700 1767263190 10+ more elided)
(def harmonic-series (->Series i/harmonic nil))
(1 3/2 11/6 25/12 137/60 49/20 363/140 761/280 7129/2520 7381/2520 83711/27720 86021/27720 1145993/360360 1171733/360360 1195757/360360 2436559/720720 42142223/12252240 14274301/4084080 275295799/77597520 55835135/15519504 10+ more elided)
(def bell-series (->Series i/bell nil))
(1 2 5 15 52 203 877 4140 21147 115975 678570 4213597 27644437 190899322 1382958545 10480142147 82864869804 682076806159 5832742205057 51724158235372 10+ more elided)

Generic Implementations

This section installs Series and PowerSeries into the Emmy generic arithmetic system.

A key idea here is that all "coefficients" of a series must be some kind derived from ::coseries. This is not true in the Scheme scmutils library; in that library, anything that responds false to series? is game for interaction with series objects.

NOTE This might be the right way to go. Feel free to experiment.

(derive ::v/scalar ::coseries)
nil
(derive ::v/function ::coseries)
nil

All generic methods:

  • unwrap the supplied series instances,
  • operate directly on their backing sequences
  • repackage the result using the appropriate constructor

This section does not define methods that coerce Series => PowerSeries or vice versa. Users should do this explicitly.

(doseq [[ctor kind] [[->Series ::series]
[->PowerSeries ::power-series]]]
(defmethod g/zero? [kind] [_] false)
(defmethod g/one? [kind] [_] false)
(defmethod g/identity? [kind] [_] false)
(defmethod g/add [kind kind] [s t]
(ctor (i/seq:+ (seq s) (seq t)) nil))
(defmethod g/add [::coseries kind] [c s]
(ctor (i/c+seq c (seq s))
(meta s)))
(defmethod g/add [kind ::coseries] [s c]
(ctor (i/seq+c (seq s) c)
(meta s)))
(defmethod g/negate [kind] [s]
(ctor (i/negate (seq s))
(meta s)))
(defmethod g/sub [kind kind] [s t]
(ctor (i/seq:- (seq s) (seq t)) nil))
(defmethod g/sub [::coseries kind] [c s]
(ctor (i/c-seq c (seq s))
(meta s)))
(defmethod g/sub [kind ::coseries] [s c]
(ctor (i/seq-c (seq s) c)
(meta s)))
(defmethod g/mul [kind kind] [s t]
(ctor (i/seq:* (seq s) (seq t)) nil))
(defmethod g/mul [::coseries kind] [c s]
(ctor (i/c*seq c (seq s))
(meta s)))
(defmethod g/mul [kind ::coseries] [s c]
(ctor (i/seq*c (seq s) c)
(meta s)))
(defmethod g/square [kind] [s]
(let [xs (seq s)]
(ctor (i/seq:* xs xs)
(meta s))))
(defmethod g/cube [kind] [s]
(let [xs (seq s)]
(ctor (i/seq:* (i/seq:* xs xs) xs)
(meta s))))
(defmethod g/expt [kind ::v/native-integral] [s e]
(ctor (i/expt (seq s) e)
(meta s)))
(defmethod g/invert [kind] [s]
(ctor (i/invert (seq s))
(meta s)))
(defmethod g/div [::coseries kind] [c s]
(ctor (i/c-div-seq c (seq s))
(meta s)))
(defmethod g/div [kind ::coseries] [s c]
(ctor (i/seq-div-c (seq s) c)
(meta s)))
(defmethod g/div [kind kind] [s t]
(ctor (i/div (seq s) (seq t)) nil))
(defmethod g/solve-linear-right [::coseries kind] [c s] (g/div c s))
(defmethod g/solve-linear-right [kind ::coseries] [s c] (g/div s c))
(defmethod g/solve-linear-right [kind kind] [s t] (g/div s t))
;; `g/solve-linear` acts identically to `g/div` with arguments reversed.
(defmethod g/solve-linear [::coseries kind] [c s] (g/div s c))
(defmethod g/solve-linear [kind ::coseries] [s c] (g/div c s))
(defmethod g/solve-linear [kind kind] [s t] (g/div t s))
(defmethod g/sqrt [kind] [s]
(ctor (i/sqrt (seq s))
(meta s)))
(defmethod g/simplify [kind] [s]
(fmap g/simplify s)))
nil

Power Series Generic Extensions

A PowerSeries is a single variable function; we extend the following methods to PowerSeries in the same style as they're extended for functions. Each of the following act like function composition, and compose their operation with the function represented by the PowerSeries.

If s is a PowerSeries that applies as (s x), (g/exp s) returns a series that represents (g/exp (s x)).

(defmethod g/exp [::power-series] [s]
(->PowerSeries (i/compose i/expx (seq s))
(meta s)))
#object[clojure.lang.MultiFn 0x7168c357 "
clojure.lang.MultiFn@7168c357"
]
(defmethod g/cos [::power-series] [s]
(->PowerSeries (i/compose i/cosx (seq s))
(meta s)))
#object[clojure.lang.MultiFn 0x3be6105f "
clojure.lang.MultiFn@3be6105f"
]
(defmethod g/sin [::power-series] [s]
(->PowerSeries (i/compose i/sinx (seq s))
(meta s)))
#object[clojure.lang.MultiFn 0x163caf8e "
clojure.lang.MultiFn@163caf8e"
]
(defmethod g/tan [::power-series] [s]
(->PowerSeries (i/compose i/tanx (seq s))
(meta s)))
#object[clojure.lang.MultiFn 0x428252da "
clojure.lang.MultiFn@428252da"
]
(defmethod g/sec [::power-series] [s]
(->PowerSeries (i/compose i/secx (seq s))
(meta s)))
#object[clojure.lang.MultiFn 0x3b1846b4 "
clojure.lang.MultiFn@3b1846b4"
]
(defmethod g/asin [::power-series] [s]
(->PowerSeries (i/compose i/asinx (seq s))
(meta s)))
#object[clojure.lang.MultiFn 0x43c45b34 "
clojure.lang.MultiFn@43c45b34"
]
(defmethod g/acos [::power-series] [s]
(->PowerSeries (i/compose i/acosx (seq s))
(meta s)))
#object[clojure.lang.MultiFn 0x1380d514 "
clojure.lang.MultiFn@1380d514"
]
(defmethod g/atan [::power-series] [s]
(->PowerSeries (i/compose i/atanx (seq s))
(meta s)))
#object[clojure.lang.MultiFn 0x6a540c36 "
clojure.lang.MultiFn@6a540c36"
]
(defmethod g/acot [::power-series] [s]
(->PowerSeries (i/compose i/acotx (seq s))
(meta s)))
#object[clojure.lang.MultiFn 0x7c9d3014 "
clojure.lang.MultiFn@7c9d3014"
]
(defmethod g/cosh [::power-series] [s]
(->PowerSeries (i/compose i/coshx (seq s))
(meta s)))
#object[clojure.lang.MultiFn 0x8dc9785 "
clojure.lang.MultiFn@8dc9785"
]
(defmethod g/sinh [::power-series] [s]
(->PowerSeries (i/compose i/sinhx (seq s))
(meta s)))
#object[clojure.lang.MultiFn 0x1d61d948 "
clojure.lang.MultiFn@1d61d948"
]
(defmethod g/tanh [::power-series] [s]
(->PowerSeries (i/compose i/tanhx (seq s))
(meta s)))
#object[clojure.lang.MultiFn 0x42231f14 "
clojure.lang.MultiFn@42231f14"
]
(defmethod g/asinh [::power-series] [s]
(->PowerSeries (i/compose i/asinhx (seq s))
(meta s)))
#object[clojure.lang.MultiFn 0x2d117535 "
clojure.lang.MultiFn@2d117535"
]
(defmethod g/atanh [::power-series] [s]
(->PowerSeries (i/compose i/atanhx (seq s))
(meta s)))
#object[clojure.lang.MultiFn 0x38f36f68 "
clojure.lang.MultiFn@38f36f68"
]

Derivatives

For a Series, the derivative operation assumes that the series contains applicables that can take their own partial derivatives.

(defmethod g/partial-derivative [::series v/seqtype] [^Series s selectors]
(->Series (map #(g/partial-derivative % selectors)
(.-xs s))
(.-m s)))
#object[clojure.lang.MultiFn 0x3d9e83a6 "
clojure.lang.MultiFn@3d9e83a6"
]

A PowerSeries is itself a single-variable function, so g/partial-derivative simply takes the series derivative of the contained sequence.

(defmethod g/partial-derivative [::power-series v/seqtype] [^PowerSeries s selectors]
(if (empty? selectors)
(->PowerSeries (i/deriv (.-xs s))
(.-m s))
(u/illegal
(str "Cannot yet take partial derivatives of a power series: " s selectors))))
#object[clojure.lang.MultiFn 0x3d9e83a6 "
clojure.lang.MultiFn@3d9e83a6"
]
(defmethod g/zero-like [::power-series] [_] zero)
#object[clojure.lang.MultiFn 0x79337717 "
clojure.lang.MultiFn@79337717"
]
(defmethod g/one-like [::power-series] [_] one)
#object[clojure.lang.MultiFn 0x46a24adf "
clojure.lang.MultiFn@46a24adf"
]
(defmethod g/identity-like [::power-series] [_] identity)
#object[clojure.lang.MultiFn 0x787144b0 "
clojure.lang.MultiFn@787144b0"
]
(defmethod g/zero-like [::series] [_] s-zero)
#object[clojure.lang.MultiFn 0x79337717 "
clojure.lang.MultiFn@79337717"
]
(defmethod g/one-like [::series] [_] s-one)
#object[clojure.lang.MultiFn 0x46a24adf "
clojure.lang.MultiFn@46a24adf"
]

This is suspect, since [[Series]], unlike [[PowerSeries]], are general infinite sequences and not necessarily interpreted as polynomials. This decision follows scmutils convention.

(defmethod g/identity-like [::series] [_] s-identity)
#object[clojure.lang.MultiFn 0x787144b0 "
clojure.lang.MultiFn@787144b0"
]
(defmethod g/exact? [::series] [_] false)
#object[clojure.lang.MultiFn 0x1864e04e "
clojure.lang.MultiFn@1864e04e"
]
(defmethod g/exact? [::power-series] [_] false)
#object[clojure.lang.MultiFn 0x1864e04e "
clojure.lang.MultiFn@1864e04e"
]
(defmethod g/freeze [::power-series] [^PowerSeries s]
(let [prefix (->> (g/simplify (take 4 (.-xs s)))
(g/freeze)
(filter (complement g/zero?))
(map-indexed
(fn [n a]
(if (g/one? a)
`(~'expt ~'_ ~n)
`(~'* ~a (~'expt ~'_ ~n))))))]
`(~'+ ~@prefix ~'...)))
#object[clojure.lang.MultiFn 0x7a0db79a "
clojure.lang.MultiFn@7a0db79a"
]
(defmethod g/freeze [::series] [^Series s]
(let [prefix (g/freeze
(g/simplify (take 4 (.-xs s))))]
`(~'+ ~@prefix ~'...)))
#object[clojure.lang.MultiFn 0x7a0db79a "
clojure.lang.MultiFn@7a0db79a"
]