Bellard-Pi benchmark (#3225)

performance/build.gradle

@@ -163,6 +163,11 @@ task ring {
     dependsOn 'ring:konanRun'
 }
 
+task numerical {
+    dependsOn 'clean'
+    dependsOn 'numerical:konanRun'
+}
+
 task swiftinterop {
     dependsOn 'clean'
     dependsOn 'swiftinterop:konanRun'

performance/numerical/build.gradle.kts

@@ -0,0 +1,47 @@
+/*
+ * Copyright 2010-2019 JetBrains s.r.o. Use of this source code is governed by the Apache 2.0 license
+ * that can be found in the LICENSE file.
+ */
+
+import org.jetbrains.kotlin.benchmark.BenchmarkingPlugin
+import org.jetbrains.kotlin.ExecClang
+import org.jetbrains.kotlin.gradle.plugin.mpp.KotlinNativeTarget
+import org.jetbrains.kotlin.konan.target.HostManager
+
+plugins {
+    id("benchmarking")
+}
+
+benchmark {
+    applicationName = "Numerical"
+    commonSrcDirs = listOf("src/main/kotlin", "../../tools/benchmarks/shared/src", "../shared/src/main/kotlin", "../../endorsedLibraries/kliopt/src/main/kotlin")
+    jvmSrcDirs = listOf("src/main/kotlin-jvm", "../shared/src/main/kotlin-jvm", "../../endorsedLibraries/kliopt/src/main/kotlin-jvm")
+    nativeSrcDirs = listOf("src/main/kotlin-native", "../shared/src/main/kotlin-native/common", "../../endorsedLibraries/kliopt/src/main/kotlin-native")
+    mingwSrcDirs = listOf("../shared/src/main/kotlin-native/mingw")
+    posixSrcDirs = listOf("../shared/src/main/kotlin-native/posix")
+    linkerOpts = listOf("$buildDir/pi.o")
+}
+
+val compileLibary by tasks.creating {
+    doFirst {
+        mkdir(buildDir)
+
+        project.withConvention(ExecClang::class) {
+            execKonanClang(HostManager.host) {
+                args("-O3")
+                args("-c", "$projectDir/src/nativeInterop/cinterop/pi.c")
+                args("-o", "$buildDir/pi.o")
+            }
+        }
+    }
+}
+
+val native = kotlin.targets.getByName("native") as KotlinNativeTarget
+native.apply {
+    compilations["main"].cinterops {
+        create("cinterop") {
+            headers("$projectDir/src/nativeInterop/cinterop/pi.h")
+        }
+    }
+    binaries.getExecutable(BenchmarkingPlugin.NATIVE_EXECUTABLE_NAME, "RELEASE").linkTask.dependsOn(compileLibary)
+}

performance/numerical/gradle.properties

@@ -0,0 +1 @@
+org.jetbrains.kotlin.native.home=../../dist

performance/numerical/src/main/kotlin-jvm/launcher.kt

@@ -0,0 +1,19 @@
+/*
+ * Copyright 2010-2019 JetBrains s.r.o. and Kotlin Programming Language contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the licenses/LICENSE.txt file.
+ */
+
+import org.jetbrains.benchmarksLauncher.*
+
+actual class NumericalLauncher : Launcher() {
+    override val benchmarks = BenchmarksCollection(
+            mutableMapOf(
+                    "BellardPi" to BenchmarkEntry(::jvmBellardPi)
+            )
+    )
+}
+
+fun jvmBellardPi() {
+    for (n in 1 .. 1000 step 9)
+            pi_nth_digit(n)
+}

performance/numerical/src/main/kotlin-native/launcher.kt

@@ -0,0 +1,25 @@
+/*
+ * Copyright 2010-2019 JetBrains s.r.o. and Kotlin Programming Language contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the licenses/LICENSE.txt file.
+ */
+
+import org.jetbrains.benchmarksLauncher.*
+
+actual class NumericalLauncher : Launcher() {
+    override val benchmarks = BenchmarksCollection(
+            mutableMapOf(
+                    "BellardPi" to BenchmarkEntry(::konanBellardPi),
+                    "BellardPiCinterop" to BenchmarkEntry(::clangBellardPi)
+            )
+    )
+}
+
+fun konanBellardPi() {
+    for (n in 1 .. 1000 step 9)
+            pi_nth_digit(n)
+}
+
+fun clangBellardPi() {
+    for (n in 1 .. 1000 step 9)
+            cinterop.pi_nth_digit(n)
+}

performance/numerical/src/main/kotlin/main.kt

@@ -0,0 +1,20 @@
+/*
+ * Copyright 2010-2019 JetBrains s.r.o. and Kotlin Programming Language contributors.
+ * Use of this source code is governed by the Apache 2.0 license that can be found in the licenses/LICENSE.txt file.
+ */
+
+import org.jetbrains.benchmarksLauncher.*
+import org.jetbrains.kliopt.*
+
+expect class NumericalLauncher() : Launcher {
+}
+
+fun main(args: Array<String>) {
+    val launcher = NumericalLauncher()
+    BenchmarksRunner.runBenchmarks(args, { arguments: BenchmarkArguments ->
+        if (arguments is BaseBenchmarkArguments) {
+            launcher.launch(arguments.warmup, arguments.repeat, arguments.prefix,
+                    arguments.filter, arguments.filterRegex, arguments.verbose)
+        } else emptyList()
+    }, benchmarksListAction = launcher::benchmarksListAction)
+}

performance/numerical/src/main/kotlin/pi.kt

@@ -0,0 +1,152 @@
+/*
+ * Computation of the n'th decimal digit of \pi with very little memory.
+ * Written by Fabrice Bellard on January 8, 1997.
+ *
+ * We use a slightly modified version of the method described by Simon
+ * Plouffe in "On the Computation of the n'th decimal digit of various
+ * transcendental numbers" (November 1996). We have modified the algorithm
+ * to get a running time of O(n^2) instead of O(n^3log(n)^3).
+ */
+
+import kotlin.math.ln
+import kotlin.math.sqrt
+
+private fun mul_mod(a: Int, b: Int, m: Int)
+        = ((a.toLong() * b.toLong()) % m).toInt()
+
+/* return the inverse of x mod y */
+private fun inv_mod(x: Int, y: Int): Int {
+    var u = x
+    var v = y
+    var c = 1
+    var a = 0
+    do {
+        val q = v / u
+        var t = c
+        c = a - q * c
+        a = t
+        t = u
+        u = v - q * u
+        v = t
+    } while (u != 0)
+    a = a % y
+    if (a < 0)
+        a = y + a
+    return a
+}
+
+/* return (a^b) mod m */
+private fun pow_mod(a: Int, b: Int, m: Int): Int {
+    var b = b
+    var r = 1
+    var aa = a
+    while (true) {
+        if (b and 1 != 0)
+            r = mul_mod(r, aa, m)
+        b = b shr 1
+        if (b == 0)
+            break
+        aa = mul_mod(aa, aa, m)
+    }
+    return r
+}
+
+/* return true if n is prime */
+private fun is_prime(n: Int): Boolean {
+    if (n % 2 == 0)
+        return false
+    val r = sqrt(n.toDouble()).toInt()
+    var i = 3
+    while (i <= r) {
+        if (n % i == 0)
+            return false
+        i += 2
+    }
+    return true
+}
+
+/* return the prime number immediatly after n */
+private fun next_prime(n: Int): Int {
+    var n = n
+    do {
+        n++
+    } while (!is_prime(n))
+    return n
+}
+
+fun pi_nth_digit(n: Int): Int {
+
+    val N = ((n + 20) * ln(10.0) / ln(2.0)).toInt()
+    var sum = 0.0
+    var a = 3
+    var t: Int
+
+    while (a <= 2 * N) {
+
+        val vmax = (ln((2 * N).toDouble()) / ln(a.toDouble())).toInt()
+        var av = 1
+        var i = 0
+        while (i < vmax) {
+            av = av * a
+            i++
+        }
+
+        var s = 0
+        var num = 1
+        var den = 1
+        var v = 0
+        var kq = 1
+        var kq2 = 1
+
+        var k = 1
+        while (k <= N) {
+
+            t = k
+            if (kq >= a) {
+                do {
+                    t = t / a
+                    v--
+                } while (t % a == 0)
+                kq = 0
+            }
+            kq++
+            num = mul_mod(num, t, av)
+
+            t = 2 * k - 1
+            if (kq2 >= a) {
+                if (kq2 == a) {
+                    do {
+                        t = t / a
+                        v++
+                    } while (t % a == 0)
+                }
+                kq2 -= a
+            }
+            den = mul_mod(den, t, av)
+            kq2 += 2
+
+            if (v > 0) {
+                t = inv_mod(den, av)
+                t = mul_mod(t, num, av)
+                t = mul_mod(t, k, av)
+                i = v
+                while (i < vmax) {
+                    t = mul_mod(t, a, av)
+                    i++
+                }
+                s += t
+                if (s >= av)
+                    s -= av
+            }
+            k++
+
+        }
+
+        t = pow_mod(10, n - 1, av)
+        s = mul_mod(s, t, av)
+        sum = (sum + s.toDouble() / av.toDouble()) % 1.0
+        a = next_prime(a)
+    }
+
+    return (sum * 1e9).toInt()
+}

performance/numerical/src/nativeInterop/cinterop/cinterop.def

@@ -0,0 +1 @@
+package = cinterop

performance/numerical/src/nativeInterop/cinterop/pi.c

@@ -0,0 +1,163 @@
+/*
+ * Computation of the n'th decimal digit of \pi with very little memory.
+ * Written by Fabrice Bellard on January 8, 1997.
+ *
+ * We use a slightly modified version of the method described by Simon
+ * Plouffe in "On the Computation of the n'th decimal digit of various
+ * transcendental numbers" (November 1996). We have modified the algorithm
+ * to get a running time of O(n^2) instead of O(n^3log(n)^3).
+ *
+ * This program uses mostly integer arithmetic. It may be slow on some
+ * hardwares where integer multiplications and divisons must be done
+ * by software. We have supposed that 'int' has a size of 32 bits. If
+ * your compiler supports 'long long' integers of 64 bits, you may use
+ * the integer version of 'mul_mod' (see HAS_LONG_LONG).
+ */
+
+#include <math.h>
+
+/* uncomment the following line to use 'long long' integers */
+#define HAS_LONG_LONG
+
+#ifdef HAS_LONG_LONG
+#define mul_mod(a,b,m) (( (long long) (a) * (long long) (b) ) % (m))
+#else
+#define mul_mod(a,b,m) fmod( (double) a * (double) b, m)
+#endif
+
+/* return the inverse of x mod y */
+static int inv_mod(int x, int y)
+{
+    int q, u, v, a, c, t;
+
+    u = x;
+    v = y;
+    c = 1;
+    a = 0;
+    do {
+        q = v / u;
+
+        t = c;
+        c = a - q * c;
+        a = t;
+
+        t = u;
+        u = v - q * u;
+        v = t;
+    } while (u != 0);
+    a = a % y;
+    if (a < 0)
+        a = y + a;
+    return a;
+}
+
+/* return (a^b) mod m */
+static int pow_mod(int a, int b, int m)
+{
+    int r, aa;
+
+    r = 1;
+    aa = a;
+    while (1) {
+        if (b & 1)
+            r = mul_mod(r, aa, m);
+        b = b >> 1;
+        if (b == 0)
+            break;
+        aa = mul_mod(aa, aa, m);
+    }
+    return r;
+}
+
+/* return true if n is prime */
+static int is_prime(int n)
+{
+    int r, i;
+    if ((n % 2) == 0)
+        return 0;
+
+    r = (int) (sqrt(n));
+    for (i = 3; i <= r; i += 2)
+        if ((n % i) == 0)
+            return 0;
+    return 1;
+}
+
+/* return the prime number immediatly after n */
+static int next_prime(int n)
+{
+    do {
+        n++;
+    } while (!is_prime(n));
+    return n;
+}
+
+int pi_nth_digit(int n)
+{
+    int av, a, vmax, N, num, den, k, kq, kq2, t, v, s, i;
+    double sum;
+
+    N = (int) ((n + 20) * log(10) / log(2));
+
+    sum = 0;
+
+    for (a = 3; a <= (2 * N); a = next_prime(a)) {
+
+        vmax = (int) (log(2 * N) / log(a));
+        av = 1;
+        for (i = 0; i < vmax; i++)
+            av = av * a;
+
+        s = 0;
+        num = 1;
+        den = 1;
+        v = 0;
+        kq = 1;
+        kq2 = 1;
+
+        for (k = 1; k <= N; k++) {
+
+            t = k;
+            if (kq >= a) {
+                do {
+                    t = t / a;
+                    v--;
+                } while ((t % a) == 0);
+                kq = 0;
+            }
+            kq++;
+            num = mul_mod(num, t, av);
+
+            t = (2 * k - 1);
+            if (kq2 >= a) {
+                if (kq2 == a) {
+                    do {
+                        t = t / a;
+                        v++;
+                    } while ((t % a) == 0);
+                }
+                kq2 -= a;
+            }
+            den = mul_mod(den, t, av);
+            kq2 += 2;
+
+            if (v > 0) {
+                t = inv_mod(den, av);
+                t = mul_mod(t, num, av);
+                t = mul_mod(t, k, av);
+                for (i = v; i < vmax; i++)
+                    t = mul_mod(t, a, av);
+                s += t;
+                if (s >= av)
+                    s -= av;
+            }
+
+        }
+
+        t = pow_mod(10, n - 1, av);
+        s = mul_mod(s, t, av);
+        sum = fmod(sum + (double) s / (double) av, 1.0);
+    }
+
+    return (int) (sum * 1e9);
+}

performance/numerical/src/nativeInterop/cinterop/pi.h

@@ -0,0 +1,16 @@
+/*
+ * Computation of the n'th decimal digit of \pi with very little memory.
+ * Written by Fabrice Bellard on January 8, 1997.
+ *
+ * We use a slightly modified version of the method described by Simon
+ * Plouffe in "On the Computation of the n'th decimal digit of various
+ * transcendental numbers" (November 1996). We have modified the algorithm
+ * to get a running time of O(n^2) instead of O(n^3log(n)^3).
+ */
+
+#ifndef _BELLARD_PI_H
+#define _BELLARD_PI_H
+
+int pi_nth_digit(int n);
+
+#endif /*_BELLARD_PI_H*/

settings.gradle

@@ -35,6 +35,7 @@ include ':performance'
 include ':performance:ring'
 include ':performance:cinterop'
 include ':performance:helloworld'
+include ':performance:numerical'
 include ':performance:videoplayer'
 include ':performance:framework'
 if (System.getProperty("os.name") == "Mac OS X") {