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
toRational(sal_Int32 n,sal_Int32 d)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.
Fraction(sal_Int64 nNum,sal_Int64 nDen)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 */
Fraction(double nNum,double nDen)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
Fraction(double dVal)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
operator double() const99 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
operator +=(const Fraction & rVal)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
operator -=(const Fraction & rVal)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 {
checked_multiply_by(boost::rational<T> & i,const boost::rational<T> & r)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
operator *=(const Fraction & rVal)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
operator /=(const Fraction & rVal)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 */
ReduceInaccurate(unsigned nSignificantBits)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
GetNumerator() const263 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
GetDenominator() const273 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
operator sal_Int32() const283 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
operator +(const Fraction & rVal1,const Fraction & rVal2)293 Fraction operator+( const Fraction& rVal1, const Fraction& rVal2 )
294 {
295 Fraction aErg( rVal1 );
296 aErg += rVal2;
297 return aErg;
298 }
299
operator -(const Fraction & rVal1,const Fraction & rVal2)300 Fraction operator-( const Fraction& rVal1, const Fraction& rVal2 )
301 {
302 Fraction aErg( rVal1 );
303 aErg -= rVal2;
304 return aErg;
305 }
306
operator *(const Fraction & rVal1,const Fraction & rVal2)307 Fraction operator*( const Fraction& rVal1, const Fraction& rVal2 )
308 {
309 Fraction aErg( rVal1 );
310 aErg *= rVal2;
311 return aErg;
312 }
313
operator /(const Fraction & rVal1,const Fraction & rVal2)314 Fraction operator/( const Fraction& rVal1, const Fraction& rVal2 )
315 {
316 Fraction aErg( rVal1 );
317 aErg /= rVal2;
318 return aErg;
319 }
320
operator !=(const Fraction & rVal1,const Fraction & rVal2)321 bool operator !=( const Fraction& rVal1, const Fraction& rVal2 )
322 {
323 return !(rVal1 == rVal2);
324 }
325
operator <=(const Fraction & rVal1,const Fraction & rVal2)326 bool operator <=( const Fraction& rVal1, const Fraction& rVal2 )
327 {
328 return !(rVal1 > rVal2);
329 }
330
operator >=(const Fraction & rVal1,const Fraction & rVal2)331 bool operator >=( const Fraction& rVal1, const Fraction& rVal2 )
332 {
333 return !(rVal1 < rVal2);
334 }
335
operator ==(const Fraction & rVal1,const Fraction & rVal2)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
operator <(const Fraction & rVal1,const Fraction & rVal2)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
operator >(const Fraction & rVal1,const Fraction & rVal2)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
ReadFraction(SvStream & rIStream,Fraction & rFract)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
WriteFraction(SvStream & rOStream,const Fraction & rFract)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.
rational_FromDouble(double dVal)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
impl_NumberOfBits(sal_uInt32 nNum)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 */
rational_ReduceInaccurate(boost::rational<sal_Int32> & rRational,unsigned nSignificantBits)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: */
520