multi_precision.h

Go to the documentation of this file.

/*
 *  Copyright (c) 2000-2022 Inria
 *  All rights reserved.
 *
 *  Redistribution and use in source and binary forms, with or without
 *  modification, are permitted provided that the following conditions are met:
 *
 *  * Redistributions of source code must retain the above copyright notice,
 *  this list of conditions and the following disclaimer.
 *  * Redistributions in binary form must reproduce the above copyright notice,
 *  this list of conditions and the following disclaimer in the documentation
 *  and/or other materials provided with the distribution.
 *  * Neither the name of the ALICE Project-Team nor the names of its
 *  contributors may be used to endorse or promote products derived from this
 *  software without specific prior written permission.
 *
 *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 *  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 *  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 *  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 *  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 *  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 *  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 *  POSSIBILITY OF SUCH DAMAGE.
 *
 *  Contact: Bruno Levy
 *
 *     https://www.inria.fr/fr/bruno-levy
 *
 *     Inria,
 *     Domaine de Voluceau,
 *     78150 Le Chesnay - Rocquencourt
 *     FRANCE
 *
 */
 
#ifndef GEOGRAM_NUMERICS_MULTI_PRECISION
#define GEOGRAM_NUMERICS_MULTI_PRECISION
 
#include <geogram/basic/common.h>
#include <geogram/basic/numeric.h>
#include <geogram/basic/memory.h>
#include <geogram/basic/assert.h>
#include <iostream>
#include <sstream>
#include <new>
#include <math.h>
 
namespace GEO {
 
    extern double expansion_splitter_;
    extern double expansion_epsilon_;
 
    inline void two_sum(double a, double b, double& x, double& y) {
        x = a + b;
        double bvirt = x - a;
        double avirt = x - bvirt;
        double bround = b - bvirt;
        double around = a - avirt;
        y = around + bround;
    }
 
    inline void two_diff(double a, double b, double& x, double& y) {
        x = a - b;
        double bvirt = a - x;
        double avirt = x + bvirt;
        double bround = bvirt - b;
        double around = a - avirt;
        y = around + bround;
    }
 
    inline void split(double a, double& ahi, double& alo) {
        double c = expansion_splitter_ * a;
        double abig = c - a;
        ahi = c - abig;
        alo = a - ahi;
    }
 
    inline void two_product(double a, double b, double& x, double& y) {
#ifdef FP_FAST_FMA
        // If the target processor supports the FMA (Fused Multiply Add)
        // instruction, then the product of two doubles into a length-2
        // expansion can be implemented as follows. Thanks to Marc Glisse
        // for the information.
        // Note: under gcc, automatic generations of fma() for a*b+c needs
        // to be deactivated, using -ffp-contract=off, else it may break
        // other functions such as fast_expansion_sum_zeroelim().
        x = a*b;
        y = fma(a,b,-x);
#else
        x = a * b;
        double ahi, alo;
        split(a, ahi, alo);
        double bhi, blo;
        split(b, bhi, blo);
        double err1 = x - (ahi * bhi);
        double err2 = err1 - (alo * bhi);
        double err3 = err2 - (ahi * blo);
        y = (alo * blo) - err3;
#endif
    }
 
    inline void square(double a, double& x, double& y) {
#ifdef FP_FAST_FMA
        // If the target processor supports the FMA (Fused Multiply Add)
        // instruction, then the product of two doubles into a length-2
        // expansion can be implemented as follows. Thanks to Marc Glisse
        // for the information.
        // Note: under gcc, automatic generations of fma() for a*b+c needs
        // to be deactivated, using -ffp-contract=off, else it may break
        // other functions such as fast_expansion_sum_zeroelim().
        x = a*a;
        y = fma(a,a,-x);
#else
        x = a * a;
        double ahi, alo;
        split(a, ahi, alo);
        double err1 = x - (ahi * ahi);
        double err3 = err1 - ((ahi + ahi) * alo);
        y = (alo * alo) - err3;
#endif
    }
 
    /************************************************************************/
 
    class GEOGRAM_API expansion {
    public:
    index_t length() const {
        return length_;
    }
 
    index_t capacity() const {
        return capacity_;
    }
 
    void set_length(index_t new_length) {
        geo_debug_assert(new_length <= capacity());
        length_ = new_length;
    }
 
    const double& operator[] (index_t i) const {
        // Note: we allocate capacity+1 storage
        // systematically, since basic functions
        // may access one additional value (without
        // using it)
        geo_debug_assert(i <= capacity_);
        return x_[i];
    }
 
    double& operator[] (index_t i) {
        // Note: we allocate capacity+1 storage
        // systematically, since basic functions
        // may access one additional value (without
        // using it)
        geo_debug_assert(i <= capacity_);
        return x_[i];
    }
 
    double* data() {
        return x_;
    }
 
    const double* data() const {
        return x_;
    }
 
    static size_t bytes(index_t capa) {
        // --> 2*sizeof(double) because x_ is declared of size [2]
        // to avoid compiler's warning.
        // --> capa+1 to have an additional 'sentry' at the end
        // because fast_expansion_sum_zeroelim() may access
        // an entry past the end (without using it).
        return
            sizeof(expansion) - 2 * sizeof(double) +
            (capa + 1) * sizeof(double);
    }
 
    static size_t bytes_on_stack(index_t capa) {
#ifndef GEO_HAS_BIG_STACK
        // Note: standard predicates need at least 512, hence the min.
        // index_t(MAX_CAPACITY_ON_STACK) is necessary, else with
        // MAX_CAPACITY_ON_STACK alone the compiler tries to generate a
        // reference to NOT_IN_LIST resulting in a link error.
        // (weird, even with constexpr, I do not understand...)
        // Probably when the function excepts a *reference*
        geo_debug_assert(
            capa <= std::max(index_t(MAX_CAPACITY_ON_STACK),index_t(512))
        );
#endif
        return bytes(capa);
    }
 
    expansion(index_t capa) :
    length_(0),
    capacity_(capa) {
    }
 
#ifdef CPPCHECK
    // cppcheck does not understand that the result
    // of alloca() is passed to the placement syntax
    // of operator new.
    expansion& new_expansion_on_stack(index_t capa);
#else
#define new_expansion_on_stack(capa)                                    \
    (new (alloca(expansion::bytes_on_stack(capa)))expansion(capa))
#endif
 
    static expansion* new_expansion_on_heap(index_t capa);
 
    static void delete_expansion_on_heap(expansion* e);
 
    // ========================== Initialization from doubles
 
    expansion& assign(double a) {
        set_length(1);
        x_[0] = a;
        return *this;
    }
 
    expansion& assign(const expansion& rhs) {
        geo_debug_assert(capacity() >= rhs.length());
        set_length(rhs.length());
        for(index_t i=0; i<rhs.length(); ++i) {
            x_[i] = rhs.x_[i];
        }
        return *this;
    }
 
    expansion& assign_abs(const expansion& rhs) {
        assign(rhs);
        if(sign() == NEGATIVE) {
            negate();
        }
        return *this;
    }
 
    static index_t sum_capacity(double a, double b) {
        geo_argused(a);
        geo_argused(b);
        return 2;
    }
 
    expansion& assign_sum(double a, double b) {
        set_length(2);
        two_sum(a, b, x_[1], x_[0]);
        return *this;
    }
 
    static index_t diff_capacity(double a, double b) {
        geo_argused(a);
        geo_argused(b);
        return 2;
    }
 
    expansion& assign_diff(double a, double b) {
        set_length(2);
        two_diff(a, b, x_[1], x_[0]);
        return *this;
    }
 
    static index_t product_capacity(double a, double b) {
        geo_argused(a);
        geo_argused(b);
        return 2;
    }
 
    expansion& assign_product(double a, double b) {
        set_length(2);
        two_product(a, b, x_[1], x_[0]);
        return *this;
    }
 
    static index_t square_capacity(double a) {
        geo_argused(a);
        return 2;
    }
 
    expansion& assign_square(double a) {
        set_length(2);
        square(a, x_[1], x_[0]);
        return *this;
    }
 
    // ====== Initialization from expansion and double
 
    static index_t sum_capacity(const expansion& a, double b) {
        geo_argused(b);
        return a.length() + 1;
    }
 
    expansion& assign_sum(const expansion& a, double b);
 
    static index_t diff_capacity(const expansion& a, double b) {
        geo_argused(b);
        return a.length() + 1;
    }
 
    expansion& assign_diff(const expansion& a, double b);
 
    static index_t product_capacity(const expansion& a, double b) {
        geo_argused(b);
        // TODO: implement special case where the double argument
        // is a power of two.
        return a.length() * 2;
    }
 
    expansion& assign_product(const expansion& a, double b);
 
    // ========================== Initialization from expansions
 
    static index_t sum_capacity(const expansion& a, const expansion& b) {
        return a.length() + b.length();
    }
 
    expansion& assign_sum(const expansion& a, const expansion& b);
 
    static index_t sum_capacity(
        const expansion& a, const expansion& b, const expansion& c
    ) {
        return a.length() + b.length() + c.length();
    }
 
    expansion& assign_sum(
        const expansion& a, const expansion& b, const expansion& c
    );
 
    static index_t sum_capacity(
        const expansion& a, const expansion& b,
        const expansion& c, const expansion& d
    ) {
        return a.length() + b.length() + c.length() + d.length();
    }
 
    expansion& assign_sum(
        const expansion& a, const expansion& b,
        const expansion& c, const expansion& d
    );
 
    static index_t diff_capacity(const expansion& a, const expansion& b) {
        return a.length() + b.length();
    }
 
    expansion& assign_diff(const expansion& a, const expansion& b);
 
    static index_t product_capacity(
        const expansion& a, const expansion& b
    ) {
        return a.length() * b.length() * 2;
    }
 
    expansion& assign_product(const expansion& a, const expansion& b);
 
    static index_t product_capacity(
        const expansion& a, const expansion& b, const expansion& c
    ) {
        return a.length() * b.length() * c.length() * 4;
    }
 
    expansion& assign_product(
        const expansion& a, const expansion& b, const expansion& c
    );
 
    static index_t square_capacity(const expansion& a) {
        if(a.length() == 2) {
            return 6;
        }                                  // see two_square()
        return a.length() * a.length() * 2;
    }
 
    expansion& assign_square(const expansion& a);
 
    // ====== Determinants =============================
 
    static index_t det2x2_capacity(
        const expansion& a11, const expansion& a12,
        const expansion& a21, const expansion& a22
    ) {
        return
            product_capacity(a11, a22) +
            product_capacity(a21, a12);
    }
 
    expansion& assign_det2x2(
        const expansion& a11, const expansion& a12,
        const expansion& a21, const expansion& a22
    );
 
    static index_t det3x3_capacity(
        const expansion& a11, const expansion& a12, const expansion& a13,
        const expansion& a21, const expansion& a22, const expansion& a23,
        const expansion& a31, const expansion& a32, const expansion& a33
    ) {
        // Development w.r.t. first row
        index_t c11_capa = det2x2_capacity(a22, a23, a32, a33);
        index_t c12_capa = det2x2_capacity(a21, a23, a31, a33);
        index_t c13_capa = det2x2_capacity(a21, a22, a31, a32);
        return 2 * (
            a11.length() * c11_capa +
            a12.length() * c12_capa +
            a13.length() * c13_capa
        );
    }
 
    expansion& assign_det3x3(
        const expansion& a11, const expansion& a12, const expansion& a13,
        const expansion& a21, const expansion& a22, const expansion& a23,
        const expansion& a31, const expansion& a32, const expansion& a33
    );
 
    static index_t det_111_2x3_capacity(
        const expansion& a21, const expansion& a22, const expansion& a23,
        const expansion& a31, const expansion& a32, const expansion& a33
    ) {
        return
            det2x2_capacity(a22, a23, a32, a33) +
            det2x2_capacity(a23, a21, a33, a31) +
            det2x2_capacity(a21, a22, a31, a32);
    }
 
    expansion& assign_det_111_2x3(
        const expansion& a21, const expansion& a22, const expansion& a23,
        const expansion& a31, const expansion& a32, const expansion& a33
    );
 
    // ======= Geometry-specific initializations =======
 
    static index_t sq_dist_capacity(coord_index_t dim) {
        return index_t(dim) * 6;
    }
 
    expansion& assign_sq_dist(
        const double* p1, const double* p2, coord_index_t dim
    );
 
    static index_t dot_at_capacity(coord_index_t dim) {
        return index_t(dim) * 8;
    }
 
    expansion& assign_dot_at(
        const double* p1, const double* p2, const double* p0,
        coord_index_t dim
    );
 
 
    static index_t length2_capacity(
        const expansion& x, const expansion& y, const expansion& z
    ) {
        return square_capacity(x) + square_capacity(y) + square_capacity(z);
    }
 
    expansion& assign_length2(
        const expansion& x, const expansion& y, const expansion& z
    );
 
    // =============== some general purpose functions =========
 
    static void initialize();
 
    expansion& negate() {
        for(index_t i = 0; i < length_; ++i) {
            x_[i] = -x_[i];
        }
        return *this;
    }
 
    expansion& scale_fast(double s) {
        // TODO: debug assert is_power_of_two(s)
        for(index_t i = 0; i < length_; ++i) {
            x_[i] *= s;
        }
        return *this;
    }
 
    double estimate() const {
        double result = 0.0;
        for(index_t i = 0; i < length(); ++i) {
            result += x_[i];
        }
        return result;
    }
 
    Sign sign() const {
        if(length() == 0) {
            return ZERO;
        }
        return geo_sgn(x_[length() - 1]);
    }
 
    bool is_same_as(const expansion& rhs) const;
 
    bool is_same_as(double rhs) const;
 
 
    Sign compare(const expansion& rhs) const;
 
    Sign compare(double rhs) const;
 
    bool equals(const expansion& rhs) const {
        return (compare(rhs) == ZERO);
    }
 
    bool equals(double rhs) const {
        return (compare(rhs) == ZERO);
    }
 
    std::ostream& show(std::ostream& out) const {
        out << "expansion[" << length() << "] = [";
        for(index_t i=0; i<length(); ++i) {
            out << (*this)[i] << " ";
        }
        out << "]";
        return out;
    }
 
    std::string to_string() const {
        std::ostringstream out;
        show(out);
        return out.str();
    }
 
    void optimize();
 
    static void show_all_stats();
 
    protected:
    static index_t sub_product_capacity(
        index_t a_length, index_t b_length
    ) {
        return a_length * b_length * 2;
    }
 
    expansion& assign_sub_product(
        const double* a, index_t a_length, const expansion& b
    );
 
    expansion(const expansion& rhs) = delete;
 
    expansion& operator= (const expansion& rhs) = delete;
 
    private:
 
#ifdef GEO_OS_APPLE
    static constexpr index_t MAX_CAPACITY_ON_STACK = 256;
#else
    static constexpr index_t MAX_CAPACITY_ON_STACK = 1024;
#endif
    index_t length_;
    index_t capacity_;
    double x_[2];  // x_ is in fact of size [capacity_]
 
    friend class expansion_nt;
    };
 
    // =============== arithmetic operations ===========================
 
#define expansion_create(a)                     \
    new_expansion_on_stack(1)->assign(a)
 
 
#define expansion_abs(e)                                \
    new_expansion_on_stack(e.length())->assign_abs(e)
 
#define expansion_sum(a, b)                     \
    new_expansion_on_stack(                     \
        expansion::sum_capacity(a, b)           \
    )->assign_sum(a, b)
 
#define expansion_sum3(a, b, c)                 \
    new_expansion_on_stack(                     \
        expansion::sum_capacity(a, b, c)        \
    )->assign_sum(a, b, c)
 
#define expansion_sum4(a, b, c, d)              \
    new_expansion_on_stack(                     \
        expansion::sum_capacity(a, b, c, d)     \
    )->assign_sum(a, b, c, d)
 
#define expansion_diff(a, b)                    \
    new_expansion_on_stack(                     \
        expansion::diff_capacity(a, b)          \
    )->assign_diff(a, b)
 
#define expansion_product(a, b)                 \
    new_expansion_on_stack(                     \
        expansion::product_capacity(a, b)       \
    )->assign_product(a, b)
 
#define expansion_product3(a, b, c)             \
    new_expansion_on_stack(                     \
        expansion::product_capacity(a, b, c)    \
    )->assign_product(a, b, c)
 
#define expansion_square(a)                     \
    new_expansion_on_stack(                     \
        expansion::square_capacity(a)           \
    )->assign_square(a)
 
    // =============== determinants =====================================
 
#define expansion_det2x2(a11, a12, a21, a22)            \
    new_expansion_on_stack(                             \
        expansion::det2x2_capacity(a11, a12, a21, a22)  \
    )->assign_det2x2(a11, a12, a21, a22)
 
#define expansion_det3x3(a11, a12, a13, a21, a22, a23, a31, a32, a33)   \
    new_expansion_on_stack(                                             \
        expansion::det3x3_capacity(a11,a12,a13,a21,a22,a23,a31,a32,a33) \
    )->assign_det3x3(a11, a12, a13, a21, a22, a23, a31, a32, a33)
 
#define expansion_det_111_2x3(a21, a22, a23, a31, a32, a33)             \
    new_expansion_on_stack(                                             \
        expansion::det_111_2x3_capacity(a21, a22, a23, a31, a32, a33)   \
    )->assign_det_111_2x3(a21, a22, a23, a31, a32, a33)
 
    // =============== geometric functions ==============================
 
#define expansion_sq_dist(a, b, dim)            \
    new_expansion_on_stack(                     \
        expansion::sq_dist_capacity(dim)        \
    )->assign_sq_dist(a, b, dim)
 
#define expansion_dot_at(a, b, c, dim)          \
    new_expansion_on_stack(                     \
        expansion::dot_at_capacity(dim)         \
    )->assign_dot_at(a, b, c, dim)
 
 
#define expansion_length2(x,y,z)                \
    new_expansion_on_stack(                     \
        expansion::length2_capacity(x,y,z)      \
    )->assign_length2(x,y,z)
 
    /************************************************************************/
 
    Sign GEOGRAM_API sign_of_expansion_determinant(
        const expansion& a00,const expansion& a01,
        const expansion& a10,const expansion& a11
    );
 
    Sign GEOGRAM_API sign_of_expansion_determinant(
        const expansion& a00,const expansion& a01,const expansion& a02,
        const expansion& a10,const expansion& a11,const expansion& a12,
        const expansion& a20,const expansion& a21,const expansion& a22
    );
 
    Sign GEOGRAM_API sign_of_expansion_determinant(
        const expansion& a00,const expansion& a01,
        const expansion& a02,const expansion& a03,
        const expansion& a10,const expansion& a11,
        const expansion& a12,const expansion& a13,
        const expansion& a20,const expansion& a21,
        const expansion& a22,const expansion& a23,
        const expansion& a30,const expansion& a31,
        const expansion& a32,const expansion& a33
    );
 
    /************************************************************************/
 
    void GEOGRAM_API grow_expansion_zeroelim(
        const expansion& e, double b, expansion& h
    );
 
    void GEOGRAM_API scale_expansion_zeroelim(
        const expansion& e, double b, expansion& h
    );
 
    void GEOGRAM_API fast_expansion_sum_zeroelim(
        const expansion& e, const expansion& f, expansion& h
    );
 
 
    void GEOGRAM_API fast_expansion_diff_zeroelim(
        const expansion& e, const expansion& f, expansion& h
    );
 
    /************************************************************************/
}
 
#endif