Skip to content

Commit a09e3b0

Browse files
authored
Taylor series using streams (#974)
A nice example of streams: computing using infinite Taylor series. ```scala def sine() = integrate(0.0) {cosine} def cosine() = for[Double] { integrate(-1.0) {sine} } { x => do emit(x.neg) } >> println(eval(100, 2.0) {sine}) // this example 0.9092974268256817 >> println(sin(2.0)) // stdlib 0.9092974268256817 ```
1 parent f175258 commit a09e3b0

File tree

2 files changed

+72
-0
lines changed

2 files changed

+72
-0
lines changed
Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
Cons(1, Cons(2, Cons(3, Cons(4, Cons(5, Nil())))))
2+
Cons(1, Cons(1, Cons(0.5, Cons(0.16667, Cons(0.04167, Nil())))))
3+
7.38906
4+
0.9093
5+
-0.41615
6+
0.84147
7+
0.5403
8+
0
9+
1
10+
0
11+
-1
Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
import stream
2+
3+
/// Warning: This currently has the wrong time complexity.
4+
def ints(): Unit / emit[Int] = {
5+
do emit(1)
6+
for[Int] {ints} { x => do emit(x + 1) }
7+
}
8+
9+
/// Take the `n` first elements of a `stream`, put them in a list.
10+
def take[A](n: Int) { stream: () => Unit / emit[A] }: List[A] = {
11+
with collectList[A]
12+
with boundary
13+
with limit[A](n + 1)
14+
stream()
15+
}
16+
17+
// ∫x^n := 1/(n+1) * x^(n + 1) + C
18+
// => integral = C + all original Taylor series terms shifted by one and divided by n + 1
19+
def integrate(c: Double) { taylor: () => Unit / emit[Double] } = {
20+
do emit(c)
21+
zip[Double, Int] {taylor} {ints} { (f, n) => do emit(f / n.toDouble) }
22+
}
23+
24+
// e^x := ∫_{0}^x e^y dy + 1
25+
def exponential(): Unit / emit[Double] =
26+
integrate(1.0) {exponential}
27+
28+
// evaluate a taylor series, taking `steps`
29+
def eval(steps: Int, x: Double) { taylor: () => Unit / emit[Double] } = {
30+
var acc = 0.0
31+
for[Double] { take[Double](steps) {taylor}.reverse.each } { a =>
32+
acc = a + acc * x
33+
}
34+
acc
35+
}
36+
37+
// sin x := ∫_{0}^x cos y dy
38+
def sine(): Unit / emit[Double] = integrate(0.0) {cosine}
39+
40+
// cos x := - ∫_{0}^x sin y dy + 1
41+
// = -(∫_{0}^x sin y dy - 1)
42+
def cosine(): Unit / emit[Double] =
43+
for[Double] { integrate(-1.0) {sine} } { x => do emit(x.neg) }
44+
45+
def main() = {
46+
def rnd(one: Double) = one.round(5)
47+
def rnd(many: List[Double]) = many.map { one => rnd(one) }
48+
49+
50+
println(take[Int](5) {ints})
51+
println(take[Double](5) {exponential}.rnd)
52+
println(eval(20, 2.0) {exponential}.rnd)
53+
println(eval(20, 2.0) {sine}.rnd)
54+
println(eval(20, 2.0) {cosine}.rnd)
55+
println(eval(20, 1.0) {sine}.rnd)
56+
println(eval(20, 1.0) {cosine}.rnd)
57+
println(eval(20, 0.0) {sine}.rnd)
58+
println(eval(20, 0.0) {cosine}.rnd)
59+
println(eval(20, PI) {sine}.rnd)
60+
println(eval(20, PI) {cosine}.rnd)
61+
}

0 commit comments

Comments
 (0)