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 <tools/stream.hxx>
  23 #include <o3tl/safeint.hxx>
  24 #include <sal/log.hxx>
  25 #include <osl/diagnose.h>
  26
  27 #include <algorithm>
  28 #include <cmath>
  29
  30 #include <boost/version.hpp>
  31 #if BOOST_VERSION >= 106700
  32 #include <boost/integer/common_factor_rt.hpp>
  33 #else
  34 #include <boost/math/common_factor_rt.hpp>
  35 #endif
  36 #include <boost/rational.hpp>
  37
  38 static boost::rational<sal_Int32> rational_FromDouble(double dVal);
  39
  40 static void rational_ReduceInaccurate(boost::rational<sal_Int32>& rRational, unsigned nSignificantBits);
  41
  42 static boost::rational<sal_Int32> toRational(sal_Int32 n, sal_Int32 d)
  43 {
  44     return boost::rational<sal_Int32>(n, d);
  45 }
  46
  47 // Initialized by setting nNum as nominator and nDen as denominator
  48 // Negative values in the denominator are invalid and cause the
  49 // inversion of both nominator and denominator signs
  50 // in order to return the correct value.
  51 Fraction::Fraction( sal_Int64 nNum, sal_Int64 nDen ) : mnNumerator(nNum), mnDenominator(nDen)
  52 {
  53     assert( nNum >= std::numeric_limits<sal_Int32>::min() );
  54     assert( nNum <= std::numeric_limits<sal_Int32>::max( ));
  55     assert( nDen >= std::numeric_limits<sal_Int32>::min() );
  56     assert( nDen <= std::numeric_limits<sal_Int32>::max( ));
  57     if ( nDen == 0 )
  58     {
  59         mbValid = false;
  60         SAL_WARN( "tools.fraction", "'Fraction(" << nNum << ",0)' invalid fraction created" );
  61         return;
  62     }
  63 }
  64
  65 /**
  66  * only here to prevent passing of NaN
  67  */
  68 Fraction::Fraction( double nNum, double nDen ) : mnNumerator(sal_Int64(nNum)), mnDenominator(sal_Int64(nDen))
  69 {
  70     assert( !std::isnan(nNum) );
  71     assert( !std::isnan(nDen) );
  72     assert( nNum >= std::numeric_limits<sal_Int32>::min() );
  73     assert( nNum <= std::numeric_limits<sal_Int32>::max( ));
  74     assert( nDen >= std::numeric_limits<sal_Int32>::min() );
  75     assert( nDen <= std::numeric_limits<sal_Int32>::max( ));
  76     if ( nDen == 0 )
  77     {
  78         mbValid = false;
  79         SAL_WARN( "tools.fraction", "'Fraction(" << nNum << ",0)' invalid fraction created" );
  80         return;
  81     }
  82 }
  83
  84 Fraction::Fraction( double dVal )
  85 {
  86     try
  87     {
  88         boost::rational<sal_Int32> v = rational_FromDouble( dVal );
  89         mnNumerator = v.numerator();
  90         mnDenominator = v.denominator();
  91     }
  92     catch (const boost::bad_rational&)
  93     {
  94         mbValid = false;
  95         SAL_WARN( "tools.fraction", "'Fraction(" << dVal << ")' invalid fraction created" );
  96     }
  97 }
  98
  99 Fraction::operator double() const
 100 {
 101     if (!mbValid)
 102     {
 103         SAL_WARN( "tools.fraction", "'double()' on invalid fraction" );
 104         return 0.0;
 105     }
 106
 107     // https://github.com/boostorg/boost/issues/335 when these are std::numeric_limits<sal_Int32>::min
 108     if (mnNumerator == mnDenominator)
 109         return 1.0;
 110
 111     return boost::rational_cast<double>(toRational(mnNumerator, mnDenominator));
 112 }
 113
 114 // This methods first validates both values.
 115 // If one of the arguments is invalid, the whole operation is invalid.
 116 // After computation detect if result overflows a sal_Int32 value
 117 // which cause the operation to be marked as invalid
 118 Fraction& Fraction::operator += ( const Fraction& rVal )
 119 {
 120     if ( !rVal.mbValid )
 121         mbValid = false;
 122
 123     if ( !mbValid )
 124     {
 125         SAL_WARN( "tools.fraction", "'operator +=' with invalid fraction" );
 126         return *this;
 127     }
 128
 129     boost::rational<sal_Int32> a = toRational(mnNumerator, mnDenominator);
 130     a += toRational(rVal.mnNumerator, rVal.mnDenominator);
 131     mnNumerator = a.numerator();
 132     mnDenominator = a.denominator();
 133
 134     return *this;
 135 }
 136
 137 Fraction& Fraction::operator -= ( const Fraction& rVal )
 138 {
 139     if ( !rVal.mbValid )
 140         mbValid = false;
 141
 142     if ( !mbValid )
 143     {
 144         SAL_WARN( "tools.fraction", "'operator -=' with invalid fraction" );
 145         return *this;
 146     }
 147
 148     boost::rational<sal_Int32> a = toRational(mnNumerator, mnDenominator);
 149     a -= toRational(rVal.mnNumerator, rVal.mnDenominator);
 150     mnNumerator = a.numerator();
 151     mnDenominator = a.denominator();
 152
 153     return *this;
 154 }
 155
 156 namespace
 157 {
 158     template<typename T> bool checked_multiply_by(boost::rational<T>& i, const boost::rational<T>& r)
 159     {
 160         // Protect against self-modification
 161         T num = r.numerator();
 162         T den = r.denominator();
 163
 164         // Avoid overflow and preserve normalization
 165 #if BOOST_VERSION >= 106700
 166         T gcd1 = boost::integer::gcd(i.numerator(), den);
 167         T gcd2 = boost::integer::gcd(num, i.denominator());
 168 #else
 169         T gcd1 = boost::math::gcd(i.numerator(), den);
 170         T gcd2 = boost::math::gcd(num, i.denominator());
 171 #endif
 172
 173         bool fail = false;
 174         fail |= o3tl::checked_multiply(i.numerator() / gcd1, num / gcd2, num);
 175         fail |= o3tl::checked_multiply(i.denominator() / gcd2, den / gcd1, den);
 176
 177         if (!fail)
 178             i.assign(num, den);
 179
 180         return fail;
 181     }
 182 }
 183
 184 Fraction& Fraction::operator *= ( const Fraction& rVal )
 185 {
 186     if ( !rVal.mbValid )
 187         mbValid = false;
 188
 189     if ( !mbValid )
 190     {
 191         SAL_WARN( "tools.fraction", "'operator *=' with invalid fraction" );
 192         return *this;
 193     }
 194
 195     boost::rational<sal_Int32> a = toRational(mnNumerator, mnDenominator);
 196     boost::rational<sal_Int32> b = toRational(rVal.mnNumerator, rVal.mnDenominator);
 197     bool bFail = checked_multiply_by(a, b);
 198     mnNumerator = a.numerator();
 199     mnDenominator = a.denominator();
 200
 201     if (bFail)
 202     {
 203         mbValid = false;
 204     }
 205
 206     return *this;
 207 }
 208
 209 Fraction& Fraction::operator /= ( const Fraction& rVal )
 210 {
 211     if ( !rVal.mbValid )
 212         mbValid = false;
 213
 214     if ( !mbValid )
 215     {
 216         SAL_WARN( "tools.fraction", "'operator /=' with invalid fraction" );
 217         return *this;
 218     }
 219
 220     boost::rational<sal_Int32> a = toRational(mnNumerator, mnDenominator);
 221     a /= toRational(rVal.mnNumerator, rVal.mnDenominator);
 222     mnNumerator = a.numerator();
 223     mnDenominator = a.denominator();
 224
 225     return *this;
 226 }
 227
 228 /** Inaccurate cancellation for a fraction.
 229
 230     Clip both nominator and denominator to said number of bits. If
 231     either of those already have equal or less number of bits used,
 232     this method does nothing.
 233
 234     @param nSignificantBits denotes, how many significant binary
 235     digits to maintain, in both nominator and denominator.
 236
 237     @example ReduceInaccurate(8) has an error <1% [1/2^(8-1)] - the
 238     largest error occurs with the following pair of values:
 239
 240     binary    1000000011111111111111111111111b/1000000000000000000000000000000b
 241     =         1082130431/1073741824
 242     = approx. 1.007812499
 243
 244     A ReduceInaccurate(8) yields 1/1.
 245 */
 246 void Fraction::ReduceInaccurate( unsigned nSignificantBits )
 247 {
 248     if ( !mbValid )
 249     {
 250         SAL_WARN( "tools.fraction", "'ReduceInaccurate' on invalid fraction" );
 251         return;
 252     }
 253
 254     if ( !mnNumerator )
 255         return;
 256
 257     auto a = toRational(mnNumerator, mnDenominator);
 258     rational_ReduceInaccurate(a, nSignificantBits);
 259     mnNumerator = a.numerator();
 260     mnDenominator = a.denominator();
 261 }
 262
 263 sal_Int32 Fraction::GetNumerator() const
 264 {
 265     if ( !mbValid )
 266     {
 267         SAL_WARN( "tools.fraction", "'GetNumerator()' on invalid fraction" );
 268         return 0;
 269     }
 270     return mnNumerator;
 271 }
 272
 273 sal_Int32 Fraction::GetDenominator() const
 274 {
 275     if ( !mbValid )
 276     {
 277         SAL_WARN( "tools.fraction", "'GetDenominator()' on invalid fraction" );
 278         return -1;
 279     }
 280     return mnDenominator;
 281 }
 282
 283 Fraction::operator sal_Int32() const
 284 {
 285     if ( !mbValid )
 286     {
 287         SAL_WARN( "tools.fraction", "'operator sal_Int32()' on invalid fraction" );
 288         return 0;
 289     }
 290     return boost::rational_cast<sal_Int32>(toRational(mnNumerator, mnDenominator));
 291 }
 292
 293 Fraction operator+( const Fraction& rVal1, const Fraction& rVal2 )
 294 {
 295     Fraction aErg( rVal1 );
 296     aErg += rVal2;
 297     return aErg;
 298 }
 299
 300 Fraction operator-( const Fraction& rVal1, const Fraction& rVal2 )
 301 {
 302     Fraction aErg( rVal1 );
 303     aErg -= rVal2;
 304     return aErg;
 305 }
 306
 307 Fraction operator*( const Fraction& rVal1, const Fraction& rVal2 )
 308 {
 309     Fraction aErg( rVal1 );
 310     aErg *= rVal2;
 311     return aErg;
 312 }
 313
 314 Fraction operator/( const Fraction& rVal1, const Fraction& rVal2 )
 315 {
 316     Fraction aErg( rVal1 );
 317     aErg /= rVal2;
 318     return aErg;
 319 }
 320
 321 bool operator !=( const Fraction& rVal1, const Fraction& rVal2 )
 322 {
 323     return !(rVal1 == rVal2);
 324 }
 325
 326 bool operator <=( const Fraction& rVal1, const Fraction& rVal2 )
 327 {
 328     return !(rVal1 > rVal2);
 329 }
 330
 331 bool operator >=( const Fraction& rVal1, const Fraction& rVal2 )
 332 {
 333     return !(rVal1 < rVal2);
 334 }
 335
 336 bool operator == ( const Fraction& rVal1, const Fraction& rVal2 )
 337 {
 338     if ( !rVal1.mbValid || !rVal2.mbValid )
 339     {
 340         SAL_WARN( "tools.fraction", "'operator ==' with an invalid fraction" );
 341         return false;
 342     }
 343
 344     return toRational(rVal1.mnNumerator, rVal1.mnDenominator) == toRational(rVal2.mnNumerator, rVal2.mnDenominator);
 345 }
 346
 347 bool operator < ( const Fraction& rVal1, const Fraction& rVal2 )
 348 {
 349     if ( !rVal1.mbValid || !rVal2.mbValid )
 350     {
 351         SAL_WARN( "tools.fraction", "'operator <' with an invalid fraction" );
 352         return false;
 353     }
 354
 355     return toRational(rVal1.mnNumerator, rVal1.mnDenominator) < toRational(rVal2.mnNumerator, rVal2.mnDenominator);
 356 }
 357
 358 bool operator > ( const Fraction& rVal1, const Fraction& rVal2 )
 359 {
 360     if ( !rVal1.mbValid || !rVal2.mbValid )
 361     {
 362         SAL_WARN( "tools.fraction", "'operator >' with an invalid fraction" );
 363         return false;
 364     }
 365
 366     return toRational(rVal1.mnNumerator, rVal1.mnDenominator) > toRational(rVal2.mnNumerator, rVal2.mnDenominator);
 367 }
 368
 369 SvStream& ReadFraction( SvStream& rIStream, Fraction & rFract )
 370 {
 371     sal_Int32 num(0), den(0);
 372     rIStream.ReadInt32( num );
 373     rIStream.ReadInt32( den );
 374     if ( den <= 0 )
 375     {
 376         SAL_WARN( "tools.fraction", "'ReadFraction()' read an invalid fraction" );
 377         rFract.mbValid = false;
 378     }
 379     else
 380     {
 381         rFract.mnNumerator = num;
 382         rFract.mnDenominator = den;
 383         rFract.mbValid = true;
 384     }
 385     return rIStream;
 386 }
 387
 388 SvStream& WriteFraction( SvStream& rOStream, const Fraction& rFract )
 389 {
 390     if ( !rFract.mbValid )
 391     {
 392         SAL_WARN( "tools.fraction", "'WriteFraction()' write an invalid fraction" );
 393         rOStream.WriteInt32( 0 );
 394         rOStream.WriteInt32( -1 );
 395     } else {
 396         rOStream.WriteInt32( rFract.mnNumerator );
 397         rOStream.WriteInt32( rFract.mnDenominator );
 398     }
 399     return rOStream;
 400 }
 401
 402 // If dVal > LONG_MAX or dVal < LONG_MIN, the rational throws a boost::bad_rational.
 403 // Otherwise, dVal and denominator are multiplied by 10, until one of them
 404 // is larger than (LONG_MAX / 10).
 405 //
 406 // NOTE: here we use 'sal_Int32' due that only values in sal_Int32 range are valid.
 407 static boost::rational<sal_Int32> rational_FromDouble(double dVal)
 408 {
 409     if ( dVal > std::numeric_limits<sal_Int32>::max() ||
 410          dVal < std::numeric_limits<sal_Int32>::min() ||
 411          std::isnan(dVal) )
 412         throw boost::bad_rational();
 413
 414     const sal_Int32 nMAX = std::numeric_limits<sal_Int32>::max() / 10;
 415     sal_Int32 nDen = 1;
 416     while ( std::abs( dVal ) < nMAX && nDen < nMAX ) {
 417         dVal *= 10;
 418         nDen *= 10;
 419     }
 420     return boost::rational<sal_Int32>( sal_Int32(dVal), nDen );
 421 }
 422
 423 // Similar to clz_table that can be googled
 424 const char nbits_table[32] =
 425 {
 426     32,  1, 23,  2, 29, 24, 14,  3,
 427     30, 27, 25, 18, 20, 15, 10,  4,
 428     31, 22, 28, 13, 26, 17, 19,  9,
 429     21, 12, 16,  8, 11,  7,  6,  5
 430 };
 431
 432 static int impl_NumberOfBits( sal_uInt32 nNum )
 433 {
 434     // http://en.wikipedia.org/wiki/De_Bruijn_sequence
 435     // background paper: Using de Bruijn Sequences to Index a 1 in a
 436     // Computer Word (1998) Charles E. Leiserson,
 437     // Harald Prokop, Keith H. Randall
 438     // (e.g. http://citeseer.ist.psu.edu/leiserson98using.html)
 439     const sal_uInt32 nDeBruijn = 0x7DCD629;
 440
 441     if ( nNum == 0 )
 442         return 0;
 443
 444     // Get it to form like 0000001111111111b
 445     nNum |= ( nNum >>  1 );
 446     nNum |= ( nNum >>  2 );
 447     nNum |= ( nNum >>  4 );
 448     nNum |= ( nNum >>  8 );
 449     nNum |= ( nNum >> 16 );
 450
 451     sal_uInt32 nNumber;
 452     int nBonus;
 453
 454     nNumber = nNum;
 455     nBonus = 0;
 456
 457     // De facto shift left of nDeBruijn using multiplication (nNumber
 458     // is all ones from topmost bit, thus nDeBruijn + (nDeBruijn *
 459     // nNumber) => nDeBruijn * (nNumber+1) clears all those bits to
 460     // zero, sets the next bit to one, and thus effectively shift-left
 461     // nDeBruijn by lg2(nNumber+1). This generates a distinct 5bit
 462     // sequence in the msb for each distinct position of the last
 463     // leading 0 bit - that's the property of a de Bruijn number.
 464     nNumber = nDeBruijn + ( nDeBruijn * nNumber );
 465
 466     // 5-bit window indexes the result
 467     return ( nbits_table[nNumber >> 27] ) + nBonus;
 468 }
 469
 470 /** Inaccurate cancellation for a fraction.
 471
 472     Clip both nominator and denominator to said number of bits. If
 473     either of those already have equal or less number of bits used,
 474     this method does nothing.
 475
 476     @param nSignificantBits denotes, how many significant binary
 477     digits to maintain, in both nominator and denominator.
 478
 479     @example ReduceInaccurate(8) has an error <1% [1/2^(8-1)] - the
 480     largest error occurs with the following pair of values:
 481
 482     binary    1000000011111111111111111111111b/1000000000000000000000000000000b
 483     =         1082130431/1073741824
 484     = approx. 1.007812499
 485
 486     A ReduceInaccurate(8) yields 1/1.
 487 */
 488 static void rational_ReduceInaccurate(boost::rational<sal_Int32>& rRational, unsigned nSignificantBits)
 489 {
 490     if ( !rRational )
 491         return;
 492
 493     // http://www.boost.org/doc/libs/release/libs/rational/rational.html#Internal%20representation
 494     const bool bNeg = ( rRational.numerator() < 0 );
 495     sal_Int32 nMul = bNeg? -rRational.numerator(): rRational.numerator();
 496     sal_Int32 nDiv = rRational.denominator();
 497
 498     DBG_ASSERT(nSignificantBits<65, "More than 64 bit of significance is overkill!");
 499
 500     // How much bits can we lose?
 501     const int nMulBitsToLose = std::max( ( impl_NumberOfBits( nMul ) - int( nSignificantBits ) ), 0 );
 502     const int nDivBitsToLose = std::max( ( impl_NumberOfBits( nDiv ) - int( nSignificantBits ) ), 0 );
 503
 504     const int nToLose = std::min( nMulBitsToLose, nDivBitsToLose );
 505
 506     // Remove the bits
 507     nMul >>= nToLose;
 508     nDiv >>= nToLose;
 509
 510     if ( !nMul || !nDiv ) {
 511         // Return without reduction
 512         OSL_FAIL( "Oops, we reduced too much..." );
 513         return;
 514     }
 515
 516     rRational.assign( bNeg ? -nMul : nMul, nDiv );
 517 }
 518
 519 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */