tools/source/generic/fract.cxx

   1 /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
   2 /*
   3  * This file is part of the LibreOffice project.
   4  *
   5  * This Source Code Form is subject to the terms of the Mozilla Public
   6  * License, v. 2.0. If a copy of the MPL was not distributed with this
   7  * file, You can obtain one at http://mozilla.org/MPL/2.0/.
   8  *
   9  * This file incorporates work covered by the following license notice:
  10  *
  11  *   Licensed to the Apache Software Foundation (ASF) under one or more
  12  *   contributor license agreements. See the NOTICE file distributed
  13  *   with this work for additional information regarding copyright
  14  *   ownership. The ASF licenses this file to you under the Apache
  15  *   License, Version 2.0 (the "License"); you may not use this file
  16  *   except in compliance with the License. You may obtain a copy of
  17  *   the License at http://www.apache.org/licenses/LICENSE-2.0 .
  18  */
  19
  20 #include <tools/fract.hxx>
  21 #include <tools/debug.hxx>
  22 #include <o3tl/hash_combine.hxx>
  23 #include <o3tl/safeint.hxx>
  24 #include <sal/log.hxx>
  25 #include <osl/diagnose.h>
  26
  27 #include <algorithm>
  28 #include <cmath>
  29 #include <numeric>
  30
  31 #include <boost/rational.hpp>
  32
  33 #ifdef _MSC_VER
  34 #include <intrin.h>
  35 #endif
  36
  37 static boost::rational<sal_Int32> rational_FromDouble(double dVal);
  38 static void rational_ReduceInaccurate(boost::rational<sal_Int32>& rRational, unsigned nSignificantBits);
  39 static int impl_NumberOfBits( sal_uInt32 nNum );
  40
  41 static boost::rational<sal_Int32> toRational(sal_Int32 n, sal_Int32 d)
  42 {
  43     // https://github.com/boostorg/boost/issues/335 when these are std::numeric_limits<sal_Int32>::min
  44     if (n == d)
  45         return 1;
  46     // tdf#144319 avoid boost::bad_rational e.g. if numerator=-476741369, denominator=-2147483648
  47     if (d < -std::numeric_limits<sal_Int32>::max())
  48         return 0;
  49     return boost::rational<sal_Int32>(n, d);
  50 }
  51
  52 static constexpr bool isOutOfRange(sal_Int64 nNum)
  53 {
  54     return nNum < std::numeric_limits<sal_Int32>::min()
  55             || nNum > std::numeric_limits<sal_Int32>::max();
  56 }
  57
  58 Fraction::Fraction( sal_Int64 nNum, sal_Int64 nDen ) : mnNumerator(nNum), mnDenominator(nDen)
  59 {
  60     if ( isOutOfRange(nNum) || isOutOfRange(nDen) )
  61     {
  62         // tdf#143200
  63         if (const auto gcd = std::gcd(nNum, nDen); gcd > 1)
  64         {
  65             nNum /= gcd;
  66             nDen /= gcd;
  67         }
  68         SAL_WARN_IF(isOutOfRange(nNum) || isOutOfRange(nDen),
  69                     "tools.fraction", "values outside of range we can represent, doing reduction, which will reduce precision");
  70         while (isOutOfRange(nNum) || isOutOfRange(nDen))
  71         {
  72             nNum /= 2;
  73             nDen /= 2;
  74         }
  75         mnNumerator = nNum;
  76         mnDenominator = nDen;
  77     }
  78     if ( mnDenominator == 0 )
  79     {
  80         mbValid = false;
  81         SAL_WARN( "tools.fraction", "'Fraction(" << nNum << ",0)' invalid fraction created" );
  82         return;
  83     }
  84     else if ((nDen == -1 && nNum == std::numeric_limits<sal_Int32>::min()) ||
  85             (nNum == -1 && nDen == std::numeric_limits<sal_Int32>::min()))
  86     {
  87         mbValid = false;
  88         SAL_WARN("tools.fraction", "'Fraction(" << nNum << "," << nDen << ")' invalid fraction created");
  89         return;
  90     }
  91 }
  92
  93 /**
  94  * only here to prevent passing of NaN
  95  */
  96 Fraction::Fraction( double nNum, double nDen )
  97     : Fraction(sal_Int64(nNum), sal_Int64(nDen))
  98 {}
  99
 100 Fraction::Fraction( double dVal )
 101 {
 102     try
 103     {
 104         boost::rational<sal_Int32> v = rational_FromDouble( dVal );
 105         mnNumerator = v.numerator();
 106         mnDenominator = v.denominator();
 107     }
 108     catch (const boost::bad_rational&)
 109     {
 110         mbValid = false;
 111         SAL_WARN( "tools.fraction", "'Fraction(" << dVal << ")' invalid fraction created" );
 112     }
 113 }
 114
 115 Fraction::operator double() const
 116 {
 117     if (!mbValid)
 118     {
 119         SAL_WARN( "tools.fraction", "'double()' on invalid fraction" );
 120         return 0.0;
 121     }
 122
 123     return boost::rational_cast<double>(toRational(mnNumerator, mnDenominator));
 124 }
 125
 126 // This methods first validates both values.
 127 // If one of the arguments is invalid, the whole operation is invalid.
 128 // After computation detect if result overflows a sal_Int32 value
 129 // which cause the operation to be marked as invalid
 130 Fraction& Fraction::operator += ( const Fraction& rVal )
 131 {
 132     if ( !rVal.mbValid )
 133         mbValid = false;
 134
 135     if ( !mbValid )
 136     {
 137         SAL_WARN( "tools.fraction", "'operator +=' with invalid fraction" );
 138         return *this;
 139     }
 140
 141     boost::rational<sal_Int32> a = toRational(mnNumerator, mnDenominator);
 142     a += toRational(rVal.mnNumerator, rVal.mnDenominator);
 143     mnNumerator = a.numerator();
 144     mnDenominator = a.denominator();
 145
 146     return *this;
 147 }
 148
 149 Fraction& Fraction::operator -= ( const Fraction& rVal )
 150 {
 151     if ( !rVal.mbValid )
 152         mbValid = false;
 153
 154     if ( !mbValid )
 155     {
 156         SAL_WARN( "tools.fraction", "'operator -=' with invalid fraction" );
 157         return *this;
 158     }
 159
 160     boost::rational<sal_Int32> a = toRational(mnNumerator, mnDenominator);
 161     a -= toRational(rVal.mnNumerator, rVal.mnDenominator);
 162     mnNumerator = a.numerator();
 163     mnDenominator = a.denominator();
 164
 165     return *this;
 166 }
 167
 168 namespace
 169 {
 170     bool checked_multiply_by(boost::rational<sal_Int32>& i, const boost::rational<sal_Int32>& r)
 171     {
 172         // Protect against self-modification
 173         sal_Int32 num = r.numerator();
 174         sal_Int32 den = r.denominator();
 175
 176         // Fast-path if the number of bits in input is < the number of bits in the output, overflow cannot happen
 177         // This is considerably faster than repeated std::gcd() operations
 178         if ((impl_NumberOfBits(std::abs(i.numerator())) + impl_NumberOfBits(std::abs(r.numerator()))) < 32 &&
 179             (impl_NumberOfBits(std::abs(i.denominator())) + impl_NumberOfBits(std::abs(r.denominator()))) < 32)
 180         {
 181             i *= r;
 182             return false;
 183         }
 184
 185         // Avoid overflow and preserve normalization
 186         sal_Int32 gcd1 = std::gcd(i.numerator(), den);
 187         sal_Int32 gcd2 = std::gcd(num, i.denominator());
 188
 189         if (!gcd1 || !gcd2)
 190             return true;
 191
 192         bool fail = false;
 193         fail |= o3tl::checked_multiply(i.numerator() / gcd1, num / gcd2, num);
 194         fail |= o3tl::checked_multiply(i.denominator() / gcd2, den / gcd1, den);
 195
 196         if (!fail)
 197             i.assign(num, den);
 198
 199         return fail;
 200     }
 201 }
 202
 203 Fraction& Fraction::operator *= ( const Fraction& rVal )
 204 {
 205     if ( !rVal.mbValid )
 206         mbValid = false;
 207
 208     if ( !mbValid )
 209     {
 210         SAL_WARN( "tools.fraction", "'operator *=' with invalid fraction" );
 211         return *this;
 212     }
 213
 214     boost::rational<sal_Int32> a = toRational(mnNumerator, mnDenominator);
 215     boost::rational<sal_Int32> b = toRational(rVal.mnNumerator, rVal.mnDenominator);
 216     bool bFail = checked_multiply_by(a, b);
 217     mnNumerator = a.numerator();
 218     mnDenominator = a.denominator();
 219
 220     if (bFail)
 221     {
 222         mbValid = false;
 223     }
 224
 225     return *this;
 226 }
 227
 228 Fraction& Fraction::operator /= ( const Fraction& rVal )
 229 {
 230     if ( !rVal.mbValid )
 231         mbValid = false;
 232
 233     if ( !mbValid )
 234     {
 235         SAL_WARN( "tools.fraction", "'operator /=' with invalid fraction" );
 236         return *this;
 237     }
 238
 239     boost::rational<sal_Int32> a = toRational(mnNumerator, mnDenominator);
 240     a /= toRational(rVal.mnNumerator, rVal.mnDenominator);
 241     mnNumerator = a.numerator();
 242     mnDenominator = a.denominator();
 243
 244     return *this;
 245 }
 246
 247 /** Inaccurate cancellation for a fraction.
 248
 249     Clip both nominator and denominator to said number of bits. If
 250     either of those already have equal or less number of bits used,
 251     this method does nothing.
 252
 253     @param nSignificantBits denotes, how many significant binary
 254     digits to maintain, in both nominator and denominator.
 255
 256     @example ReduceInaccurate(8) has an error <1% [1/2^(8-1)] - the
 257     largest error occurs with the following pair of values:
 258
 259     binary    1000000011111111111111111111111b/1000000000000000000000000000000b
 260     =         1082130431/1073741824
 261     = approx. 1.007812499
 262
 263     A ReduceInaccurate(8) yields 1/1.
 264 */
 265 void Fraction::ReduceInaccurate( unsigned nSignificantBits )
 266 {
 267     if ( !mbValid )
 268     {
 269         SAL_WARN( "tools.fraction", "'ReduceInaccurate' on invalid fraction" );
 270         return;
 271     }
 272
 273     if ( !mnNumerator )
 274         return;
 275
 276     auto a = toRational(mnNumerator, mnDenominator);
 277     rational_ReduceInaccurate(a, nSignificantBits);
 278     mnNumerator = a.numerator();
 279     mnDenominator = a.denominator();
 280 }
 281
 282 sal_Int32 Fraction::GetNumerator() const
 283 {
 284     if ( !mbValid )
 285     {
 286         SAL_WARN( "tools.fraction", "'GetNumerator()' on invalid fraction" );
 287         return 0;
 288     }
 289     return mnNumerator;
 290 }
 291
 292 sal_Int32 Fraction::GetDenominator() const
 293 {
 294     if ( !mbValid )
 295     {
 296         SAL_WARN( "tools.fraction", "'GetDenominator()' on invalid fraction" );
 297         return -1;
 298     }
 299     return mnDenominator;
 300 }
 301
 302 Fraction::operator sal_Int32() const
 303 {
 304     if ( !mbValid )
 305     {
 306         SAL_WARN( "tools.fraction", "'operator sal_Int32()' on invalid fraction" );
 307         return 0;
 308     }
 309     return boost::rational_cast<sal_Int32>(toRational(mnNumerator, mnDenominator));
 310 }
 311
 312 Fraction operator+( const Fraction& rVal1, const Fraction& rVal2 )
 313 {
 314     Fraction aErg( rVal1 );
 315     aErg += rVal2;
 316     return aErg;
 317 }
 318
 319 Fraction operator-( const Fraction& rVal1, const Fraction& rVal2 )
 320 {
 321     Fraction aErg( rVal1 );
 322     aErg -= rVal2;
 323     return aErg;
 324 }
 325
 326 Fraction operator*( const Fraction& rVal1, const Fraction& rVal2 )
 327 {
 328     Fraction aErg( rVal1 );
 329     aErg *= rVal2;
 330     return aErg;
 331 }
 332
 333 Fraction operator/( const Fraction& rVal1, const Fraction& rVal2 )
 334 {
 335     Fraction aErg( rVal1 );
 336     aErg /= rVal2;
 337     return aErg;
 338 }
 339
 340 bool operator !=( const Fraction& rVal1, const Fraction& rVal2 )
 341 {
 342     return !(rVal1 == rVal2);
 343 }
 344
 345 bool operator <=( const Fraction& rVal1, const Fraction& rVal2 )
 346 {
 347     return !(rVal1 > rVal2);
 348 }
 349
 350 bool operator >=( const Fraction& rVal1, const Fraction& rVal2 )
 351 {
 352     return !(rVal1 < rVal2);
 353 }
 354
 355 bool operator == ( const Fraction& rVal1, const Fraction& rVal2 )
 356 {
 357     if ( !rVal1.mbValid || !rVal2.mbValid )
 358     {
 359         SAL_WARN( "tools.fraction", "'operator ==' with an invalid fraction" );
 360         return false;
 361     }
 362
 363     return toRational(rVal1.mnNumerator, rVal1.mnDenominator) == toRational(rVal2.mnNumerator, rVal2.mnDenominator);
 364 }
 365
 366 bool operator < ( const Fraction& rVal1, const Fraction& rVal2 )
 367 {
 368     if ( !rVal1.mbValid || !rVal2.mbValid )
 369     {
 370         SAL_WARN( "tools.fraction", "'operator <' with an invalid fraction" );
 371         return false;
 372     }
 373
 374     return toRational(rVal1.mnNumerator, rVal1.mnDenominator) < toRational(rVal2.mnNumerator, rVal2.mnDenominator);
 375 }
 376
 377 bool operator > ( const Fraction& rVal1, const Fraction& rVal2 )
 378 {
 379     if ( !rVal1.mbValid || !rVal2.mbValid )
 380     {
 381         SAL_WARN( "tools.fraction", "'operator >' with an invalid fraction" );
 382         return false;
 383     }
 384
 385     return toRational(rVal1.mnNumerator, rVal1.mnDenominator) > toRational(rVal2.mnNumerator, rVal2.mnDenominator);
 386 }
 387
 388 // If dVal > LONG_MAX or dVal < LONG_MIN, the rational throws a boost::bad_rational.
 389 // Otherwise, dVal and denominator are multiplied by 10, until one of them
 390 // is larger than (LONG_MAX / 10).
 391 //
 392 // NOTE: here we use 'sal_Int32' due that only values in sal_Int32 range are valid.
 393 static boost::rational<sal_Int32> rational_FromDouble(double dVal)
 394 {
 395     if ( dVal > std::numeric_limits<sal_Int32>::max() ||
 396          dVal < std::numeric_limits<sal_Int32>::min() ||
 397          std::isnan(dVal) )
 398         throw boost::bad_rational();
 399
 400     const sal_Int32 nMAX = std::numeric_limits<sal_Int32>::max() / 10;
 401     sal_Int32 nDen = 1;
 402     while ( std::abs( dVal ) < nMAX && nDen < nMAX ) {
 403         dVal *= 10;
 404         nDen *= 10;
 405     }
 406     return boost::rational<sal_Int32>( sal_Int32(dVal), nDen );
 407 }
 408
 409 /**
 410  * Find the number of bits required to represent this number, using the CLZ intrinsic
 411  */
 412 static int impl_NumberOfBits( sal_uInt32 nNum )
 413 {
 414     if (nNum == 0)
 415         return 0;
 416 #ifdef _MSC_VER
 417     unsigned long r = 0;
 418     _BitScanReverse(&r, nNum);
 419     return r + 1;
 420 #else
 421     return 32 - __builtin_clz(nNum);
 422 #endif
 423 }
 424
 425 /** Inaccurate cancellation for a fraction.
 426
 427     Clip both nominator and denominator to said number of bits. If
 428     either of those already have equal or less number of bits used,
 429     this method does nothing.
 430
 431     @param nSignificantBits denotes, how many significant binary
 432     digits to maintain, in both nominator and denominator.
 433
 434     @example ReduceInaccurate(8) has an error <1% [1/2^(8-1)] - the
 435     largest error occurs with the following pair of values:
 436
 437     binary    1000000011111111111111111111111b/1000000000000000000000000000000b
 438     =         1082130431/1073741824
 439     = approx. 1.007812499
 440
 441     A ReduceInaccurate(8) yields 1/1.
 442 */
 443 static void rational_ReduceInaccurate(boost::rational<sal_Int32>& rRational, unsigned nSignificantBits)
 444 {
 445     if ( !rRational )
 446         return;
 447
 448     // http://www.boost.org/doc/libs/release/libs/rational/rational.html#Internal%20representation
 449     sal_Int32 nMul = rRational.numerator();
 450     if (nMul == std::numeric_limits<sal_Int32>::min())
 451     {
 452         // ofz#32973 Integer-overflow
 453         return;
 454     }
 455     const bool bNeg = nMul < 0;
 456     if (bNeg)
 457         nMul = -nMul;
 458     sal_Int32 nDiv = rRational.denominator();
 459
 460     DBG_ASSERT(nSignificantBits<65, "More than 64 bit of significance is overkill!");
 461
 462     // How much bits can we lose?
 463     const int nMulBitsToLose = std::max( ( impl_NumberOfBits( nMul ) - int( nSignificantBits ) ), 0 );
 464     const int nDivBitsToLose = std::max( ( impl_NumberOfBits( nDiv ) - int( nSignificantBits ) ), 0 );
 465
 466     const int nToLose = std::min( nMulBitsToLose, nDivBitsToLose );
 467
 468     // Remove the bits
 469     nMul >>= nToLose;
 470     nDiv >>= nToLose;
 471
 472     if ( !nMul || !nDiv ) {
 473         // Return without reduction
 474         OSL_FAIL( "Oops, we reduced too much..." );
 475         return;
 476     }
 477
 478     rRational.assign( bNeg ? -nMul : nMul, nDiv );
 479 }
 480
 481 size_t Fraction::GetHashValue() const
 482 {
 483     size_t hash = 0;
 484     o3tl::hash_combine( hash, mnNumerator );
 485     o3tl::hash_combine( hash, mnDenominator );
 486     o3tl::hash_combine( hash, mbValid );
 487     return hash;
 488 }
 489
 490 Fraction Fraction::MakeFraction( tools::Long nN1, tools::Long nN2, tools::Long nD1, tools::Long nD2 )
 491 {
 492     if( nD1 == 0 || nD2 == 0 ) //under these bad circumstances the following while loop will be endless
 493     {
 494         SAL_WARN("tools.fraction", "Invalid parameter for ImplMakeFraction");
 495         return Fraction( 1, 1 );
 496     }
 497
 498     tools::Long i = 1;
 499
 500     if ( nN1 < 0 ) { i = -i; nN1 = -nN1; }
 501     if ( nN2 < 0 ) { i = -i; nN2 = -nN2; }
 502     if ( nD1 < 0 ) { i = -i; nD1 = -nD1; }
 503     if ( nD2 < 0 ) { i = -i; nD2 = -nD2; }
 504     // all positive; i sign
 505
 506     assert( nN1 >= std::numeric_limits<sal_Int32>::min() );
 507     assert( nN1 <= std::numeric_limits<sal_Int32>::max( ));
 508     assert( nD1 >= std::numeric_limits<sal_Int32>::min() );
 509     assert( nD1 <= std::numeric_limits<sal_Int32>::max( ));
 510     assert( nN2 >= std::numeric_limits<sal_Int32>::min() );
 511     assert( nN2 <= std::numeric_limits<sal_Int32>::max( ));
 512     assert( nD2 >= std::numeric_limits<sal_Int32>::min() );
 513     assert( nD2 <= std::numeric_limits<sal_Int32>::max( ));
 514
 515     boost::rational<sal_Int32> a = toRational(i*nN1, nD1);
 516     boost::rational<sal_Int32> b = toRational(nN2, nD2);
 517     bool bFail = checked_multiply_by(a, b);
 518
 519     while ( bFail ) {
 520         if ( nN1 > nN2 )
 521             nN1 = (nN1 + 1) / 2;
 522         else
 523             nN2 = (nN2 + 1) / 2;
 524         if ( nD1 > nD2 )
 525             nD1 = (nD1 + 1) / 2;
 526         else
 527             nD2 = (nD2 + 1) / 2;
 528
 529         a = toRational(i*nN1, nD1);
 530         b = toRational(nN2, nD2);
 531         bFail = checked_multiply_by(a, b);
 532     }
 533
 534     return Fraction(a.numerator(), a.denominator());
 535 }
 536
 537 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */