Skip to content

Commit

Permalink
MathR LCM, IDiv, IMod, DivRem, benchmarks
Browse files Browse the repository at this point in the history
  • Loading branch information
@c-ohle
c-ohle committed Jul 9, 2022
1 parent 67ac466 commit 920464f
Show file tree
Hide file tree
Showing 127 changed files with 5,223 additions and 1,103 deletions.
32 changes: 22 additions & 10 deletions System.Numerics.Rational/BigRational.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1008,7 +1008,7 @@ public sealed partial class CPU
/// <param name="capacity">The initial stack capacity that must be greater than 0 (zero).</param>
public CPU(uint capacity = 32)
{
p = new uint[capacity][]; Debug.Assert(capacity > 0);
p = new uint[capacity][];
}
/// <summary>
/// Returns a temporary absolute index of the current stack top
Expand Down Expand Up @@ -2647,6 +2647,8 @@ public void free()
{
cpu = null;
}

#if false //todo: remove, div(a, b); mod(); swp(); pop(); must!!! have the same performance
/// <summary>
/// Performs an integer division of the numerators of the first two values on top of the stack<br/>
/// and replaces them with the integral result.
Expand Down Expand Up @@ -2690,6 +2692,7 @@ public void imod()
}
pop();
}
#endif

uint[] rent(uint n)
{
Expand Down Expand Up @@ -2812,19 +2815,28 @@ static void mul(uint* a, uint* b, uint* r)
uint f = b[1];
if (nb == 1)
{
if (f == 0) { *(ulong*)r = 1; return; }
if (f == 1) { r[0] = na; for (uint i = 1; i <= na; i++) r[i] = a[i]; return; }
switch (f) //jump table
{
case 0: *(ulong*)r = 1; return;
case 1: copy(r, a, na + 1); r[0] = na; return;
//case 2: //todo: opt. shl(1)
//case 3: //todo: opt. + +
//case 4: //todo: opt. shl(2)
//...
}
//if (f == 0) { *(ulong*)r = 1; return; }
//if (f == 1) { copy(r, a, na + 1); r[0] = na; /* r[0] = na; for (uint i = 1; i <= na; i++) r[i] = a[i]; */ return; }
}
//if (na == 1) *(ulong*)(r + 1) = (ulong)a[1] * b[1];
//else
if (na == 1) *(ulong*)(r + 1) = (ulong)a[1] * b[1]; //todo: jt
else
if (na == 2 && Bmi2.X64.IsSupported)
{
*(ulong*)(r + 3) = Bmi2.X64.MultiplyNoFlags(*(ulong*)(a + 1), nb == 2 ? *(ulong*)(b + 1) : b[1], (ulong*)(r + 1));
}
else if (na >= 0020) //nb <= na
{
var n = na + nb; for (uint i = 1; i <= n; i++) r[i] = 0; //todo: remove, kmu opt will make it needles
kmu(a + 1, na, b + 1, nb, r + 1, n); //r[0] = n; while (r[r[0]] == 0) r[0]--; return;
kmu(a + 1, na, b + 1, nb, r + 1, n);
}
else
{
Expand All @@ -2849,7 +2861,7 @@ static void mul(uint* a, uint* b, uint* r)
}
if (r[r[0] = na + nb] == 0 && r[0] > 1) { r[0]--; Debug.Assert(!(r[r[0]] == 0 && r[0] > 1)); }
}
static void sqr(uint* a, uint* r)
static void sqr(uint* a, uint* r) //todo: opt kmu for sqr
{
uint n = *a++ & 0x3fffffff; var v = r + 1;
if (n == 1)
Expand Down Expand Up @@ -3001,8 +3013,8 @@ internal static void shr(uint* p, int c)
}
static uint* gcd(uint* u, uint* v)
{
//var su = clz(u); if (su != 0) shr(u, su);
//var sv = clz(v); if (sv != 0) shr(v, sv); if (su > sv) su = sv;
var su = clz(u); if (su != 0) shr(u, su);
var sv = clz(v); if (sv != 0) shr(v, sv); if (su > sv) su = sv;
if (cms(u, v) < 0) { var t = u; u = v; v = t; }
while (v[0] > 2)
{
Expand Down Expand Up @@ -3060,7 +3072,7 @@ internal static void shr(uint* p, int c)
while (y != 0) { var t = unchecked((uint)x) % unchecked((uint)y); x = y; y = t; }
*(ulong*)(u + 1) = x; u[0] = u[2] != 0 ? 2u : 1u; m1:;
}
//if (su != 0) shl(u, su);
if (su != 0) shl(u, su);
return u;
}
static int clz(uint* p)
Expand Down
97 changes: 97 additions & 0 deletions System.Numerics.Rational/MathR.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -493,6 +493,103 @@ public static BigRational Atan2(BigRational y, BigRational x, int digits = 0)
return default(BigRational) / 0; //NaN
}
/// <summary>
/// Finds the greatest common divisor (GCD) of two <see cref="BigRational"/> integer values.
/// </summary>
/// <remarks>
/// This operation makes only sense for integer values.
/// </remarks>
/// <param name="a">The first value.</param>
/// <param name="b">The second value.</param>
/// <returns>The greatest common divisor of <paramref name="a"/> and <paramref name="b"/>.</returns>
public static BigRational GreatestCommonDivisor(BigRational a, BigRational b)
{
var cpu = BigRational.task_cpu; cpu.push(a); cpu.push(b);
cpu.gcd(); //cpu.pop(); return default;
return cpu.popr();
}
/// <summary>
/// Finds the least common multiple (LCM) of two <see cref="BigRational"/> integer values.
/// </summary>
/// <remarks>
/// This operation makes only sense for integer values.
/// </remarks>
/// <param name="a">The first value.</param>
/// <param name="b">The second value.</param>
/// <returns>The least common multiple of <paramref name="a"/> and <paramref name="b"/>.</returns>
public static BigRational LeastCommonMultiple(BigRational a, BigRational b)
{
//|a * b| / gcd(a, b) == |a / gcd(a, b) * b|
var cpu = BigRational.task_cpu; cpu.push(a); cpu.push(b);
cpu.dup(); cpu.dup(2); cpu.gcd(); cpu.div(); cpu.mul(); cpu.abs();
return cpu.popr();
}
/// <summary>
/// Performes an integer division <paramref name="a"/> / <paramref name="b"/>
/// </summary>
/// <remarks>
/// For integer values <paramref name="a"/> and <paramref name="b"/>, the result equals a <see cref="BigInteger.Divide(BigInteger, BigInteger)"/> division.<br/>
/// This in contrast to <paramref name="a"/> / <paramref name="b"/>, where a corresponding fraction results.
/// </remarks>
/// <param name="a">The value to be divided. (dividend)</param>
/// <param name="b">The value to divide by. (devisor)</param>
/// <returns>A <see cref="BigRational "/> integer value.</returns>
/// <exception cref="DivideByZeroException"><see cref="DivideByZeroException"/>: <paramref name="b"/> is zero.</exception>
public static BigRational IDiv(BigRational a, BigRational b)
{
if (BigRational.Sign(b) == 0) throw new DivideByZeroException(nameof(b)); // b.p == null
var cpu = rat.task_cpu; //cpu.push(a); cpu.push(b); cpu.idiv(); return cpu.popr();
cpu.div(a, b); cpu.mod(); cpu.swp(); cpu.pop();
return cpu.popr();
}
/// <summary>
/// Performes an integer modulo operation <paramref name="a"/> % <paramref name="b"/> what is the remainder that results from a division.
/// </summary>
/// <param name="a">The value to be divided. (dividend)</param>
/// <param name="b">The value to divide by. (divisor)</param>
/// <remarks>
/// For integer values <paramref name="a"/> and <paramref name="b"/>, the result equals a <see cref="BigInteger"/> modulo (%) operation.<br/>
/// This in contrast to <paramref name="a"/> % <paramref name="b"/>, where for <see cref="BigInteger"/> a corresponding fraction results.
/// </remarks>
/// <returns>A <see cref="BigRational "/> integer value.</returns>
/// <exception cref="DivideByZeroException"><see cref="DivideByZeroException"/>: <paramref name="b"/> is zero.</exception>
public static BigRational IMod(BigRational a, BigRational b)
{
if (BigRational.Sign(b) == 0) throw new DivideByZeroException(nameof(b)); // b.p == null
var cpu = rat.task_cpu; //cpu.push(a); cpu.push(b); cpu.imod(); var c = cpu.popr(); return c;
cpu.div(b, b); cpu.mod(); cpu.pop();
return cpu.popr();
}
/// <summary>
/// Calculates the quotient of two <see cref="BigRational"/> signed values and also returns the remainder in an output parameter.
/// </summary>
/// <remarks>
/// This function is for compatibility to <see cref="Math.DivRem(int, int, out int)"/> like functions.<br/>
/// For integer values <paramref name="a"/> and <paramref name="b"/>, the result equals a <see cref="BigInteger.DivRem(BigInteger, BigInteger, out BigInteger)"/> operation..<br/>
/// </remarks>
/// <param name="a">The dividend.</param>
/// <param name="b">The divisor.</param>
/// <param name="r">The remainder.</param>
/// <returns>The quotient of the specified numbers.</returns>
/// <exception cref="DivideByZeroException"><see cref="DivideByZeroException"/>: <paramref name="b"/> is zero.</exception>
public static BigRational DivRem(BigRational a, BigRational b, out BigRational r)
{
if (BigRational.Sign(b) == 0) throw new DivideByZeroException(nameof(b)); // b.p == null
var cpu = rat.task_cpu; cpu.div(b, b); cpu.mod();
r = cpu.popr(); return cpu.popr();
}
/// <summary>
/// Produces the quotient and the remainder of two signed <see cref="BigRational"/> numbers.
/// </summary>
/// This function is for compatibility to <see cref="Math.DivRem(int, int)"/> like functions.<br/>
/// For integer values <paramref name="a"/> and <paramref name="b"/>, the result equals a <see cref="BigInteger.DivRem(BigInteger, BigInteger, out BigInteger)"/> operation..<br/>
/// <param name="a">The dividend.</param>
/// <param name="b">The divisor.</param>
/// <returns>The quotient and the remainder of the specified numbers as integer values.</returns>
public static (BigRational Quotient, BigRational Remainder) DivRem(BigRational a, BigRational b)
{
var d = DivRem(a, b, out var r); return (d, r);
}
/// <summary>
/// Gets or sets the default number of digits used by <see cref="MathR"/> functions
/// with an optional digits parameter with default value = 0.<br/>
/// This allows for a flat interface, easily interchangeable with <see cref="Math"/> or <see cref="MathF"/>.
Expand Down
11 changes: 6 additions & 5 deletions Test/BenchmarksPage.Designer.cs
View file

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

82 changes: 33 additions & 49 deletions Test/BenchmarksPage.cs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ void textBox1_Leave(object sender, EventArgs e)
void wb_DocumentCompleted(object? sender, WebBrowserDocumentCompletedEventArgs? e)
{
var doc = wb.Document; var body = doc.Body;
var a = typeof(BigRational).Assembly; var name = a.GetName();
var a = typeof(BigRational).Assembly; var name = a.GetName();
var t1 = a.GetCustomAttributes(typeof(System.Runtime.Versioning.TargetFrameworkAttribute), false).FirstOrDefault() as System.Runtime.Versioning.TargetFrameworkAttribute;
doc.GetElementById("1").InnerText = name.Version.ToString();
doc.GetElementById("2").InnerText = a.ImageRuntimeVersion + " " + t1?.FrameworkName;
Expand All @@ -56,7 +56,7 @@ void wb_DocumentCompleted(object? sender, WebBrowserDocumentCompletedEventArgs?
esvg = XElement.Parse(doc.GetElementById("svg").InnerHtml);
}

(float[] times, int count)[]? passes;
(float[] times, int count)[]? passes; float lastav;

void runtest(HtmlElement svg, Action<int, Func<bool>> test, int ex)
{
Expand All @@ -78,21 +78,25 @@ void runtest(HtmlElement svg, Action<int, Func<bool>> test, int ex)
});
}
test(8, () => true);

if (ex == 0x00003fbf && (lastav = lastav * 0.01f) < 1) //to get realist times for cpu internal bigint div it makes sense to sub the av from mul as IDiv changed for comapibillity and has an additional mul
for (int i = 1; i < 256; i++) passes[0].times[i] *= lastav;

update_svg(svg, ex);

var a1 = passes[0].times.Take(255);//.Average();
var a2 = passes[1].times.Take(255);//.Average();
var av = MathF.Round(a1.Average() * 100 / a2.Average()).ToString();
var av = MathF.Round(a1.Average() * 100 / a2.Average()); this.lastav = av;
var mi = Math.Min(a1.Max(), a2.Max());
var ma = Math.Max(a1.Max(), a2.Max());

var info = svg.NextSibling;
var s = info.InnerHtml;
var x = s.LastIndexOf(':'); if (x > 0) s = s.Substring(0, x+1);
s += $" {av} % (min {mi} ms; max {ma} ms)";
var x = s.LastIndexOf(':'); if (x > 0) s = s.Substring(0, x + 1);
s += $" {av} % (min {mi} ms; max {ma} ms)";
info.InnerHtml = s;
s = svg.Parent.Style;
if(s.Contains("none")) svg.Parent.Style = s.Replace("none","block");// "display: block; margin-left: 16px";
if (s.Contains("none")) svg.Parent.Style = s.Replace("none", "block"); // "display: block; margin-left: 16px";
}

void button_run_Click(object sender, EventArgs e)
Expand Down Expand Up @@ -158,39 +162,19 @@ void update_svg(HtmlElement? psvg, int ex)
psvg.InnerHtml = esvg.ToString();
}

static Action<int, Func<bool>> _test_gcd(out int e)
{
var rr = new BigRational[256]; var ii = new BigInteger[256];
var cr = new BigRational[256]; var ci = new BigInteger[256];
var ra = new Random(13); //E+2040
for (int i = 0; i < 256; i++) ii[i] = (BigInteger)(rr[i] = random_rat(ra, i << 3));
e = MathR.ILog10(rr[255]); return test;
void test(int pass, Func<bool> f)
{
if (pass == 0)
for (int i = 1; f() && i < 256; i++)
cr[i] = MathRE.GreatestCommonDivisor(rr[i - 1], rr[i]);

if (pass == 1)
for (int i = 1; f() && i < 256; i++)
ci[i] = BigInteger.GreatestCommonDivisor(ii[i - 1], ii[i]);

if (pass == 8) Debug.Assert(cr.Zip(ci).All(p => p.First == p.Second));
}
}
static Action<int, Func<bool>> _test_mul(out int e)
{
var rr = new BigRational[256]; var ii = new BigInteger[256];
var cr = new BigRational[256]; var ci = new BigInteger[256];
var ra = new Random(13); //E+4080
for (int i = 0; i < 256; i++) ii[i] = (BigInteger)(rr[i] = random_rat(ra, i << 4));//<< 3
e = MathR.ILog10(rr[255]); return test;
void test(int pass, Func<bool> f)
var ra = new Random(13);
for (int i = 0; i < 256; i++) ii[i] = (BigInteger)(rr[i] = random_rat(ra, i << 4)); //E+4080
e = MathR.ILog10(rr[255]);
return (pass, f) =>
{
if (pass == 0) for (int i = 1; f() && i < 256; i++) cr[i] = rr[i - 1] * rr[i];
if (pass == 1) for (int i = 1; f() && i < 256; i++) ci[i] = ii[i - 1] * ii[i];
if (pass == 8) Debug.Assert(cr.Zip(ci).All(p => p.First == p.Second));
}
};
}
static Action<int, Func<bool>> _test_div(out int e)
{
Expand All @@ -199,29 +183,29 @@ static Action<int, Func<bool>> _test_div(out int e)
var ra = new Random(13);
for (int i = 0; i < 256; i++) ii[i] = (BigInteger)(rr[i] = random_rat(ra, i << 6));
e = MathR.ILog10(rr[255]);
return (pass, f) => //void test(int pass, Func<bool> f)
return (pass, f) =>
{
//var cpu = rat.task_cpu;
if (pass == 0) for (int i = 1; f() && i < 256; i++) cr[i] = MathRE.IntegerDivide(rr[i], rr[i - 1]);
if (pass == 0) for (int i = 1; f() && i < 256; i++) cr[i] = MathR.IDiv(rr[i], rr[i - 1]);
if (pass == 1) for (int i = 1; f() && i < 256; i++) ci[i] = ii[i] / ii[i - 1];
if (pass == 8) Debug.Assert(cr.Zip(ci).All(p => p.First == p.Second));
};
}
//static BigRational intdiv(BigRational a, BigRational b)
//{
// var cpu = rat.task_cpu;
// cpu.push(a); cpu.push(b); cpu.idiv();
// var c = cpu.popr(); return c;
//}
//static BigRational intdiv(rat.CPU cpu, BigRational a, BigRational b)
//{
// //cpu.push(a); cpu.push(b); cpu.idiv();
// cpu.push(a); cpu.idiv(b);
// //var c = cpu.popr(); return c;
// cpu.pop(); return default;
//}

static BigRational random_rat(Random rnd, int digits)
static Action<int, Func<bool>> _test_gcd(out int e)
{
var rr = new BigRational[256]; var ii = new BigInteger[256];
var cr = new BigRational[256]; var ci = new BigInteger[256];
var ra = new Random(13); //E+2040
for (int i = 0; i < 256; i++) ii[i] = (BigInteger)(rr[i] = random_rat(ra, i << 3));
e = MathR.ILog10(rr[255]);
return (pass, f) =>
{
if (pass == 0) for (int i = 1; f() && i < 256; i++) cr[i] = MathR.GreatestCommonDivisor(rr[i - 1], rr[i]);
if (pass == 1) for (int i = 1; f() && i < 256; i++) ci[i] = BigInteger.GreatestCommonDivisor(ii[i - 1], ii[i]);
if (pass == 8) Debug.Assert(cr.Zip(ci).All(p => p.First == p.Second));
};
}

static BigRational random_rat(Random rnd, int digits) //todo: make more precise
{
var cpu = rat.task_cpu;
cpu.pow(10, digits);
Expand Down
Loading

0 comments on commit 920464f

Please sign in to comment.