MkComplex and Rat structures

lib/github.com/mpllang/mpllib/MkComplex.sml

@@ -0,0 +1,111 @@
+signature COMPLEX =
+sig
+  structure R: REAL
+  type r = R.real
+
+  type t
+
+  val toString: t -> string
+
+  val make: (r * r) -> t
+  val view: t -> (r * r)
+
+  val defaultReal: Real.real -> t
+  val defaultImag: Real.real -> t
+
+  val real: r -> t
+  val imag: r -> t
+  val rotateBy: r -> t (* rotateBy x = e^(ix) *)
+
+  val zeroThreshold: r
+  val realIsZero: r -> bool
+
+  val isZero: t -> bool
+  val isNonZero: t -> bool
+
+  val zero: t
+  val i: t
+
+  val magnitude: t -> r
+
+  val ~ : t -> t
+  val - : t * t -> t
+  val + : t * t -> t
+  val * : t * t -> t
+
+  val scale: r * t -> t
+end
+
+
+functor MkComplex(R: REAL): COMPLEX =
+struct
+  structure R = R
+  open R
+  type r = real
+
+  val fromLarge = fromLarge IEEEReal.TO_NEAREST
+
+  datatype t = C of {re: real, im: real}
+
+  val rtos = fmt (StringCvt.FIX (SOME 8))
+
+  fun toString (C {re, im}) =
+    let
+      val (front, re) = if Int.< (R.sign re, 0) then ("-", R.~ re) else ("", re)
+      val (middle, im) =
+        if Int.< (R.sign im, 0) then ("-", R.~ im) else ("+", im)
+    in
+      front ^ rtos re ^ middle ^ rtos im ^ "i"
+    end
+
+  fun make (re, im) = C {re = re, im = im}
+
+  fun view (C {re, im}) = (re, im)
+
+  val zeroThreshold = fromLarge 0.00000001
+  fun realIsZero x = R.abs x < zeroThreshold
+
+  fun magnitude (C {re, im}) =
+    R.Math.sqrt (R.+ (R.* (re, re), R.* (im, im)))
+
+  fun isZero (C {re, im}) = realIsZero re andalso realIsZero im
+
+  fun isNonZero c =
+    not (isZero c)
+
+  fun rotateBy r =
+    C {re = Math.cos r, im = Math.sin r}
+
+  fun real r =
+    C {re = r, im = fromLarge 0.0}
+  fun imag i =
+    C {re = fromLarge 0.0, im = i}
+
+  fun defaultReal r =
+    real (fromLarge r)
+  fun defaultImag r =
+    imag (fromLarge r)
+
+  val zero = C {re = fromLarge 0.0, im = fromLarge 0.0}
+  val i = C {re = fromLarge 0.0, im = fromLarge 1.0}
+
+  fun neg (C {re, im}) =
+    C {re = ~re, im = ~im}
+
+  fun add (C x, C y) =
+    C {re = #re x + #re y, im = #im x + #im y}
+
+  fun sub (C x, C y) =
+    C {re = #re x - #re y, im = #im x - #im y}
+
+  fun mul (C {re = a, im = b}, C {re = c, im = d}) =
+    C {re = a * c - b * d, im = a * d + b * c}
+
+  fun scale (r, C {re, im}) =
+    C {re = r * re, im = r * im}
+
+  val ~ = neg
+  val op- = sub
+  val op+ = add
+  val op* = mul
+end

lib/github.com/mpllang/mpllib/Rat.sml

@@ -0,0 +1,145 @@
+structure Rat :>
+sig
+  type t
+  type i = IntInf.int
+
+  (* make(n, d) ~> n/d *)
+  val make: i * i -> t
+  val view: t -> i * i
+
+  val normalize: t -> t
+
+  val * : t * t -> t
+  val - : t * t -> t
+  val + : t * t -> t
+  val div: t * t -> t
+
+  val max: t * t -> t
+
+  val sign: t -> int
+  val compare: t * t -> order
+
+  val approx: t -> Real64.real
+
+  val toString: t -> string
+end =
+struct
+
+  type i = IntInf.int
+  type t = i * i
+
+  fun make (n, d) = (n, d)
+  fun view (n, d) = (n, d)
+
+  fun gcd (a, b) =
+    if b = 0 then a else gcd (b, IntInf.mod (a, b))
+
+  fun normalize (n, d) =
+    if n = 0 then
+      (0, 1)
+    else
+      let
+        val same = IntInf.sameSign (n, d)
+
+        val na = IntInf.abs n
+        val da = IntInf.abs d
+
+        val g = gcd (na, da)
+
+        val n' = IntInf.div (na, g)
+        val d' = IntInf.div (da, g)
+      in
+        if same then (n', d') else (IntInf.~ n', d')
+      end
+
+  fun mul ((a, b): t, (c, d)) = (a * c, b * d)
+
+  fun add ((a, b): t, (c, d)) =
+    (a * d + b * c, b * d)
+
+  fun sub ((a, b): t, (c, d)) =
+    (a * d - b * c, b * d)
+
+  fun divv ((a, b): t, (c, d)) = (a * d, b * c)
+
+
+  fun sign (n, d) =
+    if n = 0 then 0 else if IntInf.sameSign (n, d) then 1 else ~1
+
+
+  fun compare (r1, r2) =
+    let
+      val diff = sub (r1, r2)
+      val s = sign diff
+    in
+      if s < 0 then LESS else if s = 0 then EQUAL else GREATER
+    end
+
+
+  fun max (r1, r2) =
+    case compare (r1, r2) of
+      LESS => r2
+    | _ => r1
+
+
+  fun itor x =
+    Real64.fromLargeInt (IntInf.toLarge x)
+
+
+  (* =========================================================================
+   * approximate a rational with Real64
+   *
+   * TODO: not sure what the best way to do this is. Kinda just threw something
+   * together. It's probably kinda messed up in a subtle way.
+   *)
+
+  local
+    fun loopApprox acc (r, d) =
+      let
+        val r' = itor r
+        val d' = itor d
+      in
+        if Real64.isFinite r' andalso Real64.isFinite d' then
+          acc + r' / d'
+        else
+          let
+            val d2 = d div 2
+          in
+            if r > d2 then
+              loopApprox (acc + 0.5) (normalize (r - d2, d))
+            else
+              (* no idea how good this is... *)
+              loopApprox acc (normalize (r div 2, d div 2))
+          end
+      end
+  in
+    fun approx (n, d) =
+      let
+        val (n, d) = normalize (n, d)
+        val s = Real64.fromInt (IntInf.sign n)
+        val n = IntInf.abs n
+
+        (* abs(n/d) = m + r/d
+         * where: m is a natural number
+         *   and: r/d is a proper fraction
+         *)
+        val (m, r) = IntInf.divMod (n, d)
+        val m = itor m
+      in
+        if not (Real64.isFinite m) then s * m
+        else s * (m + loopApprox 0.0 (normalize (r, d)))
+      end
+  end
+
+  (* ======================================================================= *)
+
+  fun toString (n, d) =
+    IntInf.toString n ^ "/" ^ IntInf.toString d
+
+  (* ======================================================================= *)
+
+  val op* = mul
+  val op+ = add
+  val op- = sub
+  val op div = divv
+end

lib/github.com/mpllang/mpllib/sources.mlton.mlb

@@ -52,6 +52,8 @@ MatCOO.sml
 
 CheckSort.sml
 
+MkComplex.sml
+Rat.sml
 Geometry3D.sml
 Geometry2D.sml
 Topology2D.sml

lib/github.com/mpllang/mpllib/sources.mpl.mlb

@@ -52,6 +52,8 @@ MatCOO.sml
 
 CheckSort.sml
 
+MkComplex.sml
+Rat.sml
 Geometry3D.sml
 Geometry2D.sml
 Topology2D.sml