Implement .Rat

TODO

@@ -27,7 +27,7 @@ EASY
    Fudge and run your favorite spectest file.
 
    Stuff spectests are blocking on: "Block", "&hash", "closure for",
-   "ranges of chars", "gather for", constants in signatures, ::T, Int, ...
+   "ranges of chars", "gather for", constants in signatures, ::T, ...
 
 MEDIUM
 
@@ -92,9 +92,6 @@ HARD
    the optimized builtins to use nominal type checking.  Bonus points
    if the optimizer can turn if $x == any(1,3,5) into a loop.
 
-  *Create multiple numeric types with sensible overloads between them.  Add
-   more numeric operators.
-
    Niecza is probably the ideal platform for prototyping a concurrency
    system which will whirlpool back into the synopses.
 

lib/Builtins.cs

@@ -214,6 +214,46 @@ class SubstrLValue: Variable {
         else return Kernel.BoxAnyMO<BigInteger>(v, Kernel.IntMO);
     }
 
+    public static void GetAsRational(Variable v,
+            out BigInteger num, out BigInteger den) {
+        P6any n = v.Fetch().mo.mro_Numeric.Get(v).Fetch();
+        int rk = GetNumRank(n);
+
+        if (rk == NR_COMPLEX || rk == NR_FLOAT) {
+            double dbl = 0;
+            if (rk == NR_COMPLEX) {
+                Complex c = Kernel.UnboxAny<Complex>(n);
+                if (c.im != 0)
+                    throw new NieczaException("Complex cannot be used here");
+                dbl = c.re;
+            } else {
+                dbl = Kernel.UnboxAny<double>(n);
+            }
+            ulong bits = (ulong)BitConverter.DoubleToInt64Bits(dbl);
+            num = (bits & ((1UL << 52) - 1)) + (1UL << 52);
+            den = (1UL << 52);
+            if ((bits & (1UL << 63)) != 0) num = -num;
+            int power = ((int)((bits >> 52) & 0x7FF)) - 0x3FF;
+            if (power > 0) num <<= power;
+            else den <<= -power;
+            SimplifyFrac(ref num, ref den);
+        }
+        else if (rk == NR_FATRAT) {
+            FatRat r = Kernel.UnboxAny<FatRat>(n);
+            num = r.num; den = r.den;
+        }
+        else if (rk == NR_FIXRAT) {
+            Rat r = Kernel.UnboxAny<Rat>(n);
+            num = r.num; den = r.den;
+        }
+        else if (rk == NR_BIGINT) {
+            num = Kernel.UnboxAny<BigInteger>(n); den = BigInteger.One;
+        }
+        else {
+            num = Kernel.UnboxAny<int>(n); den = BigInteger.One;
+        }
+    }
+
     public static void SimplifyFrac(ref BigInteger num, ref BigInteger den) {
         if (den.Sign < 0) {
             den = -den;
@@ -664,6 +704,27 @@ class SubstrLValue: Variable {
         return MakeInt(~r1);
     }
 
+    public static Variable RatApprox(Variable v1, Variable v2) {
+        NominalCheck("$x", Kernel.AnyMO, v1);
+        NominalCheck("$y", Kernel.AnyMO, v2);
+
+        BigInteger nc, dc, ne, de, na, da;
+        GetAsRational(v1, out nc, out dc);
+        GetAsRational(v2, out ne, out de);
+
+        RatApproxer.Simplest(nc*de-ne*dc,dc*de,nc*de+ne*dc,dc*de,out na,out da);
+        SimplifyFrac(ref na, ref da);
+
+        // since the user controls the denominator size here, use FatRat freely
+        // XXX: is it appropriate to return FatRat from a method named Rat?
+        ulong sda;
+        if (da.AsUInt64(out sda)) {
+            return MakeFixRat(na,da);
+        } else {
+            return MakeFatRat(na,da);
+        }
+    }
+
     public static Variable PostIncrement(Variable v) {
         P6any o1 = NominalCheck("$x", Kernel.AnyMO, v);
         AssignV(v, o1.mo.mro_succ.Get(v));

lib/CLRBackend.cs

@@ -3388,6 +3388,7 @@ class NamProcessor {
             thandlers["bif_chars"] = SimpleB("Chars");
             thandlers["bif_substr3"] = SimpleB("Substr3");
             thandlers["bif_simple_eval"] = SimpleB("SimpleEval");
+            thandlers["bif_rat_approx"] = SimpleB("RatApprox");
 
             thandlers["bif_defined"] = Contexty("mro_defined");
             thandlers["bif_bool"] = Contexty("mro_Bool");

lib/CORE.setting

@@ -64,6 +64,8 @@ my class Regex  { ... }
 my class Num    { ... }
 my class Str    { ... }
 my class Cool {
+    method Rat($eps = 1e-6) { Q:CgOp { (bif_rat_approx {self} {$eps}) } }
+
     method grep(Mu $sm) { grep $sm, @(self) }
     method map($func)   { map $func, @(self) }
     method for (&cb) {

lib/Utils.cs

@@ -528,4 +528,99 @@ public sealed class FatRat {
             this.num = num; this.den = den;
         }
     }
+
+    public sealed class RatApproxer {
+        // This could be optimized a fair amount if it becomes necessary.
+        // In particular, due to common prefixes the code is doing about
+        // four times as much work as it has to; also, it's probably not
+        // necessary to actually store the continued fractions in memory.
+        // http://en.wikipedia.org/w/index.php?title=Continued_fraction&oldid=429751049#Best_rational_within_an_interval
+
+        static List<BigInteger> GetContinuedFraction(BigInteger n, BigInteger d,
+                bool extend) {
+            List<BigInteger> r = new List<BigInteger>();
+            while (d.Sign != 0) {
+                BigInteger nn;
+                r.Add(BigInteger.DivRem(n, d, out nn));
+                n = d;
+                d = nn;
+            }
+            if (extend) {
+                r[r.Count-1]--;
+                r.Add(BigInteger.One);
+            }
+            return r;
+        }
+
+        // entry invariant: 0 < l < h
+        static void CandidateSimplest(BigInteger numl, BigInteger denl, bool xl,
+                BigInteger numh, BigInteger denh, bool xh,
+                ref BigInteger snum, ref BigInteger sden) {
+            List<BigInteger> cfl = GetContinuedFraction(numl, denl, xl);
+            List<BigInteger> cfh = GetContinuedFraction(numh, denh, xh);
+            List<BigInteger> cfo = new List<BigInteger>();
+
+            int i = 0;
+            while (i != cfl.Count && i != cfh.Count && cfl[i] == cfh[i])
+                cfo.Add(cfl[i++]);
+
+            if (i != cfh.Count && (i == cfl.Count || cfl[i] > cfh[i]))
+                cfo.Add(cfh[i] + BigInteger.One);
+            else
+                cfo.Add(cfl[i] + BigInteger.One);
+
+            BigInteger onum = cfo[i];
+            BigInteger oden = BigInteger.One;
+
+            for (i--; i >= 0; i--) {
+                BigInteger t = onum; onum = oden; oden = t;
+                onum += oden * cfo[i];
+            }
+
+            if (oden < sden) {
+                sden = oden;
+                snum = onum;
+            }
+        }
+
+        public static void Simplest(BigInteger numl, BigInteger denl,
+                BigInteger numh, BigInteger denh,
+                out BigInteger numo, out BigInteger deno) {
+            int signl = numl.Sign * denl.Sign;
+            int signh = numh.Sign * denh.Sign;
+
+            if (signl * signh <= 0) {
+                // interval includes 0, automatically simplest
+                numo = BigInteger.Zero;
+                deno = BigInteger.One;
+                return;
+            }
+
+            if (signl < 0) {
+                numl = -numl; numh = -numh;
+            }
+
+            BigInteger gl = BigInteger.GreatestCommonDivisor(numl, denl);
+            numl /= gl; denl /= gl;
+            BigInteger gh = BigInteger.GreatestCommonDivisor(numh, denh);
+            numh /= gh; denh /= gh;
+
+            BigInteger snum = numl;
+            BigInteger sden = denl;
+
+            if (denh < denl) {
+                snum = numh;
+                sden = denh;
+            }
+
+            for (int k = 0; k < 4; k++)
+                CandidateSimplest(numl, denl, (k&1)!=0, numh, denh, (k&2)!=0, ref snum, ref sden);
+
+            if (signl < 0)
+                snum = -snum;
+
+            numo = snum;
+            deno = sden;
+        }
+    }
 }

src/niecza

@@ -75,6 +75,7 @@ method outerlex ($n) { self._cgop("outerlex", $n) }
 method newtypedscalar ($t) { self._cgop("newtypedscalar", $t) }
 method bif_rand () { self._cgop("bif_rand") }
 method exactnum ($base, $digits) { self._cgop("exactnum", $base, $digits) }
+method bif_rat_approx ($center, $eps) { self._cgop("bif_rat_approx", $center, $eps) }
 }
 
 my $usage = q:to/EOM/;