From 27ada6566b248882c1c625ae961353d1f6032c8d Mon Sep 17 00:00:00 2001
From: Aaron Fenyes <aaron.fenyes@fareycircles.ooo>
Date: Sat, 10 Aug 2024 19:11:55 -0700
Subject: [PATCH] Scala benchmark: step rotation by multiplying

This makes the algorithm more consistent with the Rust benchmark.
---
 .../src/main/scala/CircularLawApp.scala       | 52 +++++++++++--------
 1 file changed, 31 insertions(+), 21 deletions(-)

diff --git a/lang-trials/scala-benchmark/src/main/scala/CircularLawApp.scala b/lang-trials/scala-benchmark/src/main/scala/CircularLawApp.scala
index 950e661..12a0a00 100644
--- a/lang-trials/scala-benchmark/src/main/scala/CircularLawApp.scala
+++ b/lang-trials/scala-benchmark/src/main/scala/CircularLawApp.scala
@@ -2,6 +2,7 @@ import com.raquo.laminar.api.L.{*, given}
 import narr.*
 import org.scalajs.dom
 import org.scalajs.dom.document
+import scala.collection.mutable.ArrayBuffer
 import scala.math.{cos, sin}
 import slash.matrix.Matrix
 import slash.matrix.decomposition.Eigen
@@ -38,36 +39,45 @@ object CircularLawApp:
       )
       ctx.fill()
   
-  def eigvalsRotated[N <: Int](A: Matrix[N, N], time: Double)(using ValueOf[N]): (NArray[Double], NArray[Double]) =
-    // create transformation
-    val maxFreq = 4
-    val T = Matrix.identity[N, N]
-    val dim: Int = valueOf[N]
-    for n <- 0 to dim by 2 do
-      val a = cos(math.Pi * time * (n % maxFreq))
-      val b = sin(math.Pi * time * (n % maxFreq))
-      T(n, n) = a
-      T(n+1, n) = b
-      T(n, n+1) = -b
-      T(n+1, n+1) = a
-    
-    // find eigenvalues
-    val eigen = Eigen(T*A)
+  def complexEigenvalues[N <: Int](mat: Matrix[N, N])(using ValueOf[N]): (NArray[Double], NArray[Double]) =
+    val eigen = Eigen(mat)
     (
       eigen.realEigenvalues.asInstanceOf[NArray[Double]],
       eigen.imaginaryEigenvalues.asInstanceOf[NArray[Double]]
     )
   
-  def randEigvalSeries[N <: Int]()(using ValueOf[N]): (List[(NArray[Double], NArray[Double])], String) =
-    val timeRes = 100
-    val dim: Int = valueOf[N]
+  def randEigvalSeries[N <: Int]()(using ValueOf[N]): (ArrayBuffer[(NArray[Double], NArray[Double])], String) =
+    // start timing
     val startTime = System.currentTimeMillis()
-    val A = new Matrix[N, N](
+    
+    // initialize the random matrix step
+    val dim: Int = valueOf[N]
+    var randMat = new Matrix[N, N](
       NArray.tabulate(dim*dim)(k => (math.E*k*k) % 2 - 1)
     ).times(math.sqrt(3d / dim))
-    val series = List.tabulate(timeRes)(t => eigvalsRotated(A, t.toDouble / timeRes))
+    
+    // initialize the rotation step
+    val timeRes = 100
+    val maxFreq = 4
+    val rotStep = Matrix.identity[N, N]
+    for n <- 0 to dim by 2 do
+      val ang = math.Pi * (n % maxFreq) / timeRes
+      val cos_ang = cos(ang)
+      val sin_ang = sin(ang)
+      rotStep(n, n) = cos_ang
+      rotStep(n+1, n) = sin_ang
+      rotStep(n, n+1) = -sin_ang
+      rotStep(n+1, n+1) = cos_ang
+    
+    // find the eigenvalues
+    val eigvalSeries = ArrayBuffer(complexEigenvalues(randMat))
+    for _ <- 1 to timeRes-1 do
+      randMat = rotStep * randMat
+      eigvalSeries += complexEigenvalues(randMat)
+    
+    // finish timing
     val runTime = System.currentTimeMillis() - startTime
-    (series, runTime.toString() + " ms")
+    (eigvalSeries, runTime.toString() + " ms")
   
   def main(args: Array[String]): Unit =
     ctx.fillStyle = "white"