dbldblarith.h

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the program and library             */
/*         SCIP --- Solving Constraint Integer Programs                      */
/*                                                                           */
/*    Copyright (C) 2002-2020 Konrad-Zuse-Zentrum                            */
/*                            fuer Informationstechnik Berlin                */
/*                                                                           */
/*  SCIP is distributed under the terms of the ZIB Academic License.         */
/*                                                                           */
/*  You should have received a copy of the ZIB Academic License              */
/*  along with SCIP; see the file COPYING. If not visit scipopt.org.         */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/**@file   dbldblarith.h
 * @brief  defines macros for basic operations in double-double arithmetic giving roughly twice the precision of a double
 * @author Leona Gottwald
 *
 */

/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/

#ifndef _SCIP_DBLDBL_ARITH_
#define _SCIP_DBLDBL_ARITH_

#include "math.h"


#ifndef DISABLE_QUADPREC

/* smaller epsilon value for use with quadprecision */
#define QUAD_EPSILON 1e-12

/* convenience macros for nicer usage of double double arithmetic */
#define QUAD_HI(x)  x ## hi
#define QUAD_LO(x)  x ## lo
#define QUAD(x) QUAD_HI(x), QUAD_LO(x)
#define QUAD_MEMBER(x) QUAD_HI(x); QUAD_LO(x)
#define QUAD_TO_DBL(x) ( QUAD_HI(x) + QUAD_LO(x) )
#define QUAD_SCALE(x, a) do { QUAD_HI(x) *= (a); QUAD_LO(x) *= (a); } while(0)
#define QUAD_ASSIGN(a, constant)  do { QUAD_HI(a) = (constant); QUAD_LO(a) = 0.0; } while(0)
#define QUAD_ASSIGN_Q(a, b)  do { QUAD_HI(a) = QUAD_HI(b); QUAD_LO(a) = QUAD_LO(b); } while(0)
#define QUAD_ARRAY_SIZE(size) ((size)*2)
#define QUAD_ARRAY_LOAD(r, a, idx) do { QUAD_HI(r) = (a)[2*(idx)]; QUAD_LO(r) = (a)[2*(idx) + 1]; } while(0)
#define QUAD_ARRAY_STORE(a, idx, x) do { (a)[2*(idx)] = QUAD_HI(x); (a)[2*(idx) + 1] = QUAD_LO(x); } while(0)

/* define all the SCIPquadprec... macros such that they use the SCIPdbldbl... macros that expands the quad precision arguments using the above macros */
#define SCIPquadprecProdDD(r, a, b)  SCIPdbldblProd(QUAD_HI(r), QUAD_LO(r), a, b)
#define SCIPquadprecSquareD(r, a) SCIPdbldblSquare(QUAD_HI(r), QUAD_LO(r), a)
#define SCIPquadprecSumDD(r, a, b) SCIPdbldblSum(QUAD_HI(r), QUAD_LO(r), a, b)
#define SCIPquadprecDivDD(r, a, b) SCIPdbldblDiv(QUAD_HI(r), QUAD_LO(r), a, b)
#define SCIPquadprecSumQD(r, a, b) SCIPdbldblSum21(QUAD_HI(r), QUAD_LO(r), QUAD_HI(a), QUAD_LO(a), b)
#define SCIPquadprecProdQD(r, a, b) SCIPdbldblProd21(QUAD_HI(r), QUAD_LO(r), QUAD_HI(a), QUAD_LO(a), b)
#define SCIPquadprecDivDQ(r, a, b) SCIPdbldblDiv12(QUAD_HI(r), QUAD_LO(r), a, QUAD_HI(b), QUAD_LO(b))
#define SCIPquadprecDivQD(r, a, b) SCIPdbldblDiv21(QUAD_HI(r), QUAD_LO(r), QUAD_HI(a), QUAD_LO(a), b)
#define SCIPquadprecProdQQ(r, a, b) SCIPdbldblProd22(QUAD_HI(r), QUAD_LO(r), QUAD_HI(a), QUAD_LO(a), QUAD_HI(b), QUAD_LO(b))
#define SCIPquadprecSumQQ(r, a, b) SCIPdbldblSum22(QUAD_HI(r), QUAD_LO(r), QUAD_HI(a), QUAD_LO(a), QUAD_HI(b), QUAD_LO(b))
#define SCIPquadprecSquareQ(r, a) SCIPdbldblSquare2(QUAD_HI(r), QUAD_LO(r), QUAD_HI(a), QUAD_LO(a))
#define SCIPquadprecDivQQ(r, a, b) SCIPdbldblDiv22(QUAD_HI(r), QUAD_LO(r), QUAD_HI(a), QUAD_LO(a), QUAD_HI(b), QUAD_LO(b))
#define SCIPquadprecSqrtD(r, a) SCIPdbldblSqrt(QUAD_HI(r), QUAD_LO(r), a)
#define SCIPquadprecSqrtQ(r, a) SCIPdbldblSqrt2(QUAD_HI(r), QUAD_LO(r), QUAD_HI(a), QUAD_LO(a))
#define SCIPquadprecAbsQ(r, a) SCIPdbldblAbs2(QUAD_HI(r), QUAD_LO(r), QUAD_HI(a), QUAD_LO(a))
#define SCIPquadprecFloorQ(r, a) SCIPdbldblFloor2(QUAD_HI(r), QUAD_LO(r), QUAD_HI(a), QUAD_LO(a))
#define SCIPquadprecCeilQ(r, a) SCIPdbldblCeil2(QUAD_HI(r), QUAD_LO(r), QUAD_HI(a), QUAD_LO(a))
#define SCIPquadprecEpsFloorQ(r, a, eps) SCIPdbldblEpsFloor2(QUAD_HI(r), QUAD_LO(r), QUAD_HI(a), QUAD_LO(a), eps)
#define SCIPquadprecEpsCeilQ(r, a, eps) SCIPdbldblEpsCeil2(QUAD_HI(r), QUAD_LO(r), QUAD_HI(a), QUAD_LO(a), eps)

#else

/* normal epsilon value if quadprecision is disabled */
#define QUAD_EPSILON 1e-9

/* dummy macros that use normal arithmetic */
#define QUAD_HI(x)  x
#define QUAD_LO(x)  0.0
#define QUAD(x)     x
#define QUAD_MEMBER(x) x
#define QUAD_TO_DBL(x) (x)
#define QUAD_SCALE(x, a) do { (x) *= (a); } while(0)
#define QUAD_ASSIGN(a, constant)  do { (a) = constant; } while(0)
#define QUAD_ASSIGN_Q(a, b)  do { (a) = (b); } while(0)
#define QUAD_ARRAY_SIZE(size) (size)
#define QUAD_ARRAY_LOAD(r, a, idx) do { r = (a)[(idx)]; } while(0)
#define QUAD_ARRAY_STORE(a, idx, x) do { (a)[(idx)] = (x); } while(0)

#define SCIPquadprecProdDD(r, a, b)  do { (r) = (a) * (b); } while(0)
#define SCIPquadprecSquareD(r, a)    do { (r) = (a) * (a); } while(0)
#define SCIPquadprecSumDD(r, a, b)   do { (r) = (a) + (b); } while(0)
#define SCIPquadprecDivDD(r, a, b)   do { (r) = (a) / (b); } while(0)
#define SCIPquadprecSumQD(r, a, b)   do { (r) = (a) + (b); } while(0)
#define SCIPquadprecProdQD(r, a, b)  do { (r) = (a) * (b); } while(0)
#define SCIPquadprecDivDQ(r, a, b)   do { (r) = (a) / (b); } while(0)
#define SCIPquadprecDivQD(r, a, b)   do { (r) = (a) / (b); } while(0)
#define SCIPquadprecProdQQ(r, a, b)  do { (r) = (a) * (b); } while(0)
#define SCIPquadprecSumQQ(r, a, b)   do { (r) = (a) + (b); } while(0)
#define SCIPquadprecSquareQ(r, a)    do { (r) = (a) * (a); } while(0)
#define SCIPquadprecDivQQ(r, a, b)   do { (r) = (a) / (b); } while(0)
#define SCIPquadprecSqrtD(r, a)      do { (r) = sqrt(a); } while(0)
#define SCIPquadprecSqrtQ(r, a)      do { (r) = sqrt(a); } while(0)
#define SCIPquadprecAbsQ(r, a)       do { (r) = fabs(a); } while(0)
#define SCIPquadprecFloorQ(r, a)     do { (r) = floor(a); } while(0)
#define SCIPquadprecCeilQ(r, a)      do { (r) = ceil(a); } while(0)
#define SCIPquadprecEpsFloorQ(r, a, eps) do { (r) = floor((a) + (eps)); } while(0)
#define SCIPquadprecEpsCeilQ(r, a, eps) do { (r) = ceil((a) - (eps)); } while(0)

#endif

#define __SCIPdbldblSplit(rhi, rlo, x) \
    do { \
       const double __tmp_split_dbl = 134217729.0 * (x); \
       (rhi) = __tmp_split_dbl - (__tmp_split_dbl - (x)); \
       (rlo) = (x) - (rhi);\
    } while(0)

/** multiply two floating point numbers, both given by one double, and return the result as two doubles. */
#define SCIPdbldblProd(rhi, rlo, a, b) \
    do { \
        double __tmp_dbldbl_prod_ahi; \
        double __tmp_dbldbl_prod_alo; \
        double __tmp_dbldbl_prod_bhi; \
        double __tmp_dbldbl_prod_blo; \
        __SCIPdbldblSplit(__tmp_dbldbl_prod_ahi, __tmp_dbldbl_prod_alo, a); \
        __SCIPdbldblSplit(__tmp_dbldbl_prod_bhi, __tmp_dbldbl_prod_blo, b); \
        (rhi) = (a) * (b); \
        (rlo) = __tmp_dbldbl_prod_alo * __tmp_dbldbl_prod_blo - \
           ((((rhi) - __tmp_dbldbl_prod_ahi * __tmp_dbldbl_prod_bhi) \
           - __tmp_dbldbl_prod_alo * __tmp_dbldbl_prod_bhi) \
           - __tmp_dbldbl_prod_ahi * __tmp_dbldbl_prod_blo); \
    } while(0)

/** square a floating point number given by one double and return the result as two doubles. */
#define SCIPdbldblSquare(rhi, rlo, a) \
    do { \
        double __tmp_dbldbl_square_ahi; \
        double __tmp_dbldbl_square_alo; \
        __SCIPdbldblSplit(__tmp_dbldbl_square_ahi, __tmp_dbldbl_square_alo, a); \
        (rhi) = (a) * (a); \
        (rlo) = __tmp_dbldbl_square_alo * __tmp_dbldbl_square_alo - \
           ((((rhi) - __tmp_dbldbl_square_ahi * __tmp_dbldbl_square_ahi) \
           - 2.0 * __tmp_dbldbl_square_alo * __tmp_dbldbl_square_ahi)); \
    } while(0)

/** add two floating point numbers, both given by one double, and return the result as two doubles. */
#define SCIPdbldblSum(rhi, rlo, a, b) \
    do { \
        double __tmp1_dbldbl_sum; \
        double __tmp2_dbldbl_sum; \
        __tmp2_dbldbl_sum = (a) + (b); \
        __tmp1_dbldbl_sum = __tmp2_dbldbl_sum - (a); \
        (rlo) = ((a) - (__tmp2_dbldbl_sum - __tmp1_dbldbl_sum)) + ((b) - __tmp1_dbldbl_sum); \
        (rhi) = __tmp2_dbldbl_sum; \
    } while(0)

/** divide two floating point numbers, both given by one double, and return the result as two doubles. */
#define SCIPdbldblDiv(rhi, rlo, a, b) \
    do { \
       double __tmp_dbldbl_div_hi; \
       double __tmp_dbldbl_div_lo; \
       double __estim_dbldbl_div = (a)/(b); \
       SCIPdbldblProd(__tmp_dbldbl_div_hi, __tmp_dbldbl_div_lo, b, __estim_dbldbl_div); \
       SCIPdbldblSum21(__tmp_dbldbl_div_hi, __tmp_dbldbl_div_lo, __tmp_dbldbl_div_hi, __tmp_dbldbl_div_lo, -(a)); \
       __tmp_dbldbl_div_hi /= (b); \
       __tmp_dbldbl_div_lo /= (b); \
       SCIPdbldblSum21(rhi, rlo, -__tmp_dbldbl_div_hi, -__tmp_dbldbl_div_lo, __estim_dbldbl_div); \
    } while(0)

/** add two floating point numbers, the first is given by two doubles, the second is given by one double,
 *  and return the result as two doubles.
 */
#define SCIPdbldblSum21(rhi, rlo, ahi, alo, b) \
   do { \
      double __tmp_dbldbl_sum21_hi; \
      double __tmp_dbldbl_sum21_lo; \
      SCIPdbldblSum(__tmp_dbldbl_sum21_hi, __tmp_dbldbl_sum21_lo, ahi, b); \
      (rlo) = __tmp_dbldbl_sum21_lo + (alo); \
      (rhi) = __tmp_dbldbl_sum21_hi; \
   } while(0)


/** multiply two floating point numbers, the first is given by two doubles, the second is given by one double,
 *  and return the result as two doubles.
 */
#define SCIPdbldblProd21(rhi, rlo, ahi, alo, b) \
    do { \
       double __tmp_dbldbl_prod21_hi; \
       double __tmp_dbldbl_prod21_lo; \
       SCIPdbldblProd(__tmp_dbldbl_prod21_hi, __tmp_dbldbl_prod21_lo, ahi, b); \
       (rlo) = (alo) * (b) + __tmp_dbldbl_prod21_lo; \
       (rhi) = __tmp_dbldbl_prod21_hi; \
    } while(0)

/** divide two floating point numbers, the first is given by one double, the second is given by two doubles,
 *  and return the result as two doubles.
 */
#define SCIPdbldblDiv12(rhi, rlo, a, bhi, blo) \
    do { \
       double __tmp_dbldbl_div12_hi; \
       double __tmp_dbldbl_div12_lo; \
       double __estim_dbldbl_div12 = (a)/(bhi); \
       SCIPdbldblProd21(__tmp_dbldbl_div12_hi, __tmp_dbldbl_div12_lo, bhi, blo, __estim_dbldbl_div12); \
       SCIPdbldblSum21(__tmp_dbldbl_div12_hi, __tmp_dbldbl_div12_lo, __tmp_dbldbl_div12_hi, __tmp_dbldbl_div12_lo, -(a)); \
       __tmp_dbldbl_div12_hi /= (bhi); \
       __tmp_dbldbl_div12_lo /= (bhi); \
       SCIPdbldblSum21(rhi, rlo, -__tmp_dbldbl_div12_hi, -__tmp_dbldbl_div12_lo, __estim_dbldbl_div12); \
    } while(0)


/** divide two floating point numbers, the first is given by two doubles, the second is given by one double,
 *  and return the result as two doubles.
 */
#define SCIPdbldblDiv21(rhi, rlo, ahi, alo, b) \
   do { \
      double __tmp_dbldbl_div21_hi; \
      double __tmp_dbldbl_div21_lo; \
      double __estim_dbldbl_div21_hi; \
      double __estim_dbldbl_div21_lo; \
      __estim_dbldbl_div21_hi = (ahi)/(b); \
      __estim_dbldbl_div21_lo = (alo)/(b); \
      SCIPdbldblProd21(__tmp_dbldbl_div21_hi, __tmp_dbldbl_div21_lo, __estim_dbldbl_div21_hi, __estim_dbldbl_div21_lo, b); \
      SCIPdbldblSum22(__tmp_dbldbl_div21_hi, __tmp_dbldbl_div21_lo, __tmp_dbldbl_div21_hi, __tmp_dbldbl_div21_lo, -(ahi), -(alo)); \
      __tmp_dbldbl_div21_hi /= (b); \
      __tmp_dbldbl_div21_lo /= (b); \
      SCIPdbldblSum22(rhi, rlo, __estim_dbldbl_div21_hi, __estim_dbldbl_div21_lo, -__tmp_dbldbl_div21_hi, -__tmp_dbldbl_div21_lo); \
   } while(0)

/** multiply two floating point numbers, both given by two doubles, and return the result as two doubles. */
#define SCIPdbldblProd22(rhi, rlo, ahi, alo, bhi, blo) \
   do { \
      double __tmp_dbldbl_prod22_hi; \
      double __tmp_dbldbl_prod22_lo; \
      SCIPdbldblProd(__tmp_dbldbl_prod22_hi, __tmp_dbldbl_prod22_lo, ahi, bhi); \
      SCIPdbldblSum21(__tmp_dbldbl_prod22_hi, __tmp_dbldbl_prod22_lo, \
                   __tmp_dbldbl_prod22_hi, __tmp_dbldbl_prod22_lo, (alo) * (bhi)); \
      SCIPdbldblSum21(rhi, rlo, \
                   __tmp_dbldbl_prod22_hi, __tmp_dbldbl_prod22_lo, (ahi) * (blo)); \
   } while(0)

/** add two floating point numbers, both given by two doubles, and return the result as two doubles. */
#define SCIPdbldblSum22(rhi, rlo, ahi, alo, bhi, blo) \
   do { \
      double __tmp_dbldbl_sum22_hi; \
      double __tmp_dbldbl_sum22_lo; \
      SCIPdbldblSum21(__tmp_dbldbl_sum22_hi, __tmp_dbldbl_sum22_lo, ahi, alo, bhi); \
      SCIPdbldblSum21(rhi, rlo, __tmp_dbldbl_sum22_hi, __tmp_dbldbl_sum22_lo, blo); \
   } while(0)

/** square a floating point number given by two doubles and return the result as two doubles. */
#define SCIPdbldblSquare2(rhi, rlo, ahi, alo) \
   do { \
      double __tmp_dbldbl_square2_hi; \
      double __tmp_dbldbl_square2_lo; \
      SCIPdbldblSquare(__tmp_dbldbl_square2_hi, __tmp_dbldbl_square2_lo, (ahi)); \
      SCIPdbldblSum21(rhi, rlo, __tmp_dbldbl_square2_hi, __tmp_dbldbl_square2_lo, 2 * (ahi) * (alo)); \
   } while(0)

/** divide two floating point numbers, both given by two doubles, and return the result as two doubles. */
#define SCIPdbldblDiv22(rhi, rlo, ahi, alo, bhi, blo) \
   do { \
      double __tmp_dbldbl_div22_hi; \
      double __tmp_dbldbl_div22_lo; \
      double __estim_dbldbl_div22_hi = (ahi) / (bhi); \
      double __estim_dbldbl_div22_lo = (alo) / (bhi); \
      SCIPdbldblProd22(__tmp_dbldbl_div22_hi, __tmp_dbldbl_div22_lo, \
                    bhi, blo, __estim_dbldbl_div22_hi, __estim_dbldbl_div22_lo); \
      SCIPdbldblSum22(__tmp_dbldbl_div22_hi, __tmp_dbldbl_div22_lo, \
                   __tmp_dbldbl_div22_hi, __tmp_dbldbl_div22_lo, -(ahi), -(alo)); \
      __tmp_dbldbl_div22_hi /= (bhi); \
      __tmp_dbldbl_div22_lo /= (bhi); \
      SCIPdbldblSum22(rhi, rlo, __estim_dbldbl_div22_hi, __estim_dbldbl_div22_lo, \
                                -__tmp_dbldbl_div22_hi, -__tmp_dbldbl_div22_lo); \
   } while(0)


/** take the square root of a floating point number given by one double and return the result as two doubles. */
#define SCIPdbldblSqrt(rhi, rlo, a) \
   do { \
      double __estim_dbldbl_sqrt = sqrt(a); \
      SCIPdbldblDiv(rhi, rlo, a, __estim_dbldbl_sqrt); \
      SCIPdbldblSum21(rhi, rlo, rhi, rlo, __estim_dbldbl_sqrt); \
      (rhi) *= 0.5; \
      (rlo) *= 0.5; \
   } while(0)


/** take the square root of a floating point number given by two doubles and return the result as two doubles. */
#define SCIPdbldblSqrt2(rhi, rlo, ahi, alo) \
   do { \
      double __estim_dbldbl_sqrt2 = sqrt(ahi); \
      SCIPdbldblDiv21(rhi, rlo, ahi, alo, __estim_dbldbl_sqrt2); \
      SCIPdbldblSum21(rhi, rlo, rhi, rlo, __estim_dbldbl_sqrt2); \
      (rhi) *= 0.5; \
      (rlo) *= 0.5; \
   } while(0)

/** compute the absolute value of the floating point number given by two doubles */
#define SCIPdbldblAbs2(rhi, rlo, ahi, alo) \
   do { \
      if( ahi < 0.0 ) \
      { \
         (rhi) = -(ahi); \
         (rlo) = -(alo); \
      } \
      else \
      { \
         (rhi) = (ahi); \
         (rlo) = (alo); \
      } \
   } while(0)

/** compute the floored value of the floating point number given by two doubles */
#define SCIPdbldblFloor2(rhi, rlo, ahi, alo) \
   do { \
      double __tmp_dbldbl_floor; \
      __tmp_dbldbl_floor = floor((ahi) + (alo)); \
      SCIPdbldblSum21(rhi, rlo, ahi, alo, -__tmp_dbldbl_floor); \
      if( ((rhi) - 1.0) + (rlo) < 0.0 && (rhi) + (rlo) >= 0.0 ) \
      { \
         /* floor in double precision was fine */ \
         (rhi) = __tmp_dbldbl_floor; \
         (rlo) = 0.0; \
      } \
      else \
      { \
         /* floor in double precision needs to be corrected */ \
         double __tmp2_dbldbl_floor = floor((rhi) + (rlo)); \
         SCIPdbldblSum(rhi, rlo, __tmp_dbldbl_floor, __tmp2_dbldbl_floor); \
      } \
   } while(0)

/** compute the ceiled value of the floating point number given by two doubles */
#define SCIPdbldblCeil2(rhi, rlo, ahi, alo) \
   do { \
      double __tmp_dbldbl_ceil; \
      __tmp_dbldbl_ceil = ceil((ahi) + (alo)); \
      SCIPdbldblSum21(rhi, rlo, -(ahi), -(alo), __tmp_dbldbl_ceil); \
      if( ((rhi) - 1.0) + (rlo) < 0.0 && (rhi) + (rlo) >= 0.0 ) \
      { \
         /* ceil in double precision was fine */ \
         (rhi) = __tmp_dbldbl_ceil; \
         (rlo) = 0.0; \
      } \
      else \
      { \
         /* ceil in double precision needs to be corrected */ \
         double __tmp2_dbldbl_ceil = floor((rhi) + (rlo)); \
         SCIPdbldblSum(rhi, rlo, __tmp_dbldbl_ceil, -__tmp2_dbldbl_ceil); \
      } \
   } while(0)

/** compute the floored value of the floating point number given by two doubles, add epsilon first for safety */
#define SCIPdbldblEpsFloor2(rhi, rlo, ahi, alo, eps) SCIPdbldblFloor2(rhi, rlo, ahi, (alo) + (eps))

/** compute the ceiled value of the floating point number given by two doubles, subtract epsilon first for safety */
#define SCIPdbldblEpsCeil2(rhi, rlo, ahi, alo, eps) SCIPdbldblCeil2(rhi, rlo, ahi, (alo) - (eps))

#endif