4
* Copyright 2009 Tomi Jylhä-Ollila
6
* This file is part of Kunquat.
8
* Kunquat is free software: you can redistribute it and/or modify
9
* it under the terms of the GNU General Public License as published by
10
* the Free Software Foundation, either version 3 of the License, or
11
* (at your option) any later version.
13
* Kunquat is distributed in the hope that it will be useful,
14
* but WITHOUT ANY WARRANTY; without even the implied warranty of
15
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16
* GNU General Public License for more details.
18
* You should have received a copy of the GNU General Public License
19
* along with Kunquat. If not, see <http://www.gnu.org/licenses/>.
32
* Reduces a Real object into lowest terms.
34
* Does nothing unless the Real object is a fraction.
36
* \param real The Real object -- must not be \c NULL.
38
* \return The parameter \a real.
40
static Real* Real_normalise(Real* real);
45
* Validates a Real object.
47
* \param real The Real object.
49
* \return \c true if and only if the Real object is valid.
51
static bool Real_validate_(Real* real);
52
#define Real_validate(real) ( assert(Real_validate_(real)) )
54
#define Real_validate(real) ((void)0)
59
* Performs integer multiplication with overflow check.
61
* \param a The first term.
62
* \param b The second term.
63
* \param res Pointer to the result variable -- must not be \c NULL. If an
64
* overflow would occur, the stored value will not be changed.
66
* \return \c true if \a a * \a b was calculated successfully, or \c false
67
* if the result would become larger than INT64_MAX in magnitude.
69
static bool safe_mul(int64_t a, int64_t b, int64_t* res);
72
Real* Real_init(Real* real)
75
return Real_init_as_frac(real, 1, 1);
79
Real* Real_init_as_frac(Real* real, int64_t numerator, int64_t denominator)
82
assert(denominator > 0);
84
real->fod.frac.numerator = numerator;
85
real->fod.frac.denominator = denominator;
86
return Real_normalise(real);
90
Real* Real_init_as_double(Real* real, double val)
99
bool Real_is_frac(Real* real)
102
return real->is_frac;
106
int64_t Real_get_numerator(Real* real)
111
return (int64_t)real->fod.doub;
113
return real->fod.frac.numerator;
117
int64_t Real_get_denominator(Real* real)
124
return real->fod.frac.denominator;
128
double Real_get_double(Real* real)
133
return (double)real->fod.frac.numerator
134
/ (double)real->fod.frac.denominator;
136
return real->fod.doub;
140
Real* Real_copy(Real* dest, Real* src)
142
assert(dest != NULL);
147
dest->fod.frac.numerator = src->fod.frac.numerator;
148
dest->fod.frac.denominator = src->fod.frac.denominator;
153
dest->fod.doub = src->fod.doub;
159
Real* Real_mul(Real* ret, Real* real1, Real* real2)
162
Real_validate(real1);
163
Real_validate(real2);
164
if (real1->is_frac && real2->is_frac)
168
if (safe_mul(real1->fod.frac.numerator,
169
real2->fod.frac.numerator, &num)
170
&& safe_mul(real1->fod.frac.denominator,
171
real2->fod.frac.denominator, &den))
174
ret->fod.frac.numerator = num;
175
ret->fod.frac.denominator = den;
176
return Real_normalise(ret);
179
return Real_init_as_double(ret, Real_get_double(real1) * Real_get_double(real2));
183
Real* Real_div(Real* ret, Real* dividend, Real* divisor)
186
Real_validate(dividend);
187
Real_validate(divisor);
188
if (dividend->is_frac && divisor->is_frac)
190
assert(divisor->fod.frac.numerator != 0);
193
if (safe_mul(dividend->fod.frac.numerator,
194
divisor->fod.frac.denominator, &num)
195
&& safe_mul(dividend->fod.frac.denominator,
196
divisor->fod.frac.numerator, &den))
204
ret->fod.frac.numerator = num;
205
ret->fod.frac.denominator = den;
206
return Real_normalise(ret);
209
assert(Real_get_double(divisor) != 0);
210
return Real_init_as_double(ret, Real_get_double(dividend) / Real_get_double(divisor));
214
double Real_mul_float(Real* real, double d)
217
return Real_get_double(real) * d;
221
int Real_cmp(Real* real1, Real* real2)
223
Real_validate(real1);
224
Real_validate(real2);
225
if (real1->is_frac && real2->is_frac)
227
int64_t num1 = real1->fod.frac.numerator;
228
int64_t den1 = real1->fod.frac.denominator;
229
int64_t num2 = real2->fod.frac.numerator;
230
int64_t den2 = real2->fod.frac.denominator;
231
if (num1 == num2 && den1 == den2)
237
if (safe_mul(num1, den2, &term1)
238
&& safe_mul(num2, den1, &term2))
242
else if (term1 > term2)
247
double val1 = Real_get_double(real1);
248
double val2 = Real_get_double(real2);
253
else if (val1 > val2)
261
Real* Real_normalise(Real* real)
268
if (real->fod.frac.numerator == INT64_MIN)
270
// get rid of INT64_MIN for imaxabs()
271
if ((INT64_MIN % 2) == -1)
273
// Potentially suboptimal: real may not be in lowest terms here
274
// if INT_LEAST64_MIN is used.
277
if ((real->fod.frac.denominator & 1) == 0)
279
real->fod.frac.numerator = INT64_MIN / 2;
280
real->fod.frac.denominator >>= 1;
284
// INT64_MIN is has the form -2^n,
285
// so the fraction is irreducible here.
286
// Potentially suboptimal if INT_LEAST64_MIN is used.
291
int64_t num = imaxabs(real->fod.frac.numerator);
292
int64_t den = real->fod.frac.denominator;
295
real->fod.frac.denominator = 1;
298
while ((num & 1) == 0 && (den & 1) == 0)
308
else if ((num & 1) == 0)
311
den = (den - num) >> 1;
313
num = (num - den) >> 1;
317
real->fod.frac.numerator /= num;
318
real->fod.frac.denominator /= num;
325
bool Real_validate_(Real* real)
333
if (real->fod.frac.denominator <= 0)
340
if (isnan(real->fod.doub))
350
// Works the way described in http://www.fefe.de/intof.html
351
static bool safe_mul(int64_t a, int64_t b, int64_t* res)
354
bool negative = false;
355
if ((a == INT64_MIN || b == INT64_MIN)
357
&& INT64_MIN < -INT64_MAX)
364
negative = !negative;
369
negative = !negative;
371
uint64_t ahi = (uint64_t)a >> 32;
372
uint64_t bhi = (uint64_t)b >> 32;
373
if (ahi != 0 && bhi != 0)
377
uint64_t alo = a & 0xffffffffULL;
378
uint64_t blo = b & 0xffffffffULL;
379
uint64_t hi = (ahi * blo) + (alo * bhi);
380
if (hi > 0x7fffffffULL)
385
uint64_t lo = alo * blo;
386
if (UINT64_MAX - lo < hi || INT64_MAX < hi + lo)
added
removed
src/lib/Real.c
