Replace bigfraction.js with CLJS code (#155)

CHANGELOG.md

@@ -2,6 +2,26 @@
 
 ## [unreleased]
 
+- #155
+
+  - Replace the implementation of arbitrary-precision rational arithmetic
+    provided by `fraction.js` with Clojure code, allowing us to remove
+    a JavaScript dependency.
+
+  - We note here a subtlety with fraction reader syntax. When the CLJS
+    compiler runs and encounters input like  `#emmy/ratio 1/2`, since
+    the compilation occurs in a JVM Clojure environment the 1/2 will
+    deserialize as a Ratio, which is then transformed by the reader into
+    code which will generate a ClojureScript `Fraction`. The ClojureScript
+    reader, however, evaluates the input as q JS expression resulting
+    in 0.5, which is then rationalized with an expensive algorithm
+    which may lose exactness because of the unwanted floating-point
+    conversion. For this reason, when we emit a ratio in string form,
+    we quote the arguments: `#emmy/ratio "1/2"` to prevent the conversion.
+    Consider this a deprecation notice for the unquoted form. We also
+    now allow the form `#emmy/ratio [1 2]`, like Complex, and now no
+    longer allow an initial `+` on the numerator.
+
 - #154:
 
   - Adds `emmy.tape` with an implementation of reverse-mode automatic

src/emmy/bigfraction.cljs

@@ -0,0 +1,295 @@
+#_"SPDX-License-Identifier: GPL-3.0"
+
+(ns emmy.bigfraction
+  "A CLJS bigfraction is a coprime pair of JavaScript BigInts, with the
+   sign carried in the numerator."
+  (:refer-clojure :exclude [abs zero? neg?])
+  (:require
+   [clojure.core :as core]
+   [goog.array :as garray]))
+
+(def ^:private ZERO (js/BigInt 0))
+(def ^:private ONE (js/BigInt 1))
+(def ^:private TEN (js/BigInt 10))
+(def ^:private -ONE (- ONE))
+
+(defn- bigint?
+  "Returns true if x is a BigInt. There is a similar function in [[emmy.util]],
+   but we prefer that this library avoid that dependency."
+  [x]
+  (= "bigint" (goog/typeOf x)))
+
+(declare eq cmp ->real)
+
+(deftype Fraction [^js/BigInt n ^js/BigInt d]
+
+  Object
+  (valueOf [this] (->real this))
+  (toString [_] (str n "/" d))
+
+  IHash
+  (-hash [x]
+    (bit-xor
+     (-hash (.-n x))
+     (-hash (.-d x))))
+
+  IComparable
+  (-compare [this other]
+    (cond (instance? Fraction other)
+          (cmp this other)
+
+          (bigint? other)
+          (compare n (* d other))
+
+          :else
+          (let [o-value (.valueOf other)]
+            (garray/defaultCompare this o-value))))
+
+)
+
+(def ^:private F_ONE (Fraction. ONE ONE))
+
+(defn division-by-zero
+  "Throws JS exception used to signal an attempt to construct a fraction
+   with zero denominator."
+  []
+  (throw (js/Error "Fraction with zero denominator")))
+
+(defn numerator [^Fraction x]
+  (.-n x))
+
+(defn denominator [^Fraction x]
+  (.-d x))
+
+(defn zero?
+  "Returns true iff `x` is a zero fraction."
+  [^Fraction x]
+  (let [a (.-n x)]
+    (== ZERO a)))
+
+(defn one?
+  "Returns true iff `x` is a unit fraction."
+  [^Fraction x]
+  (let [a (.-n x)
+        b (.-d x)]
+    (== a b)))
+
+(defn eq [^Fraction x ^Fraction y]
+  (let [a (.-n x)
+        b (.-d x)
+        c (.-n y)
+        d (.-d y)]
+    (and (== a c) (== b d))))
+
+(defn bigint-gcd
+  "GCD assuming a and b are BigInts > 0"
+  [^js/BigInt a ^js/BigInt b]
+  (loop [a a
+         b b]
+    (if (== b ZERO) a
+        (recur b (js-mod a b)))))
+
+(defn- bigint-abs [a]
+  (if (< a 0) (* -ONE a) a))
+
+(defn integer->
+  "Create a fraction with unit denominator."
+  [n]
+  (->Fraction (js/BigInt n) ONE))
+
+(defn ->real
+  "Coerce a fraction to real by performing the division
+   in the floating point domain"
+  [^Fraction q]
+  (/ (js/Number (.-n q)) (js/Number (.-d q))))
+
+(defn ->normal-form
+  "We assume we are given two BigInts, with b > 0. The GCD is divided out, and the
+   sign is carried in the numerator."
+  [^js/BigInt a ^js/BigInt b]
+  (when (== ZERO b)
+    (division-by-zero))
+  (let [an (< a 0)
+        a (bigint-abs a)
+        bn (< b 0)
+        b (bigint-abs b)
+        g (bigint-gcd a b)
+        neg (not= an bn)
+        abs_c (/ a g)
+        c (if neg (* -ONE abs_c) abs_c)
+        d (/ b g)]
+    (->Fraction c d)))
+
+(defn make
+  "Produces a fraction in canonical form. Note that the canonical form of an integer is
+   an integer, so if `(one? b)` you just get a."
+  [a b]
+  (let [a (js/BigInt a)
+        b (js/BigInt b)]
+    (when (== 0 b)
+      (division-by-zero))
+    (->normal-form a b)))
+
+(defn- js-expt
+  "Use the `js*` escape clause to get access to the JavaScript `**` operator,
+   which can exponentiate two BigInts exactly."
+  [a b]
+  (js* "(~{} ** ~{})" a b))
+
+(def ^:private double-re #"(-?\d+)(\.(\d+))?([Ee]([+-]\d+))?")
+
+(defn real->
+  "Clojure converts the real to BigDecimal and rationalizes from that.
+   The JVM documentation explains that the BigDecimal value will correspond
+   to what would be printed for the double value. We attempt to do the
+   same thing here by converting to a string and converting from there."
+  [x]
+  (let [s (.toString x)]
+    (if-let [[_ int _ frac _ exp] (re-matches double-re s)]
+      (let [scale (- (js/parseInt (or exp "0"))
+                     (count (or frac "")))
+            scale-neg (< scale 0)
+            mantissa (js/BigInt (str int frac))
+            exponent (js/BigInt (js-expt TEN (js/BigInt (Math/abs scale))))]
+        (if scale-neg
+          (->normal-form mantissa exponent)
+          (Fraction. (* mantissa exponent) ONE)))
+      (throw (js/Error (str "Cannot convert " x " to ratio."))))))
+
+(defn abs
+  "Absolute value of the fraction `x`."
+  [^Fraction x]
+  (let [n (.-n x)
+        d (.-d x)]
+    (if (< n 0) (->Fraction (* -ONE n) d) x)))
+
+(defn neg
+  "Negation of the fraction `x`."
+  [^Fraction x]
+  (let [n (.-n x)
+        d (.-d x)]
+    (->Fraction (* -ONE n) d)))
+
+(defn neg?
+  "True if $x<0$."
+  [^Fraction x]
+  (< (.-n x) 0))
+
+(defn add
+  "Returns the sum of `x` and `y`."
+  [^Fraction x ^Fraction y]
+  (let [a (.-n x)
+        b (.-d x)
+        c (.-n y)
+        d (.-d y)]
+    (->normal-form (+ (* a d) (* b c)) (* b d))))
+
+(defn sub
+  "Returns the difference of `x` and `y`."
+  [^Fraction x ^Fraction y]
+  (let [a (.-n x)
+        b (.-d x)
+        c (.-n y)
+        d (.-d y)]
+    (->normal-form (core/- (core/* a d) (core/* b c)) (core/* b d))))
+
+(defn mul
+  "Returns the product of `x` and `y`."
+  [^Fraction x ^Fraction y]
+  (let [a (.-n x)
+        b (.-d x)
+        c (.-n y)
+        d (.-d y)]
+    (->normal-form (core/* a c) (core/* b d))))
+
+(defn div
+  "Returns the quotient of `x` and `y`."
+  [^Fraction x ^Fraction y]
+  (let [a (.-n x)
+        b (.-d x)
+        c (.-n y)
+        d (.-d y)]
+    (->normal-form (core/* a d) (core/* b c))))
+
+(defn invert
+  "Returns the reciprocal of `x`, but throws if $x=0$."
+  [^Fraction x]
+  (let [a (.-n x)
+        b (.-d x)
+        neg (< a 0)]
+    (when (== ZERO a)
+      (division-by-zero))
+    (if neg
+      (->Fraction (* -ONE b) (* -ONE a))
+      (->Fraction b a))))
+
+(defn cmp
+  "Compares the fractions `x` and `y`, returning -1, 0, or 1."
+  [^Fraction x ^Fraction y]
+  (let [a (.-n x)
+        b (.-d x)
+        c (.-n y)
+        d (.-d y)
+        s (- (* a d) (* b c))]
+    (cond (< s 0) -1
+          (> s 0) 1
+          :else 0)))
+
+(defn integer-power
+  "Raises the fraction `x` to the integer power `n`."
+  [^Fraction x n]
+  (let [a (.-n x)
+        b (.-d x)
+        N (js/BigInt n)]
+    (cond
+      (= n ZERO) F_ONE
+      (= n ONE) x
+      (> N ZERO) (->normal-form (js-expt a N) (js-expt b N))
+      :else (->normal-form (js-expt b (- N)) (js-expt a (- N))))))
+
+(defn promote
+  "If the fraction has a unit denominator, return the numerator, else the fraction."
+  [^Fraction x]
+  (if (== ONE (.-d x)) (.-n x) x))
+
+(defn ceil
+  "Ceiling. Result is a BigInt."
+  [^Fraction x]
+  (let [a (.-n x)
+        b (.-d x)]
+    (+ (/ a b) (if (or (< a 0)
+                       (== (js-mod a b) ZERO))
+                 ZERO
+                 ONE))))
+
+(defn floor
+  "Floor. Result is a BigInt."
+  [^Fraction x]
+  (let [a (.-n x)
+        b (.-d x)]
+    (- (/ a b) (if (or (> a 0)
+                       (== (js-mod a b) ZERO))
+                 ZERO
+                 ONE))))
+
+(defn quotient
+  "Fractions form a field, so somewhat dubiously the function returns
+   the largest integer N for which $Ny\\le x$."
+  [x y]
+  (let [z (div x y)]
+    (if (neg? z)
+      (ceil z)
+      (floor z))))
+
+(defn remainder
+  "If $q$ is `(quotient x y)`, returns $x-qy$."
+  [x y]
+  (sub x (mul (Fraction. (quotient x y) ONE) y)))
+
+(defn gcd
+  [^Fraction x ^Fraction y]
+  (let [a (.-n x)
+        b (.-d x)
+        c (.-n y)
+        d (.-d y)]
+    (abs (->normal-form (* (bigint-gcd a c) (bigint-gcd b d)) (* b d)))))

src/emmy/complex/impl.cljc

@@ -20,7 +20,7 @@
             [emmy.util :as u]
             [emmy.value :as v]))
 
-(declare equal?)
+(declare equal? ->string)
 
 (deftype Complex [re im]
   v/IKind
@@ -32,23 +32,27 @@
   (numerical? [_] true)
 
   #?@(:clj [Object
-            (equals [a b] (equal? a b))]
-      :cljs [IEquiv
+            (equals [a b] (equal? a b))
+            (toString [a] (->string a))]
+
+      :cljs [Object
+             (toString [a] (->string a))
+
+             IEquiv
              (-equiv [a b] (equal? a b))
 
              IPrintWithWriter
              (-pr-writer
               [z writer _]
-              (write-all
-               writer
-               "#emmy/complex "
-               (str [(.-re z) (.-im z)])))]))
+              (write-all writer (->string z)))]))
 
 #?(:clj
-   (defmethod print-method Complex [^Complex v ^java.io.Writer w]
-     (.write w (str "#emmy/complex "
-                    [(.-re v)
-                     (.-im v)]))))
+   (defmethod print-method Complex [^Complex z ^java.io.Writer w]
+     (.write w (->string z))))
+
+(defn- ->string
+  [^Complex c]
+  (str "#emmy/complex " [(.-re c) (.-im c)]))
 
 (def ZERO (Complex. 0 0))
 (def ONE (Complex. 1 0))