alordash / FastExponentiation

Algorithms and functions for fast exponentiation with small error
https://alordash.github.io/FastExponentiation/publish/wwwroot/
9 stars 1 forks source link
approximation approximation-algorithms exponent exponentiation

Fast exponentiation

Description

This repository contains realizations of various algorithms for fast exponentiation with small error in C++, C# or Java.
More info: "Ускоряем pow"

List of algorithms

  1. "Old" approximation
  2. Binary power
  3. Dividing fast power
  4. Fractional fast power
  5. "Another" approximation

Repository structure

There are 9 projects in visual studio solution and 5 java projects used to test custom algorithms, including web page for precision tests.

You can check results of performance measures in this excel table.
Tests ran on: i5-10300H, 19.8 DDR4 GB of usable RAM, 64bit, single threaded
C++: MSVC + /O2 + /Oi + /Ot
C#: "Optimize code" option

Algorithms realizations

"Old" approximation

In C++

double OldApproximatePower(double b, double e) {
    union {
        double d;
        long long i;
    } u = { b };
    u.i = (long long)(4606853616395542500L + e * (u.i - 4606853616395542500L));
    return u.d;
}
In C# ```C# double OldApproximatePower(double b, double e) { long i = BitConverter.DoubleToInt64Bits(b); i = (long)(FastMath.doubleApproximator + e * (i - FastMath.doubleApproximator)); return BitConverter.Int64BitsToDouble(i); } ```
In Java ```java double OldApproximatePower(double b, double e) { long i = Double.doubleToLongBits(b); i = (long) (FastMath.doubleApproximator + e * (i - FastMath.doubleApproximator)); return Double.longBitsToDouble(i); } ```

Binary power

In C++

double BinaryPower(double b, unsigned long long e) {
    double v = 1.0;
    while(e != 0) {
        if((e & 1) != 0) {
            v *= b;
        }
        b *= b;
        e >>= 1;
    }
    return v;
}
In C# ```c# double BinaryPower(double b, UInt64 e) { double v = 1d; while(e != 0) { if((e & 1) != 0) { v *= b; } b *= b; e >>= 1; } return v; } ```
In Java ```java double BinaryPower(double b, long e) { double v = 1d; while (e > 0) { if ((e & 1) != 0) { v *= b; } b *= b; e >>= 1; } return v; } ```

Dividing fast power

In C++

double FastPowerDividing(double b, double e) {
    // To avoid undefined behaviour near key points,
    // we can hardcode results for them, but this
    // will make function slightly slower.
    if(b == 1.0 || e == 0.0) {
        return 1.0;
    }

    double eAbs = fabs(e);
    double el = ceil(eAbs);
    double basePart = OldApproximatePower(b, eAbs / el);

    double result = BinaryPower(basePart, (unsigned long long)el);
    // Because OldApproximatePower gives inaccurate results
    // with negative exponent, we can increase precision
    // by calculating exponent of a number in positive power
    // and then dividing 1 by result of calculation
    if(e < 0.0) {
        return 1.0 / result;
    }
    return result;
}
In C# ```c# double FastPowerDividing(double b, double e) { if(b == 1d || e == 0d) { return 1d; } var eAbs = Math.Abs(e); var el = Math.Ceiling(eAbs); var basePart = OldApproximatePower(b, eAbs / el); var result = BinaryPower(basePart, (UInt64)el); if(e < 0d) { return 1d / result; } return result; } ```
In Java ```java double FastPowerDividing(double b, double e) { if (b == 1d || e == 0d) { return 1d; } var eAbs = Math.abs(e); var el = Math.ceil(eAbs); var basePart = OldApproximatePower(b, eAbs / el); var result = BinaryPower(basePart, (long) el); if (e < 0d) { return 1d / result; } return result; } ```

Fractional fast power

In C++

double FastPowerFractional(double b, double e) {
    if(b == 1.0 || e == 0.0) {
        return 1.0;
    }

    double absExp = fabs(e);
    unsigned long long eIntPart = (unsigned long long)absExp;
    double eFractPart = absExp - eIntPart;
    double result = OldApproximatePower(b, eFractPart) * BinaryPower(b, eIntPart);

    if(e < 0.0) {
        return 1.0 / result;
    }
    return result;
}
In C# ```c# double FastPowerFractional(double b, double e) { if(b == 1d || e == 0d) { return 1d; } double absExp = Math.Abs(e); UInt64 eIntPart = (UInt64)absExp; double eFractPart = absExp - eIntPart; double result = OldApproximatePower(b, eFractPart) * BinaryPower(b, eIntPart); if(e < 0d) { return 1d / result; } return result; } ```
In Java ```java double FastPowerFractional(double b, double e) { if (b == 1d || e == 0d) { return 1d; } double absExp = Math.abs(e); long eIntPart = (long)absExp; double eFractPart = absExp - eIntPart; double result = OldApproximatePower(b, eFractPart) * BinaryPower(b, eIntPart); if(e < 0d) { return 1d / result; } return result; } ```

Another approximation

In C++

double AnotherApproximatePower(double a, double b) {
    union {
        double d;
        int x[2];
    } u = { a };
    u.x[1] = (int)(b * (u.x[1] - 1072632447) + 1072632447);
    u.x[0] = 0;
    return u.d;
}
In C# ```c# double AnotherApproxPower(double a, double b) { int tmp = (int)(BitConverter.DoubleToInt64Bits(a) >> 32); int tmp2 = (int)(b * (tmp - 1072632447) + 1072632447); return BitConverter.Int64BitsToDouble(((long)tmp2) << 32); } ```
In Java ```java double AnotherApproxPower(double a, double b) { int tmp = (int)(Double.doubleToLongBits(a) >> 32); int tmp2 = (int)(b * (tmp - 1072632447) + 1072632447); return Double.longBitsToDouble(((long)tmp2) << 32); } ```