OMCompiler/SimulationRuntime/c/util/rational.c
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /* | ||
| 2 | * This file belongs to the OpenModelica Run-Time System | ||
| 3 | * | ||
| 4 | * Copyright (c) 1998-2026, Open Source Modelica Consortium (OSMC), c/o Linköpings | ||
| 5 | * universitet, Department of Computer and Information Science, SE-58183 Linköping, Sweden. All rights | ||
| 6 | * reserved. | ||
| 7 | * | ||
| 8 | * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF THE BSD NEW LICENSE OR THE | ||
| 9 | * AGPL VERSION 3 LICENSE OR THE OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.8. ANY | ||
| 10 | * USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES RECIPIENT'S | ||
| 11 | * ACCEPTANCE OF THE BSD NEW LICENSE OR THE OSMC PUBLIC LICENSE OR THE AGPL | ||
| 12 | * VERSION 3, ACCORDING TO RECIPIENTS CHOICE. | ||
| 13 | * | ||
| 14 | * The OpenModelica software and the OSMC (Open Source Modelica Consortium) Public License | ||
| 15 | * (OSMC-PL) are obtained from OSMC, either from the above address, from the URLs: | ||
| 16 | * http://www.openmodelica.org or https://github.com/OpenModelica/ or | ||
| 17 | * http://www.ida.liu.se/projects/OpenModelica, and in the OpenModelica distribution. GNU | ||
| 18 | * AGPL version 3 is obtained from: https://www.gnu.org/licenses/licenses.html#GPL. The BSD NEW | ||
| 19 | * License is obtained from: http://www.opensource.org/licenses/BSD-3-Clause. | ||
| 20 | * | ||
| 21 | * This program is distributed WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| 22 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY | ||
| 23 | * SET FORTH IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS OF | ||
| 24 | * OSMC-PL. | ||
| 25 | * | ||
| 26 | */ | ||
| 27 | |||
| 28 | #include "rational.h" | ||
| 29 | #include "omc_msvc.h" | ||
| 30 | #include "omc_error.h" | ||
| 31 | #include <assert.h> | ||
| 32 | #include <stdlib.h> | ||
| 33 | #include <limits.h> | ||
| 34 | |||
| 35 | |||
| 36 | /* | ||
| 37 | * Rational arithmetic is particularly prone to overflow because numerator | ||
| 38 | * and/or denominator can grow quickly during computations. Therefore overflow | ||
| 39 | * is checked during critical steps in the calculations below. | ||
| 40 | */ | ||
| 41 | #if defined __has_builtin | ||
| 42 | #if __has_builtin(__builtin_add_overflow) | ||
| 43 | #define RAT_INT_ADD(a, b, c, op) \ | ||
| 44 | assertStreamPrint(NULL, !__builtin_add_overflow((a), (b), &(c)), \ | ||
| 45 | "RATIONAL overflow. Unable to store result of " \ | ||
| 46 | "("RAT_FMT"/"RAT_FMT") %c ("RAT_FMT"/"RAT_FMT")", \ | ||
| 47 | r1.num, r1.den, op, r2.num, r2.den) | ||
| 48 | #define RAT_INT_MUL(a, b, c, op) \ | ||
| 49 | assertStreamPrint(NULL, !__builtin_mul_overflow((a), (b), &(c)), \ | ||
| 50 | "RATIONAL overflow. Unable to store result of " \ | ||
| 51 | "("RAT_FMT"/"RAT_FMT") %c ("RAT_FMT"/"RAT_FMT")", \ | ||
| 52 | r1.num, r1.den, op, r2.num, r2.den) | ||
| 53 | #endif | ||
| 54 | #endif | ||
| 55 | |||
| 56 | #if !(defined RAT_INT_ADD) /* no overflow checks available */ | ||
| 57 | #define RAT_INT_ADD(a, b, c, op) (c) = (a) + (b) | ||
| 58 | #define RAT_INT_MUL(a, b, c, op) (c) = (a) * (b) | ||
| 59 | #endif | ||
| 60 | |||
| 61 | |||
| 62 | /** | ||
| 63 | * @brief Greatest common divisor. | ||
| 64 | * | ||
| 65 | * Largest positive integer that divides a and b. | ||
| 66 | * gcd(a,b) | ||
| 67 | * | ||
| 68 | * @param a First integer a. | ||
| 69 | * @param b Second integer b. | ||
| 70 | * @return long long Greatest common divisor of a and b. | ||
| 71 | */ | ||
| 72 | static rat_int_t gcd(rat_int_t a, rat_int_t b) | ||
| 73 | { | ||
| 74 |
2/10✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 60 times.
✓ Branch 9 taken 92 times.
|
152 | while(a != 0) { |
| 75 | rat_int_t tmp = a; | ||
| 76 | 60 | a = b % a; | |
| 77 | b = tmp; | ||
| 78 | } | ||
| 79 | ✗ | return RAT_INT_ABS(b); | |
| 80 | } | ||
| 81 | |||
| 82 | |||
| 83 | /** | ||
| 84 | * @brief Simplify rational number. | ||
| 85 | * | ||
| 86 | * Divide numerator a and denominator b by gcd(a,b). | ||
| 87 | * | ||
| 88 | * @param a Numerator. | ||
| 89 | * @param b Denominator. | ||
| 90 | */ | ||
| 91 | static void simplifyRat(rat_int_t *a, rat_int_t *b) | ||
| 92 | { | ||
| 93 | rat_int_t tmp = gcd(*a, *b); | ||
| 94 | ✗ | if(tmp != 0) { | |
| 95 | 92 | *a /= tmp; | |
| 96 | ✗ | *b /= tmp; | |
| 97 | } | ||
| 98 | } | ||
| 99 | |||
| 100 | |||
| 101 | /** | ||
| 102 | * @brief Create rational number from numerator and denominator. | ||
| 103 | * | ||
| 104 | * Asserts denominator is non-zero and simplifies rational. | ||
| 105 | * | ||
| 106 | * @param numerator | ||
| 107 | * @param denominator | ||
| 108 | * @return RATIONAL numerator/denominator | ||
| 109 | */ | ||
| 110 | 92 | RATIONAL makeRATIONAL(rat_int_t numerator, rat_int_t denominator) | |
| 111 | { | ||
| 112 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 92 times.
|
92 | assertStreamPrint(NULL, denominator != 0, "RATIONAL zero denominator."); |
| 113 | simplifyRat(&numerator, &denominator); | ||
| 114 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 92 times.
|
92 | if(denominator < 0) { |
| 115 | ✗ | assertStreamPrint(NULL, numerator != RAT_INT_MIN, "RATIONAL numerator overflow."); | |
| 116 | ✗ | assertStreamPrint(NULL, denominator != RAT_INT_MIN, "RATIONAL denominator overflow."); | |
| 117 | ✗ | return (RATIONAL){-numerator, -denominator}; | |
| 118 | } | ||
| 119 | 92 | return (RATIONAL){numerator, denominator}; | |
| 120 | } | ||
| 121 | |||
| 122 | |||
| 123 | /** | ||
| 124 | * @brief Addition of two rational numbers. | ||
| 125 | * | ||
| 126 | * a/b + c/d = (ad + bc)/(bd) = (a(d/g) + (b/g)c)/((b/g)d) | ||
| 127 | * | ||
| 128 | * @param r1 First rational number. | ||
| 129 | * @param r2 Second rational number. | ||
| 130 | * @return RATIONAL Sum of r1 and r2. | ||
| 131 | */ | ||
| 132 | ✗ | RATIONAL addRat(RATIONAL r1, RATIONAL r2) | |
| 133 | { | ||
| 134 | rat_int_t num, den; | ||
| 135 | rat_int_t g = gcd(r1.den, r2.den); | ||
| 136 | ✗ | rat_int_t num1 = r2.den/g; | |
| 137 | ✗ | rat_int_t num2 = r1.den/g; | |
| 138 | ✗ | RAT_INT_MUL(num2, r2.den, den, '+'); /* den = (r1.den/g)*r2.den */ | |
| 139 | ✗ | RAT_INT_MUL(num1, r1.num, num1, '+'); /* num1 *= r1.num */ | |
| 140 | ✗ | RAT_INT_MUL(num2, r2.num, num2, '+'); /* num2 *= r2.num */ | |
| 141 | ✗ | RAT_INT_ADD(num1, num2, num, '+'); /* num = num1 + num2 */ | |
| 142 | simplifyRat(&num, &den); | ||
| 143 | ✗ | return (RATIONAL){num, den}; | |
| 144 | } | ||
| 145 | |||
| 146 | |||
| 147 | /** | ||
| 148 | * @brief Negation of a rational number. | ||
| 149 | * | ||
| 150 | * -(a/b) = (-a)/b | ||
| 151 | * | ||
| 152 | * @param r Rational number. | ||
| 153 | * @return RATIONAL Negative of r. | ||
| 154 | */ | ||
| 155 | ✗ | RATIONAL negRat(RATIONAL r) | |
| 156 | { | ||
| 157 | ✗ | assertStreamPrint(NULL, r.num != RAT_INT_MIN, | |
| 158 | "RATIONAL overflow. Unable to store result of -("RAT_FMT"/"RAT_FMT")", | ||
| 159 | r.num, r.den); | ||
| 160 | ✗ | return (RATIONAL){-r.num, r.den}; | |
| 161 | } | ||
| 162 | |||
| 163 | |||
| 164 | /** | ||
| 165 | * @brief Subtraction of two rational numbers. | ||
| 166 | * | ||
| 167 | * a - b = a + (-b) | ||
| 168 | * | ||
| 169 | * @param r1 First rational number. | ||
| 170 | * @param r2 Second rational number. | ||
| 171 | * @return RATIONAL Difference of r1 and r2. | ||
| 172 | */ | ||
| 173 | ✗ | RATIONAL subRat(RATIONAL r1, RATIONAL r2) { | |
| 174 | ✗ | return addRat(r1, negRat(r2)); | |
| 175 | } | ||
| 176 | |||
| 177 | |||
| 178 | /** | ||
| 179 | * @brief Multiplication of two rational numbers. | ||
| 180 | * | ||
| 181 | * a/b * c/d = (ac)/(bd) = ((a/g1)(c/g2))/((b/g2)(d/g1)) | ||
| 182 | * | ||
| 183 | * @param r1 First rational number. | ||
| 184 | * @param r2 Second rational number. | ||
| 185 | * @return RATIONAL Product of r1 and r2. | ||
| 186 | */ | ||
| 187 | ✗ | RATIONAL mulRat(RATIONAL r1, RATIONAL r2) | |
| 188 | { | ||
| 189 | rat_int_t num, den; | ||
| 190 | rat_int_t g1 = gcd(r1.num, r2.den); | ||
| 191 | rat_int_t g2 = gcd(r2.num, r1.den); | ||
| 192 | ✗ | RAT_INT_MUL(r1.num/g1, r2.num/g2, num, '*'); /* num = (a/g1)*(c/g2) */ | |
| 193 | ✗ | RAT_INT_MUL(r1.den/g2, r2.den/g1, den, '*'); /* den = (b/g2)*(d/g1) */ | |
| 194 | ✗ | return (RATIONAL){num, den}; | |
| 195 | } | ||
| 196 | |||
| 197 | |||
| 198 | /** | ||
| 199 | * @brief Reciprocal of a rational number. | ||
| 200 | * | ||
| 201 | * (a/b)^(-1) = b/a | ||
| 202 | * | ||
| 203 | * @param r Rational number. | ||
| 204 | * @return RATIONAL Reciprocal of r. | ||
| 205 | */ | ||
| 206 | ✗ | RATIONAL invRat(RATIONAL r) | |
| 207 | { | ||
| 208 | ✗ | assertStreamPrint(NULL, r.num != 0, "RATIONAL division by zero."); | |
| 209 | ✗ | if(r.num < 0) { | |
| 210 | ✗ | assertStreamPrint(NULL, r.num != RAT_INT_MIN, | |
| 211 | "RATIONAL overflow. Unable to store result of ("RAT_FMT"/"RAT_FMT")^(-1)", | ||
| 212 | r.num, r.den); | ||
| 213 | ✗ | return (RATIONAL){-r.den, -r.num}; | |
| 214 | } | ||
| 215 | ✗ | return (RATIONAL){r.den, r.num}; | |
| 216 | } | ||
| 217 | |||
| 218 | |||
| 219 | /** | ||
| 220 | * @brief Division of two rational numers. | ||
| 221 | * A.k.a multiplication with multiplicative inverse. | ||
| 222 | * | ||
| 223 | * (a/b) / (c/d) = a/b * d/c | ||
| 224 | * | ||
| 225 | * @param r1 First rational number. | ||
| 226 | * @param r2 Second rational number. | ||
| 227 | * @return RATIONAL Quotient of r1 and r2. | ||
| 228 | */ | ||
| 229 | ✗ | RATIONAL divRat(RATIONAL r1, RATIONAL r2) { | |
| 230 | ✗ | return mulRat(r1, invRat(r2)); | |
| 231 | } | ||
| 232 | |||
| 233 | |||
| 234 | /** | ||
| 235 | * @brief Get real approximation of rational number. | ||
| 236 | * | ||
| 237 | * @param a Rational number. | ||
| 238 | * @return double Real approximation. | ||
| 239 | */ | ||
| 240 | ✗ | double rat2Real(RATIONAL a) { | |
| 241 | ✗ | return (double)a.num / a.den; | |
| 242 | } | ||
| 243 | |||
| 244 | |||
| 245 | /** | ||
| 246 | * @brief Convert integer to rational number. | ||
| 247 | * | ||
| 248 | * @param n Integer | ||
| 249 | * @return RATIONAL Rational representation of n. | ||
| 250 | */ | ||
| 251 | ✗ | RATIONAL int2Rat(rat_int_t n) { | |
| 252 | ✗ | return (RATIONAL){n, 1}; | |
| 253 | } | ||
| 254 | |||
| 255 | |||
| 256 | /** | ||
| 257 | * @brief Ceil rational number. | ||
| 258 | * | ||
| 259 | * Return minimum integer a, for which a >= m / n | ||
| 260 | * | ||
| 261 | * @param a Rational number. | ||
| 262 | * @return long Smallest integer number greater or equal rational number. | ||
| 263 | */ | ||
| 264 | ✗ | rat_int_t ceilRat(RATIONAL a) { | |
| 265 | ✗ | return a.num / a.den + (a.num > 0 && a.num % a.den ? 1 : 0); | |
| 266 | } | ||
| 267 | |||
| 268 | |||
| 269 | /** | ||
| 270 | * @brief Strict ceil rational number. | ||
| 271 | * | ||
| 272 | * Return minimum integer a, for which a > m / n | ||
| 273 | * | ||
| 274 | * @param a Rational number. | ||
| 275 | * @return long Smallest integer number greater rational number. | ||
| 276 | */ | ||
| 277 | ✗ | rat_int_t ceilRatStrict(RATIONAL a) { | |
| 278 | ✗ | return a.num / a.den + (a.num < 0 && a.num % a.den ? 0 : 1); | |
| 279 | } | ||
| 280 | |||
| 281 | |||
| 282 | /** | ||
| 283 | * @brief Floor rational number | ||
| 284 | * | ||
| 285 | * Return maximum a, for which a <= m / n | ||
| 286 | * | ||
| 287 | * @param a Rational number. | ||
| 288 | * @return long Biggest integer number smaller or equal to rational number. | ||
| 289 | */ | ||
| 290 | ✗ | rat_int_t floorRat(RATIONAL a) { | |
| 291 | ✗ | return a.num / a.den - (a.num < 0 && a.num % a.den ? 1 : 0); | |
| 292 | } | ||
| 293 | |||
| 294 | |||
| 295 | /** | ||
| 296 | * @brief Strict floor rational number. | ||
| 297 | * | ||
| 298 | * Return maximum a, for which a < m / n | ||
| 299 | * | ||
| 300 | * @param a Rational number. | ||
| 301 | * @return long Biggest integer number smaller then rational number. | ||
| 302 | */ | ||
| 303 | ✗ | rat_int_t floorRatStrict(RATIONAL a) { | |
| 304 | ✗ | return a.num / a.den - (a.num > 0 && a.num % a.den ? 0 : 1); | |
| 305 | } | ||
| 306 |