diff --git a/examples/stdlib/stream/fastexp.check b/examples/stdlib/stream/fastexp.check
new file mode 100644
index 000000000..44d2d4cc8
--- /dev/null
+++ b/examples/stdlib/stream/fastexp.check
@@ -0,0 +1,52 @@
+0
+0b0
+0b0
+1
+0b1
+0b1
+2
+0b01
+0b10
+3
+0b11
+0b11
+4
+0b001
+0b100
+123456789
+0b101010001011001111011010111
+0b111010110111100110100010101
+
+2 ^ 0
+1
+1
+2 ^ 5
+32
+32
+2 ^ 20
+1048576
+1048576
+10 ^ 10
+10000000000
+10000000000
+0 ^ 1
+0
+0
+1 ^ 0
+1
+1
+0 ^ 0
+1
+1
+200 ^ 2
+40000
+40000
+128 ^ 4
+268435456
+268435456
+7 ^ 13
+96889010407
+96889010407
+13 ^ 7
+62748517
+62748517
\ No newline at end of file
diff --git a/examples/stdlib/stream/fastexp.effekt b/examples/stdlib/stream/fastexp.effekt
new file mode 100644
index 000000000..f7d36ece8
--- /dev/null
+++ b/examples/stdlib/stream/fastexp.effekt
@@ -0,0 +1,95 @@
+import stream
+
+/// Least significant bits
+def eachLSB(n: Int): Unit / emit[Bool] = {
+  if (n <= 0) {
+    do emit(false)
+  } else {
+    var tmp = n
+    while (tmp > 0) {
+      do emit(if (tmp.mod(2) == 1) true else false)
+      tmp = tmp / 2
+    }
+  }
+}
+
+/// Most significant bits
+def eachMSB(n: Int): Unit / emit[Bool] = {
+  if (n <= 0) {
+    do emit(false)
+  } else {
+    // First, calculate the highest bit position
+    var highest = 0
+    var m = n
+    while (m > 0) {
+      highest = highest + 1
+      m = m / 2
+    }
+
+    // Then emit from MSB to LSB
+    for[Int] { range(0, highest) } { i =>
+      val index = highest - i
+      val bit = n.bitwiseShr(index - 1).mod(2) == 1
+      do emit(bit)
+    }
+  }
+}
+
+def fastexp(n: Int, k: Int) = product {
+  stream::zip[Int, Bool] {n.iterate { x => x * x }} {k.eachLSB} {
+    case res, true => do emit(res)
+    case res, false => ()
+  }
+}
+
+def main() = {
+  /// Computes n^exp for integers. (from #923!)
+  def pow(n: Int, exp: Int): Int = {
+    def go(n: Int, exp: Int, acc: Int): Int = {
+      if (exp == 0) {
+        acc
+      } else if (mod(exp, 2) == 0) {
+        go(n * n, exp / 2, acc)
+      } else {
+        go(n * n, exp / 2, acc * n)
+      }
+    }
+    go(n, exp, 1)
+  }
+
+  def prettyBits(bits: List[Bool]): String =
+    "0b" ++ bits.map { b => if (b) "1" else "0"}.join("")
+
+  def testBits(n: Int) = {
+    println(n)
+    println(collectList[Bool] {n.eachLSB}.prettyBits)
+    println(collectList[Bool] {n.eachMSB}.prettyBits)
+  }
+
+  testBits(0)
+  testBits(1)
+  testBits(2)
+  testBits(3)
+  testBits(4)
+  testBits(123456789)
+
+  println("")
+
+  def testExp(n: Int, k: Int) = {
+    println(show(n) ++ " ^ " ++ show(k))
+    println(n.pow(k))
+    println(n.fastexp(k))
+  }
+
+  testExp(2, 0)
+  testExp(2, 5)
+  testExp(2, 20)
+  testExp(10, 10)
+  testExp(0, 1)
+  testExp(1, 0)
+  testExp(0, 0)
+  testExp(200, 2)
+  testExp(128, 4)
+  testExp(7, 13)
+  testExp(13, 7)
+}
diff --git a/libraries/common/stream.effekt b/libraries/common/stream.effekt
index 72b32c6ee..dc26eefdf 100644
--- a/libraries/common/stream.effekt
+++ b/libraries/common/stream.effekt
@@ -173,10 +173,23 @@ def replicate[A](number: Int) { action: () => A }: Unit / emit[A] =
     replicate(number - 1) {action}
   }
 
+/// Creates an infinite iterated stream given by an `initial` seed and a `step` function:
+/// iterate(a){f} ~> a, f(a), f(f(a)), f(f(f(a))), ...
+def iterate[A](initial: A) { step: A => A }: Unit / emit[A] = {
+  var current = initial
+  while (true) {
+    do emit(current)
+    current = step(current)
+  }
+}
+
 
-def sum { stream : () => Unit / emit[Int] }: Int =
+def sum { stream: () => Unit / emit[Int] }: Int =
   returning::sum[Unit]{stream}.second
 
+def product { stream: () => Unit / emit[Int] }: Int =
+  returning::product[Unit]{stream}.second
+
 def collectList[A] { stream: () => Unit / emit[A] }: List[A] =
   returning::collectList[A, Unit]{stream}.second
 
@@ -461,7 +474,7 @@ def limit[A, R](number: Int) { stream: () => R / emit[A] }: R / { emit[A], stop
   }
 }
 
-def sum[R] { stream : () => R / emit[Int] }: (R, Int) = {
+def sum[R] { stream: () => R / emit[Int] }: (R, Int) = {
   var s = 0;
   try {
     (stream(), s)
@@ -471,6 +484,16 @@ def sum[R] { stream : () => R / emit[Int] }: (R, Int) = {
   }
 }
 
+def product[R] { stream: () => R / emit[Int] }: (R, Int) = {
+  var s = 1;
+  try {
+    (stream(), s)
+  } with emit[Int] { v =>
+    s = s * v;
+    resume(())
+  }
+}
+
 def collectList[A, R] { stream: () => R / emit[A] }: (R, List[A]) =
   try {
     (stream(), Nil())