RoundingVariants performance comparison

[whirlicote]

Dear CGAL contributors,

I would like to prepare a performance comparison of native RoundingVariants versus a software emulation in WebAssembly.

What would make a good CGAL demo example?
How would this be approached best?

Paul

[sloriot]

Hi,

If I understand correctly, you'd like to add proper rounding in WebAssembly and you are looking for a test case in CGAL that is failing with the default to-nearest rounding?

[whirlicote]

If I understand correctly, you'd like to add proper rounding in WebAssembly and you are looking for a test case in CGAL that is failing with the default to-nearest rounding?

That would certainly be quite insightful to investigate these cases. But for starters, I think a demo would suffice that executes lots of RoundingVariants and gets a speedup when native RoundingVariants are used.

I have a software reference implementation of the RoundingVariants from the WebAssembly RoundingVariants proposal.

I am not familiar with the CGAL part, but I was thinking of trying one of these two approaches:

1. High-level approach

   • Adding my C or C++ version of the reference implementation as a backend for the used intervals
   • Marking the new function symbols as noinline and verifying that noinline is honored all the way down to the .wasm WebAssembly module by the toolchain
   • A demo assembled with an artifact like this should run on any current browser.
   • In a plain Chromium browser, I would take the baseline measurement.
   • Then I do some low-level overloading in the module by replacing the function calls to the RoundingVariants by native instructions.
   • This patched version then would only run on my browser fork that implements native RoundingVariants as a reference
   • From this second run I would take a second measurement, which I hope to be faster.
   • Then I would take the ratio between these two numbers so that I can report the performance improvement numbers to the WebAssembly Community Group.

2. Low-level approach

   • CGAL would not implement the RoundingVariants but just reference them
   • Only at the very end in the browser in JavaScript do these get resolved.
   • For the first measurement, the final .wasm module is linked in JavaScript with another WebAssembly module I have prepared that implements the RoundingVariants in plain software.
   • For the second measurement, I replace my module with a trivial module that just simply forwards the parameters to the native instruction.

[mglisse]

In CGAL/FPU.h, we currently have #define CGAL_IA_ADD(a,b) IA_up((a)+CGAL_IA_STOP_CPROP(b)), and similar for other operations. It would make sense to have a macro that enables an alternative definition of those using rounded arithmetic (with subcases for WASM, RISCV, AVX512/AVX10, nvidia, as long as we don't have a standard way to refer to such instructions), with a few extra tweaks elsewhere (similar to CGAL_ALWAYS_ROUND_TO_NEAREST). It may be easier to do this first for AVX512, temporarily disabling the SIMD path, so we can more easily run the whole testsuite.

There are suggestions of benchmarks that are sensitive to interval arithmetic around #5510 (comment) .

[whirlicote]

In CGAL/FPU.h

$ find -iname "FPU.h"
./Number_types/doc/Number_types/CGAL/FPU.h
./Number_types/include/CGAL/FPU.h

The second file contains the macro CGAL_IA_ADD. The similar macros seem to be:

CGAL_IA_ADD
CGAL_IA_SUB
CGAL_IA_MUL
CGAL_IA_DIV
CGAL_IA_SQUARE

Assuming f64 precicion the reference implementation for these instructions in C++ would be:

// MIT License
//
// Copyright (c) 2022 Paul Dennis
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.


// Requires IEEE 754 round-to-nearest-even semantics.
// Do not compile with -ffast-math or equivalent.
#include <cstdint>
#include <cmath>
#include <bit>

double f64_add_ceil(double a, double b) {
    double s = a + b;
    if (s == -INFINITY) {
        if (a != -INFINITY && b != -INFINITY)
            return -0x1.fffffffffffffp+1023;
        return s;
    }
    double v = s - a;
    double t = (a - (s - v)) + (b - v);
    if (t > 0.0) {
        int64_t bits = std::bit_cast<int64_t>(s);
        return std::bit_cast<double>(bits >= 0 ? bits + 1 : bits - 1);
    }
    return s;
}

double f64_sub_ceil(double a, double b) {
    double s = a - b;
    if (s == -INFINITY) {
        if (a != -INFINITY && b != INFINITY)
            return -0x1.fffffffffffffp+1023;
        return s;
    }
    double v = s - a;
    double t = (a - (s - v)) - (b + v);
    if (t > 0.0) {
        int64_t bits = std::bit_cast<int64_t>(s);
        return std::bit_cast<double>(bits >= 0 ? bits + 1 : bits - 1);
    }
    return s;
}

double f64_mul_ceil(double a, double b) {
    double s = a * b;
    if (s == -INFINITY) {
        if (std::fabs(a) != INFINITY && std::fabs(b) != INFINITY)
            return -0x1.fffffffffffffp+1023;
        return s;
    }
    if (std::bit_cast<uint64_t>(s) == 0 && a != 0.0 && b != 0.0)
        return 0x0.0000000000001p-1022;
    double a1 = a, b1 = b, p, t, g;
    if (std::fabs(s) < 0x1p-1022) {
        if (std::fabs(a) <= std::fabs(b))
            a1 = a * 0x1p+600;
        else
            b1 = b * 0x1p+600;
        p = a1 * b1;
        double gs = s * 0x1p+600;
        if (p > gs) {
            int64_t bits = std::bit_cast<int64_t>(s);
            return std::bit_cast<double>(bits >= 0 ? bits + 1 : bits - 1);
        }
        if (p < gs)
            return s;
    }
    if (std::fabs(a1) >= 0x1p+996) {
        a1 *= 0x1p-28; b1 *= 0x1p+28;
    } else if (std::fabs(b1) >= 0x1p+996) {
        a1 *= 0x1p+28; b1 *= 0x1p-28;
    }
    g = a1 * 134217729.0;
    double ah = g - (g - a1), al = a1 - ah;
    g = b1 * 134217729.0;
    double bh = g - (g - b1), bl = b1 - bh;
    p = a1 * b1;
    t = ((ah * bh - p) + ah * bl + al * bh) + al * bl;
    if (t > 0.0) {
        int64_t bits = std::bit_cast<int64_t>(s);
        return std::bit_cast<double>(bits >= 0 ? bits + 1 : bits - 1);
    }
    return s;
}

double f64_div_ceil(double a, double b) {
    double s = a / b;
    if (s == -INFINITY) {
        if (std::fabs(a) != INFINITY && std::fabs(b) != INFINITY && b != 0.0)
            return -0x1.fffffffffffffp+1023;
        return s;
    }
    if (std::bit_cast<uint64_t>(s) == 0 && a != 0.0)
        return 0x0.0000000000001p-1022;
    if (s == 0.0 || !(std::fabs(s) < INFINITY))
        return s;
    double pc = s * b;
    double s1 = s, b1 = b, a1 = a;
    if (std::fabs(pc) < 0x1p-1022) {
        s1 = s * 0x1p+600;
        b1 = b * 0x1p+600;
        a1 = a * 0x1p+600 * 0x1p+600;
    } else if (std::fabs(b) < 0x1p-1022) {
        s1 = s * 0x1p-600;
        b1 = b * 0x1p+600;
    } else if (std::fabs(s) >= 0x1p+996) {
        s1 = s * 0x1p-28;
        b1 = b * 0x1p+28;
    } else if (std::fabs(b) >= 0x1p+996) {
        s1 = s * 0x1p+28;
        b1 = b * 0x1p-28;
    }
    double g = s1 * 134217729.0;
    double sh = g - (g - s1), sl = s1 - sh;
    g = b1 * 134217729.0;
    double bh = g - (g - b1), bl = b1 - bh;
    double p = s1 * b1;
    double t = ((sh * bh - p) + sh * bl + sl * bh) + sl * bl;
    double r1 = a1 - p;
    if ((b > 0.0 && r1 > t) || (b < 0.0 && r1 < t)) {
        int64_t bits = std::bit_cast<int64_t>(s);
        return std::bit_cast<double>(bits >= 0 ? bits + 1 : bits - 1);
    }
    return s;
}

double f64_sqrt_ceil(double radicand) {
    double root_nearest = std::sqrt(radicand);
    double cr = root_nearest, cd = radicand;
    if (root_nearest < 0x1p-470) {
        cr = root_nearest * 0x1p511;
        cd = radicand * 0x1p1022;
    }
    if (root_nearest > 0x1p510) {
        cr = root_nearest * 0x1p-511;
        cd = radicand * 0x1p-1022;
    }
    double hi = cr * cr;
    double c = 134217729.0 * cr;
    double a_hi = c - (c - cr);
    double a_lo = cr - a_hi;
    double cross = a_hi * a_lo;
    double e1 = hi - a_hi * a_hi;
    double e2 = e1 - cross;
    double e3 = e2 - cross;
    double lo = a_lo * a_lo - e3;
    if (hi < cd || (hi == cd && lo < 0.0))
        return std::bit_cast<double>(std::bit_cast<int64_t>(root_nearest) + 1);
    return root_nearest;
}

What is the used precicion of CGAL_IA_ADD and the similar macros?

[mglisse]

In CGAL/FPU.h

./Number_types/doc/Number_types/CGAL/FPU.h

That's the doc.

./Number_types/include/CGAL/FPU.h

The second file contains the macro CGAL_IA_ADD. The similar macros seem to be:

CGAL_IA_ADD
CGAL_IA_SUB
CGAL_IA_MUL
CGAL_IA_DIV
CGAL_IA_SQUARE

Yes, plus IA_sqrt_*.

Assuming f64 precicion the reference implementation for these instructions in C++ would be:
[...]

Thanks. There may be some cases that cannot happen, for instance the upper bound of an interval is probably not allowed to be -infinity, so I think we don't have to handle the case where a or b is -inf in up(a+b), but that's a microoptimization, not worth worrying about.

What is the used precicion of CGAL_IA_ADD and the similar macros?

Not sure what you mean, but they work with double, which we expect to be 64 bits. We normally want them as precise as possible for performance reasons, but the validity of the program only requires that they be at least as large as the true result, for instance some double rounding is ok as long as they are both upwards. The current workaround in the absence of upward rounding is to perform the operation with some arbitrary rounding then call nextafter, which is neither fast nor tight, but at least correct.

[whirlicote]

IA_sqrt_*.

Here is the adaption for sqrt_floor:

double f64_sqrt_floor(double radicand) {
    // [...]
    if (hi > cd || (hi == cd && lo > 0.0))
        return std::bit_cast<double>(std::bit_cast<int64_t>(root_nearest) - 1);
    return root_nearest;
}

It would make sense to have a macro that enables an alternative definition of those using rounded arithmetic (with subcases for WASM, RISCV, AVX512/AVX10, nvidia, as long as we don't have a standard way to refer to such instructions),

I believe this to be a good idea.