expansion_nt.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_EXPANSION_NT
#define GEOGRAM_NUMERICS_EXPANSION_NT
 
#include <geogram/basic/common.h>
#include <geogram/numerics/multi_precision.h>
#include <geogram/basic/matrix.h>
#include <geogram/basic/rationalg.h>
 
namespace GEO {
 
    class expansion_nt;
 
    class GEOGRAM_API expansion_nt {
    public:
 
    enum Operation {
        SUM, DIFF, PRODUCT
    };
 
    expansion_nt() : rep_(nullptr) {
    }
 
    explicit expansion_nt(double x) {
        rep_ = expansion::new_expansion_on_heap(1);
        rep()[0] = x;
        rep().set_length(1);
    }
 
    explicit expansion_nt(const expansion& rhs) {
        rep_ = expansion::new_expansion_on_heap(rhs.length());
        rep().assign(rhs);
    }
 
    explicit expansion_nt(
        Operation op, const expansion& x, const expansion& y
    ) {
        switch(op) {
        case SUM:
            rep_ = expansion::new_expansion_on_heap(
                expansion::sum_capacity(x,y)
            );
            rep_->assign_sum(x,y);
            break;
        case DIFF:
            rep_ = expansion::new_expansion_on_heap(
                expansion::diff_capacity(x,y)
            );
            rep_->assign_diff(x,y);
            break;
        case PRODUCT:
            rep_ = expansion::new_expansion_on_heap(
                expansion::product_capacity(x,y)
            );
            rep_->assign_product(x,y);
            break;
        }
    }
 
    explicit expansion_nt(
        Operation op,
        const expansion& x, const expansion& y, const expansion& z
    ) {
        switch(op) {
        case SUM:
            rep_ = expansion::new_expansion_on_heap(
                expansion::sum_capacity(x,y,z)
            );
            rep_->assign_sum(x,y,z);
            break;
        case DIFF:
            geo_assert_not_reached;
            break;
        case PRODUCT:
            rep_ = expansion::new_expansion_on_heap(
                expansion::product_capacity(x,y,z)
            );
            rep_->assign_product(x,y,z);
            break;
        }
    }
 
    explicit expansion_nt(
        Operation op,
        const expansion& x, const expansion& y,
        const expansion& z, const expansion& t
    ) {
        switch(op) {
        case SUM:
            rep_ = expansion::new_expansion_on_heap(
                expansion::sum_capacity(x,y,z,t)
            );
            rep_->assign_sum(x,y,z,t);
            break;
        case DIFF:
            geo_assert_not_reached;
            break;
        case PRODUCT:
            // HERE: TODO CHECK SIZE
            const expansion& p1 = expansion_product(x,y);
            const expansion& p2 = expansion_product(z,t);
            rep_ = expansion::new_expansion_on_heap(
                expansion::product_capacity(p1,p2)
            );
            rep_->assign_sum(p1,p2);
            break;
        }
    }
 
    explicit expansion_nt(Operation op, double x, double y) {
        switch(op) {
        case SUM:
            rep_ = expansion::new_expansion_on_heap(
                expansion::sum_capacity(x,y)
            );
            rep_->assign_sum(x,y);
            break;
        case DIFF:
            rep_ = expansion::new_expansion_on_heap(
                expansion::diff_capacity(x,y)
            );
            rep_->assign_diff(x,y);
            break;
        case PRODUCT:
            rep_ = expansion::new_expansion_on_heap(
                expansion::product_capacity(x,y)
            );
            rep_->assign_product(x,y);
            break;
        }
    }
 
    expansion_nt(const expansion_nt& rhs) {
        copy(rhs);
    }
 
    expansion_nt(expansion_nt&& rhs) {
        rep_ = nullptr;
        std::swap(rep_, rhs.rep_);
    }
 
    expansion_nt& operator= (const expansion_nt& rhs) {
        if(&rhs != this) {
            cleanup();
            copy(rhs);
        }
        return *this;
    }
 
    expansion_nt& operator= (expansion_nt&& rhs) {
        if(&rhs != this) {
            cleanup();
            std::swap(rep_, rhs.rep_);
        }
        return *this;
    }
 
    ~expansion_nt() {
        cleanup();
    }
 
    void optimize() {
        rep().optimize();
    }
 
    void negate() {
        rep().negate();
    }
 
    /********************************************************************/
 
    expansion_nt& operator+= (const expansion_nt& rhs);
 
    expansion_nt& operator-= (const expansion_nt& rhs);
 
    expansion_nt& operator*= (const expansion_nt& rhs);
 
    expansion_nt& operator+= (double rhs);
 
    expansion_nt& operator-= (double rhs);
 
    expansion_nt& operator*= (double rhs);
 
    /********************************************************************/
 
    expansion_nt operator+ (const expansion_nt& rhs) const;
 
    expansion_nt operator- (const expansion_nt& rhs) const;
 
    expansion_nt operator* (const expansion_nt& rhs) const;
 
    expansion_nt operator+ (double rhs) const;
 
    expansion_nt operator- (double rhs) const;
 
    expansion_nt operator* (double rhs) const;
 
    /********************************************************************/
 
    expansion_nt operator- () const;
 
    /********************************************************************/
 
    Sign compare(const expansion_nt& rhs) const {
        return rep().compare(rhs.rep());
    }
 
    Sign compare(double rhs) const {
        return rep().compare(rhs);
    }
 
    bool operator> (const expansion_nt& rhs) const {
        return (int(compare(rhs))>0);
    }
 
    bool operator>= (const expansion_nt& rhs) const {
        return (int(compare(rhs))>=0);
    }
 
    bool operator< (const expansion_nt& rhs) const {
        return (int(compare(rhs))<0);
    }
 
    bool operator<= (const expansion_nt& rhs) const {
        return (int(compare(rhs))<=0);
    }
 
    bool operator> (double rhs) const {
        return (int(compare(rhs))>0);
    }
 
    bool operator>= (double rhs) const {
        return (int(compare(rhs))>=0);
    }
 
    bool operator< (double rhs) const {
        return (int(compare(rhs))<0);
    }
 
    bool operator<= (double rhs) const {
        return (int(compare(rhs))<=0);
    }
 
    /********************************************************************/
 
    double estimate() const {
        return rep().estimate();
    }
 
    Sign sign() const {
        return rep().sign();
    }
 
    index_t length() const {
        return rep().length();
    }
 
    double component(index_t i) const {
        geo_debug_assert(i < length());
        return rep()[i];
    }
 
    expansion_nt(expansion* rep) :
    rep_(rep) {
    }
 
    expansion& rep() {
        return *rep_;
    }
 
    const expansion& rep() const {
        return *rep_;
    }
 
    std::string to_string() const {
        return (rep_ == nullptr) ?
            std::string("null") :
            rep_->to_string()   ;
    }
 
    bool is_one() const {
        return rep().equals(1.0);
    }
 
    bool is_zero() const {
        return sign() == ZERO;
    }
 
    protected:
 
    void copy(const expansion_nt& rhs) {
        if(rhs.rep_ == nullptr) {
            rep_ = nullptr;
        } else {
            rep_ = expansion::new_expansion_on_heap(rhs.rep().capacity());
            rep_->set_length(rhs.rep().length());
            for(index_t i=0; i<rep_->length(); ++i) {
                (*rep_)[i] = rhs.rep()[i];
            }
        }
    }
 
    void cleanup() {
        if(rep_ != nullptr) {
            expansion::delete_expansion_on_heap(rep_);
            rep_ = nullptr;
        }
    }
 
    private:
    expansion* rep_;
    friend expansion_nt operator- (double a, const expansion_nt& b);
 
    friend expansion_nt expansion_nt_sq_dist(
        const double* a, const double* b, coord_index_t dim
    );
 
    friend expansion_nt expansion_nt_dot_at(
        const double* a, const double* b, const double* c,
        coord_index_t dim
    );
    // friend class rational_nt;
    };
 
    inline expansion_nt operator+ (double a, const expansion_nt& b) {
        return b + a;
    }
 
    inline expansion_nt operator- (double a, const expansion_nt& b) {
        expansion_nt result = b - a;
        result.rep().negate();
        return result;
    }
 
    inline expansion_nt operator* (double a, const expansion_nt& b) {
        return b * a;
    }
 
    inline bool operator== (const expansion_nt& a, const expansion_nt& b) {
        return a.rep().equals(b.rep());
    }
 
    inline bool operator== (const expansion_nt& a, double b) {
        return a.rep().equals(b);
    }
 
    inline bool operator== (double a, const expansion_nt& b) {
        return b.rep().equals(a);
    }
 
    inline bool operator!= (const expansion_nt& a, const expansion_nt& b) {
        return !a.rep().equals(b.rep());
    }
 
    inline bool operator!= (const expansion_nt& a, double b) {
        return !a.rep().equals(b);
    }
 
    inline bool operator!= (double a, const expansion_nt& b) {
        return !b.rep().equals(a);
    }
 
    inline expansion_nt expansion_nt_sq_dist(
        const double* a, const double* b, coord_index_t dim
    ) {
        expansion* result = expansion::new_expansion_on_heap(
            expansion::sq_dist_capacity(dim)
        );
        result->assign_sq_dist(a, b, dim);
        return expansion_nt(result);
    }
 
    inline expansion_nt expansion_nt_dot_at(
        const double* a, const double* b, const double* c, coord_index_t dim
    ) {
        expansion* result = expansion::new_expansion_on_heap(
            expansion::dot_at_capacity(dim)
        );
        result->assign_dot_at(a, b, c, dim);
        return expansion_nt(result);
    }
 
    /************************************************************************/
 
    template <> inline Sign geo_sgn(const expansion_nt& x) {
        return x.sign();
    }
 
    template <> inline Sign geo_cmp(
        const expansion_nt& x, const expansion_nt& y
    ) {
        return x.compare(y);
    }
 
    /************************************************************************/
 
    inline bool expansion_nt_is_zero(const expansion_nt& x) {
        return (x.sign() == GEO::ZERO);
    }
 
    inline bool expansion_nt_is_one(const expansion_nt& x) {
        return x.rep().equals(1.0);
    }
 
 
    inline Sign expansion_nt_compare(
        const expansion_nt& x, const expansion_nt& y
    ) {
        const expansion& diff = expansion_diff(x.rep(), y.rep());
        return diff.sign();
    }
 
    inline expansion_nt expansion_nt_square(const expansion_nt& x) {
        expansion_nt result(
            expansion::new_expansion_on_heap(
                expansion::square_capacity(x.rep()
                                          ))
        );
        result.rep().assign_square(x.rep());
        return result;
    }
 
 
    expansion_nt GEOGRAM_API expansion_nt_determinant(
        const expansion_nt& a00,const expansion_nt& a01,
        const expansion_nt& a10,const expansion_nt& a11
    );
 
    expansion_nt GEOGRAM_API expansion_nt_determinant(
        const expansion_nt& a00,const expansion_nt& a01,const expansion_nt& a02,
        const expansion_nt& a10,const expansion_nt& a11,const expansion_nt& a12,
        const expansion_nt& a20,const expansion_nt& a21,const expansion_nt& a22
    );
 
    expansion_nt GEOGRAM_API expansion_nt_determinant(
        const expansion_nt& a00,const expansion_nt& a01,
        const expansion_nt& a02,const expansion_nt& a03,
        const expansion_nt& a10,const expansion_nt& a11,
        const expansion_nt& a12,const expansion_nt& a13,
        const expansion_nt& a20,const expansion_nt& a21,
        const expansion_nt& a22,const expansion_nt& a23,
        const expansion_nt& a30,const expansion_nt& a31,
        const expansion_nt& a32,const expansion_nt& a33
    );
 
// Make things a bit faster if target OS has large stack size
#ifdef GEO_HAS_BIG_STACK
 
    template <> inline expansion_nt det2x2(
        const expansion_nt& a11, const expansion_nt& a12,
        const expansion_nt& a21, const expansion_nt& a22
    ) {
        return expansion_nt_determinant(
            a11,a12,
            a21,a22
        );
    }
 
    template <> inline expansion_nt det3x3(
        const expansion_nt& a11, const expansion_nt& a12,
        const expansion_nt& a13,
        const expansion_nt& a21, const expansion_nt& a22,
        const expansion_nt& a23,
        const expansion_nt& a31, const expansion_nt& a32,
        const expansion_nt& a33
    ) {
        return expansion_nt_determinant(
            a11,a12,a13,
            a21,a22,a23,
            a31,a32,a33
        );
    }
 
    template <> inline expansion_nt det4x4(
        const expansion_nt& a11, const expansion_nt& a12,
        const expansion_nt& a13, const expansion_nt& a14,
        const expansion_nt& a21, const expansion_nt& a22,
        const expansion_nt& a23, const expansion_nt& a24,
        const expansion_nt& a31, const expansion_nt& a32,
        const expansion_nt& a33, const expansion_nt& a34,
        const expansion_nt& a41, const expansion_nt& a42,
        const expansion_nt& a43, const expansion_nt& a44
    ) {
        return expansion_nt_determinant(
            a11,a12,a13,a14,
            a21,a22,a23,a24,
            a31,a32,a33,a34,
            a41,a42,a43,a44
        );
    }
 
#endif
 
    /************************************************************************/
}
 
inline std::ostream& operator<< (
    std::ostream& os, const GEO::expansion_nt& a
) {
    return os << a.estimate();
}
 
inline std::istream& operator>> ( std::istream& is, GEO::expansion_nt& a) {
    double d;
    is >> d;
    if (is) {
        a = GEO::expansion_nt(d);
    }
    return is;
}
 
/*****************************************************************************/
 
namespace GEO {
 
    /**************************************************************************/
 
    namespace Numeric {
 
        template<> inline void optimize_number_representation(expansion_nt& x) {
            x.optimize();
        }
 
        template<> Sign GEOGRAM_API ratio_compare(
            const expansion_nt& a_num, const expansion_nt& a_denom,
            const expansion_nt& b_num, const expansion_nt& b_denom
        );
    }
 
    /**************************************************************************/
 
    typedef rationalg<expansion_nt> rational_nt;
 
    /**************************************************************************/
 
    template <> struct is_scalar<expansion_nt> {
        typedef expansion_nt type;
        static constexpr bool value = true;
    };
 
    /**************************************************************************/
}
 
#endif