From cb8e1c2f921741391e4d743a272a38657d37ffa4 Mon Sep 17 00:00:00 2001 From: Stephane Letz Date: Fri, 9 Aug 2024 17:09:29 +0200 Subject: [PATCH] Formatting. --- FaustAlgebra/FaustAlgebra.hh | 3 +- format-files | 6 +- interval/bitwiseOperations.cpp | 112 +++++--- interval/bitwiseOperations.hh | 29 +- interval/check.cpp | 482 +++++++++++++++++---------------- interval/check.hh | 18 +- interval/intervalAbs.cpp | 20 +- interval/intervalAcos.cpp | 11 +- interval/intervalAcosh.cpp | 14 +- interval/intervalAdd.cpp | 53 ++-- interval/intervalAnd.cpp | 29 +- interval/intervalAsin.cpp | 5 +- interval/intervalAsinh.cpp | 15 +- interval/intervalAtan.cpp | 11 +- interval/intervalAtan2.cpp | 45 +-- interval/intervalAtanh.cpp | 14 +- interval/intervalCeil.cpp | 7 +- interval/intervalCos.cpp | 38 ++- interval/intervalCosh.cpp | 10 +- interval/intervalDiv.cpp | 39 ++- interval/intervalEq.cpp | 21 +- interval/intervalFloor.cpp | 7 +- interval/intervalGe.cpp | 25 +- interval/intervalGt.cpp | 15 +- interval/intervalHSlider.cpp | 22 +- interval/intervalIntCast.cpp | 3 +- interval/intervalIntNum.cpp | 9 +- interval/intervalInv.cpp | 6 +- interval/intervalLe.cpp | 9 +- interval/intervalLog.cpp | 6 +- interval/intervalLog10.cpp | 8 +- interval/intervalLsh.cpp | 21 +- interval/intervalLt.cpp | 9 +- interval/intervalMax.cpp | 6 +- interval/intervalMin.cpp | 14 +- interval/intervalMissing.cpp | 31 ++- interval/intervalMod.cpp | 44 +-- interval/intervalMul.cpp | 54 ++-- interval/intervalNumEntry.cpp | 20 +- interval/intervalOr.cpp | 24 +- interval/intervalPow.cpp | 75 +++-- interval/intervalRint.cpp | 6 +- interval/intervalRound.cpp | 7 +- interval/intervalRsh.cpp | 25 +- interval/intervalSin.cpp | 54 ++-- interval/intervalSinh.cpp | 3 +- interval/intervalSqrt.cpp | 12 +- interval/intervalSub.cpp | 24 +- interval/intervalTan.cpp | 16 +- interval/intervalTanh.cpp | 8 +- interval/intervalVSlider.cpp | 20 +- interval/intervalXor.cpp | 54 ++-- interval/interval_algebra.hh | 27 +- interval/interval_def.hh | 25 +- interval/precision_utils.hh | 47 ++-- interval/utils.hh | 28 +- main.cpp | 43 ++- tlib/exception.hh | 8 +- tlib/garbageable.hh | 6 +- tlib/list.cpp | 3 +- tlib/node.hh | 39 ++- tlib/occur.hh | 2 +- tlib/recursive-tree.cpp | 12 +- tlib/shlysis.cpp | 4 +- tlib/symbol.cpp | 20 +- tlib/symbol.hh | 29 +- tlib/tree.cpp | 4 +- tlib/tree.hh | 80 +++--- 68 files changed, 1167 insertions(+), 829 deletions(-) diff --git a/FaustAlgebra/FaustAlgebra.hh b/FaustAlgebra/FaustAlgebra.hh index 5b2e095..2013e09 100644 --- a/FaustAlgebra/FaustAlgebra.hh +++ b/FaustAlgebra/FaustAlgebra.hh @@ -22,8 +22,7 @@ //===================================================================================================================== template -class FaustAlgebra -{ +class FaustAlgebra { public: //-------------------------------------------------------------------------------- // Dispatch tables diff --git a/format-files b/format-files index 5b12c71..4e08cc1 100755 --- a/format-files +++ b/format-files @@ -1,4 +1,4 @@ -find interval -path -prune -o -iname '*.cpp' -execdir clang-format -i -style=file {} \; -find interval -path -prune -o -iname '*.hh' -execdir clang-format -i -style=file {} \; -find interval -path -prune -o -iname '*.h' -execdir clang-format -i -style=file {} \; +find interval tlib FaustAlgebra main.cpp -path -prune -o -iname '*.cpp' -execdir clang-format -i -style=file {} \; +find interval tlib FaustAlgebra -path -prune -o -iname '*.hh' -execdir clang-format -i -style=file {} \; +find interval tlib FaustAlgebra -path -prune -o -iname '*.h' -execdir clang-format -i -style=file {} \; diff --git a/interval/bitwiseOperations.cpp b/interval/bitwiseOperations.cpp index 3a16de2..84b9a1e 100644 --- a/interval/bitwiseOperations.cpp +++ b/interval/bitwiseOperations.cpp @@ -1,11 +1,11 @@ /* Copyright 2023 Yann ORLAREY - * + * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at - * + * * http://www.apache.org/licenses/LICENSE-2.0 - * + * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. @@ -18,17 +18,22 @@ #include "bitwiseOperations.hh" -namespace itv -{ +namespace itv { /** * Split a signed interval into two unsigned intervals, for the negative and the positive part */ std::pair signSplit(const SInterval& x) { - if (isEmpty(x)) return {UEMPTY, UEMPTY}; - if (x.hi < 0) return {{(unsigned int)(x.lo), (unsigned int)(x.hi)}, UEMPTY}; - if (x.lo >= 0) return {UEMPTY, {(unsigned int)(x.lo), (unsigned int)(x.hi)}}; + if (isEmpty(x)) { + return {UEMPTY, UEMPTY}; + } + if (x.hi < 0) { + return {{(unsigned int)(x.lo), (unsigned int)(x.hi)}, UEMPTY}; + } + if (x.lo >= 0) { + return {UEMPTY, {(unsigned int)(x.lo), (unsigned int)(x.hi)}}; + } return {{(unsigned int)(x.lo), (unsigned int)(-1)}, {(unsigned int)(0), (unsigned int)(x.hi)}}; } @@ -38,7 +43,9 @@ std::pair signSplit(const SInterval& x) SInterval signMerge(const UInterval& np, const UInterval& pp) { if (isEmpty(np)) { - if (isEmpty(pp)) return SEMPTY; + if (isEmpty(pp)) { + return SEMPTY; + } return {(int)(pp.lo), (int)(pp.hi)}; } if (isEmpty(pp)) { @@ -69,10 +76,18 @@ static bool contains(const UInterval& i, unsi UInterval bitwiseUnsignedOr(const UInterval& a, const UInterval& b) { - if (a == UInterval{0, 0}) return b; - if (b == UInterval{0, 0}) return a; - if (isEmpty(a)) return a; - if (isEmpty(b)) return b; + if (a == UInterval{0, 0}) { + return b; + } + if (b == UInterval{0, 0}) { + return a; + } + if (isEmpty(a)) { + return a; + } + if (isEmpty(b)) { + return b; + } UInterval r{loOr2(a, b), hiOr2(a, b)}; return r; } @@ -86,8 +101,12 @@ static bool contains(const UInterval& i, unsigned int x) unsigned int hiOr2(UInterval a, UInterval b) { // simple cases - if (a.lo == 0 && a.hi == 0) return b.hi; - if (b.lo == 0 && b.hi == 0) return a.hi; + if (a.lo == 0 && a.hi == 0) { + return b.hi; + } + if (b.lo == 0 && b.hi == 0) { + return a.hi; + } // analyze and split the intervals auto [ma, a0, a1] = splitInterval(a); @@ -99,28 +118,44 @@ unsigned int hiOr2(UInterval a, UInterval b) } if (mb > ma) { - if (contains(a, mb - 1)) return 2 * mb - 1; + if (contains(a, mb - 1)) { + return 2 * mb - 1; + } return hiOr2(b1 - mb, a) + mb; } if (ma > mb) { - if (contains(b, ma - 1)) return 2 * ma - 1; + if (contains(b, ma - 1)) { + return 2 * ma - 1; + } return hiOr2(a1 - ma, b) + ma; } // ma == mb != 0 - if (isEmpty(a0) && isEmpty(b0)) return hiOr2(a1 - ma, b1 - ma) + ma; - if (isEmpty(a0)) return std::max(hiOr2(a1 - ma, b1 - ma), hiOr2(a1 - ma, b0)) + ma; - if (isEmpty(b0)) return std::max(hiOr2(a1 - ma, b1 - ma), hiOr2(a0, b1 - ma)) + ma; + if (isEmpty(a0) && isEmpty(b0)) { + return hiOr2(a1 - ma, b1 - ma) + ma; + } + if (isEmpty(a0)) { + return std::max(hiOr2(a1 - ma, b1 - ma), hiOr2(a1 - ma, b0)) + ma; + } + if (isEmpty(b0)) { + return std::max(hiOr2(a1 - ma, b1 - ma), hiOr2(a0, b1 - ma)) + ma; + } return std::max(hiOr2(a1 - ma, b1 - ma), std::max(hiOr2(a1 - ma, b0), hiOr2(a0, b1 - ma))) + ma; } unsigned int loOr2(UInterval a, UInterval b) { // isEmpty case - if (isEmpty(a) || isEmpty(b)) return 0; + if (isEmpty(a) || isEmpty(b)) { + return 0; + } // zero cases - if (a.lo == 0) return b.lo; - if (b.lo == 0) return a.lo; + if (a.lo == 0) { + return b.lo; + } + if (b.lo == 0) { + return a.lo; + } // non zero cases auto [ma, a0, a1] = splitInterval(a); @@ -129,17 +164,27 @@ unsigned int loOr2(UInterval a, UInterval b) // obvious cases if (ma > mb) { - if (isEmpty(a0)) return loOr2(a1 - ma, b) | ma; // ma bit unavoidable ! + if (isEmpty(a0)) { + return loOr2(a1 - ma, b) | ma; // ma bit unavoidable ! + } return loOr2(a0, b); } if (mb > ma) { - if (isEmpty(b0)) return loOr2(a, b1 - mb) | mb; + if (isEmpty(b0)) { + return loOr2(a, b1 - mb) | mb; + } return loOr2(a, b0); } // ma == mb != 0 - if (!isEmpty(a0) && !isEmpty(b0)) return loOr2(a0, b0); // obvious case ! - if (isEmpty(a0) && isEmpty(b0)) return loOr2(a1 - ma, b1 - ma) | ma; // ma bit unavoidable ! - if (isEmpty(a0)) return std::min(loOr2(a1 - ma, b0) | ma, loOr2(a1 - ma, b1 - ma) | ma); + if (!isEmpty(a0) && !isEmpty(b0)) { + return loOr2(a0, b0); // obvious case ! + } + if (isEmpty(a0) && isEmpty(b0)) { + return loOr2(a1 - ma, b1 - ma) | ma; // ma bit unavoidable ! + } + if (isEmpty(a0)) { + return std::min(loOr2(a1 - ma, b0) | ma, loOr2(a1 - ma, b1 - ma) | ma); + } return std::min(loOr2(a0, b1 - mb) | mb, loOr2(a1 - ma, b1 - ma) | ma); } @@ -169,11 +214,15 @@ unsigned int msb32(unsigned int x) // split interval according to its msb std::tuple splitInterval(UInterval x) { - if (x.lo == 0 && x.hi == 0) return {0, {1, 0}, x}; // special case, no msb + if (x.lo == 0 && x.hi == 0) { + return {0, {1, 0}, x}; // special case, no msb + } unsigned int m = msb32(x.hi); assert(m > 0); - if (m <= x.lo) return {m, {1, 0}, x}; // no msb in the interval + if (m <= x.lo) { + return {m, {1, 0}, x}; // no msb in the interval + } return {m, {x.lo, (unsigned int)(m - 1)}, {m, x.hi}}; } @@ -210,7 +259,8 @@ SInterval bitwiseSignedAnd(const SInterval& a, const SInterval& b) UInterval bitwiseUnsignedXOr(const UInterval& a, const UInterval& b) { - return bitwiseUnsignedAnd(bitwiseUnsignedOr(a, b), bitwiseUnsignedNot(bitwiseUnsignedAnd(a, b))); + return bitwiseUnsignedAnd(bitwiseUnsignedOr(a, b), + bitwiseUnsignedNot(bitwiseUnsignedAnd(a, b))); } SInterval bitwiseSignedXOr(const SInterval& a, const SInterval& b) diff --git a/interval/bitwiseOperations.hh b/interval/bitwiseOperations.hh index 00f8acf..962bcf2 100644 --- a/interval/bitwiseOperations.hh +++ b/interval/bitwiseOperations.hh @@ -1,12 +1,11 @@ #pragma once #include +#include #include #include -#include -namespace itv -{ +namespace itv { //============================================================================== // Definitions of signed and unsigned intervals for bitwise operations //============================================================================== @@ -41,7 +40,9 @@ bool isEmpty(const BitwiseInterval& i) template bool operator==(const BitwiseInterval& a, const BitwiseInterval& b) { - if (isEmpty(a)) return isEmpty(b); + if (isEmpty(a)) { + return isEmpty(b); + } return (a.lo == b.lo) && (a.hi == b.hi); } @@ -49,7 +50,9 @@ bool operator==(const BitwiseInterval& a, const BitwiseInterval& b) template bool operator!=(const BitwiseInterval& a, const BitwiseInterval& b) { - if (isEmpty(a)) return !isEmpty(b); + if (isEmpty(a)) { + return !isEmpty(b); + } return (a.lo != b.lo) || (a.hi != b.hi); } @@ -59,13 +62,17 @@ bool operator!=(const BitwiseInterval& a, const BitwiseInterval& b) inline std::ostream& operator<<(std::ostream& os, const UInterval& x) { - if (isEmpty(x)) return os << "U[]"; + if (isEmpty(x)) { + return os << "U[]"; + } return os << "U[" << x.lo << ".." << x.hi << "]"; } inline std::ostream& operator<<(std::ostream& os, const SInterval& x) { - if (isEmpty(x)) return os << "S[]"; + if (isEmpty(x)) { + return os << "S[]"; + } return os << "S[" << x.lo << ".." << x.hi << "]"; } @@ -77,8 +84,12 @@ inline std::ostream& operator<<(std::ostream& os, const SInterval& x) template BitwiseInterval operator+(const BitwiseInterval& a, const BitwiseInterval& b) { - if (isEmpty(a)) return b; - if (isEmpty(b)) return a; + if (isEmpty(a)) { + return b; + } + if (isEmpty(b)) { + return a; + } return {std::min(a.lo, b.lo), std::max(a.hi, b.hi)}; } diff --git a/interval/check.cpp b/interval/check.cpp index f5eac14..1e7f2f7 100644 --- a/interval/check.cpp +++ b/interval/check.cpp @@ -1,11 +1,11 @@ /* Copyright 2023 Yann ORLAREY - * + * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at - * + * * http://www.apache.org/licenses/LICENSE-2.0 - * + * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. @@ -35,12 +35,11 @@ void check(const std::string& expected, const itv::interval& exp) ss << exp; if (ss.str().compare(expected) == 0) { std::cout << "\033[32m" - << "OK: " << expected - << "\033[0m" << std::endl; + << "OK: " << expected << "\033[0m" << std::endl; } else { std::cout << "\033[31m" - << "ERR: We got " << ss.str() << " instead of " << expected - << "\033[0m" << std::endl; + << "ERR: We got " << ss.str() << " instead of " << expected << "\033[0m" + << std::endl; } } @@ -55,16 +54,16 @@ void check(const std::string& testname, const itv::interval& exp, const itv::int { if (exp == res) { std::cout << "\033[32m" - << "OK: " << testname << " " << exp << " = " << res - << "\033[0m" << std::endl; - if (exp.lsb() != res.lsb()) + << "OK: " << testname << " " << exp << " = " << res << "\033[0m" << std::endl; + if (exp.lsb() != res.lsb()) { std::cout << "\033[33m" - << "\t But precisions differ by " << exp.lsb() - res.lsb() - << "\033[0m" << std::endl; + << "\t But precisions differ by " << exp.lsb() - res.lsb() << "\033[0m" + << std::endl; + } } else { std::cout << "\033[31m" - << "ERR:" << testname << " FAILED. We got " << exp << " instead of " << res - << "\033[0m" << std::endl; + << "ERR:" << testname << " FAILED. We got " << exp << " instead of " << res + << "\033[0m" << std::endl; } } @@ -79,12 +78,11 @@ void check(const std::string& testname, bool exp, bool res) { if (exp == res) { std::cout << "\033[32m" - << "OK: " << testname - << "\033[0m" << std::endl; + << "OK: " << testname << "\033[0m" << std::endl; } else { std::cout << "\033[31m" - << "ERR:" << testname << " FAILED. We got " << exp << " instead of " << res - << "\033[0m" << std::endl; + << "ERR:" << testname << " FAILED. We got " << exp << " instead of " << res + << "\033[0m" << std::endl; } } @@ -99,7 +97,7 @@ void check(const std::string& testname, bool exp, bool res) */ itv::interval testfun(int N, bfun f, const itv::interval& x, const itv::interval& y) { - std::random_device R; // used to generate a random seed, based on some hardware randomness + std::random_device R; // used to generate a random seed, based on some hardware randomness std::default_random_engine generator((gRandom) ? R() : 12345); std::uniform_real_distribution rx(x.lo(), x.hi()); std::uniform_real_distribution ry(y.lo(), y.hi()); @@ -137,23 +135,24 @@ itv::interval testfun(int N, bfun f, const itv::interval& x, const itv::interval void analyzemod(itv::interval x, itv::interval y) { itv::interval_algebra A; - std::cout << "simulated fmod(" << x << "," << y << ") = " << testfun(10000, fmod, x, y) << std::endl; + std::cout << "simulated fmod(" << x << "," << y << ") = " << testfun(10000, fmod, x, y) + << std::endl; std::cout << "computed fmod(" << x << "," << y << ") = " << A.Mod(x, y) << std::endl; std::cout << std::endl; } /** * @brief Compute numerically the output interval of a function - * + * * @param E number of intervals/experiments * @param M number of measurements used to estimate the resulting interval * @param title name of the tested function * @param D maximal interval for the argument * @param f the numerical function of reference -*/ + */ void analyzeUnaryFunction(int E, int M, const char* title, const itv::interval& D, ufun f) { - std::random_device R; // used to generate a random seed, based on some hardware randomness + std::random_device R; // used to generate a random seed, based on some hardware randomness std::default_random_engine generator((gRandom) ? R() : 12345); std::uniform_real_distribution rd(D.lo(), D.hi()); @@ -191,31 +190,31 @@ void analyzeUnaryFunction(int E, int M, const char* title, const itv::interval& std::cout << std::endl; } - /** * @brief Check the unary interval function gives a good approximation of the numerical function - * + * * @param E number of intervals/experiments * @param M number of measurements used to estimate the resulting interval * @param title name of the tested function * @param D maximal interval for the argument * @param f the numerical function of reference * @param mp the interval method corresponding to f -*/ + */ void analyzeUnaryMethod(int E, int M, const char* title, const itv::interval& D, ufun f, umth mp) { - std::random_device R; // used to generate a random seed, based on some hardware randomness + std::random_device R; // used to generate a random seed, based on some hardware randomness std::default_random_engine generator((gRandom) ? R() : 12345); std::uniform_real_distribution rd(D.lo(), D.hi()); itv::interval_algebra A; - std::cout << "Analysis of " << title << " in domain " << D << " (u = " << pow(2, D.lsb()) << ")" << std::endl; + std::cout << "Analysis of " << title << " in domain " << D << " (u = " << pow(2, D.lsb()) << ")" + << std::endl; for (int e = 0; e < E; e++) { // E experiments // X: random input interval X < I - double a = truncate(rd(generator), D.lsb()); - double b = truncate(rd(generator), D.lsb()); + double a = truncate(rd(generator), D.lsb()); + double b = truncate(rd(generator), D.lsb()); itv::interval X(std::min(a, b), std::max(a, b), D.lsb()); // boundaries of the resulting interval @@ -232,7 +231,8 @@ void analyzeUnaryMethod(int E, int M, const char* title, const itv::interval& D, std::set measurements; // the loop has almost no chance of drawing X.hi(): we manually add it - double sample = X.hi(); // not truncated since morally the interval boundaries should already have the right precision + double sample = X.hi(); // not truncated since morally the interval boundaries should + // already have the right precision double y = f(sample); // y = truncate(y, -30); @@ -247,12 +247,13 @@ void analyzeUnaryMethod(int E, int M, const char* title, const itv::interval& D, } } - for (int m = 0; m < M; m++) { // M measurements - double presample = rx(generator); // non-truncated sample - sample = truncate(presample, D.lsb()); // truncated sample - double pre_y = f(presample); - y = f(sample); - // y = truncate(y, -30); // workaround to avoid artefacts in trigonometric functions + for (int m = 0; m < M; m++) { // M measurements + double presample = rx(generator); // non-truncated sample + sample = truncate(presample, D.lsb()); // truncated sample + double pre_y = f(presample); + y = f(sample); + // y = truncate(y, -30); // workaround to avoid artefacts in trigonometric + // functions measurements.insert(y); @@ -275,12 +276,10 @@ void analyzeUnaryMethod(int E, int M, const char* title, const itv::interval& D, double meas = *(measurements.begin()); - for (auto it = std::next(measurements.begin()); it != measurements.end(); ++it) - { + for (auto it = std::next(measurements.begin()); it != measurements.end(); ++it) { double next = *it; - double l = log2(next - meas); - if (l < lsb) - { + double l = log2(next - meas); + if (l < lsb) { lsb = floor(l); } @@ -288,19 +287,23 @@ void analyzeUnaryMethod(int E, int M, const char* title, const itv::interval& D, } itv::interval Y(y0, y1, lsb); - if (y0 > y1) Y = itv::empty(); // if we didn't manage to draw any samples + if (y0 > y1) { + Y = itv::empty(); // if we didn't manage to draw any samples + } itv::interval Z = (A.*mp)(X); if (Z >= Y and Z.lsb() <= Y.lsb()) { double precision = (Z.size() == 0) ? 1 : Y.size() / Z.size(); std::cout << "\033[32m" - << "OK " << e << ": " << title << "(" << X << ") = \n" << Z << "(c)\t >= \t" << Y - << "(m)\t (precision " << precision << ", LSB diff = " << Y.lsb() - Z.lsb() << ")" + << "OK " << e << ": " << title << "(" << X << ") = \n" + << Z << "(c)\t >= \t" << Y << "(m)\t (precision " << precision + << ", LSB diff = " << Y.lsb() - Z.lsb() << ")" << "\033[0m" << std::endl; } else { std::cout << "\033[31m" - << "ERROR " << e << ": " << title << "(" << X << ") = \n" << Z << "(c)\t INSTEAD OF \t" << Y + << "ERROR " << e << ": " << title << "(" << X << ") = \n" + << Z << "(c)\t INSTEAD OF \t" << Y << "(m), \t LSB diff = " << Y.lsb() - Z.lsb() << "\033[0m" << std::endl; } std::cout << std::endl; @@ -320,10 +323,10 @@ void analyzeUnaryMethod(int E, int M, const char* title, const itv::interval& D, * @param bm the interval method corresponding to f */ -void analyzeBinaryMethod(int E, int M, const char* title, const itv::interval& Dx, const itv::interval& Dy, bfun f, - bmth bm) +void analyzeBinaryMethod(int E, int M, const char* title, const itv::interval& Dx, + const itv::interval& Dy, bfun f, bmth bm) { - std::random_device R; // used to generate a random seed, based on some hardware randomness + std::random_device R; // used to generate a random seed, based on some hardware randomness std::default_random_engine generator((gRandom) ? R() : 12345); std::uniform_real_distribution rdx(Dx.lo(), Dx.hi()); std::uniform_real_distribution rdy(Dy.lo(), Dy.hi()); @@ -335,46 +338,32 @@ void analyzeBinaryMethod(int E, int M, const char* title, const itv::interval& D // store output values in order to measure the output precision std::set measurements; - if (Dx.lsb() < 0 or Dy.lsb() < 0) // if we're not doing an integer operation + if (Dx.lsb() < 0 or Dy.lsb() < 0) // if we're not doing an integer operation { - // X: random input interval X < Dx - double x0 = truncate(rdx(generator), Dx.lsb()); - double x1 = truncate(rdx(generator), Dx.lsb()); - itv::interval X(std::min(x0, x1), std::max(x0, x1), Dx.lsb()); - - // Y: random input interval Y < Dy - double y0 = truncate(rdy(generator), Dy.lsb()); - double y1 = truncate(rdy(generator), Dy.lsb()); - itv::interval Y(std::min(y0, y1), std::max(y0, y1), Dy.lsb()); - - // boundaries of the resulting interval Z - double zlo = HUGE_VAL; // std::min(t0, t1); - double zhi = -HUGE_VAL; // std::max(t0, t1); - - // precision of the resulting interval - int lsb = INT_MAX; - - // random values in X - std::uniform_real_distribution rvx(X.lo(), X.hi()); - std::uniform_real_distribution rvy(Y.lo(), Y.hi()); - - // draw the upper bounds manually - double z = f(X.hi(), Y.hi()); // no need to truncate: interval boundaries are already truncated - measurements.insert(z); - - if (!std::isnan(z)) { - if (z < zlo) { - zlo = z; - } - if (z > zhi) { - zhi = z; - } - } - - // measure the interval Z using the numerical function f - for (int m = 0; m < M; m++) { // M measurements - z = f(truncate(rvx(generator), Dx.lsb()), truncate(rvy(generator), Dy.lsb())); - + // X: random input interval X < Dx + double x0 = truncate(rdx(generator), Dx.lsb()); + double x1 = truncate(rdx(generator), Dx.lsb()); + itv::interval X(std::min(x0, x1), std::max(x0, x1), Dx.lsb()); + + // Y: random input interval Y < Dy + double y0 = truncate(rdy(generator), Dy.lsb()); + double y1 = truncate(rdy(generator), Dy.lsb()); + itv::interval Y(std::min(y0, y1), std::max(y0, y1), Dy.lsb()); + + // boundaries of the resulting interval Z + double zlo = HUGE_VAL; // std::min(t0, t1); + double zhi = -HUGE_VAL; // std::max(t0, t1); + + // precision of the resulting interval + int lsb = INT_MAX; + + // random values in X + std::uniform_real_distribution rvx(X.lo(), X.hi()); + std::uniform_real_distribution rvy(Y.lo(), Y.hi()); + + // draw the upper bounds manually + double z = f(X.hi(), + Y.hi()); // no need to truncate: interval boundaries are already truncated measurements.insert(z); if (!std::isnan(z)) { @@ -385,190 +374,209 @@ void analyzeBinaryMethod(int E, int M, const char* title, const itv::interval& D zhi = z; } } - } - double meas = *(measurements.begin()); + // measure the interval Z using the numerical function f + for (int m = 0; m < M; m++) { // M measurements + z = f(truncate(rvx(generator), Dx.lsb()), truncate(rvy(generator), Dy.lsb())); - for (auto it = std::next(measurements.begin()); it != measurements.end(); ++it) - { - double next = *it; - double l = log2(next - meas); - if (l < lsb) - { - lsb = floor(l); + measurements.insert(z); + + if (!std::isnan(z)) { + if (z < zlo) { + zlo = z; + } + if (z > zhi) { + zhi = z; + } + } } - meas = next; - } + double meas = *(measurements.begin()); - itv::interval Zm(zlo, zhi, lsb); // the measured Z - if (zlo > zhi) Zm = itv::empty(); // if we didn't manage to draw any samples - itv::interval Zc = (A.*bm)(X, Y); // the computed Z - double precision = (Zm.size() == Zc.size()) ? 1 : Zm.size() / Zc.size(); - - if (Zc >= Zm and Zc.lsb() <= Zm.lsb()) { - std::string color = "\033[32m"; - if (precision < 0.8 or Zm.lsb() - Zc.lsb() >= 10) color = "\033[36m"; // cyan instead of green if approximation is technically correct but of poor quality - std::cout << color - << "OK " << e << ": " << title << "(" << X << ",\t" << Y << ")\n =c=> " << Zc << "(c) >= " << Zm << "(m)" - << "\t (precision " << precision << "), \t LSB diff = " << Zm.lsb() - Zc.lsb() - << "\033[0m" << std::endl; - } else { - std::cout << "\033[31m" - << "ERROR " << e << ": " << title << "(" << X << ",\t" << Y << ")\n =c=> " << Zc << "(c) != " << Zm << "(m)" - << "\t LSB diff = " << Zm.lsb() - Zc.lsb() - << "\033[0m" << std::endl; - } - - } else { // integer operation - // std::cout << "Testing integer version of " << title << std::endl; - // X: random input interval X < Dx - double x0 = truncate(rdx(generator), Dx.lsb()); - double x1 = truncate(rdx(generator), Dx.lsb()); - itv::interval X(std::min(x0, x1), std::max(x0, x1), Dx.lsb()); - - // Y: random input interval Y < Dy - double y0 = truncate(rdy(generator), Dy.lsb()); - double y1 = truncate(rdy(generator), Dy.lsb()); - itv::interval Y(std::min(y0, y1), std::max(y0, y1), Dy.lsb()); - - // boundaries of the resulting interval Z - int zlo = INT_MAX; // std::min(t0, t1); - int zhi = INT_MIN; // std::max(t0, t1); + for (auto it = std::next(measurements.begin()); it != measurements.end(); ++it) { + double next = *it; + double l = log2(next - meas); + if (l < lsb) { + lsb = floor(l); + } - // precision of the resulting interval - int lsb = INT_MAX; + meas = next; + } - // random values in X and Y - std::uniform_int_distribution ivx((int)X.lo(), (int)X.hi()); - std::uniform_int_distribution ivy((int)Y.lo(), (int)Y.hi()); + itv::interval Zm(zlo, zhi, lsb); // the measured Z + if (zlo > zhi) { + Zm = itv::empty(); // if we didn't manage to draw any samples + } + itv::interval Zc = (A.*bm)(X, Y); // the computed Z + double precision = (Zm.size() == Zc.size()) ? 1 : Zm.size() / Zc.size(); + + if (Zc >= Zm and Zc.lsb() <= Zm.lsb()) { + std::string color = "\033[32m"; + if (precision < 0.8 or Zm.lsb() - Zc.lsb() >= 10) { + color = "\033[36m"; // cyan instead of green if approximation is technically + // correct but of poor quality + } + std::cout << color << "OK " << e << ": " << title << "(" << X << ",\t" << Y + << ")\n =c=> " << Zc << "(c) >= " << Zm << "(m)" + << "\t (precision " << precision + << "), \t LSB diff = " << Zm.lsb() - Zc.lsb() << "\033[0m" << std::endl; + } else { + std::cout << "\033[31m" + << "ERROR " << e << ": " << title << "(" << X << ",\t" << Y << ")\n =c=> " + << Zc << "(c) != " << Zm << "(m)" + << "\t LSB diff = " << Zm.lsb() - Zc.lsb() << "\033[0m" << std::endl; + } - // draw the upper bounds manually - int z = f(X.hi(), Y.hi()); // no need to truncate: interval boundaries are already truncated - measurements.insert((double)z); + } else { // integer operation + // std::cout << "Testing integer version of " << title << std::endl; + // X: random input interval X < Dx + double x0 = truncate(rdx(generator), Dx.lsb()); + double x1 = truncate(rdx(generator), Dx.lsb()); + itv::interval X(std::min(x0, x1), std::max(x0, x1), Dx.lsb()); - if (z < zlo) { - zlo = z; - } - if (z > zhi) { - zhi = z; - } + // Y: random input interval Y < Dy + double y0 = truncate(rdy(generator), Dy.lsb()); + double y1 = truncate(rdy(generator), Dy.lsb()); + itv::interval Y(std::min(y0, y1), std::max(y0, y1), Dy.lsb()); - // measure the interval Z using the numerical function f - for (int m = 0; m < M; m++) { // M measurements - int pre_x = ivx(generator); - int x = truncate(pre_x, Dx.lsb()); - int pre_y = ivy(generator); - int y = truncate(pre_y, Dy.lsb()); - z = f(x, y); - int pre_z = f(pre_x, pre_y); + // boundaries of the resulting interval Z + int zlo = INT_MAX; // std::min(t0, t1); + int zhi = INT_MIN; // std::max(t0, t1); - measurements.insert(z); + // precision of the resulting interval + int lsb = INT_MAX; - if (!std::isnan(pre_z)){ - if (z < zlo) { - zlo = pre_z; - } - if (z > zhi) { - zhi = pre_z; + // random values in X and Y + std::uniform_int_distribution ivx((int)X.lo(), (int)X.hi()); + std::uniform_int_distribution ivy((int)Y.lo(), (int)Y.hi()); + + // draw the upper bounds manually + int z = f(X.hi(), + Y.hi()); // no need to truncate: interval boundaries are already truncated + measurements.insert((double)z); + + if (z < zlo) { + zlo = z; + } + if (z > zhi) { + zhi = z; + } + + // measure the interval Z using the numerical function f + for (int m = 0; m < M; m++) { // M measurements + int pre_x = ivx(generator); + int x = truncate(pre_x, Dx.lsb()); + int pre_y = ivy(generator); + int y = truncate(pre_y, Dy.lsb()); + z = f(x, y); + int pre_z = f(pre_x, pre_y); + + measurements.insert(z); + + if (!std::isnan(pre_z)) { + if (z < zlo) { + zlo = pre_z; + } + if (z > zhi) { + zhi = pre_z; + } } } - } - double meas = *(measurements.begin()); + double meas = *(measurements.begin()); - for (auto it = std::next(measurements.begin()); it != measurements.end(); ++it) - { - double next = *it; - double l = log2(next - meas); - if (l < lsb) - { - lsb = floor(l); - } + for (auto it = std::next(measurements.begin()); it != measurements.end(); ++it) { + double next = *it; + double l = log2(next - meas); + if (l < lsb) { + lsb = floor(l); + } - meas = next; - } + meas = next; + } - itv::interval Zm(zlo, zhi, lsb); // the measured Z - itv::interval Zc = (A.*bm)(X, Y); // the computed Z - double precision = (Zm.size() == Zc.size()) ? 1 : Zm.size() / Zc.size(); + itv::interval Zm(zlo, zhi, lsb); // the measured Z + itv::interval Zc = (A.*bm)(X, Y); // the computed Z + double precision = (Zm.size() == Zc.size()) ? 1 : Zm.size() / Zc.size(); - if (Zc >= Zm and Zc.lsb() <= Zm.lsb()) { - std::string color = "\033[32m"; - if (precision < 0.8 or Zm.lsb() - Zc.lsb() >= 10) color = "\033[36m"; // cyan instead of green if approximation is technically correct but of poor quality - std::cout << color - << "OK " << e << ": " << title << "(" << X << ",\t" << Y << ")\n =c=> " << Zc << "(c) >= " << Zm << "(m)" - << "\t (precision " << precision << "), \t LSB diff = " << Zm.lsb() - Zc.lsb() - << "\033[0m" << std::endl; - } else { - std::cout << "\033[31m" - << "ERROR " << e << ": " << title << "(" << X << ",\t" << Y << ")\n =c=> " << Zc << "(c) != " << Zm << "(m)" - << "\t LSB diff = " << Zm.lsb() - Zc.lsb() - << "\033[0m" << std::endl; - } + if (Zc >= Zm and Zc.lsb() <= Zm.lsb()) { + std::string color = "\033[32m"; + if (precision < 0.8 or Zm.lsb() - Zc.lsb() >= 10) { + color = "\033[36m"; // cyan instead of green if approximation is technically + // correct but of poor quality + } + std::cout << color << "OK " << e << ": " << title << "(" << X << ",\t" << Y + << ")\n =c=> " << Zc << "(c) >= " << Zm << "(m)" + << "\t (precision " << precision + << "), \t LSB diff = " << Zm.lsb() - Zc.lsb() << "\033[0m" << std::endl; + } else { + std::cout << "\033[31m" + << "ERROR " << e << ": " << title << "(" << X << ",\t" << Y << ")\n =c=> " + << Zc << "(c) != " << Zm << "(m)" + << "\t LSB diff = " << Zm.lsb() - Zc.lsb() << "\033[0m" << std::endl; + } } - } std::cout << std::endl; } /** * @brief Adjusts the lsb of an input interval to match a target output lsb - * + * * @param title name of the tested function * @param mp the interval method of the studied function * @param X the input interval * @param l the target lsb for the output -*/ + */ void propagateBackwardsUnaryMethod(const char* title, umth mp, itv::interval& X, int l) { - std::cout << "Shaving input " << X << " of " << title << " to achieve an output lsb of " << l << std::endl; + std::cout << "Shaving input " << X << " of " << title << " to achieve an output lsb of " << l + << std::endl; itv::interval_algebra A; // itv::interval X = itv::interval(D.lo(), D.hi(), D.lsb()); itv::interval Z = (A.*mp)(X); - while (Z.lsb() < l) // the lsb of Z is more precise than l + while (Z.lsb() < l) // the lsb of Z is more precise than l { X = itv::interval(X.lo(), X.hi(), X.lsb() + 1); Z = (A.*mp)(X); std::cout << X.lsb() << " -> " << Z.lsb() << std::endl; } - if (Z.lsb() > l) // if we've overshot the target lsb + if (Z.lsb() > l) { // if we've overshot the target lsb X = itv::interval(X.lo(), X.hi(), X.lsb() - 1); + } std::cout << "Input interval " << X << " is sufficient" << std::endl; } /** * @brief Adjusts the lsbs of two input intervals to a binary function to match a target output lbs - * + * * @param title name of the tested function * @param bm the interval method of the studied function * @param X the first input interval * @param Y the second input interval * @param l the target lsb for the output -*/ -void propagateBackwardsBinaryMethod(const char* title, bmth bm, itv::interval& X, itv::interval& Y, int l) + */ +void propagateBackwardsBinaryMethod(const char* title, bmth bm, itv::interval& X, itv::interval& Y, + int l) { - std::cout << "Shaving inputs " << X << " and " << Y << " of " << title << " to achieve an output lsb of " << l << std::endl; + std::cout << "Shaving inputs " << X << " and " << Y << " of " << title + << " to achieve an output lsb of " << l << std::endl; itv::interval_algebra A; itv::interval Z = (A.*bm)(X, Y); - while (Z.lsb() < l) - { + while (Z.lsb() < l) { // std::cout << "X = " << X << "; Y = " << Y << std::endl; - if (X.lsb() < Y.lsb()) - { - X = itv::interval(X.lo(), X.hi(), X.lsb()+1); + if (X.lsb() < Y.lsb()) { + X = itv::interval(X.lo(), X.hi(), X.lsb() + 1); std::cout << "Shaving interval X = " << X << std::endl; - } - else - { - Y = itv::interval(Y.lo(), Y.hi(), Y.lsb()+1); + } else { + Y = itv::interval(Y.lo(), Y.hi(), Y.lsb() + 1); std::cout << "Shaving interval Y = " << Y << std::endl; } Z = (A.*bm)(X, Y); @@ -577,19 +585,19 @@ void propagateBackwardsBinaryMethod(const char* title, bmth bm, itv::interval& X std::cout << "Input intervals " << X << " and " << Y << " are sufficient" << std::endl; } - /** - * @brief Adjusts the lsb of an input iterval to a list of composed functions to match a target output lsb - * + * @brief Adjusts the lsb of an input iterval to a list of composed functions to match a target + * output lsb + * * @param titles names of the tested functions, from outermost to innermost * @param mps the interval methods of the functions, from outermost to innermost * @param X the input interval * @param l the target lsb for the output -*/ -void propagateBackwardsComposition(std::vector titles, std::vector mps, itv::interval& X, int l) + */ +void propagateBackwardsComposition(std::vector titles, std::vector mps, + itv::interval& X, int l) { - if (titles.size() != mps.size()) - { + if (titles.size() != mps.size()) { std::cout << "Incompatible vector sizes" << std::endl; return; } @@ -597,35 +605,37 @@ void propagateBackwardsComposition(std::vector titles, std::vector< unsigned long n = titles.size(); std::cout << "Shaving input " << X << " of "; - for(const auto *t: titles) + for (const auto* t : titles) { std::cout << t << " ○ "; + } std::cout << "\b\b\b"; std::cout << " to achieve an output lsb of " << l << std::endl << std::endl; - itv::interval_algebra A; - std::vector intermediate_intervals{X}; // should be one element bigger than titles and mps + itv::interval_algebra A; + std::vector intermediate_intervals{ + X}; // should be one element bigger than titles and mps - for(int i=0; i < n; i++) - { - intermediate_intervals.push_back((A.*(mps[n-i-1]))(intermediate_intervals[i])); - std::cout << titles[n-i-1] << "(" << intermediate_intervals[i] << ") = " << intermediate_intervals[i+1] << std::endl; + for (int i = 0; i < n; i++) { + intermediate_intervals.push_back((A.*(mps[n - i - 1]))(intermediate_intervals[i])); + std::cout << titles[n - i - 1] << "(" << intermediate_intervals[i] + << ") = " << intermediate_intervals[i + 1] << std::endl; } std::cout << std::endl << "Intermediate intervals before shaving:" << std::endl; - for(auto Y: intermediate_intervals) + for (auto Y : intermediate_intervals) { std::cout << Y << std::endl; + } std::cout << std::endl; int li = l; - - for(int i=0; i titles, std::vector< itv::interval Y = X; - for(int i=0; i titles, std::vector mps, itv::interval& X, int l); +void propagateBackwardsBinaryMethod(const char* title, bmth bm, itv::interval& X, itv::interval& Y, + int l); +void propagateBackwardsComposition(std::vector titles, std::vector mps, + itv::interval& X, int l); diff --git a/interval/intervalAbs.cpp b/interval/intervalAbs.cpp index 1da4f87..4ca523c 100644 --- a/interval/intervalAbs.cpp +++ b/interval/intervalAbs.cpp @@ -28,20 +28,19 @@ namespace itv { interval interval_algebra::Abs(const interval& x) { - if (x.isEmpty()) + if (x.isEmpty()) { return empty(); + } // precision stays the same if (x.lo() >= 0) { return x; } - + // integer overflowing - if (x.lsb() >= 0 and x.lo() <= (double)INT_MIN){ + if (x.lsb() >= 0 and x.lo() <= (double)INT_MIN) { double lo = (x.hi() >= 0) ? 0 : std::min(std::abs(x.hi()), (double)INT_MAX); - return {lo, - (double)INT_MAX, - x.lsb()}; + return {lo, (double)INT_MAX, x.lsb()}; } if (x.hi() <= 0) { return {-x.hi(), -x.lo(), x.lsb()}; @@ -54,9 +53,11 @@ void interval_algebra::testAbs() { std::cout << "Testing abs on finite intervals" << std::endl; analyzeUnaryMethod( - 10, 10000, "abs", interval(-10, 10, 0), [](double x) { return std::abs(x); }, &interval_algebra::Abs); + 10, 10000, "abs", interval(-10, 10, 0), [](double x) { return std::abs(x); }, + &interval_algebra::Abs); analyzeUnaryMethod( - 10, 10000, "abs", interval(-10, 10, -15), [](double x) { return std::abs(x); }, &interval_algebra::Abs); + 10, 10000, "abs", interval(-10, 10, -15), [](double x) { return std::abs(x); }, + &interval_algebra::Abs); std::cout << "Testing abs on intervals with limit bounds" << std::endl; check("abs", Abs(interval(INT_MIN, INT_MIN, 0)), interval(INT_MAX, INT_MAX, 0)); @@ -65,6 +66,7 @@ void interval_algebra::testAbs() check("abs", Abs(interval(-10, INT_MAX, 0)), interval(0, INT_MAX, 0)); check("abs", Abs(interval(INT_MIN, INT_MIN, 5)), interval(INT_MAX, INT_MAX, 5)); - check("abs", Abs(interval(INT_MIN, INT_MIN, -5)), interval(-(float)INT_MIN, -(float)INT_MIN, -5)); + check("abs", Abs(interval(INT_MIN, INT_MIN, -5)), + interval(-(float)INT_MIN, -(float)INT_MIN, -5)); } } // namespace itv \ No newline at end of file diff --git a/interval/intervalAcos.cpp b/interval/intervalAcos.cpp index 1ae13dd..072f4b7 100644 --- a/interval/intervalAcos.cpp +++ b/interval/intervalAcos.cpp @@ -30,14 +30,14 @@ static const interval AcosDomain(-1, 1, 0); // this interval needs 0 digits of interval interval_algebra::Acos(const interval& x) { - interval i = intersection(AcosDomain, x); // TODO: warn about interval violations + interval i = intersection(AcosDomain, x); // TODO: warn about interval violations if (i.isEmpty()) { return empty(); } double v = 0; // value at which the min slope is attained, zero if it is present - int sign = 1; // whether we compute the difference between f(x) and f(x+ε) or f(x-ε), chosing the point that lies - // in the interval + int sign = 1; // whether we compute the difference between f(x) and f(x+ε) or f(x-ε), chosing + // the point that lies in the interval if (not i.has(0)) // if zero is not present, it's the bound closer to zero { v = minValAbs(i); @@ -46,8 +46,9 @@ interval interval_algebra::Acos(const interval& x) int precision = exactPrecisionUnary(acos, v, sign * pow(2, i.lsb())); - if (precision == INT_MIN or taylor_lsb) - precision = floor(i.lsb() - (double)log2(1 - v*v)/2); + if (precision == INT_MIN or taylor_lsb) { + precision = floor(i.lsb() - (double)log2(1 - v * v) / 2); + } return {acos(i.hi()), acos(i.lo()), precision}; } diff --git a/interval/intervalAcosh.cpp b/interval/intervalAcosh.cpp index 27418a1..c1612d0 100644 --- a/interval/intervalAcosh.cpp +++ b/interval/intervalAcosh.cpp @@ -39,8 +39,9 @@ interval interval_algebra::Acosh(const interval& x) // we thus compute the gap between f(hi) and f(hi-ε), to remain in the interval int precision = exactPrecisionUnary(acosh, x.hi(), -pow(2, x.lsb())); - if (precision == INT_MIN or taylor_lsb) - precision = floor(x.lsb() - (double)log2(x.hi()*x.hi() - 1)/2); + if (precision == INT_MIN or taylor_lsb) { + precision = floor(x.lsb() - (double)log2(x.hi() * x.hi() - 1) / 2); + } return {acosh(i.lo()), acosh(i.hi()), precision}; } @@ -49,9 +50,12 @@ void interval_algebra::testAcosh() { analyzeUnaryMethod(10, 1000, "acosh", interval(950, 1000, 0), acosh, &interval_algebra::Acosh); analyzeUnaryMethod(10, 1000, "acosh", interval(950, 1000, -5), acosh, &interval_algebra::Acosh); - analyzeUnaryMethod(10, 1000, "acosh", interval(950, 1000, -10), acosh, &interval_algebra::Acosh); - analyzeUnaryMethod(10, 1000, "acosh", interval(950, 1000, -15), acosh, &interval_algebra::Acosh); - analyzeUnaryMethod(10, 1000, "acosh", interval(950, 1000, -20), acosh, &interval_algebra::Acosh); + analyzeUnaryMethod(10, 1000, "acosh", interval(950, 1000, -10), acosh, + &interval_algebra::Acosh); + analyzeUnaryMethod(10, 1000, "acosh", interval(950, 1000, -15), acosh, + &interval_algebra::Acosh); + analyzeUnaryMethod(10, 1000, "acosh", interval(950, 1000, -20), acosh, + &interval_algebra::Acosh); /* analyzeUnaryMethod(10, 1000, "acosh", interval(0, 10, 0), acosh, &interval_algebra::Acosh); analyzeUnaryMethod(10, 1000, "acosh", interval(0, 10, -5), acosh, &interval_algebra::Acosh); diff --git a/interval/intervalAdd.cpp b/interval/intervalAdd.cpp index fcbde07..b2e250a 100644 --- a/interval/intervalAdd.cpp +++ b/interval/intervalAdd.cpp @@ -41,7 +41,7 @@ interval interval_algebra::Add(const interval& x, const interval& y) return empty(); } - if (x.lsb() >= 0 and y.lsb() >= 0) { // if both intervals are integers + if (x.lsb() >= 0 and y.lsb() >= 0) { // if both intervals are integers // if we're dealing with integers the interval has to wrap around 0 const int xlo = (int)x.lo(); const int xhi = (int)x.hi(); @@ -49,24 +49,22 @@ interval interval_algebra::Add(const interval& x, const interval& y) const int yhi = (int)y.hi(); // detect wrapping - /* if (std::abs((double)xhi + (double)yhi) >= (double) INT_MAX - or std::abs((double)xhi + (double)yhi) <= (double) INT_MIN + /* if (std::abs((double)xhi + (double)yhi) >= (double) INT_MAX + or std::abs((double)xhi + (double)yhi) <= (double) INT_MIN or std::abs((double)xlo + (double)ylo) >= (double) INT_MAX or std::abs((double)xlo + (double)ylo) <= (double) INT_MIN) return {(double) INT_MIN, (double) INT_MAX, std::min(x.lsb(), y.lsb())};*/ - + double lo = x.lo() + y.lo(); double hi = x.hi() + y.hi(); // if there is a discontinuity by the lower end of integers - if (lo <= (double)INT_MIN -1 and hi >= (double)INT_MIN) - { + if (lo <= (double)INT_MIN - 1 and hi >= (double)INT_MIN) { return {(double)INT_MIN, (double)INT_MAX, std::min(x.lsb(), y.lsb())}; } - + // if there is a discontinuity by the higher end of integers - if (lo <= (double)INT_MAX and hi >= (double)INT_MAX+1) - { + if (lo <= (double)INT_MAX and hi >= (double)INT_MAX + 1) { return {(double)INT_MIN, (double)INT_MAX, std::min(x.lsb(), y.lsb())}; } @@ -75,8 +73,8 @@ interval interval_algebra::Add(const interval& x, const interval& y) } return {x.lo() + y.lo(), x.hi() + y.hi(), - std::min(x.lsb(), - y.lsb())}; // the result of an addition needs to be only as precise as the most precise of the operands + std::min(x.lsb(), y.lsb())}; // the result of an addition needs to be only as precise + // as the most precise of the operands } void interval_algebra::testAdd() @@ -85,20 +83,33 @@ void interval_algebra::testAdd() check("test algebra Add", Add(interval(0, 100), interval(10, 500)), interval(10, 600)); std::cout << "Testing add on positive intervals" << std::endl; - analyzeBinaryMethod(5, 2000, "add", interval(0, 10, 0), interval(0, 10, 0), add, &interval_algebra::Add); - analyzeBinaryMethod(5, 2000, "add", interval(0, 10, -10), interval(0, 10, 0), add, &interval_algebra::Add); - analyzeBinaryMethod(5, 2000, "add", interval(0, 10, 0), interval(0, 10, -10), add, &interval_algebra::Add); - analyzeBinaryMethod(5, 2000, "add", interval(0, 10, -10), interval(0, 10, -10), add, &interval_algebra::Add); + analyzeBinaryMethod(5, 2000, "add", interval(0, 10, 0), interval(0, 10, 0), add, + &interval_algebra::Add); + analyzeBinaryMethod(5, 2000, "add", interval(0, 10, -10), interval(0, 10, 0), add, + &interval_algebra::Add); + analyzeBinaryMethod(5, 2000, "add", interval(0, 10, 0), interval(0, 10, -10), add, + &interval_algebra::Add); + analyzeBinaryMethod(5, 2000, "add", interval(0, 10, -10), interval(0, 10, -10), add, + &interval_algebra::Add); std::cout << "Testing add on negative intervals" << std::endl; - analyzeBinaryMethod(5, 2000, "add", interval(-10, 10, -5), interval(-10, 0, -5), add, &interval_algebra::Add); - analyzeBinaryMethod(5, 2000, "add", interval(-10, 0, -5), interval(-10, 0, -5), add, &interval_algebra::Add); + analyzeBinaryMethod(5, 2000, "add", interval(-10, 10, -5), interval(-10, 0, -5), add, + &interval_algebra::Add); + analyzeBinaryMethod(5, 2000, "add", interval(-10, 0, -5), interval(-10, 0, -5), add, + &interval_algebra::Add); std::cout << "Testing add with wrapping" << std::endl; - analyzeBinaryMethod(10, 2000, "add", interval(0, pow(2, 31), 0), interval(0, pow(2, 31), 0), addint, &interval_algebra::Add); - analyzeBinaryMethod(10, 2000, "add", interval((double)INT_MAX/2, (double)INT_MAX, 0), interval((double)INT_MAX/2, (double)INT_MAX, 0), addint, &interval_algebra::Add); - analyzeBinaryMethod(10, 2000, "add", interval((double)INT_MIN, (double)INT_MIN/2, 0), interval((double)INT_MAX/2, (double)INT_MAX, 0), addint, &interval_algebra::Add); - analyzeBinaryMethod(10, 2000, "add", interval((double)INT_MIN, (double)INT_MIN/2, 0), interval((double)INT_MIN, (double)INT_MIN/2, 0), addint, &interval_algebra::Add); + analyzeBinaryMethod(10, 2000, "add", interval(0, pow(2, 31), 0), interval(0, pow(2, 31), 0), + addint, &interval_algebra::Add); + analyzeBinaryMethod(10, 2000, "add", interval((double)INT_MAX / 2, (double)INT_MAX, 0), + interval((double)INT_MAX / 2, (double)INT_MAX, 0), addint, + &interval_algebra::Add); + analyzeBinaryMethod(10, 2000, "add", interval((double)INT_MIN, (double)INT_MIN / 2, 0), + interval((double)INT_MAX / 2, (double)INT_MAX, 0), addint, + &interval_algebra::Add); + analyzeBinaryMethod(10, 2000, "add", interval((double)INT_MIN, (double)INT_MIN / 2, 0), + interval((double)INT_MIN, (double)INT_MIN / 2, 0), addint, + &interval_algebra::Add); } } // namespace itv diff --git a/interval/intervalAnd.cpp b/interval/intervalAnd.cpp index ad600f6..8c417ad 100644 --- a/interval/intervalAnd.cpp +++ b/interval/intervalAnd.cpp @@ -141,7 +141,8 @@ interval interval_algebra::And(const interval& x, const interval& y) int precision = std::min(x.lsb(), y.lsb()); - /* int precision = std::max(x.lsb(), y.lsb()); // output precision cannot be finer than that of the input intervals + /* int precision = std::max(x.lsb(), y.lsb()); // output precision cannot be finer than that of + the input intervals // however, if one of the intervals is reduced to one element, the mask can make it so int precisionx = 0; @@ -166,21 +167,29 @@ interval interval_algebra::And(const interval& x, const interval& y) } }*/ - return {double(z.lo), double(z.hi), + return {double(z.lo), double(z.hi), // std::max(precision, std::max(precisionx, precisiony)) precision}; } void interval_algebra::testAnd() { - analyzeBinaryMethod(10, 2000, "And", interval(0, 257, 0), singleton(12), myAnd, &interval_algebra::And); - analyzeBinaryMethod(10, 2000, "And", interval(-1000, -800, 0), singleton(12), myAnd, &interval_algebra::And); - analyzeBinaryMethod(10, 2000, "And", interval(-1000, -800, 0), singleton(12), myAnd, &interval_algebra::And); - analyzeBinaryMethod(10, 2000, "And", interval(-128, 128, 0), singleton(127), myAnd, &interval_algebra::And); - analyzeBinaryMethod(10, 2000, "And", interval(0, 1000, 0), interval(63, 127), myAnd, &interval_algebra::And); - analyzeBinaryMethod(10, 2000, "And", interval(-1000, 1000, 0), interval(63, 127), myAnd, &interval_algebra::And); - analyzeBinaryMethod(10, 2000, "And", interval(10, 20, 0), singleton(0), myAnd, &interval_algebra::And); - analyzeBinaryMethod(10, 2000, "And", singleton(0), interval(15, 25, 0), myAnd, &interval_algebra::And); + analyzeBinaryMethod(10, 2000, "And", interval(0, 257, 0), singleton(12), myAnd, + &interval_algebra::And); + analyzeBinaryMethod(10, 2000, "And", interval(-1000, -800, 0), singleton(12), myAnd, + &interval_algebra::And); + analyzeBinaryMethod(10, 2000, "And", interval(-1000, -800, 0), singleton(12), myAnd, + &interval_algebra::And); + analyzeBinaryMethod(10, 2000, "And", interval(-128, 128, 0), singleton(127), myAnd, + &interval_algebra::And); + analyzeBinaryMethod(10, 2000, "And", interval(0, 1000, 0), interval(63, 127), myAnd, + &interval_algebra::And); + analyzeBinaryMethod(10, 2000, "And", interval(-1000, 1000, 0), interval(63, 127), myAnd, + &interval_algebra::And); + analyzeBinaryMethod(10, 2000, "And", interval(10, 20, 0), singleton(0), myAnd, + &interval_algebra::And); + analyzeBinaryMethod(10, 2000, "And", singleton(0), interval(15, 25, 0), myAnd, + &interval_algebra::And); analyzeBinaryMethod(10, 2000, "And", singleton(0), singleton(0), myAnd, &interval_algebra::And); } } // namespace itv diff --git a/interval/intervalAsin.cpp b/interval/intervalAsin.cpp index 0e37e23..a67eda2 100644 --- a/interval/intervalAsin.cpp +++ b/interval/intervalAsin.cpp @@ -47,8 +47,9 @@ interval interval_algebra::Asin(const interval& x) } int precision = exactPrecisionUnary(asin, v, sign * pow(2, i.lsb())); - if (precision == INT_MIN or taylor_lsb) - precision = floor(x.lsb() - (double)log2(1 - v*v)/2); + if (precision == INT_MIN or taylor_lsb) { + precision = floor(x.lsb() - (double)log2(1 - v * v) / 2); + } return {asin(i.lo()), asin(i.hi()), precision}; } diff --git a/interval/intervalAsinh.cpp b/interval/intervalAsinh.cpp index be978b7..e363504 100644 --- a/interval/intervalAsinh.cpp +++ b/interval/intervalAsinh.cpp @@ -29,17 +29,20 @@ static const interval domain(-HUGE_VAL, HUGE_VAL); interval interval_algebra::Asinh(const interval& x) { - if (x.isEmpty()) + if (x.isEmpty()) { return empty(); + } - double v = maxValAbs(x); // value at which the min slope is attained, here the bound of highest absolute value - int sign = signMaxValAbs(x); // whether we compute the difference between f(v) and f(v+ε) or f(v-ε) + double v = maxValAbs( + x); // value at which the min slope is attained, here the bound of highest absolute value + int sign = + signMaxValAbs(x); // whether we compute the difference between f(v) and f(v+ε) or f(v-ε) int precision = exactPrecisionUnary(asinh, v, sign * pow(2, x.lsb())); - if (precision == INT_MIN or taylor_lsb) - precision = floor(x.lsb() - (double)log2(1 + v*v)/2); - + if (precision == INT_MIN or taylor_lsb) { + precision = floor(x.lsb() - (double)log2(1 + v * v) / 2); + } return {asinh(x.lo()), asinh(x.hi()), precision}; } diff --git a/interval/intervalAtan.cpp b/interval/intervalAtan.cpp index 6839e05..0a2a60c 100644 --- a/interval/intervalAtan.cpp +++ b/interval/intervalAtan.cpp @@ -35,13 +35,16 @@ interval interval_algebra::Atan(const interval& x) return empty(); } - double v = maxValAbs(x); // value at which the min slope is attained, here the bound of highest absolute value - int sign = signMaxValAbs(x); // whether we compute the difference between f(v) and f(v+ε) or f(v-ε) + double v = maxValAbs( + x); // value at which the min slope is attained, here the bound of highest absolute value + int sign = + signMaxValAbs(x); // whether we compute the difference between f(v) and f(v+ε) or f(v-ε) int precision = exactPrecisionUnary(atan, v, sign * pow(2, x.lsb())); - if (precision == INT_MIN or taylor_lsb) - precision = floor(x.lsb() - (double)log2(1 + v*v)); + if (precision == INT_MIN or taylor_lsb) { + precision = floor(x.lsb() - (double)log2(1 + v * v)); + } return {atan(x.lo()), atan(x.hi()), precision}; } diff --git a/interval/intervalAtan2.cpp b/interval/intervalAtan2.cpp index 0f9a24f..58fb004 100644 --- a/interval/intervalAtan2.cpp +++ b/interval/intervalAtan2.cpp @@ -31,21 +31,22 @@ namespace itv { // (where (x,y) are the cartesian coordinates of the point we wish to retrieve the angle of) interval interval_algebra::Atan2(const interval& y, const interval& x) { - if (x.isEmpty() || y.isEmpty()) + if (x.isEmpty() || y.isEmpty()) { return empty(); + } double lo = -M_PI; double hi = M_PI; - // atan2 is continuous on the plane except on Ox- = {(x,y)| x<=0 and y=0} where the angle gap happens - // if the domain spans the discontinuity, we split in along the Ox axis - // in order to have a continuous function on each domain - // we study it on each of the sub-domains and then combine the results + // atan2 is continuous on the plane except on Ox- = {(x,y)| x<=0 and y=0} where the angle gap + // happens if the domain spans the discontinuity, we split in along the Ox axis in order to have + // a continuous function on each domain we study it on each of the sub-domains and then combine + // the results // atan2(y, x) = atan(y/x) + constant: precision is that of y/x compounded with that of atan // cf https://en.wikipedia.org/wiki/Atan2#Definition_and_computation - if (y.lo() <= 0 and x.hasZero()) // if we intersect the Ox- axis + if (y.lo() <= 0 and x.hasZero()) // if we intersect the Ox- axis { /* interval yp = {0, y.hi(), y.lsb()}; // positive part of y interval yn = {y.lo(), 0, y.lsb()}; // negative part of y*/ @@ -56,19 +57,22 @@ interval interval_algebra::Atan2(const interval& y, const interval& x) interval dp = interval_algebra::Div(y, xp); interval dn = interval_algebra::Div(y, xn); - int precisionp = exactPrecisionUnary(atan, maxValAbs(dp), signMaxValAbs(dp) * pow(2, dp.lsb())); - int precisionn = exactPrecisionUnary(atan, maxValAbs(dn), signMaxValAbs(dn) * pow(2, dn.lsb())); + int precisionp = + exactPrecisionUnary(atan, maxValAbs(dp), signMaxValAbs(dp) * pow(2, dp.lsb())); + int precisionn = + exactPrecisionUnary(atan, maxValAbs(dn), signMaxValAbs(dn) * pow(2, dn.lsb())); return {lo, hi, - std::min(precisionp, - precisionn)}; // final precision is the finest precision attained on either of the domains + std::min(precisionp, precisionn)}; // final precision is the finest precision + // attained on either of the domains } - interval d = interval_algebra::Div(y, x); - int precision = exactPrecisionUnary(atan, maxValAbs(d), signMaxValAbs(d) * pow(2, d.lsb())); + interval d = interval_algebra::Div(y, x); + int precision = exactPrecisionUnary(atan, maxValAbs(d), signMaxValAbs(d) * pow(2, d.lsb())); // highest angle between a point of XxY and the x-axis - if (y.lo() >= 0) // the domain XxY is entirely included in the higher half of the plane, where the angle is highest + if (y.lo() >= 0) // the domain XxY is entirely included in the higher half of the plane, where + // the angle is highest { if (x.lo() <= 0) { // we intersect the quadrant in which atan2 takes the highest values hi = atan2(y.lo(), x.lo()); @@ -88,7 +92,8 @@ interval interval_algebra::Atan2(const interval& y, const interval& x) } // lowest angle between a point of XxY and the x-axis - if (y.hi() <= 0) // the domain XxY is entirely included in the lower half of the plane, where the angle is highest + if (y.hi() <= 0) // the domain XxY is entirely included in the lower half of the plane, where + // the angle is highest { if (x.lo() <= 0) { lo = atan2(y.hi(), x.lo()); @@ -113,11 +118,13 @@ interval interval_algebra::Atan2(const interval& y, const interval& x) void interval_algebra::testAtan2() { // std::cout << "Atan2 not implemented" << std::endl; - /* analyzeBinaryMethod(10, 1000000, "atan2", interval(1, 2, -24), interval(1, 2, -24), atan2, &interval_algebra::Atan2); - analyzeBinaryMethod(10, 1000000, "atan2", interval(-1, 2, -24), interval(1, 2, -24), atan2, &interval_algebra::Atan2); - analyzeBinaryMethod(10, 1000000, "atan2", interval(-2, -1, -24), interval(1, 2, -24), atan2, &interval_algebra::Atan2); - analyzeBinaryMethod(10, 1000000, "atan2", interval(-2, -1, -24), interval(-1, 2, -24), atan2, &interval_algebra::Atan2); - analyzeBinaryMethod(10, 1000000, "atan2", interval(-2, -1, -24), interval(-2, -1, -24), atan2, &interval_algebra::Atan2);*/ + /* analyzeBinaryMethod(10, 1000000, "atan2", interval(1, 2, -24), interval(1, 2, -24), atan2, + &interval_algebra::Atan2); analyzeBinaryMethod(10, 1000000, "atan2", interval(-1, 2, -24), + interval(1, 2, -24), atan2, &interval_algebra::Atan2); analyzeBinaryMethod(10, 1000000, "atan2", + interval(-2, -1, -24), interval(1, 2, -24), atan2, &interval_algebra::Atan2); + analyzeBinaryMethod(10, 1000000, "atan2", interval(-2, -1, -24), interval(-1, 2, -24), atan2, + &interval_algebra::Atan2); analyzeBinaryMethod(10, 1000000, "atan2", interval(-2, -1, -24), + interval(-2, -1, -24), atan2, &interval_algebra::Atan2);*/ analyzeBinaryMethod(10, 1000000, "atan2", interval(-1, 2, -24), interval(-1, 2, -24), atan2, &interval_algebra::Atan2); diff --git a/interval/intervalAtanh.cpp b/interval/intervalAtanh.cpp index 32eaa06..9164041 100644 --- a/interval/intervalAtanh.cpp +++ b/interval/intervalAtanh.cpp @@ -27,7 +27,8 @@ namespace itv { // interval Atanh(const interval& x); // void testAtanh(); -static const interval domain(std::nexttoward(-1, 0), std::nexttoward(1, 0), 0); // interval ]-1,1[, precision 0 +static const interval domain(std::nexttoward(-1, 0), std::nexttoward(1, 0), + 0); // interval ]-1,1[, precision 0 interval interval_algebra::Atanh(const interval& x) { @@ -36,12 +37,13 @@ interval interval_algebra::Atanh(const interval& x) return empty(); } - double v = minValAbs(x); - double sign = signMinValAbs(x); - int precision = exactPrecisionUnary(atanh, v, sign*pow(2, x.lsb())); + double v = minValAbs(x); + double sign = signMinValAbs(x); + int precision = exactPrecisionUnary(atanh, v, sign * pow(2, x.lsb())); - if (precision == INT_MIN or taylor_lsb) - precision = floor(x.lsb() - (double)log2(1 - v*v)); + if (precision == INT_MIN or taylor_lsb) { + precision = floor(x.lsb() - (double)log2(1 - v * v)); + } return {atanh(i.lo()), atanh(i.hi()), precision}; } diff --git a/interval/intervalCeil.cpp b/interval/intervalCeil.cpp index e2e1187..ed84b50 100644 --- a/interval/intervalCeil.cpp +++ b/interval/intervalCeil.cpp @@ -31,9 +31,10 @@ interval interval_algebra::Ceil(const interval& x) if (x.isEmpty()) { return empty(); } - return {ceil(x.lo()), ceil(x.hi()), -1 }; // even though the output of floor are mathematical integers, - // they are implemented as floats and thus should not be given precision 0, - // lest it be cast as an int + return {ceil(x.lo()), ceil(x.hi()), + -1}; // even though the output of floor are mathematical integers, + // they are implemented as floats and thus should not be given precision 0, + // lest it be cast as an int } void interval_algebra::testCeil() diff --git a/interval/intervalCos.cpp b/interval/intervalCos.cpp index 1e188ee..ebfcaa9 100644 --- a/interval/intervalCos.cpp +++ b/interval/intervalCos.cpp @@ -39,16 +39,19 @@ interval interval_algebra::Cos(const interval& x) } int precision = exactPrecisionUnary(cos, 0, pow(2, x.lsb())); - if (precision == INT_MIN or taylor_lsb) precision = 2*x.lsb() - 1; // if x.lsb() is so small that the automatic computation doesn't work + if (precision == INT_MIN or taylor_lsb) { + precision = + 2 * x.lsb() - 1; // if x.lsb() is so small that the automatic computation doesn't work + } - if (x.size() >= 2*M_PI) { + if (x.size() >= 2 * M_PI) { return {-1, 1, precision}; } // normalize input interval between 0..2 (corresponding to 0..2PI) - double l = fmod(x.lo(), 2*M_PI); + double l = fmod(x.lo(), 2 * M_PI); if (l < 0) { - l += 2*M_PI; + l += 2 * M_PI; } interval i(l, l + x.size(), x.lsb()); @@ -59,20 +62,21 @@ interval interval_algebra::Cos(const interval& x) double hi = std::max(a, b); // check if integers are included - if (i.has(0) || i.has(2*M_PI)) { + if (i.has(0) || i.has(2 * M_PI)) { hi = 1; } - if (i.has(1*M_PI) || i.has(3*M_PI)) { + if (i.has(1 * M_PI) || i.has(3 * M_PI)) { lo = -1; } double v = 0; // value of the interval at which the finest precision is computed - if (i.hi() < 1*M_PI or - (i.lo() > 1*M_PI and i.hi() < 2*M_PI)) // if there are no integers in i, i.e i is included in ]0;1[ or ]1;2[ + if (i.hi() < 1 * M_PI or + (i.lo() > 1 * M_PI and + i.hi() < 2 * M_PI)) // if there are no integers in i, i.e i is included in ]0;1[ or ]1;2[ { - double delta_hi = ceil(x.hi()/M_PI) - x.hi()/M_PI; - double delta_lo = x.lo()/M_PI - floor(x.lo()/M_PI); + double delta_hi = ceil(x.hi() / M_PI) - x.hi() / M_PI; + double delta_lo = x.lo() / M_PI - floor(x.lo() / M_PI); if (delta_hi < delta_lo) { // if the lowest slope is attained for the higher bound v = x.hi(); } else { // ... for the lower bound @@ -81,13 +85,17 @@ interval interval_algebra::Cos(const interval& x) } precision = exactPrecisionUnary(cos, v, pow(2, x.lsb())); - if (precision == INT_MIN or taylor_lsb) - { + if (precision == INT_MIN or taylor_lsb) { /* cos(x + u) - cos(x) = - u·sin(x) if x != 0 = - u^2/2 · cos(x) = -u^2/2 if x == 0*/ - if (v != 0) precision = x.lsb() + (int)floor(log2(abs(sin(v))));// (int)floor(log2(M_PI*abs(cos(M_PI*v)))) + x.lsb(); - else precision = 2*x.lsb() -1 ; //- (int)floor(2*log2(M_PI)); + if (v != 0) { + precision = + x.lsb() + (int)floor(log2( + abs(sin(v)))); // (int)floor(log2(M_PI*abs(cos(M_PI*v)))) + x.lsb(); + } else { + precision = 2 * x.lsb() - 1; //- (int)floor(2*log2(M_PI)); + } } return {lo, hi, precision}; @@ -95,7 +103,7 @@ interval interval_algebra::Cos(const interval& x) void interval_algebra::testCos() { - analyzeUnaryMethod(10, 40000, "cos", interval(0, 2*M_PI, -3), cos, &interval_algebra::Cos); + analyzeUnaryMethod(10, 40000, "cos", interval(0, 2 * M_PI, -3), cos, &interval_algebra::Cos); analyzeUnaryMethod(10, 40000, "cos", interval(0, M_PI, -5), cos, &interval_algebra::Cos); analyzeUnaryMethod(10, 40000, "cos", interval(0, M_PI, -10), cos, &interval_algebra::Cos); analyzeUnaryMethod(10, 40000, "cos", interval(0, M_PI, -15), cos, &interval_algebra::Cos); diff --git a/interval/intervalCosh.cpp b/interval/intervalCosh.cpp index 34802cc..ee5649f 100644 --- a/interval/intervalCosh.cpp +++ b/interval/intervalCosh.cpp @@ -36,8 +36,9 @@ interval interval_algebra::Cosh(const interval& x) if (x.hasZero()) { int precision = exactPrecisionUnary(cosh, v, pow(2, x.lsb())); - if (precision == INT_MIN or taylor_lsb) - precision = floor(2*x.lsb() - 1); // cosh(u) - cosh(0) = u^2/2 + o(u^2) + if (precision == INT_MIN or taylor_lsb) { + precision = floor(2 * x.lsb() - 1); // cosh(u) - cosh(0) = u^2/2 + o(u^2) + } return {1, std::max(cosh(x.lo()), cosh(x.hi())), precision}; } @@ -53,8 +54,9 @@ interval interval_algebra::Cosh(const interval& x) } int precision = exactPrecisionUnary(cosh, v, sign * pow(2, x.lsb())); - if (precision == INT_MIN or taylor_lsb) - precision = floor(x.lsb() + log2(abs(sinh(v)))); // cosh(x+u) - cosh(x) = u sinh(x) + o(u) + if (precision == INT_MIN or taylor_lsb) { + precision = floor(x.lsb() + log2(abs(sinh(v)))); // cosh(x+u) - cosh(x) = u sinh(x) + o(u) + } return {std::min(cosh(x.lo()), cosh(x.hi())), std::max(cosh(x.lo()), cosh(x.hi())), precision}; } diff --git a/interval/intervalDiv.cpp b/interval/intervalDiv.cpp index b4a847f..bd8d4bd 100644 --- a/interval/intervalDiv.cpp +++ b/interval/intervalDiv.cpp @@ -38,24 +38,33 @@ double div(double x, double y) void interval_algebra::testDiv() { // lots of experiments because of the quadratic size of the input - analyzeBinaryMethod(10, 5000000, "Div", interval(-100, 100, -15), interval(0.001, 1000, -15), div, + analyzeBinaryMethod(10, 5000000, "Div", interval(-100, 100, -15), interval(0.001, 1000, -15), + div, &interval_algebra::Div); + analyzeBinaryMethod(10, 5000000, "Div", interval(-100, 100, -15), interval(-1000, -0.001, -15), + div, &interval_algebra::Div); + + analyzeBinaryMethod(10, 500000, "div", interval(0, 10, 0), interval(0, 10, 0), div, &interval_algebra::Div); - analyzeBinaryMethod(10, 5000000, "Div", interval(-100, 100, -15), interval(-1000, -0.001, -15), div, + analyzeBinaryMethod(10, 500000, "div", interval(0, 10, -5), interval(0, 10, 0), div, + &interval_algebra::Div); + analyzeBinaryMethod(10, 500000, "div", interval(0, 10, -10), interval(0, 10, 0), div, + &interval_algebra::Div); + analyzeBinaryMethod(10, 500000, "div", interval(0, 10, -15), interval(0, 10, 0), div, + &interval_algebra::Div); + analyzeBinaryMethod(10, 500000, "div", interval(0, 10, 0), interval(0, 10, -10), div, + &interval_algebra::Div); + analyzeBinaryMethod(10, 500000, "div", interval(0, 10, -5), interval(0, 10, -10), div, + &interval_algebra::Div); + analyzeBinaryMethod(10, 500000, "div", interval(0, 10, -10), interval(0, 10, -10), div, + &interval_algebra::Div); + analyzeBinaryMethod(10, 500000, "div", interval(0, 10, -15), interval(0, 10, -10), div, &interval_algebra::Div); - - analyzeBinaryMethod(10, 500000, "div", interval(0, 10, 0), interval(0, 10, 0), div, &interval_algebra::Div); - analyzeBinaryMethod(10, 500000, "div", interval(0, 10, -5), interval(0, 10, 0), div, &interval_algebra::Div); - analyzeBinaryMethod(10, 500000, "div", interval(0, 10, -10), interval(0, 10, 0), div, &interval_algebra::Div); - analyzeBinaryMethod(10, 500000, "div", interval(0, 10, -15), interval(0, 10, 0), div, &interval_algebra::Div); - analyzeBinaryMethod(10, 500000, "div", interval(0, 10, 0), interval(0, 10, -10), div, &interval_algebra::Div); - analyzeBinaryMethod(10, 500000, "div", interval(0, 10, -5), interval(0, 10, -10), div, &interval_algebra::Div); - analyzeBinaryMethod(10, 500000, "div", interval(0, 10, -10), interval(0, 10, -10), div, &interval_algebra::Div); - analyzeBinaryMethod(10, 500000, "div", interval(0, 10, -15), interval(0, 10, -10), div, &interval_algebra::Div); // check("test algebra Div", Div(interval(-2, 3), interval(1, 10)), interval(-2, 3)); - // check("test algebra Div", Div(interval(-2, 3), interval(-1, 10)), interval(-HUGE_VAL, +HUGE_VAL)); - // check("test algebra Div", Div(interval(-2, 3), interval(-0.1, -0.01)), interval(-300, 200)); - // check("test algebra Div", Div(interval(0), interval(0)), interval(0)); - // check("test algebra Div", Div(interval(0, 1), interval(0, 1)), interval(0, +HUGE_VAL)); + // check("test algebra Div", Div(interval(-2, 3), interval(-1, 10)), interval(-HUGE_VAL, + // +HUGE_VAL)); check("test algebra Div", Div(interval(-2, 3), interval(-0.1, -0.01)), + // interval(-300, 200)); check("test algebra Div", Div(interval(0), interval(0)), + // interval(0)); check("test algebra Div", Div(interval(0, 1), interval(0, 1)), interval(0, + // +HUGE_VAL)); } } // namespace itv diff --git a/interval/intervalEq.cpp b/interval/intervalEq.cpp index 6792d49..9af14eb 100644 --- a/interval/intervalEq.cpp +++ b/interval/intervalEq.cpp @@ -34,9 +34,15 @@ static double myEq(double x, double y) interval interval_algebra::Eq(const interval& x, const interval& y) { // boolean value => precision 0 - if (x.isEmpty() || y.isEmpty()) return empty(); - if (x.lo() == x.hi() && x.lo() == y.lo() && x.lo() == y.hi()) return interval{1,1,0}; - if (x.hi() < y.lo() || x.lo() > y.hi()) return interval{0, 0, 0}; + if (x.isEmpty() || y.isEmpty()) { + return empty(); + } + if (x.lo() == x.hi() && x.lo() == y.lo() && x.lo() == y.hi()) { + return interval{1, 1, 0}; + } + if (x.hi() < y.lo() || x.lo() > y.hi()) { + return interval{0, 0, 0}; + } return interval{0, 1, 0}; } @@ -46,8 +52,11 @@ void interval_algebra::testEq() check("test algebra Eq", Eq(interval(2, 5), interval(0, 1)), interval(0)); check("test algebra Eq", Eq(interval(-1, 1), interval(0, 10)), interval(0, 1));*/ - analyzeBinaryMethod(10, 200, "Eq", interval(-1, 1, 0), interval(-1, 1, 0), myEq, &interval_algebra::Eq); - analyzeBinaryMethod(10, 200, "Eq", interval(-10, 10, 0), interval(-10, 10, 0), myEq, &interval_algebra::Eq); - analyzeBinaryMethod(10, 20000, "Eq", interval(-1, 1, -2), interval(-1, 1, -2), myEq, &interval_algebra::Eq); + analyzeBinaryMethod(10, 200, "Eq", interval(-1, 1, 0), interval(-1, 1, 0), myEq, + &interval_algebra::Eq); + analyzeBinaryMethod(10, 200, "Eq", interval(-10, 10, 0), interval(-10, 10, 0), myEq, + &interval_algebra::Eq); + analyzeBinaryMethod(10, 20000, "Eq", interval(-1, 1, -2), interval(-1, 1, -2), myEq, + &interval_algebra::Eq); } } // namespace itv diff --git a/interval/intervalFloor.cpp b/interval/intervalFloor.cpp index 0aafdcd..768c406 100644 --- a/interval/intervalFloor.cpp +++ b/interval/intervalFloor.cpp @@ -31,9 +31,10 @@ interval interval_algebra::Floor(const interval& x) if (x.isEmpty()) { return empty(); } - return {floor(x.lo()), floor(x.hi()), -1}; // even though the output of floor are mathematical integers, - // they are implemented as floats and thus should not be given precision 0, - // lest it be cast as an int + return {floor(x.lo()), floor(x.hi()), + -1}; // even though the output of floor are mathematical integers, + // they are implemented as floats and thus should not be given precision 0, + // lest it be cast as an int } void interval_algebra::testFloor() diff --git a/interval/intervalGe.cpp b/interval/intervalGe.cpp index 8d8c729..f6fcb60 100644 --- a/interval/intervalGe.cpp +++ b/interval/intervalGe.cpp @@ -34,9 +34,15 @@ static double myGe(double x, double y) interval interval_algebra::Ge(const interval& x, const interval& y) { // boolean value => precision 0 - if (x.isEmpty() || y.isEmpty()) return empty(); - if (x.lo() >= y.hi()) return interval{1,1,0}; - if (x.hi() < y.lo()) return interval{0,0,0}; + if (x.isEmpty() || y.isEmpty()) { + return empty(); + } + if (x.lo() >= y.hi()) { + return interval{1, 1, 0}; + } + if (x.hi() < y.lo()) { + return interval{0, 0, 0}; + } return interval{0, 1, 0}; } @@ -44,10 +50,13 @@ void interval_algebra::testGe() { /* check("test algebra Ge", Ge(interval(5), interval(5)), interval(1)); check("test algebra Ge", Ge(interval(2, 5), interval(0, 1)), interval(1)); - check("test algebra Ge", Ge(interval(-1, 1), interval(0, 10)), interval(0, 1));*/ - - analyzeBinaryMethod(10, 200, "Ge", interval(-1, 1, 0), interval(-1, 1, 0), myGe, &interval_algebra::Ge); - analyzeBinaryMethod(10, 200, "Ge", interval(-10, 10, 0), interval(-10, 10, 0), myGe, &interval_algebra::Ge); - analyzeBinaryMethod(10, 2000, "Ge", interval(-10, 10), interval(-10, 10), myGe, &interval_algebra::Ge); + check("test algebra Ge", Ge(interval(-1, 1), interval(0, 10)), interval(0, 1));*/ + + analyzeBinaryMethod(10, 200, "Ge", interval(-1, 1, 0), interval(-1, 1, 0), myGe, + &interval_algebra::Ge); + analyzeBinaryMethod(10, 200, "Ge", interval(-10, 10, 0), interval(-10, 10, 0), myGe, + &interval_algebra::Ge); + analyzeBinaryMethod(10, 2000, "Ge", interval(-10, 10), interval(-10, 10), myGe, + &interval_algebra::Ge); } } // namespace itv diff --git a/interval/intervalGt.cpp b/interval/intervalGt.cpp index 5c89a29..6f09514 100644 --- a/interval/intervalGt.cpp +++ b/interval/intervalGt.cpp @@ -27,7 +27,7 @@ namespace itv { // void testGt(); static double myGt(double x, double y) -{ +{ return (x > y); } @@ -37,10 +37,10 @@ interval interval_algebra::Gt(const interval& x, const interval& y) return empty(); } if (x.lo() > y.hi()) { - return interval{1,1,0}; + return interval{1, 1, 0}; } if (x.hi() <= y.lo()) { - return interval{0,0,0}; + return interval{0, 0, 0}; } return interval{0, 1, 0}; } @@ -53,8 +53,11 @@ void interval_algebra::testGt() check("test algebra Gt", Gt(interval(2, 5), interval(5, 20)), interval(0)); check("test algebra Gt", Gt(interval(2, 5), interval(0, 20)), interval(0, 1)); */ - analyzeBinaryMethod(10, 200, "Gt", interval(-1, 1, 0), interval(-1, 1, 0), myGt, &interval_algebra::Gt); - analyzeBinaryMethod(10, 200, "Gt", interval(-10, 10, 0), interval(-10, 10, 0), myGt, &interval_algebra::Gt); - analyzeBinaryMethod(10, 2000, "Gt", interval(-10, 10), interval(-10, 10), myGt, &interval_algebra::Gt); + analyzeBinaryMethod(10, 200, "Gt", interval(-1, 1, 0), interval(-1, 1, 0), myGt, + &interval_algebra::Gt); + analyzeBinaryMethod(10, 200, "Gt", interval(-10, 10, 0), interval(-10, 10, 0), myGt, + &interval_algebra::Gt); + analyzeBinaryMethod(10, 2000, "Gt", interval(-10, 10), interval(-10, 10), myGt, + &interval_algebra::Gt); } } // namespace itv diff --git a/interval/intervalHSlider.cpp b/interval/intervalHSlider.cpp index 85c14ae..af8b183 100644 --- a/interval/intervalHSlider.cpp +++ b/interval/intervalHSlider.cpp @@ -21,19 +21,23 @@ namespace itv { //------------------------------------------------------------------------------------------ // Interval IntNum -interval interval_algebra::HSlider(const interval& name, const interval& init, const interval& lo, const interval& hi, - const interval& step) +interval interval_algebra::HSlider(const interval& name, const interval& init, const interval& lo, + const interval& hi, const interval& step) { - if (init.isEmpty() || lo.isEmpty() || hi.isEmpty() || step.isEmpty()) + if (init.isEmpty() || lo.isEmpty() || hi.isEmpty() || step.isEmpty()) { return empty(); + } - // elements of a slider with range [lo; hi] and step s are of the form lo + k·s <= hi with k an integer - // the precision needed to represent such elements is the minimum between - int lsb = std::min(step.lsb(), lo.lsb()); // the precision of the lower bound and that of the step - if (step.lo() > 0) { // if we don't have negative or zero steps - lsb = std::min(lsb, (int)log2(step.lo())); // and that associated to the smallest value the step can take + // elements of a slider with range [lo; hi] and step s are of the form lo + k·s <= hi with k an + // integer the precision needed to represent such elements is the minimum between + int lsb = + std::min(step.lsb(), lo.lsb()); // the precision of the lower bound and that of the step + if (step.lo() > 0) { // if we don't have negative or zero steps + lsb = std::min( + lsb, + (int)log2(step.lo())); // and that associated to the smallest value the step can take } - return {lo.lo(), hi.hi(), lsb}; + return {lo.lo(), hi.hi(), lsb}; } } // namespace itv diff --git a/interval/intervalIntCast.cpp b/interval/intervalIntCast.cpp index 4ad3aa2..b2a77f6 100644 --- a/interval/intervalIntCast.cpp +++ b/interval/intervalIntCast.cpp @@ -40,6 +40,7 @@ interval interval_algebra::IntCast(const interval& x) void interval_algebra::testIntCast() { check("test algebra IntCast", IntCast(interval{-3.8, 4.9}), interval{-3.0, 4.0, 0}); - check("test algebra IntCast", IntCast(interval{-HUGE_VAL, HUGE_VAL}), interval{-2147483648.0, 2147483647.0, 0}); + check("test algebra IntCast", IntCast(interval{-HUGE_VAL, HUGE_VAL}), + interval{-2147483648.0, 2147483647.0, 0}); } } // namespace itv diff --git a/interval/intervalIntNum.cpp b/interval/intervalIntNum.cpp index 558d8eb..8ae7fa1 100644 --- a/interval/intervalIntNum.cpp +++ b/interval/intervalIntNum.cpp @@ -23,9 +23,11 @@ namespace itv { interval interval_algebra::IntNum(int x) { - /* int lsb = -24; // lsb_number(x); // x is an integer so lsb is bound to be >=0, but we might be able to shave a + /* int lsb = -24; // lsb_number(x); // x is an integer so lsb is bound to be >=0, but we might + be able to shave a // couple more bits - // actually this is a terrible idea, in regard of the consequences on later inferred precisions + // actually this is a terrible idea, in regard of the consequences on later + inferred precisions // why should _,16:*:sin have a different output precision than _,15:*:sin ? // shaving a couple extra bits is not worth the later consequences @@ -39,7 +41,8 @@ interval interval_algebra::IntNum(int x) interval interval_algebra::Int64Num(int64_t x) { - /* int lsb = -24; // lsb_number(x); // x is an integer so lsb is bound to be >=0, but we might be able to shave a + /* int lsb = -24; // lsb_number(x); // x is an integer so lsb is bound to be >=0, but we might + be able to shave a // couple more bits while (floor(x * pow(2, -lsb - 1)) == x * pow(2, -lsb - 1) and x != 0) { diff --git a/interval/intervalInv.cpp b/interval/intervalInv.cpp index 8d96ef0..f7bfd31 100644 --- a/interval/intervalInv.cpp +++ b/interval/intervalInv.cpp @@ -39,7 +39,7 @@ interval interval_algebra::Inv(const interval& x) } int sign = signMaxValAbs(x); - double v = maxValAbs(x); // precision is computed at the bound with the highest absolute value + double v = maxValAbs(x); // precision is computed at the bound with the highest absolute value // if v is infinite, this means that images of interval elements will get closer and closer // but since we're writing the fixed-point numbers on at most 31 bits of MSB @@ -51,7 +51,7 @@ interval interval_algebra::Inv(const interval& x) int precision = exactPrecisionUnary(inv, v, sign * pow(2, x.lsb())); if (precision == INT_MIN or taylor_lsb) { - precision = floor(x.lsb() - 2*log2(abs(v))); // 1/(x+u) - 1/x = -u/x^2 + o(u) + precision = floor(x.lsb() - 2 * log2(abs(v))); // 1/(x+u) - 1/x = -u/x^2 + o(u) } // precision = std::max(precision, -31); @@ -63,7 +63,7 @@ interval interval_algebra::Inv(const interval& x) return {-HUGE_VAL, 1.0 / x.lo(), precision}; // normally, we don't divide by 0 // so the bounds of the image interval are not infinity but 1/ulp - // return {-pow(2, -x.lsb()), 1.0 / x.lo(), precision}; + // return {-pow(2, -x.lsb()), 1.0 / x.lo(), precision}; } if (x.lo() == 0 && x.hi() > 0) { return {1 / x.hi(), HUGE_VAL, precision}; diff --git a/interval/intervalLe.cpp b/interval/intervalLe.cpp index 396e640..78fa59a 100644 --- a/interval/intervalLe.cpp +++ b/interval/intervalLe.cpp @@ -42,8 +42,11 @@ void interval_algebra::testLe() check("test algebra Le", Le(interval(2, 5), interval(0, 1)), interval(0)); check("test algebra Le", Le(interval(-1, 1), interval(0, 10)), interval(0, 1));*/ - analyzeBinaryMethod(10, 200, "Le", interval(-1, 1, 0), interval(-1, 1, 0), myLe, &interval_algebra::Le); - analyzeBinaryMethod(10, 200, "Le", interval(-10, 10, 0), interval(-10, 10, 0), myLe, &interval_algebra::Le); - analyzeBinaryMethod(10, 2000, "Le", interval(-10, 10), interval(-10, 10), myLe, &interval_algebra::Le); + analyzeBinaryMethod(10, 200, "Le", interval(-1, 1, 0), interval(-1, 1, 0), myLe, + &interval_algebra::Le); + analyzeBinaryMethod(10, 200, "Le", interval(-10, 10, 0), interval(-10, 10, 0), myLe, + &interval_algebra::Le); + analyzeBinaryMethod(10, 2000, "Le", interval(-10, 10), interval(-10, 10), myLe, + &interval_algebra::Le); } } // namespace itv diff --git a/interval/intervalLog.cpp b/interval/intervalLog.cpp index b077d1d..94bff18 100644 --- a/interval/intervalLog.cpp +++ b/interval/intervalLog.cpp @@ -38,9 +38,9 @@ interval interval_algebra::Log(const interval& x) // lowest slope is at the highest bound of the interval int precision = exactPrecisionUnary( - log, i.hi(), -pow(2, i.lsb())); // -pow because we take the FP number right before the higher bound - if (precision == INT_MIN or taylor_lsb) - { + log, i.hi(), + -pow(2, i.lsb())); // -pow because we take the FP number right before the higher bound + if (precision == INT_MIN or taylor_lsb) { /* double delta = -log(1 - pow(2, i.lsb()) / i.hi()); precision = floor((double)log2(delta));*/ precision = floor(i.lsb() - (double)log2(abs(i.hi()))); diff --git a/interval/intervalLog10.cpp b/interval/intervalLog10.cpp index 95f645c..d2124b8 100644 --- a/interval/intervalLog10.cpp +++ b/interval/intervalLog10.cpp @@ -34,10 +34,12 @@ interval interval_algebra::Log10(const interval& x) // lowest slope is at the highest bound of the interval int precision = exactPrecisionUnary( - log10, x.hi(), -pow(2, x.lsb())); // -pow because we take the FP number right before the higher bound - if (precision == INT_MIN or taylor_lsb) + log10, x.hi(), + -pow(2, x.lsb())); // -pow because we take the FP number right before the higher bound + if (precision == INT_MIN or taylor_lsb) { precision = floor(x.lsb() - (double)log2(abs(x.hi())) - log2(log(10))); - + } + interval i = intersection(interval(0, HUGE_VAL), x); return {log10(i.lo()), log10(i.hi()), precision}; } diff --git a/interval/intervalLsh.cpp b/interval/intervalLsh.cpp index 182a723..5a5049b 100644 --- a/interval/intervalLsh.cpp +++ b/interval/intervalLsh.cpp @@ -34,14 +34,16 @@ static double lsh(double x, double k) interval interval_algebra::Lsh(const interval& x, const interval& k) { - if (x.isEmpty() || k.isEmpty()) + if (x.isEmpty() || k.isEmpty()) { return empty(); + } interval j{pow(2, k.lo()), pow(2, k.hi())}; interval z = Mul(x, j); return {z.lo(), z.hi(), - x.lsb() + (int)k.lo()}; // lshifts shave some precision off the numbers, at least y.lo() bits + x.lsb() + + (int)k.lo()}; // lshifts shave some precision off the numbers, at least y.lo() bits } void interval_algebra::testLsh() @@ -49,10 +51,15 @@ void interval_algebra::testLsh() /* check("test algebra Lsh", Lsh(interval(0, 1), interval(4)), interval(0, 16)); check("test algebra Lsh", Lsh(interval(0.5, 1), interval(-1, 4)), interval(0.25, 16)); check("test algebra Lsh", Lsh(interval(-10, 10), interval(-1, 4)), interval(-160, 160));*/ - analyzeBinaryMethod(10, 100000, "lshift", interval(0, 32, 0), interval(8, 8, 1), lsh, &interval_algebra::Lsh); - analyzeBinaryMethod(10, 100000, "lshift", interval(0, 1024, 0), interval(-10, 10, 0), lsh, &interval_algebra::Lsh); - analyzeBinaryMethod(10, 100000, "lshift", interval(0, 1024, 2), interval(-10, 10, 0), lsh, &interval_algebra::Lsh); - analyzeBinaryMethod(10, 100000, "lshift", interval(0, 1024, 0), interval(-10, 10, 1), lsh, &interval_algebra::Lsh); - analyzeBinaryMethod(10, 100000, "lshift", interval(0, 1024, 2), interval(-10, 10, 1), lsh, &interval_algebra::Lsh); + analyzeBinaryMethod(10, 100000, "lshift", interval(0, 32, 0), interval(8, 8, 1), lsh, + &interval_algebra::Lsh); + analyzeBinaryMethod(10, 100000, "lshift", interval(0, 1024, 0), interval(-10, 10, 0), lsh, + &interval_algebra::Lsh); + analyzeBinaryMethod(10, 100000, "lshift", interval(0, 1024, 2), interval(-10, 10, 0), lsh, + &interval_algebra::Lsh); + analyzeBinaryMethod(10, 100000, "lshift", interval(0, 1024, 0), interval(-10, 10, 1), lsh, + &interval_algebra::Lsh); + analyzeBinaryMethod(10, 100000, "lshift", interval(0, 1024, 2), interval(-10, 10, 1), lsh, + &interval_algebra::Lsh); } } // namespace itv diff --git a/interval/intervalLt.cpp b/interval/intervalLt.cpp index bd8c7b7..0392721 100644 --- a/interval/intervalLt.cpp +++ b/interval/intervalLt.cpp @@ -42,8 +42,11 @@ void interval_algebra::testLt() check("test algebra Lt", Lt(interval(2, 5), interval(0, 2)), interval(0)); check("test algebra Lt", Lt(interval(-1, 1), interval(0, 10)), interval(0, 1)); */ - analyzeBinaryMethod(10, 200, "Lt", interval(-1, 1, 0), interval(-1, 1, 0), myLt, &interval_algebra::Lt); - analyzeBinaryMethod(10, 200, "Lt", interval(-10, 10, 0), interval(-10, 10, 0), myLt, &interval_algebra::Lt); - analyzeBinaryMethod(10, 2000, "Lt", interval(-10, 10), interval(-10, 10), myLt, &interval_algebra::Lt); + analyzeBinaryMethod(10, 200, "Lt", interval(-1, 1, 0), interval(-1, 1, 0), myLt, + &interval_algebra::Lt); + analyzeBinaryMethod(10, 200, "Lt", interval(-10, 10, 0), interval(-10, 10, 0), myLt, + &interval_algebra::Lt); + analyzeBinaryMethod(10, 2000, "Lt", interval(-10, 10), interval(-10, 10), myLt, + &interval_algebra::Lt); } } // namespace itv diff --git a/interval/intervalMax.cpp b/interval/intervalMax.cpp index 359bb5a..03ed6b4 100644 --- a/interval/intervalMax.cpp +++ b/interval/intervalMax.cpp @@ -32,9 +32,9 @@ interval interval_algebra::Max(const interval& x, const interval& y) return empty(); } - return { - std::max(x.lo(), y.lo()), std::max(x.hi(), y.hi()), - std::min(x.lsb(), y.lsb())}; // the resulting interval should be as precise as the most precise of the operands + return {std::max(x.lo(), y.lo()), std::max(x.hi(), y.hi()), + std::min(x.lsb(), y.lsb())}; // the resulting interval should be as precise as the most + // precise of the operands } void interval_algebra::testMax() diff --git a/interval/intervalMin.cpp b/interval/intervalMin.cpp index 6c8247d..5b4cb1c 100644 --- a/interval/intervalMin.cpp +++ b/interval/intervalMin.cpp @@ -1,11 +1,11 @@ /* Copyright 2023 Yann ORLAREY - * + * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at - * + * * http://www.apache.org/licenses/LICENSE-2.0 - * + * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. @@ -28,9 +28,13 @@ namespace itv { interval interval_algebra::Min(const interval& x, const interval& y) { - if (x.isEmpty() || y.isEmpty()) return empty(); + if (x.isEmpty() || y.isEmpty()) { + return empty(); + } - return {std::min(x.lo(), y.lo()), std::min(x.hi(), y.hi()), std::min(x.lsb(), y.lsb())}; // resulting interval should be as precise as the most precise of the operands + return {std::min(x.lo(), y.lo()), std::min(x.hi(), y.hi()), + std::min(x.lsb(), y.lsb())}; // resulting interval should be as precise as the most + // precise of the operands } void interval_algebra::testMin() diff --git a/interval/intervalMissing.cpp b/interval/intervalMissing.cpp index 50725ce..d88780b 100644 --- a/interval/intervalMissing.cpp +++ b/interval/intervalMissing.cpp @@ -1,11 +1,10 @@ #include "interval_algebra.hh" #include "interval_def.hh" -namespace itv -{ +namespace itv { //------------------------------------------------------------------------------------------ -// Missing operations. A default implementation is provided for the code to compile. A real implementation has to be -// provided. +// Missing operations. A default implementation is provided for the code to compile. A real +// implementation has to be provided. interval interval_algebra::Nil() { @@ -68,7 +67,8 @@ interval interval_algebra::RDTbl(const interval& wtbl, const interval& ri) { return interval(0); } -interval interval_algebra::WRTbl(const interval& n, const interval& g, const interval& wi, const interval& ws) +interval interval_algebra::WRTbl(const interval& n, const interval& g, const interval& wi, + const interval& ws) { return interval(0); } @@ -89,7 +89,8 @@ interval interval_algebra::SoundFileLength(const interval& sf, const interval& x { return interval(0); } -interval interval_algebra::SoundFileBuffer(const interval& sf, const interval& x, const interval& y, const interval& z) +interval interval_algebra::SoundFileBuffer(const interval& sf, const interval& x, const interval& y, + const interval& z) { return interval(0); } @@ -98,18 +99,18 @@ interval interval_algebra::Waveform(const std::vector& w) return interval(0); } -interval interval_algebra::VSlider(const interval& name, const interval& init, const interval& lo, const interval& hi, - const interval& step) +interval interval_algebra::VSlider(const interval& name, const interval& init, const interval& lo, + const interval& hi, const interval& step) { return interval(0); } -interval interval_algebra::HSlider(const interval& name, const interval& init, const interval& lo, const interval& hi, - const interval& step) +interval interval_algebra::HSlider(const interval& name, const interval& init, const interval& lo, + const interval& hi, const interval& step) { return interval(0); } -interval interval_algebra::NumEntry(const interval& name, const interval& init, const interval& lo, const interval& hi, - const interval& step) +interval interval_algebra::NumEntry(const interval& name, const interval& init, const interval& lo, + const interval& hi, const interval& step) { return interval(0); } @@ -136,11 +137,13 @@ interval interval_algebra::ForeignFunction(const std::vector& ff) { return interval(0); } -interval interval_algebra::ForeignVar(const interval& type, const interval& name, const interval& file) +interval interval_algebra::ForeignVar(const interval& type, const interval& name, + const interval& file) { return interval(0); } -interval interval_algebra::ForeignConst(const interval& type, const interval& name, const interval& file) +interval interval_algebra::ForeignConst(const interval& type, const interval& name, + const interval& file) { return interval(0); } diff --git a/interval/intervalMod.cpp b/interval/intervalMod.cpp index ae0c726..684702d 100644 --- a/interval/intervalMod.cpp +++ b/interval/intervalMod.cpp @@ -76,7 +76,8 @@ std::pair splitnz(interval x) if (x.hi() < 0) { return {x, empty()}; } - return {interval{x.lo(), nexttoward(0.0, -1.0), x.lsb()}, interval{nexttoward(0.0, 1.0), x.hi(), x.lsb()}}; + return {interval{x.lo(), nexttoward(0.0, -1.0), x.lsb()}, + interval{nexttoward(0.0, 1.0), x.hi(), x.lsb()}}; } //------------------------------------------------------------------------------------------ @@ -97,7 +98,7 @@ interval positiveFMod(const interval& x, const interval& y) return empty(); } int n = int(x.lo() / y.hi()); - //std::cout << "n = " << n << std::endl; + // std::cout << "n = " << n << std::endl; int precision = std::min(x.lsb(), y.lsb()); // n == 0 obeys the same rules as the general case @@ -126,7 +127,7 @@ interval positiveFMod(const interval& x, const interval& y) } // prop : y.lo() > hi // in that case, the quotient between x and y is constant and equal - return interval{x.lo() - n*y.hi(), x.hi() - n*y.lo(), precision}; + return interval{x.lo() - n * y.hi(), x.hi() - n * y.lo(), precision}; } // fmod of two signed intervals @@ -142,10 +143,8 @@ interval interval_algebra::Mod(const interval& x, const interval& y) auto xpyp = positiveFMod(xp, yp); // Make sure these 4 values are in the resulting interval - auto bb = singleton(fmod(x.hi(), y.hi())) - + singleton(fmod(x.lo(), y.hi())) - + singleton(fmod(x.hi(), y.lo())) - + singleton(fmod(x.lo(), y.lo())); + auto bb = singleton(fmod(x.hi(), y.hi())) + singleton(fmod(x.lo(), y.hi())) + + singleton(fmod(x.hi(), y.lo())) + singleton(fmod(x.lo(), y.lo())); bb = interval{bb.lo(), bb.hi(), std::min(x.lsb(), y.lsb())}; @@ -155,21 +154,28 @@ interval interval_algebra::Mod(const interval& x, const interval& y) void interval_algebra::testMod() { - // check("test algebra Mod", Mod(interval(-100, 100), 1.0), interval(nextafter(-1.0, 0.0), nextafter(1.0, 0.0))); - // check("test algebra Mod", Mod(interval(0, 100), 2), interval(0, nextafter(2.0, 0))); - // check("test algebra Mod", Mod(interval(0, 100), -1.0), interval(0, nextafter(1.0, 0))); + // check("test algebra Mod", Mod(interval(-100, 100), 1.0), interval(nextafter(-1.0, 0.0), + // nextafter(1.0, 0.0))); check("test algebra Mod", Mod(interval(0, 100), 2), interval(0, + // nextafter(2.0, 0))); check("test algebra Mod", Mod(interval(0, 100), -1.0), interval(0, + // nextafter(1.0, 0))); /* check("test algebra Mod", Mod(interval(5, 7), interval(4, 4.5)), interval(0.5, 3)); check("test algebra Mod", Mod(interval(5, 7), interval(8, 10)), interval(5, 7)); check("test algebra Mod", Mod(interval(-7, 7), interval(8, 10)), interval(-7, 7)); - check("test algebra Mod", Mod(interval(0, 100), interval(7, 7)), interval(0, nextafter(7.0, 0.0)));*/ - - // analyzeBinaryMethod(10, 10000, "mod", interval(0, 10, -5), interval(0, 10, -5), fmod, &interval_algebra::Mod); - analyzeBinaryMethod(10, 10000, "mod", interval(0, 10, 1), interval(0, 10, 0), fmod, &interval_algebra::Mod); - /* analyzeBinaryMethod(10, 10000, "mod", interval(0, 10, 0), interval(0, 10, 0), fmod, &interval_algebra::Mod); - analyzeBinaryMethod(10, 10000, "mod", interval(0, 10, 0), interval(0, 10, -5), fmod, &interval_algebra::Mod); - - analyzeBinaryMethod(10, 100000000, "mod", interval(3, 4, -3), interval(1.2, 1.4, -3), fmod, &interval_algebra::Mod);*/ - // analyzeBinaryMethod(10, 10000, "mod", interval(-10, 10, -5), interval(-10, 10, -5), fmod, &interval_algebra::Mod); + check("test algebra Mod", Mod(interval(0, 100), interval(7, 7)), interval(0, nextafter(7.0, + 0.0)));*/ + + // analyzeBinaryMethod(10, 10000, "mod", interval(0, 10, -5), interval(0, 10, -5), fmod, + // &interval_algebra::Mod); + analyzeBinaryMethod(10, 10000, "mod", interval(0, 10, 1), interval(0, 10, 0), fmod, + &interval_algebra::Mod); + /* analyzeBinaryMethod(10, 10000, "mod", interval(0, 10, 0), interval(0, 10, 0), fmod, + &interval_algebra::Mod); analyzeBinaryMethod(10, 10000, "mod", interval(0, 10, 0), interval(0, + 10, -5), fmod, &interval_algebra::Mod); + + analyzeBinaryMethod(10, 100000000, "mod", interval(3, 4, -3), interval(1.2, 1.4, -3), fmod, + &interval_algebra::Mod);*/ + // analyzeBinaryMethod(10, 10000, "mod", interval(-10, 10, -5), interval(-10, 10, -5), fmod, + // &interval_algebra::Mod); } } // namespace itv diff --git a/interval/intervalMul.cpp b/interval/intervalMul.cpp index fff5017..c860b7d 100644 --- a/interval/intervalMul.cpp +++ b/interval/intervalMul.cpp @@ -53,12 +53,15 @@ interval interval_algebra::Mul(const interval& x, const interval& y) double lo = min4(a, b, c, d); double hi = max4(a, b, c, d); - if (x.lsb() >= 0 and y.lsb() >= 0) // operation between integers + if (x.lsb() >= 0 and y.lsb() >= 0) // operation between integers { - // if the quotient of an INT limit by an interval limit is below a limit of the other interval - // ie, if there is something big enough in the other interval to make the interval limit go beyond an INT limit - if (std::max(std::abs(x.lo()), std::abs(x.hi()))*std::max(std::abs(y.lo()), std::abs(y.hi())) >= (double) INT_MAX) { - return {(double) INT_MIN, (double) INT_MAX, x.lsb() + y.lsb()}; + // if the quotient of an INT limit by an interval limit is below a limit of the other + // interval ie, if there is something big enough in the other interval to make the interval + // limit go beyond an INT limit + if (std::max(std::abs(x.lo()), std::abs(x.hi())) * + std::max(std::abs(y.lo()), std::abs(y.hi())) >= + (double)INT_MAX) { + return {(double)INT_MIN, (double)INT_MAX, x.lsb() + y.lsb()}; } /* interval z{lo, hi, x.lsb()+y.lsb()}; interval shift{pow(2, 31), pow(2, 31), 31}; @@ -66,9 +69,8 @@ interval interval_algebra::Mul(const interval& x, const interval& y) return Sub(Mod(Add(z, shift), m), shift);*/ - /* if ((lo <= (double)INT_MIN - 1 and hi >= (double)INT_MIN) // discontinuity at the lower end - or - (lo <= (double)INT_MAX and hi >= (double)INT_MAX+1)) + /* if ((lo <= (double)INT_MIN - 1 and hi >= (double)INT_MIN) // discontinuity at the lower + end or (lo <= (double)INT_MAX and hi >= (double)INT_MAX+1)) { return {(double)INT_MIN, (double)INT_MAX, std::min(x.lsb(), y.lsb())}; } @@ -81,8 +83,10 @@ interval interval_algebra::Mul(const interval& x, const interval& y) return {min4(aint, bint, cint, dint), max4(aint, bint, cint, dint), x.lsb() + y.lsb()};*/ } - return {lo, hi, - x.lsb() + y.lsb()}; // the worst case, we need all the precision digits from both the operands + return { + lo, hi, + x.lsb() + + y.lsb()}; // the worst case, we need all the precision digits from both the operands } void interval_algebra::testMul() @@ -93,11 +97,13 @@ void interval_algebra::testMul() check("test algebra Mul", Mul(interval(0), interval(-HUGE_VAL, HUGE_VAL)), interval(0)); check("test algebra Mul", Mul(interval(-HUGE_VAL, HUGE_VAL), interval(0)), interval(0));*/ - // analyzeBinaryMethod(10, 2000, "mul", interval(0, 10, 0), interval(0, 10, 0), specialmult, &interval_algebra::Mul); - /* analyzeBinaryMethod(10, 20000, "mul", interval(0, 10, -5), interval(0, 10, 0), specialmult, &interval_algebra::Mul); - analyzeBinaryMethod(10, 20000, "mul", interval(0, 10, -10), interval(0, 10, 0), specialmult, &interval_algebra::Mul); - analyzeBinaryMethod(10, 20000, "mul", interval(0, 10, -15), interval(0, 10, 0), specialmult, &interval_algebra::Mul); - analyzeBinaryMethod(10, 20000, "mul", interval(0, 10, 0), interval(0, 10, -10), specialmult, &interval_algebra::Mul); + // analyzeBinaryMethod(10, 2000, "mul", interval(0, 10, 0), interval(0, 10, 0), specialmult, + // &interval_algebra::Mul); + /* analyzeBinaryMethod(10, 20000, "mul", interval(0, 10, -5), interval(0, 10, 0), specialmult, + &interval_algebra::Mul); analyzeBinaryMethod(10, 20000, "mul", interval(0, 10, -10), interval(0, + 10, 0), specialmult, &interval_algebra::Mul); analyzeBinaryMethod(10, 20000, "mul", interval(0, + 10, -15), interval(0, 10, 0), specialmult, &interval_algebra::Mul); analyzeBinaryMethod(10, + 20000, "mul", interval(0, 10, 0), interval(0, 10, -10), specialmult, &interval_algebra::Mul); analyzeBinaryMethod(10, 20000, "mul", interval(0, 10, -5), interval(0, 10, -10), specialmult, &interval_algebra::Mul); analyzeBinaryMethod(10, 20000, "mul", interval(0, 10, -10), interval(0, 10, -10), specialmult, @@ -105,11 +111,17 @@ void interval_algebra::testMul() analyzeBinaryMethod(10, 20000, "mul", interval(0, 10, -15), interval(0, 10, -10), specialmult, &interval_algebra::Mul); - analyzeBinaryMethod(10, 200000, "mul", interval(-pow(2, 31), pow(2, 31), 0), interval(-pow(2, 31), pow(2, 31), 0), specialmultint, &interval_algebra::Mul); - analyzeBinaryMethod(10, 200000, "mul", interval((double)INT_MAX/2, (double)INT_MAX, 0), interval((double)INT_MAX/2, (double)INT_MAX, 0), specialmultint, &interval_algebra::Mul); - analyzeBinaryMethod(10, 200000, "mul", interval((double)INT_MIN, (double)INT_MIN/2, 0), interval((double)INT_MAX/2, (double)INT_MAX, 0), specialmultint, &interval_algebra::Mul); - analyzeBinaryMethod(10, 2000000, "mul", interval((double)2*INT_MAX/3, (double)INT_MAX, 0), interval(0, 10, 0), specialmultint, &interval_algebra::Mul);*/ - check("Test integer Mult", Mul(interval(pow(2, 30), pow(2, 30)+2, 2), interval(1, 2, 0)), interval((double)INT_MIN, pow(2, 30)+2, 0)); - // analyzeBinaryMethod(10, 2000, "mul", interval((double)INT_MAX-1, (double)INT_MAX, 0), interval((double)INT_MAX-1, (double)INT_MAX, 0), specialmultint, &interval_algebra::Mul); + analyzeBinaryMethod(10, 200000, "mul", interval(-pow(2, 31), pow(2, 31), 0), interval(-pow(2, + 31), pow(2, 31), 0), specialmultint, &interval_algebra::Mul); analyzeBinaryMethod(10, 200000, + "mul", interval((double)INT_MAX/2, (double)INT_MAX, 0), interval((double)INT_MAX/2, + (double)INT_MAX, 0), specialmultint, &interval_algebra::Mul); analyzeBinaryMethod(10, 200000, + "mul", interval((double)INT_MIN, (double)INT_MIN/2, 0), interval((double)INT_MAX/2, + (double)INT_MAX, 0), specialmultint, &interval_algebra::Mul); analyzeBinaryMethod(10, 2000000, + "mul", interval((double)2*INT_MAX/3, (double)INT_MAX, 0), interval(0, 10, 0), specialmultint, + &interval_algebra::Mul);*/ + check("Test integer Mult", Mul(interval(pow(2, 30), pow(2, 30) + 2, 2), interval(1, 2, 0)), + interval((double)INT_MIN, pow(2, 30) + 2, 0)); + // analyzeBinaryMethod(10, 2000, "mul", interval((double)INT_MAX-1, (double)INT_MAX, 0), + // interval((double)INT_MAX-1, (double)INT_MAX, 0), specialmultint, &interval_algebra::Mul); } } // namespace itv diff --git a/interval/intervalNumEntry.cpp b/interval/intervalNumEntry.cpp index 94a2bce..85f5d97 100644 --- a/interval/intervalNumEntry.cpp +++ b/interval/intervalNumEntry.cpp @@ -19,17 +19,21 @@ namespace itv { //------------------------------------------------------------------------------------------ // Interval IntNum -interval interval_algebra::NumEntry(const interval& name, const interval& init, const interval& lo, const interval& hi, - const interval& step) +interval interval_algebra::NumEntry(const interval& name, const interval& init, const interval& lo, + const interval& hi, const interval& step) { - if (init.isEmpty() || lo.isEmpty() || hi.isEmpty() || step.isEmpty()) + if (init.isEmpty() || lo.isEmpty() || hi.isEmpty() || step.isEmpty()) { return empty(); + } - // elements of a slider with range [lo; hi] and step step are of the form lo + k·step <= hi with k an integer - // the precision needed to represent such elements is the minimum between - int lsb = std::min(step.lsb(), lo.lsb()); // the precision of the lower bound and that of the step - if (step.lo() > 0) { // if we don't have negative or zero steps - lsb = std::min(lsb, (int)log2(step.lo())); // and that associated to the smallest value the step can take + // elements of a slider with range [lo; hi] and step step are of the form lo + k·step <= hi with + // k an integer the precision needed to represent such elements is the minimum between + int lsb = + std::min(step.lsb(), lo.lsb()); // the precision of the lower bound and that of the step + if (step.lo() > 0) { // if we don't have negative or zero steps + lsb = std::min( + lsb, + (int)log2(step.lo())); // and that associated to the smallest value the step can take } return {lo.lo(), hi.hi(), lsb}; diff --git a/interval/intervalOr.cpp b/interval/intervalOr.cpp index 99608b8..0468a40 100644 --- a/interval/intervalOr.cpp +++ b/interval/intervalOr.cpp @@ -47,11 +47,12 @@ interval interval_algebra::Or(const interval& x, const interval& y) SInterval z = bitwiseSignedOr({x0, x1}, {y0, y1}); - int precision = std::min(x.lsb(), y.lsb()); // all input bits are significant + int precision = std::min(x.lsb(), y.lsb()); // all input bits are significant + + /* int precision = std::max(x.lsb(), y.lsb()); // output precision cannot be finer than that of + the input intervals - /* int precision = std::max(x.lsb(), y.lsb()); // output precision cannot be finer than that of the input intervals - // however, if one of the intervals is reduced to one element, the mask can make it so int precisionx = 0; @@ -74,18 +75,23 @@ interval interval_algebra::Or(const interval& x, const interval& y) precisiony++; } }*/ - return {double(z.lo), double(z.hi), + return {double(z.lo), double(z.hi), // std::max(precision, std::max(precisionx, precisiony)) precision}; } void interval_algebra::testOr() { - analyzeBinaryMethod(10, 20000, "Or", interval(-1000, -800), interval(127), myOr, &interval_algebra::Or); - analyzeBinaryMethod(10, 20000, "Or", interval(-1000, -800), interval(123), myOr, &interval_algebra::Or); - analyzeBinaryMethod(10, 20000, "Or", interval(-128, 128), interval(127), myOr, &interval_algebra::Or); - analyzeBinaryMethod(10, 20000, "Or", interval(0, 1000), interval(63, 127), myOr, &interval_algebra::Or); - analyzeBinaryMethod(10, 20000, "Or", interval(-1000, 1000), interval(63, 127), myOr, &interval_algebra::Or); + analyzeBinaryMethod(10, 20000, "Or", interval(-1000, -800), interval(127), myOr, + &interval_algebra::Or); + analyzeBinaryMethod(10, 20000, "Or", interval(-1000, -800), interval(123), myOr, + &interval_algebra::Or); + analyzeBinaryMethod(10, 20000, "Or", interval(-128, 128), interval(127), myOr, + &interval_algebra::Or); + analyzeBinaryMethod(10, 20000, "Or", interval(0, 1000), interval(63, 127), myOr, + &interval_algebra::Or); + analyzeBinaryMethod(10, 20000, "Or", interval(-1000, 1000), interval(63, 127), myOr, + &interval_algebra::Or); analyzeBinaryMethod(10, 2000, "Or", interval(10, 20), interval(0), myOr, &interval_algebra::Or); analyzeBinaryMethod(10, 2000, "Or", interval(0), interval(15, 25), myOr, &interval_algebra::Or); analyzeBinaryMethod(10, 2000, "Or", interval(0), interval(0), myOr, &interval_algebra::Or); diff --git a/interval/intervalPow.cpp b/interval/intervalPow.cpp index 8abc58c..844f009 100644 --- a/interval/intervalPow.cpp +++ b/interval/intervalPow.cpp @@ -36,10 +36,11 @@ static interval ipow(const interval& x, int k) { assert(k >= 0); if (k == 0) { - return interval{1,1,0}; + return interval{1, 1, 0}; } - // explicit expression because passing an anonymous function to exactPrecisionUnary is complicated + // explicit expression because passing an anonymous function to exactPrecisionUnary is + // complicated int precision = x.lsb() * k; // if x contains 0: finest precision is attained in 0 if (not x.hasZero()) { double v = minValAbs(x); @@ -47,8 +48,8 @@ static interval ipow(const interval& x, int k) int p1 = k * (int)log2(abs(v)); int p2 = 0; - double u = pow(2, x.lsb()); // ulp - double delta = abs(pow(1 + sign * u / v, k) - 1); + double u = pow(2, x.lsb()); // ulp + double delta = abs(pow(1 + sign * u / v, k) - 1); if (delta == 0) { // in case of u << x p2 = floor((double)log2(k) + x.lsb() - (double)log2(abs(v))); // (1 + u/v)^k - 1 ≃ k*u/v if u/v very small @@ -63,8 +64,12 @@ static interval ipow(const interval& x, int k) // k is even double z0 = std::pow(x.lo(), k); double z1 = std::pow(x.hi(), k); - return {x.hasZero() ? 0 : std::min(z0, z1), // 0 is in the output interval only if it is in the input interval - std::max(z0, z1), precision}; + return { + x.hasZero() + ? 0 + : std::min(z0, + z1), // 0 is in the output interval only if it is in the input interval + std::max(z0, z1), precision}; } // k is odd @@ -76,8 +81,9 @@ static interval ipow(const interval& x, int k) */ interval interval_algebra::fPow(const interval& x, const interval& y) { - if (x.isEmpty() || y.isEmpty()) + if (x.isEmpty() || y.isEmpty()) { return empty(); + } assert(x.lo() > 0); // x all positive @@ -86,8 +92,9 @@ interval interval_algebra::fPow(const interval& x, const interval& y) interval interval_algebra::iPow(const interval& x, const interval& y) { - if (x.isEmpty() || y.isEmpty()) + if (x.isEmpty() || y.isEmpty()) { return empty(); + } int y0 = std::max(0, saturatedIntCast(y.lo())); int y1 = std::max(0, saturatedIntCast(y.hi())); @@ -121,26 +128,36 @@ static double myiPow(double x, double y) void interval_algebra::testPow() { - /* analyzeBinaryMethod(10, 2000000, "iPow^2", interval(-100, 100), interval(2), myiPow, &interval_algebra::iPow); - analyzeBinaryMethod(10, 2000000, "iPow^3", interval(-100, 100), interval(3), myiPow, &interval_algebra::iPow); - analyzeBinaryMethod(10, 2000000, "iPow^2", interval(-1, 1), interval(2), myiPow, &interval_algebra::iPow); - analyzeBinaryMethod(10, 2000000, "iPow^3", interval(-1, 1), interval(3), myiPow, &interval_algebra::iPow);*/ - - analyzeBinaryMethod(5, 2000000, "iPow2", interval(-100, 100, 0), interval(0, 200, 0), myiPow, &interval_algebra::iPow); - analyzeBinaryMethod(5, 2000000, "iPow2", interval(-100, 100, -5), interval(0, 200, 0), myiPow, &interval_algebra::iPow); - analyzeBinaryMethod(5, 2000000, "iPow2", interval(-100, 100, 0), interval(0, 200, -5), myiPow, &interval_algebra::iPow); - - analyzeBinaryMethod(5, 2000000, "iPow2", interval(-1, 1, 0), interval(1, 3, 0), myiPow, &interval_algebra::iPow); - analyzeBinaryMethod(5, 2000000, "iPow2", interval(-1, 1, -5), interval(1, 3, 0), myiPow, &interval_algebra::iPow); - analyzeBinaryMethod(5, 2000000, "iPow2", interval(-1, 1, 0), interval(1, 3, -5), myiPow, &interval_algebra::iPow); - - analyzeBinaryMethod(5, 2000000, "iPow2", interval(-1, 1), interval(1, 10), myiPow, &interval_algebra::iPow); - analyzeBinaryMethod(5, 2000000, "iPow2", interval(-1, 1), interval(1, 10), myiPow, &interval_algebra::iPow); - analyzeBinaryMethod(5, 2000000, "iPow2", interval(-1, 1), interval(1, 10), myiPow, &interval_algebra::iPow); - - /* analyzeBinaryMethod(10, 2000000, "fPow2", interval(1, 1000), interval(-10, 10), myfPow, &interval_algebra::fPow); - analyzeBinaryMethod(10, 2000000, "fPow2", interval(0.001, 1), interval(-10, 10), myfPow, &interval_algebra::fPow); - analyzeBinaryMethod(10, 2000000, "fPow2", interval(0.001, 10), interval(-200, 200), myfPow, - &interval_algebra::fPow);*/ + /* analyzeBinaryMethod(10, 2000000, "iPow^2", interval(-100, 100), interval(2), myiPow, + &interval_algebra::iPow); analyzeBinaryMethod(10, 2000000, "iPow^3", interval(-100, 100), + interval(3), myiPow, &interval_algebra::iPow); analyzeBinaryMethod(10, 2000000, "iPow^2", + interval(-1, 1), interval(2), myiPow, &interval_algebra::iPow); analyzeBinaryMethod(10, 2000000, + "iPow^3", interval(-1, 1), interval(3), myiPow, &interval_algebra::iPow);*/ + + analyzeBinaryMethod(5, 2000000, "iPow2", interval(-100, 100, 0), interval(0, 200, 0), myiPow, + &interval_algebra::iPow); + analyzeBinaryMethod(5, 2000000, "iPow2", interval(-100, 100, -5), interval(0, 200, 0), myiPow, + &interval_algebra::iPow); + analyzeBinaryMethod(5, 2000000, "iPow2", interval(-100, 100, 0), interval(0, 200, -5), myiPow, + &interval_algebra::iPow); + + analyzeBinaryMethod(5, 2000000, "iPow2", interval(-1, 1, 0), interval(1, 3, 0), myiPow, + &interval_algebra::iPow); + analyzeBinaryMethod(5, 2000000, "iPow2", interval(-1, 1, -5), interval(1, 3, 0), myiPow, + &interval_algebra::iPow); + analyzeBinaryMethod(5, 2000000, "iPow2", interval(-1, 1, 0), interval(1, 3, -5), myiPow, + &interval_algebra::iPow); + + analyzeBinaryMethod(5, 2000000, "iPow2", interval(-1, 1), interval(1, 10), myiPow, + &interval_algebra::iPow); + analyzeBinaryMethod(5, 2000000, "iPow2", interval(-1, 1), interval(1, 10), myiPow, + &interval_algebra::iPow); + analyzeBinaryMethod(5, 2000000, "iPow2", interval(-1, 1), interval(1, 10), myiPow, + &interval_algebra::iPow); + + /* analyzeBinaryMethod(10, 2000000, "fPow2", interval(1, 1000), interval(-10, 10), myfPow, + &interval_algebra::fPow); analyzeBinaryMethod(10, 2000000, "fPow2", interval(0.001, 1), + interval(-10, 10), myfPow, &interval_algebra::fPow); analyzeBinaryMethod(10, 2000000, "fPow2", + interval(0.001, 10), interval(-200, 200), myfPow, &interval_algebra::fPow);*/ } } // namespace itv diff --git a/interval/intervalRint.cpp b/interval/intervalRint.cpp index 4178486..46a8cbe 100644 --- a/interval/intervalRint.cpp +++ b/interval/intervalRint.cpp @@ -28,10 +28,12 @@ namespace itv { interval interval_algebra::Rint(const interval& x) { - if (x.isEmpty()) + if (x.isEmpty()) { return empty(); + } - return {rint(x.lo()), rint(x.hi()), std::max(0, x.lsb())}; // round to nearest integral value => integer => precision 0 + return {rint(x.lo()), rint(x.hi()), + std::max(0, x.lsb())}; // round to nearest integral value => integer => precision 0 } void interval_algebra::testRint() diff --git a/interval/intervalRound.cpp b/interval/intervalRound.cpp index 609b91a..ff3456a 100644 --- a/interval/intervalRound.cpp +++ b/interval/intervalRound.cpp @@ -28,10 +28,13 @@ namespace itv { interval interval_algebra::Round(const interval& x) { - if (x.isEmpty()) + if (x.isEmpty()) { return empty(); + } - return {round(x.lo()), round(x.hi()), std::max(0, x.lsb())}; // round to integral value (regardless of rounding direction) => integer => precision 0 + return {round(x.lo()), round(x.hi()), + std::max(0, x.lsb())}; // round to integral value (regardless of rounding direction) => + // integer => precision 0 } void interval_algebra::testRound() diff --git a/interval/intervalRsh.cpp b/interval/intervalRsh.cpp index 9238690..b7bee73 100644 --- a/interval/intervalRsh.cpp +++ b/interval/intervalRsh.cpp @@ -33,23 +33,32 @@ static double rsh(double x, double k) interval interval_algebra::Rsh(const interval& x, const interval& k) { - if (x.isEmpty() || k.isEmpty()) + if (x.isEmpty() || k.isEmpty()) { return empty(); + } interval j{pow(2, -k.hi()), pow(2, -k.lo())}; interval z = Mul(x, j); - return {z.lo(), z.hi(), x.lsb() - (int)k.hi()}; // rshifts add some precision to the numbers, at most y.hi() bits + return { + z.lo(), z.hi(), + x.lsb() - (int)k.hi()}; // rshifts add some precision to the numbers, at most y.hi() bits } void interval_algebra::testRsh() { // check("test algebra Rsh", Rsh(interval(8, 16), interval(4)), interval(0.5, 1)); - analyzeBinaryMethod(10, 1000, "rshift", interval(0, 32, 0), interval(8, 8, 1), rsh, &interval_algebra::Rsh); - analyzeBinaryMethod(10, 1000, "rshift", interval(0, 1024, 0), interval(-10, 10, 0), rsh, &interval_algebra::Rsh); - analyzeBinaryMethod(10, 1000, "rshift", interval(0, 1024, 2), interval(-10, 10, 0), rsh, &interval_algebra::Rsh); - analyzeBinaryMethod(10, 1000, "rshift", interval(0, 1024, 0), interval(-10, 10, 1), rsh, &interval_algebra::Rsh); - analyzeBinaryMethod(10, 1000, "rshift", interval(0, 1024, 2), interval(-10, 10, 1), rsh, &interval_algebra::Rsh); - // analyzeBinaryMethod(10, 1000, "rshift", interval(0, 32, 0), interval(-3, 0, 0), rsh, &interval_algebra::Rsh); + analyzeBinaryMethod(10, 1000, "rshift", interval(0, 32, 0), interval(8, 8, 1), rsh, + &interval_algebra::Rsh); + analyzeBinaryMethod(10, 1000, "rshift", interval(0, 1024, 0), interval(-10, 10, 0), rsh, + &interval_algebra::Rsh); + analyzeBinaryMethod(10, 1000, "rshift", interval(0, 1024, 2), interval(-10, 10, 0), rsh, + &interval_algebra::Rsh); + analyzeBinaryMethod(10, 1000, "rshift", interval(0, 1024, 0), interval(-10, 10, 1), rsh, + &interval_algebra::Rsh); + analyzeBinaryMethod(10, 1000, "rshift", interval(0, 1024, 2), interval(-10, 10, 1), rsh, + &interval_algebra::Rsh); + // analyzeBinaryMethod(10, 1000, "rshift", interval(0, 32, 0), interval(-3, 0, 0), rsh, + // &interval_algebra::Rsh); } } // namespace itv diff --git a/interval/intervalSin.cpp b/interval/intervalSin.cpp index 38661f1..e4e1cd8 100644 --- a/interval/intervalSin.cpp +++ b/interval/intervalSin.cpp @@ -34,20 +34,24 @@ static double sinPi(double x) interval interval_algebra::Sin(const interval& x) { - if (x.isEmpty()) + if (x.isEmpty()) { return empty(); + } int precision = exactPrecisionUnary(sin, 0.5, pow(2, x.lsb())); - if (precision == INT_MIN or taylor_lsb) precision = 2*x.lsb() - 1; // if x.lsb() is so small that the automatic computation doesn't work + if (precision == INT_MIN or taylor_lsb) { + precision = + 2 * x.lsb() - 1; // if x.lsb() is so small that the automatic computation doesn't work + } - if (x.size() >= 2*M_PI) { + if (x.size() >= 2 * M_PI) { return {-1, 1, precision}; } // normalize input interval between 0..2 - double l = fmod(x.lo(), 2*M_PI); + double l = fmod(x.lo(), 2 * M_PI); if (l < 0) { - l += 2*M_PI; + l += 2 * M_PI; } interval i(l, l + x.size(), x.lsb()); @@ -58,23 +62,25 @@ interval interval_algebra::Sin(const interval& x) double hi = std::max(a, b); // check if integers are included - if (i.has(M_PI_2) || i.has(5*M_PI_2)) { + if (i.has(M_PI_2) || i.has(5 * M_PI_2)) { hi = 1; } - if (i.has(3*M_PI_2) || i.has(7*M_PI_2)) { + if (i.has(3 * M_PI_2) || i.has(7 * M_PI_2)) { lo = -1; } double v = M_PI_2; // value of the interval at which the finest precision is computed - // defaults at 0.5, interchangeable with any other half-integer + // defaults at 0.5, interchangeable with any other half-integer // precision if we don't hit the half integers if (i.hi() < M_PI_2) { v = x.hi(); - } else if ((i.lo() > M_PI_2 and i.hi() < 3*M_PI_2) or (i.lo() > 3*M_PI_2 and i.hi() < 2.5*M_PI)) { - double delta_hi = ceil(i.hi()/M_PI + 0.5) - i.hi()/M_PI; - double delta_lo = i.lo()/M_PI - floor(i.lo()/M_PI - 0.5); - if (delta_lo > delta_hi) { // if i.hi is closer to its higher half-integer than i.lo() to its lower half-integer + } else if ((i.lo() > M_PI_2 and i.hi() < 3 * M_PI_2) or + (i.lo() > 3 * M_PI_2 and i.hi() < 2.5 * M_PI)) { + double delta_hi = ceil(i.hi() / M_PI + 0.5) - i.hi() / M_PI; + double delta_lo = i.lo() / M_PI - floor(i.lo() / M_PI - 0.5); + if (delta_lo > delta_hi) { // if i.hi is closer to its higher half-integer than i.lo() to + // its lower half-integer v = x.hi(); } else { v = x.lo(); @@ -82,10 +88,14 @@ interval interval_algebra::Sin(const interval& x) } precision = exactPrecisionUnary(sin, v, pow(2, x.lsb())); - if (precision == INT_MIN or taylor_lsb) - { - if (v != 0.5*M_PI) precision = x.lsb() + (int)floor(log2(abs(cos(v))));// (int)floor(log2(M_PI*cos(M_PI*v))) + x.lsb(); - else precision = 2*x.lsb() -1; // - (int)floor(2*log2(M_PI)); + if (precision == INT_MIN or taylor_lsb) { + if (v != 0.5 * M_PI) { + precision = + x.lsb() + + (int)floor(log2(abs(cos(v)))); // (int)floor(log2(M_PI*cos(M_PI*v))) + x.lsb(); + } else { + precision = 2 * x.lsb() - 1; // - (int)floor(2*log2(M_PI)); + } } return {lo, hi, precision}; @@ -94,11 +104,11 @@ interval interval_algebra::Sin(const interval& x) void interval_algebra::testSin() { // analyzeUnaryMethod(5, 20000, "sin", interval(-1, 1, -3), sin, &interval_algebra::Sin); - analyzeUnaryMethod(10, 40000, "sin", interval(0, 2*M_PI, -3), sin, &interval_algebra::Sin); - analyzeUnaryMethod(10, 40000, "sin", interval(0, 2*M_PI, -5), sin, &interval_algebra::Sin); - analyzeUnaryMethod(10, 40000, "sin", interval(0, 2*M_PI, -10), sin, &interval_algebra::Sin); - analyzeUnaryMethod(10, 40000, "sin", interval(0, 2*M_PI, -15), sin, &interval_algebra::Sin); - analyzeUnaryMethod(10, 40000, "sin", interval(0, 2*M_PI, -20), sin, &interval_algebra::Sin); - analyzeUnaryMethod(10, 40000, "sin", interval(0, 2*M_PI, -24), sin, &interval_algebra::Sin); + analyzeUnaryMethod(10, 40000, "sin", interval(0, 2 * M_PI, -3), sin, &interval_algebra::Sin); + analyzeUnaryMethod(10, 40000, "sin", interval(0, 2 * M_PI, -5), sin, &interval_algebra::Sin); + analyzeUnaryMethod(10, 40000, "sin", interval(0, 2 * M_PI, -10), sin, &interval_algebra::Sin); + analyzeUnaryMethod(10, 40000, "sin", interval(0, 2 * M_PI, -15), sin, &interval_algebra::Sin); + analyzeUnaryMethod(10, 40000, "sin", interval(0, 2 * M_PI, -20), sin, &interval_algebra::Sin); + analyzeUnaryMethod(10, 40000, "sin", interval(0, 2 * M_PI, -24), sin, &interval_algebra::Sin); } } // namespace itv diff --git a/interval/intervalSinh.cpp b/interval/intervalSinh.cpp index 603d83c..e5d517a 100644 --- a/interval/intervalSinh.cpp +++ b/interval/intervalSinh.cpp @@ -42,8 +42,9 @@ interval interval_algebra::Sinh(const interval& x) } int precision = exactPrecisionUnary(sinh, v, sign * pow(2, x.lsb())); - if (precision == INT_MIN or taylor_lsb) + if (precision == INT_MIN or taylor_lsb) { precision = floor(x.lsb() + log2(cosh(v))); + } return {sinh(x.lo()), sinh(x.hi()), precision}; } diff --git a/interval/intervalSqrt.cpp b/interval/intervalSqrt.cpp index 3b412c6..d41fad6 100644 --- a/interval/intervalSqrt.cpp +++ b/interval/intervalSqrt.cpp @@ -41,12 +41,12 @@ interval interval_algebra::Sqrt(const interval& x) // lowest slope at the highest bound of the interval int precision = exactPrecisionUnary(sqrt, i.hi(), -pow(2, i.lsb())); - if (precision == INT_MIN or taylor_lsb) - { - if (i.hi() == 0) - precision = floor(i.lsb()/2); - else - precision = floor(i.lsb() - log2(i.hi()) - 1); // sqrt(x+u) - sqrt(x) = 1/2 u/sqrt(x) + if (precision == INT_MIN or taylor_lsb) { + if (i.hi() == 0) { + precision = floor(i.lsb() / 2); + } else { + precision = floor(i.lsb() - log2(i.hi()) - 1); // sqrt(x+u) - sqrt(x) = 1/2 u/sqrt(x) + } } return {sqrt(i.lo()), sqrt(i.hi()), precision}; diff --git a/interval/intervalSub.cpp b/interval/intervalSub.cpp index ac5dcc5..382a141 100644 --- a/interval/intervalSub.cpp +++ b/interval/intervalSub.cpp @@ -41,13 +41,21 @@ interval interval_algebra::Sub(const interval& x, const interval& y) void interval_algebra::testSub() { // check("test algebra Sub", Sub(interval(0, 100), interval(10, 500)), interval(-500, 90)); - analyzeBinaryMethod(10, 2000, "sub", interval(0, 10, 0), interval(0, 10, 0), sub, &interval_algebra::Sub); - analyzeBinaryMethod(10, 2000, "sub", interval(0, 10, -5), interval(0, 10, 0), sub, &interval_algebra::Sub); - analyzeBinaryMethod(10, 2000, "sub", interval(0, 10, -10), interval(0, 10, 0), sub, &interval_algebra::Sub); - analyzeBinaryMethod(10, 2000, "sub", interval(0, 10, -15), interval(0, 10, 0), sub, &interval_algebra::Sub); - analyzeBinaryMethod(10, 2000, "sub", interval(0, 10, 0), interval(0, 10, -10), sub, &interval_algebra::Sub); - analyzeBinaryMethod(10, 2000, "sub", interval(0, 10, -5), interval(0, 10, -10), sub, &interval_algebra::Sub); - analyzeBinaryMethod(10, 2000, "sub", interval(0, 10, -10), interval(0, 10, -10), sub, &interval_algebra::Sub); - analyzeBinaryMethod(10, 2000, "sub", interval(0, 10, -15), interval(0, 10, -10), sub, &interval_algebra::Sub); + analyzeBinaryMethod(10, 2000, "sub", interval(0, 10, 0), interval(0, 10, 0), sub, + &interval_algebra::Sub); + analyzeBinaryMethod(10, 2000, "sub", interval(0, 10, -5), interval(0, 10, 0), sub, + &interval_algebra::Sub); + analyzeBinaryMethod(10, 2000, "sub", interval(0, 10, -10), interval(0, 10, 0), sub, + &interval_algebra::Sub); + analyzeBinaryMethod(10, 2000, "sub", interval(0, 10, -15), interval(0, 10, 0), sub, + &interval_algebra::Sub); + analyzeBinaryMethod(10, 2000, "sub", interval(0, 10, 0), interval(0, 10, -10), sub, + &interval_algebra::Sub); + analyzeBinaryMethod(10, 2000, "sub", interval(0, 10, -5), interval(0, 10, -10), sub, + &interval_algebra::Sub); + analyzeBinaryMethod(10, 2000, "sub", interval(0, 10, -10), interval(0, 10, -10), sub, + &interval_algebra::Sub); + analyzeBinaryMethod(10, 2000, "sub", interval(0, 10, -15), interval(0, 10, -10), sub, + &interval_algebra::Sub); } } // namespace itv diff --git a/interval/intervalTan.cpp b/interval/intervalTan.cpp index b67042e..892494b 100644 --- a/interval/intervalTan.cpp +++ b/interval/intervalTan.cpp @@ -40,12 +40,10 @@ interval interval_algebra::Tan(const interval& x) if (x.size() >= M_PI) { // spans asymptots and contains an integer multiple of pi // std::cout << "Spanning more than one period" << std::endl; int precision = exactPrecisionUnary(tan, 0, pow(2, x.lsb())); - if (precision == INT_MIN or taylor_lsb) + if (precision == INT_MIN or taylor_lsb) { precision = x.lsb(); - return { - -HUGE_VAL, HUGE_VAL, - precision - }; + } + return {-HUGE_VAL, HUGE_VAL, precision}; } // normalize input interval between -PI/2..PI/2 @@ -58,14 +56,14 @@ interval interval_algebra::Tan(const interval& x) if (i.lo() > 0) { v = i.lo(); } else if (i.hi() < 0) { - v = i.hi(); + v = i.hi(); sign = -1; } int precision = exactPrecisionUnary(tan, v, sign * pow(2, x.lsb())); - if (precision == INT_MIN or taylor_lsb) - precision = floor(x.lsb() - 2*(double)log2(cos(v))); - + if (precision == INT_MIN or taylor_lsb) { + precision = floor(x.lsb() - 2 * (double)log2(cos(v))); + } if (i.has(-M_PI_2) or i.has(M_PI_2)) { // asymptots at PI/2 return {-HUGE_VAL, HUGE_VAL, precision}; diff --git a/interval/intervalTanh.cpp b/interval/intervalTanh.cpp index b7067a3..36b944b 100644 --- a/interval/intervalTanh.cpp +++ b/interval/intervalTanh.cpp @@ -33,13 +33,15 @@ interval interval_algebra::Tanh(const interval& x) return empty(); } - // value at which the lowest slope is attained: bound of the interval with the highest absolute value + // value at which the lowest slope is attained: bound of the interval with the highest absolute + // value double v = maxValAbs(x); int sign = signMaxValAbs(x); int precision = exactPrecisionUnary(tanh, v, sign * pow(2, x.lsb())); - if (precision == INT_MIN or taylor_lsb) - precision = floor(x.lsb() - 2*(double)log2(cosh(v))); + if (precision == INT_MIN or taylor_lsb) { + precision = floor(x.lsb() - 2 * (double)log2(cosh(v))); + } return {tanh(x.lo()), tanh(x.hi()), precision}; } diff --git a/interval/intervalVSlider.cpp b/interval/intervalVSlider.cpp index 1362532..339ab67 100644 --- a/interval/intervalVSlider.cpp +++ b/interval/intervalVSlider.cpp @@ -20,17 +20,21 @@ namespace itv { //------------------------------------------------------------------------------------------ // Interval IntNum -interval interval_algebra::VSlider(const interval& name, const interval& init, const interval& lo, const interval& hi, - const interval& step) +interval interval_algebra::VSlider(const interval& name, const interval& init, const interval& lo, + const interval& hi, const interval& step) { - if (init.isEmpty() || lo.isEmpty() || hi.isEmpty() || step.isEmpty()) + if (init.isEmpty() || lo.isEmpty() || hi.isEmpty() || step.isEmpty()) { return empty(); + } - // elements of a slider with range [lo; hi] and step step are of the form lo + k·step <= hi with k an integer - // the precision needed to represent such elements is the minimum between - int lsb = std::min(step.lsb(), lo.lsb()); // the precision of the lower bound and that of the step - if (step.lo() > 0) { // if we don't have negative or zero steps - lsb = std::min(lsb, (int)log2(step.lo())); // and that associated to the smallest value the step can take + // elements of a slider with range [lo; hi] and step step are of the form lo + k·step <= hi with + // k an integer the precision needed to represent such elements is the minimum between + int lsb = + std::min(step.lsb(), lo.lsb()); // the precision of the lower bound and that of the step + if (step.lo() > 0) { // if we don't have negative or zero steps + lsb = std::min( + lsb, + (int)log2(step.lo())); // and that associated to the smallest value the step can take } return {lo.lo(), hi.hi(), lsb}; // TODO: step, init diff --git a/interval/intervalXor.cpp b/interval/intervalXor.cpp index 400d292..0de5d46 100644 --- a/interval/intervalXor.cpp +++ b/interval/intervalXor.cpp @@ -47,9 +47,11 @@ interval interval_algebra::Xor(const interval& x, const interval& y) SInterval z = bitwiseSignedXOr({x0, x1}, {y0, y1}); - int precision = std::min(x.lsb(), y.lsb()); // output precision cannot be finer than that of the input intervals + int precision = std::min( + x.lsb(), y.lsb()); // output precision cannot be finer than that of the input intervals - // if both intervals are singletons, the lsb is the least significant bit of the only element of the interval + // if both intervals are singletons, the lsb is the least significant bit of the only element of + // the interval if (x0 == x1 and y0 == y1) { int v = x0 ^ y0; precision = 0; @@ -60,8 +62,8 @@ interval interval_algebra::Xor(const interval& x, const interval& y) } } - // if only one of the intervals is a singleton, all of the variation is due to the other interval, which transmits - // its lsb + // if only one of the intervals is a singleton, all of the variation is due to the other + // interval, which transmits its lsb if (x0 == x1) { precision = y.lsb(); } @@ -80,33 +82,35 @@ void interval_algebra::testXor() std::uniform_int_distribution lx(0, 10); std::uniform_int_distribution ly(0, 10); - analyzeBinaryMethod(10, 20000, "Xor", interval(-1000, -800, lx(generator)), interval(127, 127, ly(generator)), - myXor, &interval_algebra::Xor); - analyzeBinaryMethod(10, 20000, "Xor", interval(-1000, -800, lx(generator)), interval(127, 127, ly(generator)), - myXor, &interval_algebra::Xor); + analyzeBinaryMethod(10, 20000, "Xor", interval(-1000, -800, lx(generator)), + interval(127, 127, ly(generator)), myXor, &interval_algebra::Xor); + analyzeBinaryMethod(10, 20000, "Xor", interval(-1000, -800, lx(generator)), + interval(127, 127, ly(generator)), myXor, &interval_algebra::Xor); - analyzeBinaryMethod(10, 20000, "Xor", interval(-1000, -800, lx(generator)), interval(123, 123, ly(generator)), - myXor, &interval_algebra::Xor); - analyzeBinaryMethod(10, 20000, "Xor", interval(-1000, -800, lx(generator)), interval(123, 123, ly(generator)), - myXor, &interval_algebra::Xor); + analyzeBinaryMethod(10, 20000, "Xor", interval(-1000, -800, lx(generator)), + interval(123, 123, ly(generator)), myXor, &interval_algebra::Xor); + analyzeBinaryMethod(10, 20000, "Xor", interval(-1000, -800, lx(generator)), + interval(123, 123, ly(generator)), myXor, &interval_algebra::Xor); - analyzeBinaryMethod(10, 20000, "Xor", interval(-128, 128, lx(generator)), interval(127, 127, ly(generator)), myXor, - &interval_algebra::Xor); - analyzeBinaryMethod(10, 20000, "Xor", interval(-128, 128, lx(generator)), interval(127, 127, ly(generator)), myXor, - &interval_algebra::Xor); + analyzeBinaryMethod(10, 20000, "Xor", interval(-128, 128, lx(generator)), + interval(127, 127, ly(generator)), myXor, &interval_algebra::Xor); + analyzeBinaryMethod(10, 20000, "Xor", interval(-128, 128, lx(generator)), + interval(127, 127, ly(generator)), myXor, &interval_algebra::Xor); - analyzeBinaryMethod(10, 20000, "Xor", interval(0, 1000, lx(generator)), interval(63, 127, ly(generator)), myXor, - &interval_algebra::Xor); - analyzeBinaryMethod(10, 20000, "Xor", interval(0, 1000, lx(generator)), interval(63, 127, ly(generator)), myXor, - &interval_algebra::Xor); + analyzeBinaryMethod(10, 20000, "Xor", interval(0, 1000, lx(generator)), + interval(63, 127, ly(generator)), myXor, &interval_algebra::Xor); + analyzeBinaryMethod(10, 20000, "Xor", interval(0, 1000, lx(generator)), + interval(63, 127, ly(generator)), myXor, &interval_algebra::Xor); + + analyzeBinaryMethod(10, 20000, "Xor", interval(-1000, 1000, lx(generator)), + interval(63, 127, ly(generator)), myXor, &interval_algebra::Xor); + analyzeBinaryMethod(10, 20000, "Xor", interval(-1000, 1000, lx(generator)), + interval(63, 127, ly(generator)), myXor, &interval_algebra::Xor); - analyzeBinaryMethod(10, 20000, "Xor", interval(-1000, 1000, lx(generator)), interval(63, 127, ly(generator)), myXor, + analyzeBinaryMethod(10, 2000, "Xor", interval(10, 20), interval(0), myXor, &interval_algebra::Xor); - analyzeBinaryMethod(10, 20000, "Xor", interval(-1000, 1000, lx(generator)), interval(63, 127, ly(generator)), myXor, + analyzeBinaryMethod(10, 2000, "Xor", interval(0), interval(15, 25), myXor, &interval_algebra::Xor); - - analyzeBinaryMethod(10, 2000, "Xor", interval(10, 20), interval(0), myXor, &interval_algebra::Xor); - analyzeBinaryMethod(10, 2000, "Xor", interval(0), interval(15, 25), myXor, &interval_algebra::Xor); analyzeBinaryMethod(10, 2000, "Xor", interval(0), interval(0), myXor, &interval_algebra::Xor); } } // namespace itv diff --git a/interval/interval_algebra.hh b/interval/interval_algebra.hh index 8c56b74..e4038d9 100644 --- a/interval/interval_algebra.hh +++ b/interval/interval_algebra.hh @@ -4,10 +4,8 @@ #include "FaustAlgebra.hh" -namespace itv -{ -class interval_algebra : public FaustAlgebra -{ +namespace itv { +class interval_algebra : public FaustAlgebra { private: interval iPow(const interval& x, const interval& y); // integer power, when x can be negative interval fPow(const interval& x, const interval& y); // float power, when x is positive @@ -34,27 +32,30 @@ class interval_algebra : public FaustAlgebra interval Select2(const interval& x, const interval& y, const interval& z) override; interval Prefix(const interval& x, const interval& y) override; interval RDTbl(const interval& wtbl, const interval& ri) override; - interval WRTbl(const interval& n, const interval& g, const interval& wi, const interval& ws) override; + interval WRTbl(const interval& n, const interval& g, const interval& wi, + const interval& ws) override; interval SoundFile(const interval& label) override; interval SoundFileRate(const interval& sf, const interval& x) override; interval SoundFileLength(const interval& sf, const interval& x) override; - interval SoundFileBuffer(const interval& sf, const interval& x, const interval& y, const interval& z) override; + interval SoundFileBuffer(const interval& sf, const interval& x, const interval& y, + const interval& z) override; interval Waveform(const std::vector& w) override; // Foreign functions interval ForeignFunction(const std::vector& ff) override; interval ForeignVar(const interval& type, const interval& name, const interval& file) override; - interval ForeignConst(const interval& type, const interval& name, const interval& file) override; + interval ForeignConst(const interval& type, const interval& name, + const interval& file) override; // User interface elements interval Button(const interval& name) override; interval Checkbox(const interval& name) override; - interval VSlider(const interval& name, const interval& init, const interval& lo, const interval& hi, - const interval& step) override; - interval HSlider(const interval& name, const interval& init, const interval& lo, const interval& hi, - const interval& step) override; - interval NumEntry(const interval& name, const interval& init, const interval& lo, const interval& hi, - const interval& step) override; + interval VSlider(const interval& name, const interval& init, const interval& lo, + const interval& hi, const interval& step) override; + interval HSlider(const interval& name, const interval& init, const interval& lo, + const interval& hi, const interval& step) override; + interval NumEntry(const interval& name, const interval& init, const interval& lo, + const interval& hi, const interval& step) override; interval Abs(const interval& x) override; void testAbs(); diff --git a/interval/interval_def.hh b/interval/interval_def.hh index 4ee4d7b..46a5cd6 100644 --- a/interval/interval_def.hh +++ b/interval/interval_def.hh @@ -60,10 +60,11 @@ class interval { interval(double n, double m, int lsb = -24) noexcept { - if (lsb == INT_MIN) + if (lsb == INT_MIN) { fLSB = -24; - else + } else { fLSB = lsb; + } if (std::isnan(n) || std::isnan(m)) { fLo = NAN; @@ -113,18 +114,20 @@ class interval { // position of the most significant bit of the value, without taking the sign bit into account int msb() const { - if (fLo == 0 and fHi == 0) + if (fLo == 0 and fHi == 0) { return 0; + } // amplitude of the interval - // can be < 1.0, in which case the msb will be negative and indicate the number of implicit leading zeroes + // can be < 1.0, in which case the msb will be negative and indicate the number of implicit + // leading zeroes double range = std::max(std::abs(fLo), std::abs(fHi)); - if (std::isinf(range)) { // if (fLSB == 0) // if we're dealing with integers: is that a good criterion? - return 31; - // return 20; // max MSB of the VHDL design; TODO: change when integrating in the compiler + return 31; + // return 20; // max MSB of the VHDL design; TODO: change when integrating in the + // compiler } int l = int(std::ceil(std::log2(range))); @@ -169,7 +172,8 @@ inline interval intersection(const interval& i, const interval& j) } else { double l = std::max(i.lo(), j.lo()); double h = std::min(i.hi(), j.hi()); - int p = std::min(i.lsb(), j.lsb()); // precision of the intersection should be the finest of the two + int p = std::min(i.lsb(), + j.lsb()); // precision of the intersection should be the finest of the two if (l > h) { return {}; } else { @@ -187,7 +191,8 @@ inline interval reunion(const interval& i, const interval& j) } else { double l = std::min(i.lo(), j.lo()); double h = std::max(i.hi(), j.hi()); - int p = std::min(i.lsb(), j.lsb()); // precision of the reunion should be the finest of the two + int p = + std::min(i.lsb(), j.lsb()); // precision of the reunion should be the finest of the two return {l, h, p}; } } @@ -206,7 +211,7 @@ inline interval singleton(double x) int m = std::floor(std::log2(std::abs(x))); - int precision = m - 32; // 32 = set width + int precision = m - 32; // 32 = set width return {x, x, precision}; } diff --git a/interval/precision_utils.hh b/interval/precision_utils.hh index 867b6ce..503692d 100644 --- a/interval/precision_utils.hh +++ b/interval/precision_utils.hh @@ -3,58 +3,60 @@ #include #include "check.hh" -/** +/** * @brief truncate x at the precision induced by lsb - * + * * @param x value to truncate * @param lsb the precision to which to truncate * @return x truncated with lsb bits of precision * */ double truncate(double x, int lsb) { - double u = pow(2, lsb); // ulp - double res = u*(double)floor(x/u); + double u = pow(2, lsb); // ulp + double res = u * (double)floor(x / u); return res; } -/** +/** * @brief truncate x at the precision induced by lsb (integer version) - * + * * @param x value to truncate * @param lsb the precision to which to truncate * @return x truncated with lsb bits of precision * */ int truncate(int x, int lsb) { - double u = pow(2, lsb); // ulp - double res = u*(double)floor(x/u); + double u = pow(2, lsb); // ulp + double res = u * (double)floor(x / u); return (int)res; } /** * @brief Computes the position of the least significant bit of a number * Floored to -24 - * + * * @param x The number * @return The lsb -*/ + */ int lsb_number(double x) { int precision = -24; - while (floor(x*pow(2, -precision-1)) == x*pow(2, -precision-1) and x != 0) + while (floor(x * pow(2, -precision - 1)) == x * pow(2, -precision - 1) and x != 0) { precision++; + } return precision; } /** - * @brief compute the precision needed in the output of a function - * + * @brief compute the precision needed in the output of a function + * * @param f the function to analyse * @param x the point at which the tightest precision is needed - * @param u the signed gap between the two consecutive numbers at which to compute the precision (ie the ulp) -*/ + * @param u the signed gap between the two consecutive numbers at which to compute the precision (ie + * the ulp) + */ int exactPrecisionUnary(ufun f, long double x, long double u) { int res = floor((double)log2(std::abs(f(x + u) - f(x)))); @@ -62,19 +64,16 @@ int exactPrecisionUnary(ufun f, long double x, long double u) } /** - * @brief compute the precision needed in the input of a function - * + * @brief compute the precision needed in the input of a function + * * @param f the function to analyse * @param finv a function such that f o finv = Id locally around x * @param x the input point at which the tightest precision is needed - * @param u the signed gap between the two consecutive numbers at which to compute the precision (ie the ulp) -*/ + * @param u the signed gap between the two consecutive numbers at which to compute the precision (ie + * the ulp) + */ int exactPrecisionUnaryBackwards(ufun f, ufun finv, double x, double u) { - int res = ceil( - (double)log2( - std::abs(finv(f(x) + u) - x) - ) - ); + int res = ceil((double)log2(std::abs(finv(f(x) + u) - x))); return res; } \ No newline at end of file diff --git a/interval/utils.hh b/interval/utils.hh index 06f119d..dc2489c 100644 --- a/interval/utils.hh +++ b/interval/utils.hh @@ -40,42 +40,48 @@ static double max4(int a, int b, int c, int d) /** * @brief Computes the value with minimum absolute value among the bounds of an interval -*/ + */ static double minValAbs(interval x) { - if (std::abs(x.lo()) < std::abs(x.hi())) + if (std::abs(x.lo()) < std::abs(x.hi())) { return x.lo(); + } return x.hi(); } /** * @brief Computes the value with maximum absolute value -*/ + */ static double maxValAbs(interval x) { - if (std::abs(x.lo()) < std::abs(x.hi())) + if (std::abs(x.lo()) < std::abs(x.hi())) { return x.hi(); + } return x.lo(); } /** - * @brief Computes the direction of the interior of the interval at the minimum absolute value of its bounds -*/ + * @brief Computes the direction of the interior of the interval at the minimum absolute value of + * its bounds + */ static int signMinValAbs(interval x) { - if (std::abs(x.lo()) < std::abs(x.hi())) + if (std::abs(x.lo()) < std::abs(x.hi())) { return 1; + } return -1; } /** - * @brief Computes the direction of the interior of the interval at the maximum absolute value of its bounds -*/ + * @brief Computes the direction of the interior of the interval at the maximum absolute value of + * its bounds + */ static int signMaxValAbs(interval x) { - if (std::abs(x.lo()) < std::abs(x.hi())) + if (std::abs(x.lo()) < std::abs(x.hi())) { return -1; + } return 1; } -} // namespace itv \ No newline at end of file +} // namespace itv \ No newline at end of file diff --git a/main.cpp b/main.cpp index dd58285..d3e3a91 100644 --- a/main.cpp +++ b/main.cpp @@ -28,9 +28,9 @@ using namespace itv; void print_help() { std::cout << "Usage: ./TestInterval [option]\n" - << "Options:\n" - << " -h Show this help message\n" - << " -norand Enable deterministic computation (non random)\n"; + << "Options:\n" + << " -h Show this help message\n" + << " -norand Enable deterministic computation (non random)\n"; } int main(int argc, char* argv[]) @@ -38,20 +38,20 @@ int main(int argc, char* argv[]) // Check if an argument is passed if (argc > 1) { std::string arg = argv[1]; - if (arg == "-h" || arg == "--help" ) { + if (arg == "-h" || arg == "--help") { print_help(); return 0; // Exit after showing help } else if (arg == "-norand") { gRandom = false; } } - + // test interval representation // check("interval()", interval()); check("interval(0,100,-24)", interval(100.0, 0.0)); check("interval(0,0,-24)", interval(0, 0)); check("interval(-10,0,-24)", interval(0, -10)); - check("interval(-10,10,-5)", interval(-10,10,-5)); + check("interval(-10,10,-5)", interval(-10, 10, -5)); // test union intersection @@ -108,10 +108,10 @@ int main(int argc, char* argv[]) interval_algebra A; A.testAll(); - + /*A.testExp(); A.testLog();*/ - /*A.testAcos(); + /*A.testAcos(); A.testAcosh(); A.testAsin();*/ // A.testAsinh(); @@ -186,33 +186,33 @@ int main(int argc, char* argv[]) /* interval X = interval(-10, 10, -24); std::vector titles{ "sin", "exp", "cos"}; std::vector mps{&interval_algebra::Sin, &interval_algebra::Exp, &interval_algebra::Cos}; - std::cout << std::endl << "sin(exp(cos("<< X << "))) = " << A.Sin(A.Exp(A.Cos(X))) << std::endl << std::endl; - propagateBackwardsComposition(titles, mps, X, -24); - std::cout << std::endl << "sin(exp(cos("<< X << "))) = " << A.Sin(A.Exp(A.Cos(X))) << std::endl; - std::cout << "----------------" << std::endl << std::endl; + std::cout << std::endl << "sin(exp(cos("<< X << "))) = " << A.Sin(A.Exp(A.Cos(X))) << std::endl + << std::endl; propagateBackwardsComposition(titles, mps, X, -24); std::cout << std::endl << + "sin(exp(cos("<< X << "))) = " << A.Sin(A.Exp(A.Cos(X))) << std::endl; std::cout << + "----------------" << std::endl << std::endl; X = interval(1, 10, -24); std::vector titles2{"log", "exp"}; std::vector mps2{ &interval_algebra::Log, &interval_algebra::Exp}; - std::cout << std::endl << "log(exp("<< X << ")) = " << A.Log(A.Exp(X)) << std::endl << std::endl; - propagateBackwardsComposition(titles2, mps2, X, -24); - std::cout << std::endl << "log(exp("<< X << ")) = " << A.Log(A.Exp(X)) << std::endl; - std::cout << "----------------" << std::endl << std::endl; + std::cout << std::endl << "log(exp("<< X << ")) = " << A.Log(A.Exp(X)) << std::endl << + std::endl; propagateBackwardsComposition(titles2, mps2, X, -24); std::cout << std::endl << + "log(exp("<< X << ")) = " << A.Log(A.Exp(X)) << std::endl; std::cout << "----------------" << + std::endl << std::endl; X = interval(1, 10, -24); std::vector titles3{"exp", "log"}; std::vector mps3{ &interval_algebra::Exp, &interval_algebra::Log}; - std::cout << std::endl << "exp(log("<< X << ")) = " << A.Exp(A.Log(X)) << std::endl << std::endl; - propagateBackwardsComposition(titles3, mps3, X, -24); - std::cout << std::endl << "exp(log("<< X << ")) = " << A.Exp(A.Log(X)) << std::endl; - std::cout << "----------------" << std::endl << std::endl; */ + std::cout << std::endl << "exp(log("<< X << ")) = " << A.Exp(A.Log(X)) << std::endl << + std::endl; propagateBackwardsComposition(titles3, mps3, X, -24); std::cout << std::endl << + "exp(log("<< X << ")) = " << A.Exp(A.Log(X)) << std::endl; std::cout << "----------------" << + std::endl << std::endl; */ /* interval X = interval(1, 2, -24); interval Y = interval(1, 10, -24); std::cout << "pow(" << X << ", " << Y << ") = " << A.Pow(X, Y) << std::endl << std::endl; propagateBackwardsBinaryMethod("pow", &interval_algebra::Pow, X, Y, -24); std::cout << std::endl; - std::cout << "pow(" << X << ", " << Y << ") = " << A.Pow(X, Y) << std::endl; + std::cout << "pow(" << X << ", " << Y << ") = " << A.Pow(X, Y) << std::endl; interval Z = A.Pow(X, Y); std::cout << "exp(" << Z << ") = " << A.Exp(Z) << std::endl; propagateBackwardsUnaryMethod("exp", &interval_algebra::Exp, Z, -24); @@ -262,4 +262,3 @@ int main(int argc, char* argv[]) interval Y = interval(23.4489,23.4489,-958); std::cout << X << "/" << Y << " = " << A.Div(X, Y) << std::endl;*/ } - diff --git a/tlib/exception.hh b/tlib/exception.hh index 1a2b8f5..f0c33ec 100644 --- a/tlib/exception.hh +++ b/tlib/exception.hh @@ -33,13 +33,15 @@ Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #include #endif -class faustexception : public std::runtime_error -{ +class faustexception : public std::runtime_error { public: #ifdef EMCC static const char* gJSExceptionMsg; - faustexception(const std::string& msg = "") : std::runtime_error(msg) { gJSExceptionMsg = strdup(msg.c_str()); } + faustexception(const std::string& msg = "") : std::runtime_error(msg) + { + gJSExceptionMsg = strdup(msg.c_str()); + } faustexception(char* msg) : std::runtime_error(msg) { gJSExceptionMsg = strdup(msg); } faustexception(const char* msg) : std::runtime_error(msg) { gJSExceptionMsg = strdup(msg); } #else diff --git a/tlib/garbageable.hh b/tlib/garbageable.hh index a358c68..a40bd85 100644 --- a/tlib/garbageable.hh +++ b/tlib/garbageable.hh @@ -32,8 +32,7 @@ // To be inherited by all garbageable classes -class LIBFAUST_API Garbageable -{ +class LIBFAUST_API Garbageable { public: Garbageable() {} virtual ~Garbageable() {} @@ -47,8 +46,7 @@ class LIBFAUST_API Garbageable }; template -class GarbageablePtr : public virtual Garbageable -{ +class GarbageablePtr : public virtual Garbageable { private: P* fPtr; diff --git a/tlib/list.cpp b/tlib/list.cpp index 4074649..a9d29f3 100644 --- a/tlib/list.cpp +++ b/tlib/list.cpp @@ -533,7 +533,8 @@ Tree tmap(Tree key, tfun f, Tree t) static Tree substkey(Tree t, Tree id, Tree val) { char name[256]; - snprintf(name, 255, "SUBST<%p,%p,%p> : ", (void*)(CTree*)t, (void*)(CTree*)id, (void*)(CTree*)val); + snprintf(name, 255, "SUBST<%p,%p,%p> : ", (void*)(CTree*)t, (void*)(CTree*)id, + (void*)(CTree*)val); return tree(unique(name)); } diff --git a/tlib/node.hh b/tlib/node.hh index 5f0c122..d1965ac 100644 --- a/tlib/node.hh +++ b/tlib/node.hh @@ -66,8 +66,7 @@ enum NodeType { kIntNode, kInt64Node, kDoubleNode, kSymNode, kPointerNode }; * Class Node = (type x (int + double + Sym + void*)) */ -class Node : public virtual Garbageable -{ +class Node : public virtual Garbageable { int fType; union { int i; @@ -122,8 +121,14 @@ class Node : public virtual Garbageable void* getPointer() const { return fData.p; } // conversions and promotion for numbers - operator int() const { return (fType == kIntNode) ? fData.i : (fType == kDoubleNode) ? int(fData.f) : 0; } - operator double() const { return (fType == kIntNode) ? double(fData.i) : (fType == kDoubleNode) ? fData.f : 0.0; } + operator int() const + { + return (fType == kIntNode) ? fData.i : (fType == kDoubleNode) ? int(fData.f) : 0; + } + operator double() const + { + return (fType == kIntNode) ? double(fData.i) : (fType == kDoubleNode) ? fData.f : 0.0; + } std::ostream& print(std::ostream& fout) const; ///< print a node on a stream }; @@ -188,27 +193,32 @@ inline bool isDouble(const Node& n, double* x) inline bool isZero(const Node& n) { - return ((n.type() == kDoubleNode) && (n.getDouble() == 0.0)) || ((n.type() == kIntNode) && (n.getInt() == 0)); + return ((n.type() == kDoubleNode) && (n.getDouble() == 0.0)) || + ((n.type() == kIntNode) && (n.getInt() == 0)); } inline bool isGEZero(const Node& n) { - return ((n.type() == kDoubleNode) && (n.getDouble() >= 0.0)) || ((n.type() == kIntNode) && (n.getInt() >= 0)); + return ((n.type() == kDoubleNode) && (n.getDouble() >= 0.0)) || + ((n.type() == kIntNode) && (n.getInt() >= 0)); } inline bool isGTZero(const Node& n) { - return ((n.type() == kDoubleNode) && (n.getDouble() > 0.0)) || ((n.type() == kIntNode) && (n.getInt() > 0)); + return ((n.type() == kDoubleNode) && (n.getDouble() > 0.0)) || + ((n.type() == kIntNode) && (n.getInt() > 0)); } inline bool isOne(const Node& n) { - return ((n.type() == kDoubleNode) && (n.getDouble() == 1.0)) || ((n.type() == kIntNode) && (n.getInt() == 1)); + return ((n.type() == kDoubleNode) && (n.getDouble() == 1.0)) || + ((n.type() == kIntNode) && (n.getInt() == 1)); } inline bool isMinusOne(const Node& n) { - return ((n.type() == kDoubleNode) && (n.getDouble() == -1.0)) || ((n.type() == kIntNode) && (n.getInt() == -1)); + return ((n.type() == kDoubleNode) && (n.getDouble() == -1.0)) || + ((n.type() == kIntNode) && (n.getInt() == -1)); } bool sameMagnitude(const Node& a, const Node& b); @@ -275,11 +285,16 @@ inline Node mulNode(const Node& x, const Node& y) inline Node divExtendedNode(const Node& x, const Node& y) { if (isDouble(x) || isDouble(y)) { - if (double(y) == 0) goto raise_exception; + if (double(y) == 0) { + goto raise_exception; + } return Node(double(x) / double(y)); } else { - if (int(y) == 0 || double(y) == 0) goto raise_exception; - return (double(int(x) / int(y)) == double(x) / double(y)) ? Node(int(x) / int(y)) : Node(double(x) / double(y)); + if (int(y) == 0 || double(y) == 0) { + goto raise_exception; + } + return (double(int(x) / int(y)) == double(x) / double(y)) ? Node(int(x) / int(y)) + : Node(double(x) / double(y)); } raise_exception: diff --git a/tlib/occur.hh b/tlib/occur.hh index 7784c55..02df8b7 100644 --- a/tlib/occur.hh +++ b/tlib/occur.hh @@ -34,7 +34,7 @@ class Occur : public virtual Garbageable { Tree fKey; // specific property key public: - Occur(Tree root); // count the occurrences of each subtree of root + Occur(Tree root); // count the occurrences of each subtree of root int getCount(Tree t); // return the number of occurrences of t in root private: diff --git a/tlib/recursive-tree.cpp b/tlib/recursive-tree.cpp index 2a110b4..86c8cf6 100644 --- a/tlib/recursive-tree.cpp +++ b/tlib/recursive-tree.cpp @@ -141,7 +141,9 @@ int CTree::calcTreeAperture(const Node& n, const tvec& br) tvec::const_iterator b = br.begin(); tvec::const_iterator z = br.end(); while (b != z) { - if ((*b)->aperture() > rc) rc = (*b)->aperture(); + if ((*b)->aperture() > rc) { + rc = (*b)->aperture(); + } ++b; } return rc; @@ -281,8 +283,12 @@ static Tree calcsubstitute(Tree t, int level, Tree id) // fprintf(stderr, "aperture %d < level %d !!\n", t->aperture(), level); return t; } - if (isRef(t, l)) return (l == level) ? id : t; - if (isRec(t, body)) return rec(substitute(body, level + 1, id)); + if (isRef(t, l)) { + return (l == level) ? id : t; + } + if (isRec(t, body)) { + return rec(substitute(body, level + 1, id)); + } int ar = t->arity(); // Tree br[4]; diff --git a/tlib/shlysis.cpp b/tlib/shlysis.cpp index 66dc5e6..05e924e 100644 --- a/tlib/shlysis.cpp +++ b/tlib/shlysis.cpp @@ -136,7 +136,9 @@ static void annotate(Tree k, Tree t, barrier foo) } else { int n = t->arity(); if (n > 0 && !foo(t)) { - for (int i = 0; i < n; i++) annotate(k, t->branch(i), foo); + for (int i = 0; i < n; i++) { + annotate(k, t->branch(i), foo); + } } } } else { diff --git a/tlib/symbol.cpp b/tlib/symbol.cpp index f971bfb..3f750da 100644 --- a/tlib/symbol.cpp +++ b/tlib/symbol.cpp @@ -67,7 +67,9 @@ Symbol* Symbol::get(const char* rawstr) int bckt = hsh % kHashTableSize; Symbol* item = gSymbolTable[bckt]; - while (item && !item->equiv(hsh, str.c_str())) item = item->fNext; + while (item && !item->equiv(hsh, str.c_str())) { + item = item->fNext; + } Symbol* r = item ? item : gSymbolTable[bckt] = new Symbol(str, hsh, gSymbolTable[bckt]); return r; @@ -85,13 +87,15 @@ bool Symbol::isnew(const char* str) int bckt = hsh % kHashTableSize; Symbol* item = gSymbolTable[bckt]; - while (item && !item->equiv(hsh, str)) item = item->fNext; + while (item && !item->equiv(hsh, str)) { + item = item->fNext; + } return item == 0; } /** - * Creates a new symbol with a name obtained by concatenating the \p str prefix with a number in order to make it - * unique. \param str the prefix of the name \return a symbol of name \p prefix++n + * Creates a new symbol with a name obtained by concatenating the \p str prefix with a number in + * order to make it unique. \param str the prefix of the name \return a symbol of name \p prefix++n */ Symbol* Symbol::prefix(const char* str) @@ -100,7 +104,9 @@ Symbol* Symbol::prefix(const char* str) for (int n = 0; n < 10000; n++) { snprintf(name, 256, "%s%d", str, gPrefixCounters[str]++); - if (isnew(name)) return get(name); + if (isnew(name)) { + return get(name); + } } faustassert(false); return get("UNIQUEOVERFLOW"); @@ -131,7 +137,9 @@ unsigned int Symbol::calcHashKey(const char* str) { unsigned int h = 0; - while (*str) h = (h << 1) ^ (h >> 20) ^ (*str++); + while (*str) { + h = (h << 1) ^ (h >> 20) ^ (*str++); + } return h; } diff --git a/tlib/symbol.hh b/tlib/symbol.hh index c660491..10a6e91 100644 --- a/tlib/symbol.hh +++ b/tlib/symbol.hh @@ -47,18 +47,19 @@ /** * Symbols are unique objects with a name stored in a hash table. */ -class Symbol : public virtual Garbageable -{ +class Symbol : public virtual Garbageable { private: - static const int kHashTableSize = 511; ///< Size of the hash table (a prime number is recommended) - static Symbol* gSymbolTable[kHashTableSize]; ///< Hash table used to store the symbols + static const int kHashTableSize = + 511; ///< Size of the hash table (a prime number is recommended) + static Symbol* gSymbolTable[kHashTableSize]; ///< Hash table used to store the symbols static std::map gPrefixCounters; // Fields - std::string fName; ///< Name of the symbol - unsigned int fHash; ///< Hash key computed from the name and used to determine the hash table entry - Symbol* fNext; ///< Next symbol in the hash table entry - void* fData; ///< Field to user disposal to store additional data + std::string fName; ///< Name of the symbol + unsigned int + fHash; ///< Hash key computed from the name and used to determine the hash table entry + Symbol* fNext; ///< Next symbol in the hash table entry + void* fData; ///< Field to user disposal to store additional data // Constructors & destructors Symbol(const std::string&, unsigned int hsh, @@ -66,15 +67,17 @@ class Symbol : public virtual Garbageable ~Symbol(); ///< The destructor is never used // Others - bool equiv(unsigned int hash, - const char* str) const; ///< Check if the name of the symbol is equal to string \p str - static unsigned int calcHashKey(const char* str); ///< Compute the 32-bits hash key of string \p str + bool equiv( + unsigned int hash, + const char* str) const; ///< Check if the name of the symbol is equal to string \p str + static unsigned int calcHashKey( + const char* str); ///< Compute the 32-bits hash key of string \p str // Static methods static Symbol* get(const std::string& str); ///< Get the symbol of name \p str static Symbol* get(const char* str); ///< Get the symbol of name \p str - static Symbol* prefix(const char* str); ///< Creates a new symbol of name prefixed by \p str - static bool isnew(const char* str); ///< Returns \b true if no symbol of name \p str exists + static Symbol* prefix(const char* str); ///< Creates a new symbol of name prefixed by \p str + static bool isnew(const char* str); ///< Returns \b true if no symbol of name \p str exists public: std::ostream& print(std::ostream& fout) const; ///< print a symbol on a stream diff --git a/tlib/tree.cpp b/tlib/tree.cpp index 954c633..c966cfe 100644 --- a/tlib/tree.cpp +++ b/tlib/tree.cpp @@ -159,7 +159,9 @@ Tree CTree::make(const Node& n, int ar, Tree* tbl) { tvec br(ar); - for (int i = 0; i < ar; i++) br[i] = tbl[i]; + for (int i = 0; i < ar; i++) { + br[i] = tbl[i]; + } size_t hk = calcTreeHash(n, br); Tree t = gHashTable[hk % kHashTableSize]; diff --git a/tlib/tree.hh b/tlib/tree.hh index 7043043..d30bd34 100644 --- a/tlib/tree.hh +++ b/tlib/tree.hh @@ -100,16 +100,16 @@ typedef std::vector tvec; * WARNING : in the current implementation CTrees are allocated but never deleted. **/ -class LIBFAUST_API CTree : public virtual Garbageable -{ +class LIBFAUST_API CTree : public virtual Garbageable { private: static const int kHashTableSize = 400009; ///< size of the hash table (prime number) static size_t gSerialCounter; ///< the serial number counter static Tree gHashTable[kHashTableSize]; ///< hash table used for "hash consing" public: - static bool gDetails; ///< Ctree::print() print with more details when true - static unsigned int gVisitTime; ///< Should be incremented for each new visit to keep track of visited tree + static bool gDetails; ///< Ctree::print() print with more details when true + static unsigned int + gVisitTime; ///< Should be incremented for each new visit to keep track of visited tree private: // fields @@ -123,32 +123,44 @@ class LIBFAUST_API CTree : public virtual Garbageable unsigned int fVisitTime; ///< keep track of visits tvec fBranch; ///< the subtrees - CTree(size_t hk, const Node& n, const tvec& br); ///< construction is private, uses tree::make instead + CTree(size_t hk, const Node& n, + const tvec& br); ///< construction is private, uses tree::make instead - bool equiv(const Node& n, const tvec& br) const; ///< used to check if an equivalent tree already exists - static size_t calcTreeHash(const Node& n, - const tvec& br); ///< compute the hash key of a tree according to its node and branches - static int calcTreeAperture(const Node& n, const tvec& br); ///< compute how open is a tree + bool equiv(const Node& n, + const tvec& br) const; ///< used to check if an equivalent tree already exists + static size_t calcTreeHash( + const Node& n, + const tvec& br); ///< compute the hash key of a tree according to its node and branches + static int calcTreeAperture(const Node& n, const tvec& br); ///< compute how open is a tree public: virtual ~CTree(); - static Tree make(const Node& n, int ar, Tree br[]); ///< return a new tree or an existing equivalent one - static Tree make(const Node& n, const tvec& br); ///< return a new tree or an existing equivalent one + static Tree make(const Node& n, int ar, + Tree br[]); ///< return a new tree or an existing equivalent one + static Tree make(const Node& n, + const tvec& br); ///< return a new tree or an existing equivalent one // Accessors - const Node& node() const { return fNode; } ///< return the content of the tree - int arity() const { return (int)fBranch.size(); } ///< return the number of branches (subtrees) of a tree - Tree branch(int i) const { return fBranch[i]; } ///< return the ith branch (subtree) of a tree - const tvec& branches() const { return fBranch; } ///< return all branches (subtrees) of a tree - size_t hashkey() const { return fHashKey; } ///< return the hashkey of the tree - size_t serial() const { return fSerial; } ///< return the serial of the tree - int aperture() const { return fAperture; } ///< return how "open" is a tree in terms of free variables - void setAperture(int a) { fAperture = a; } ///< modify the aperture of a tree + const Node& node() const { return fNode; } ///< return the content of the tree + int arity() const + { + return (int)fBranch.size(); + } ///< return the number of branches (subtrees) of a tree + Tree branch(int i) const { return fBranch[i]; } ///< return the ith branch (subtree) of a tree + const tvec& branches() const { return fBranch; } ///< return all branches (subtrees) of a tree + size_t hashkey() const { return fHashKey; } ///< return the hashkey of the tree + size_t serial() const { return fSerial; } ///< return the serial of the tree + int aperture() const + { + return fAperture; + } ///< return how "open" is a tree in terms of free variables + void setAperture(int a) { fAperture = a; } ///< modify the aperture of a tree // Print a tree and the hash table (for debugging purposes) - std::ostream& print(std::ostream& fout) const; ///< print recursively the content of a tree on a stream - static void control(); ///< print the hash table content (for debug purpose) + std::ostream& print( + std::ostream& fout) const; ///< print recursively the content of a tree on a stream + static void control(); ///< print the hash table content (for debug purpose) static void init(); @@ -208,7 +220,8 @@ inline Tree tree(const Node& n, const Tree& a, const Tree& b, const Tree& c, con return CTree::make(n, 4, br); } -inline Tree tree(const Node& n, const Tree& a, const Tree& b, const Tree& c, const Tree& d, const Tree& e) +inline Tree tree(const Node& n, const Tree& a, const Tree& b, const Tree& c, const Tree& d, + const Tree& e) { Tree br[] = {a, b, c, d, e}; return CTree::make(n, 5, br); @@ -219,13 +232,15 @@ inline Tree tree(const Node& n, const tvec& br) } // useful conversions -LIBFAUST_API int tree2int(Tree t); ///< if t has a node of type int, return it otherwise error -double tree2float(Tree t); ///< if t has a node of type float, return it otherwise error -double tree2double(Tree t); ///< if t has a node of type float, return it otherwise error -LIBFAUST_API const char* tree2str(Tree t); ///< if t has a node of type symbol, return its name otherwise error -std::string tree2quotedstr(Tree t); -void* tree2ptr(Tree t); ///< if t has a node of type ptr, return it otherwise error -LIBFAUST_API void* getUserData(Tree t); ///< if t has a node of type symbol, return the associated user data +LIBFAUST_API int tree2int(Tree t); ///< if t has a node of type int, return it otherwise error +double tree2float(Tree t); ///< if t has a node of type float, return it otherwise error +double tree2double(Tree t); ///< if t has a node of type float, return it otherwise error +LIBFAUST_API const char* tree2str( + Tree t); ///< if t has a node of type symbol, return its name otherwise error +std::string tree2quotedstr(Tree t); +void* tree2ptr(Tree t); ///< if t has a node of type ptr, return it otherwise error +LIBFAUST_API void* getUserData( + Tree t); ///< if t has a node of type symbol, return the associated user data // pattern matching bool isTree(const Tree& t, const Node& n); @@ -280,8 +295,7 @@ Tree deBruijn2Sym(Tree t); ////< transform a tree from deBruijn to symbolic not //--------------------------------------------------------------------------- -class Tabber -{ +class Tabber { int fIndent; int fPostInc; @@ -301,7 +315,9 @@ class Tabber std::ostream& print(std::ostream& fout) { - for (int i = 0; i < fIndent; i++) fout << '\t'; + for (int i = 0; i < fIndent; i++) { + fout << '\t'; + } fIndent += fPostInc; fPostInc = 0; return fout;