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