From b2c61f510f595ffdc8bfc494f8a6ca268b1e6783 Mon Sep 17 00:00:00 2001 From: Stefan O'Rear Date: Sun, 22 May 2011 03:48:04 -0700 Subject: [PATCH] Implement .Rat --- TODO | 5 +-- lib/Builtins.cs | 61 ++++++++++++++++++++++++++++++ lib/CLRBackend.cs | 1 + lib/CORE.setting | 2 + lib/Utils.cs | 95 +++++++++++++++++++++++++++++++++++++++++++++++ src/niecza | 1 + 6 files changed, 161 insertions(+), 4 deletions(-) diff --git a/TODO b/TODO index d2fc0ebf..d630c44e 100644 --- a/TODO +++ b/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. diff --git a/lib/Builtins.cs b/lib/Builtins.cs index aa07d9d0..efd5b9cf 100644 --- a/lib/Builtins.cs +++ b/lib/Builtins.cs @@ -214,6 +214,46 @@ class SubstrLValue: Variable { else return Kernel.BoxAnyMO(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(n); + if (c.im != 0) + throw new NieczaException("Complex cannot be used here"); + dbl = c.re; + } else { + dbl = Kernel.UnboxAny(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(n); + num = r.num; den = r.den; + } + else if (rk == NR_FIXRAT) { + Rat r = Kernel.UnboxAny(n); + num = r.num; den = r.den; + } + else if (rk == NR_BIGINT) { + num = Kernel.UnboxAny(n); den = BigInteger.One; + } + else { + num = Kernel.UnboxAny(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)); diff --git a/lib/CLRBackend.cs b/lib/CLRBackend.cs index eb23d7be..3ac2358c 100644 --- a/lib/CLRBackend.cs +++ b/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"); diff --git a/lib/CORE.setting b/lib/CORE.setting index 1aab0340..c58c64cf 100644 --- a/lib/CORE.setting +++ b/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) { diff --git a/lib/Utils.cs b/lib/Utils.cs index 2f543461..6e78310d 100644 --- a/lib/Utils.cs +++ b/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 GetContinuedFraction(BigInteger n, BigInteger d, + bool extend) { + List r = new List(); + 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 cfl = GetContinuedFraction(numl, denl, xl); + List cfh = GetContinuedFraction(numh, denh, xh); + List cfo = new List(); + + 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; + } + } } diff --git a/src/niecza b/src/niecza index c3493677..1cef62c0 100644 --- a/src/niecza +++ b/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/;