Why use rational numbers?

In a previous article I've written about pseudo-random floating-point numbers generation. Another solution I haven't spoken of consist of using rational numbers. These numbers are defined as a fraction \(\frac{a}{b}\) where \(a\) and \(b\) are integers. From my previous article you should remember that there was efficient ways to generate pseudo-random integers, but difficulties appear when trying to generate floating-point number due to the way they are represented in memory. Rational numbers allow to represent floating-point numbers while using only integers, hence it makes straightforward the generation of uniformly distributed pseudo-random floating-point numbers from the generation of pseudo-random integers. Unfortunately, C does not provide an implementation of rationals, hence their omission in my previous article. I had an incomplete implementation in my previous framework, and this gave me the motivation to come back to it.

Rational numbers are a subset of the real numbers, all real numbers can't be represented using only the rational numbers, but anyway we are limited by the underlying encoding we are dealing with. If we choose to represent the rationals with three integers \(b\), \(n\) and \(d\): \(b+\frac{n}{d}\), \(b\in \mathbb{N}\), \(n\in \mathbb{N}^+\), \(d\in \mathbb{N}^{*+}\), \(n\le d\), and if we choose to encode them using int64_t and uint64_t, we can represent more than \(2^{128}\) values (almost \(2^{192}\) in fact), way more than the \(2^{64}\) values available with a real number encoded using double, and even more than the \(2^{80}\) values available with a real number encoded using long double. Also, we could represent the rational with only two integers, using a fixed value for the denominator, for example the maximum value representable by the encoding type. It simplifies the implementation, at the cost of the number of representable values (for example if \(d\) is fixed to 10, 1/3 can't be represented exactly anymore) and consequently the accuracy, thus I'll use here the three integers approach.

The range of values that can be represented is the same as the one of the type we choose for \(b\) plus one on the positive side to account for the addition of \(\frac{n}{d}\). If it's int64_t it becomes: \([-2^{63}, 2^{63}]\), or [-9223372036854775808, 9223372036854775808]. This is far smaller than the range available when using, for example, double: from -1.7E+308 to 1.7E+308, but rational numbers provide an uniform distribution of values, which is a big advantage in cases like random numbers generation for example. Each value is exactly \(\frac{1}{D}\) away from its two immediate neighbours, where \(D\) is the largest value encoded using uint64_t: \(2^{64}-1\) or 18446744073709551615. In comparison, the largest value not greater or equal than 9223372036854775808.0 that can be encoded using a double is 9223372036854774784.0, a gap of 1024.0! So, double allows for an extremely large range of values, but its precision gets so bad for large values that its usefulness becomes questionable. Rational numbers instead are restricted to a smaller range but guaranty a high precision all along.

All these numbers are extremely large (or small) and difficult to apprehend. To put them in perspective lets compare them to physical values. Our galaxy is approximately 1,000,000,000,000,000,000 kilometer in diameter, that would fit inside the [-9,223,372,036,854,775,808, 9,223,372,036,854,775,808] range. The radius of a hydrogen atom is approximately 5E−14 kilometer, that's larger than 1/18,446,744,073,709,551,615. Put together, choosing kilometer as the unit, the rational numbers with the encoding introduced above is precise and large enough to represent any point in our galaxy with an error less than the size of a hydrogen atom. Wow!

At that point some may wonder why we bother with double anyway. First, rational numbers take more memory: a double takes 8 bytes, same as an int64_t but it takes 3 of them for a rational. Sure you could use smaller integer types than int64_t, but that would decrease dramatically the range and precision. Nowadays hardware have so huge quantity of memory it may seem irrelevant, but if you're dealing with very large amount of data, or tight data transfer speed, it may become a problem. Second, the double precision gets as bad for large values as it gets good for small values. The double value next to 0.0 is 2.3E-308, way smaller than 1/18,446,744,073,709,551,615. However, rationals' precision is already so huge that the need for more in real application is probably rare. Third, and probably the main reason, our computing hardware implement double, not rationals, so if you want to use the latter you have to do it at software level, meaning a huge cost in speed performance. I've looked for hardware implementing rational but couldn't find any. (If you know one, please let me know!)

In defense of floating-point numbers some may also argue that their lack of precision becomes sensible only for large values, whose need in real application is also questionable. This is however incorrect, floating-point precision problems do occur in real-life application. For example, it led to the failure of a Patriot missile during the Gulf War, or problems in videogame as shown in this video by SimonDev about procedural generation of planets. Looking for alternative to the floating-point representation of real numbers is indeed necessary. Rational numbers have also their own limitations and flaws and don't guarantee to avoid the same kind of problems encountered with floating-point encoding. However, it provides an alternative which could prove advantageous in specific cases. In the end, there is no definite good or bad way to do, as always. It only depends on what you are trying to do, need to be decided case by case, so the more tools you have in your toolbox the better you'll be prepared.

Then, what about my own use cases? Will I ever make a simulation of our galaxy at atomic level? Certainly not. A procedurally generated planet? Much more probable even if I have no plan in that direction for the near future. Haven't this already be done by others? Of course it does. Regardless, it just looks too much fun to do it myself!

The implementation.

Lets first define the specs. I want:

• A data structure representing a rational number as 3 integers (one signed integer for the base, two unsigned integers for the fraction)
• A way to represent NaN
• Raising exception on numerical errors
• Reduction (keeping the numerator smaller than the denominator)
• Conversion from/to double
• Addition, negation, substraction
• Comparison (smaller, equal, greater)
• Absolute value
• Multiplication, inversion, division

There is really nothing difficult from the math point of view, however there are many complications that arise during implementation due to overflows and types conversion. Overflow during integers artihmetic is undefined behaviour, thus the only solution is to detect it prior it's happening and make the calculation in a way that avoid it if possible. Using signed integers everywhere would simplify type conversion problems but the numerator and denominator being always positive would leave half of the representable values unused, I don't want that. Reciprocally, using unsigned integers everywhere means creating rationals only able to represent positive values, which is not acceptable. Using a flag for signedness instead raises probably as many problems as it solves, I don't want that neither.

When mixing different types in the same C statement, implicit conversion occurs as detailed in this link. I'll need only to consider int64_t and uint64_t, for which the relevant rules become:

• In the assignment operator, the value of the right-hand operand is converted to the unqualified type of the left-hand operand.
• In the arithmetic operators, if both operands have different types, int64_t is implicitly converted to uint64_t.

I'll explicitly cast when necessary anyway as I don't expect everyone, starting with me, to know by heart the full spec. Relying on implicit rules is the way to confusion.

Data structure

The data structure is as simple as:

// Ratio object typedef struct CapyRatio { // Components of the ratio int64_t base; uint64_t num; uint64_t den; } CapyRatio;

For the representation of NaN, to manage infinities and overflows, a rational with denominator equal to 0 seems to make sense, then I choose: CapyRatio const capyRatioNaN = {.base = 0, .num = 0, .den = 0};. It also leaves the possibility to use the base and num values to identify several types of NaN.

Rational number reduction.

The reduction of a rational number has two steps: correct the base to bring the numerator less than the denominator, then divide the numerator and denominator by their greatest common divisor. The base may overflow during the correction step, this must be taken care of by checking first by how much it can be increased. To find the greatest common divisor, I'll use the Stein's algorithm (not detailed here).

// Reduce a CapyRatio // Input: // that: the CapyRatio to reduce // Output: // Return a new CapyRatio equals to the reduced CapyRatio // Exception: // May raise CapyExc_NumericalOverflow CapyRatio CapyRatioReduce(CapyRatio const that) { // Trying to reduce NaN gives NaN if (CapyRatioIsNaN(that)) return capyRatioNaN; // Declare the variable to memorise the result CapyRatio res = that; // If the numerator is greater than the denominator if (res.num >= res.den) { // If the reduction triggers an overflow of the base uint64_t inc = res.num / res.den; uint64_t threshold = 0; if (res.base < 0) threshold = (uint64_t)INT64_MAX + (uint64_t)(-res.base); else threshold = (uint64_t)INT64_MAX - (uint64_t)(res.base); if (inc > threshold) { // Raise an exception, or return NaN if the exception is not // caught raiseExc(CapyExc_NumericalOverflow); return capyRatioNaN; } // Update the components to keep the numerator smaller than the // denominator res.base += inc; res.num -= inc * res.den; } // If the numerator is not null if (res.num != 0) { // Divide the numerator and denominator by their gcd uint64_t gcd = CapyGcd(res.num, res.den); res.num /= gcd; res.den /= gcd; // Else, the numerator is null, ensure the denominator is equal to 1 // by convention } else res.den = 1; // Return the result return res; }

Conversion to double.

The conversion from rational to double is as simple as:

// Convert a CapyRatio to a double // Input: // that: the CapyRatio to convert // Output: // Return a double representing approximating the CapyRatio. double CapyRatioToDouble(CapyRatio const that) { return (double)(that.base) + (double)(that.num) / (double)(that.den); }

Conversion from double.

The conversion from double to rational can be obtained using fractions from the Farey sequence and binary search. First we calculate the integral and fractional values of the double using modf. Then we search for the approximation of the fractional part. We start with the two fractions \(a=\frac{0}{1}\) and \(b=\frac{1}{1}\), then we create the median fraction whose numerator is the sum of the numerators of \(a\) and \(b\), and the denominator is the sum of the denominators. If the median is lower than the decimals it replaces \(a\), if it's greater it replaces \(b\). And so on until \(a\) and \(b\) have converged (meaning their converted values to double are equal) or the denominators can't be added anymore without overflowing. A problem with this method is that it can take a very long time to converge. For example, consider the largest double less than 1.0:

union {uint64_t i; double d;} v; v.d = 1.0; --(v.i);

The bounds and median of the binary search evolve as follow:

0 + 0 / 1 < 0 + 1 / 2 < 0 + 1 / 1 0 + 1 / 2 < 0 + 2 / 3 < 0 + 1 / 1 0 + 2 / 3 < 0 + 3 / 4 < 0 + 1 / 1 0 + 3 / 4 < 0 + 4 / 5 < 0 + 1 / 1 0 + 4 / 5 < 0 + 5 / 6 < 0 + 1 / 1 0 + 5 / 6 < 0 + 6 / 7 < 0 + 1 / 1 ...

and goes on forever (well, not forever but way too long for practical use). One solution consists of limiting the denominator to an arbitrary maximum value. The larger the maximum value the better precision you get in the conversion but the slower it may take to converge. This solution is unsatisfying because it implies that the result of the conversion won't be as precise as it could be, and choosing that maximum value is just replacing the first problem with another one. A better solution is to jump several steps at once: instead of calculating the median as \(\frac{a_n+b_n}{a_d+b_d}\) we calculate it by solving \(\frac{a_n+kb_n}{a_d+kb_d}=x\) for \(k\in\mathbb{N}\), which gives \(k=\lfloor\frac{xa_d-a_n}{b_n-xb_d}\rfloor\), where \(a_n,a_d,b_n,b_d\) are the numerator and denominator of the bounding fractions, and \(x\) is the decimal part of the converted double. The calculation of \(k\) is performed after each component is converted to double, this means we have to accept the imprecision from double here as there is not much we can do about it. \(b_n-xb_d\) can't be equal to 0: it would mean \(x=\frac{b_n}{b_d}\), i.e. the algorithm has converged and already stopped. The only thing to take care of is when \(k\) is equal to 0. In that case we default to the calculation of the median as \(\frac{a_n+b_n}{a_d+b_d}\), equivalent to forcing \(k\) to 1. Jumping over steps results in performing the conversion in the previous example in one single step: 0 + 0 / 1 <= 0 + 9007199254740991 / 9007199254740992 <= 0 + 1 / 1, and is better illustrated when converting M_PI (the number in parenthesis is the size of the jump in number of steps):

3 + 0 / 1 <= 3 + 1 / 2 <= 3 + 1 / 1 (1) 3 + 0 / 1 <= 3 + 1 / 3 <= 3 + 1 / 2 (1) 3 + 0 / 1 <= 3 + 1 / 4 <= 3 + 1 / 3 (1) 3 + 0 / 1 <= 3 + 1 / 5 <= 3 + 1 / 4 (1) 3 + 0 / 1 <= 3 + 1 / 6 <= 3 + 1 / 5 (1) 3 + 0 / 1 <= 3 + 1 / 7 <= 3 + 1 / 6 (1) 3 + 0 / 1 <= 3 + 15 / 106 <= 3 + 1 / 7 (15) 3 + 15 / 106 <= 3 + 16 / 113 <= 3 + 1 / 7 (1) 3 + 15 / 106 <= 3 + 4687 / 33102 <= 3 + 16 / 113 (292) 3 + 4687 / 33102 <= 3 + 4703 / 33215 <= 3 + 16 / 113 (1) 3 + 4687 / 33102 <= 3 + 9390 / 66317 <= 3 + 4703 / 33215 (1) 3 + 9390 / 66317 <= 3 + 14093 / 99532 <= 3 + 4703 / 33215 (1) 3 + 9390 / 66317 <= 3 + 37576 / 265381 <= 3 + 14093 / 99532 (2) 3 + 37576 / 265381 <= 3 + 51669 / 364913 <= 3 + 14093 / 99532 (1) 3 + 37576 / 265381 <= 3 + 192583 / 1360120 <= 3 + 51669 / 364913 (3) 3 + 192583 / 1360120 <= 3 + 244252 / 1725033 <= 3 + 51669 / 364913 (1) 3 + 192583 / 1360120 <= 3 + 3612111 / 25510582 <= 3 + 244252 / 1725033 (14) 3 + 3612111 / 25510582 <= 3 + 3612111 / 25510582 <= 3 + 3612111 / 25510582

Finally the function to convert from double to rational is as follow:

Edited on 2022/07/19 to correct a bug when (dec < 0.0).

// Create a CapyRatio from a double // Input: // a: the double // Output: // Return a new CapyRatio representing the nearest possible value to // the input double. If the input value is out of range, return // capyRatioNaN. // Exception: // May raise CapyExc_NumericalOverflow CapyRatio CapyRatioFromDouble(double const a) { // Get the integral and fractional values of the input value double base; double dec = modf(a, &base); // Ensure the decimal part is positive according to CapyRatio format if (dec < 0.0) { dec += 1.0; base -= 1.0; } // Create the bounding CapyRatios CapyRatio low = CapyRatioCreate((int64_t)base, 0, 1); CapyRatio high = CapyRatioCreate(low.base, 1, 1); // If the input value isn't in the supported range if (!equald(base, (double)(low.base))) { // Raise an exception or return NaN in case the exception isn't caught raiseExc(CapyExc_NumericalOverflow); return capyRatioNaN; } // If the decimal are equal to zero, nothing to do if (equald(dec, 0.0)) return low; // Use a binary search with the Farey sequence to find the nearest // ratio to the decimals CapyRatio med = low; while (high.den <= UINT64_MAX - low.den && (low.num != high.num || low.den != high.den)) { // Calculate the number of steps we can safely jump over to speed // up the conversion uint64_t k = (uint64_t)floor( (dec * (double)(low.den) - (double)(low.num)) / ((double)(high.num) - dec * (double)(high.den))); if (k == 0) k = 1; med.num = low.num + k * high.num; med.den = low.den + k * high.den; // Update the bounding fraction double dblMed = (double)(med.num) / (double)(med.den); if (equald(dec, dblMed)) low = high = med; else if (dblMed > dec) high = med; else low = med; } // Return the reduced result return CapyRatioReduce(low); }

Addition.

Given two rationals \(x\) and \(y\), their sum is equal to \(x_b+y_b+\frac{x_ny_d+y_nx_d}{x_dy_d}\). Nothing can be done for the overflow on \(x_b+y_b\), so I'll just raise an exception in that case, but we need to perform the addition of the fractions first (we'll see why later).

Let's first consider the case where \(x_dy_d\) does not overflow. Given that \(x\) and \(y\) are in reduced form we know that \(x_n\le x_d\) and \(y_n\le y_d\). Added to \(x_dy_d\le UINT64\_MAX\), we are sure that \(x_ny_d\) and \(y_nx_d\) do not overflow neither. Then, we only have to check for the eventual overflow of the sum \(x_ny_d+y_nx_d\). If it does overflow, it means the base must be increased by one, and the numerator can be brought back in range by substracting \(x_dy_d\) from it. If we call \(a\) and \(b\) the two components of the sum \(x_ny_d+y_nx_d\), the safe way to do it is as follow: \(a-(x_dy_d-b)\). From \(x_n\le x_d\) and \(y_n\le y_d\) we know that \(x_dy_d\ge a\) and \(x_dy_d\ge b\), so \(x_dy_d-b\) will stay in range, and because \(x_dy_d\le UINT64\_MAX\lt a+b\) we are sure that \(x_dy_d-b\lt a\), hence \(a-(x_dy_d-b)\) will stay in range.

Now, the case where \(x_dy_d\) does overflow. We can't calculate \(\frac{x_ny_d+y_nx_d}{x_dy_d}\) but instead we can calculate \(\frac{x_n+y_nx_d/y_d}{x_d}\) or \(\frac{x_ny_d/x_d+y_n}{y_d}\). What follows is independent of which one is chosen, so I'll abstract the choice (explained later) by using the notation: \(\frac{a+bc/d}{b}\). The calculation of \(bc/d\) is guaranted here to not overflow: by definition \(b\le UINT64\_MAX\), then \(1\le \frac{UINT64\_MAX}{b}\), then \(c\le \frac{c.UINT64\_MAX}{b}\), and as we also have by definition \(c\le d\) we get \(c\le \frac{d.UINT64\_MAX}{b}\), hence \(\frac{bc}{d}\le UINT64\_MAX\). But, how to calculate it correctly? If we do the multiplication first we may overflow, and if we do the division first we loose accuracy. It takes a short detour, very very short, just a few thousand years ago in Egypt!

The ancient Egytptian multiplication, also known as the Russian peasant multiplication, is a method to multiply two integers, only by dividing and multiplying by two and performing additions. The multiplication of \(a\) by \(b\) looks like this:

uint64_t mul(uint64_t a, uint64_t b) { uint64_t res = 0; while (b > 0) { if (b & 1) res += a; a *= 2; b /= 2; } return res; }

Why does this help? It does exactly the same thing as \(a*b\), and doesn't care of the division we are interested in. Well, the point is that with some modification during the calculation, it is possible to account for the division, which actually come in handy to avoid the overflow (as much as possible). The peasant multiplication progressively update the result by increment of updated \(a\). Instead of dividing before the multiplication, too early, or after, too late, we divide during the multiplication, each time \(a\) becomes larger than the divisor, keeping \(a\) as small as possible and minimizing the chance of overflow.

Let's see on some simple examples how it works. To simplify I'll consider the maximum value above which overflow occurs is 10. The multiplication of \(2*3/4=1.5\) is completely safe. It goes as follow:

a b c r res 2 3 4 0 0 initialisation: res=(a/c)*b 0 3 4 2 0 initialisation: r=a%c; a=0 2 3 4 2 0 b odd: a=a+r 2 1 4 4 0 step: b=b/2; r=r*2 2 1 4 0 1 r>=c: r=r-c; res=res+b 2 1 4 0 1 b odd: a=a+r 2 0 4 0 1 step: b=b/2 result: res+a/c = 1+2/4

Other example, \(6*2/4=3\). The multiplication \(6*2=12\) would overflow, and neither 6 or 2 are divisible by 4, hence dividing first would not give the correct result. Thanks to the modified peasant multiplication we can calculate the exact result without overflow from any of the temporary variables.

a b c r res 6 2 4 0 2 initialisation: res=(a/c)*b 0 2 4 2 2 initialisation: r=a%c; a=0 0 1 4 4 2 step: b=b/2; r=r*2 0 1 4 0 3 r>=c: r=r-c; res=res+b 0 1 4 0 3 b odd: a=a+r 0 0 4 0 3 step: b=b/2 result: res+a/c = 3+0/4

The code is as follow:

// Calculate ab/c for a,b,c positive integers while avoiding eventual // intermediate overflow using the peasant multiplication // Inputs: // a, b, c: the three integers // Output: // Return the value of ab/c as an uint64_t for the integer part and // a CapyRatio for the remaining fractional part, or {0, capyRatioNaN} // if the calculation overflows. CapyPeasantMulDivRes CapyPeasantMulDiv( uint64_t a, uint64_t b, uint64_t c) { // Eliminate trivial cases if (c == 0) return (CapyPeasantMulDivRes){0, capyRatioNaN}; if (a == 0 || b == 0) return (CapyPeasantMulDivRes){0, capyRatioZero}; if (a == c) return (CapyPeasantMulDivRes){b, capyRatioZero}; if (b == c) return (CapyPeasantMulDivRes){a, capyRatioZero}; if (a == 1) return (CapyPeasantMulDivRes){ 0, CapyRatioReduce((CapyRatio){0, b, c})}; if (b == 1) return (CapyPeasantMulDivRes){ 0, CapyRatioReduce((CapyRatio){0, a, c})}; // Constant used during calculation uint64_t const half_uint64_max = UINT64_MAX / 2; // Initialise the result value with the result of the integer division CapyPeasantMulDivRes res = {0, capyRatioZero}; res.frac.den = c; res.base = a / c; if (res.base > UINT64_MAX / b) return (CapyPeasantMulDivRes){0, capyRatioNaN}; res.base *= b; // Variable to memorise the remainder from the integer division which // we have to multiply by b and add to the result value uint64_t r = a % c; // Peasant multiplication algorithm to calculate res.frac.num=rb, // modified to update res.base when accounting for the division of // res.frac.num by c. // Loop until b has been consumed, or the remainder is null. while (b > 0 && r > 0) { // If b is odd, the remainder is added to the result of the // multiplication. if (b & 1) { // Add the current remainder to the numerator of the fractional // part. If by adding it the numerator overflows, we are sure // that it will be larger than c. To avoid the overflow we can // first increment the numerator up to c, equivalent to reset it // to 0 and increment the base by one, and then update the // numerator with the remainder decreased by the initial numerator // minus c (which can't overflow anymore). uint64_t t = r; if (res.frac.num > UINT64_MAX - t) { t -= c - res.frac.num; if (res.base == UINT64_MAX) return (CapyPeasantMulDivRes){0, capyRatioNaN}; res.base++; res.frac.num = 0; } res.frac.num += t; // If the numerator becomes larger than c it means the current // value of res.frac.num/c is larger than 1, increment res.base // by one and decrement res.frac.num by c. This helps avoiding // overflows of a*b. while (res.frac.num >= c) { res.frac.num -= c; if (res.base == UINT64_MAX) return (CapyPeasantMulDivRes){0, capyRatioNaN}; res.base++; } } // Update the multiplication. b /= 2; // Avoid updating the remainder if b reach 0, we don't need it // anymore anyway so it avoid useless calculation. if (b > 0) { // If updating the remainder would overflow if (r > half_uint64_max) { // Increment the base and update the remainder accordingly if (res.base > UINT64_MAX - b) return (CapyPeasantMulDivRes){0, capyRatioNaN}; res.base += b; r -= c - r; // Update the remainder normally. } else r *= 2; // If the remainder becomes larger than divisor we can jump // forward by increment of b to go faster and avoid overflow on // the remainder. while (r >= c) { r -= c; if (res.base > UINT64_MAX - b) return (CapyPeasantMulDivRes){0, capyRatioNaN}; res.base += b; } } } // Reduce the fraction res.frac = CapyRatioReduce(res.frac); // Return the result return res; }

Back to \(\frac{a+bc/d}{b}\), which we now have rewritten as \(\frac{a+(b'+n'/d')}{b}=\frac{a+b'}{b}+\frac{n'}{bd'}\). I'll call \(\frac{n'}{bd'}\) the residual. Let's ignore it one second: the whole addition has become equal to \(x_b+y_b+\frac{a+b'}{b}\). The eventual overflow of \(a+b'\) can be canceled thanks to the division by \(b\) by incrementing the base of the result. That's why we need to calculate the sum of the numerators before the one of the bases: it may also be corrected here.

This leaves only the residual. If \(bd'\) can be calculated, it's just a new ratio that we'll need to add to the result by calling recursively the addition function. We are sure that this recursion ends because at each step the residual gets smaller, given that we always choose \(\frac{a+bc/d}{b}\) such as \(d\le b\). The proof goes as follow: \(d\le b \Rightarrow\frac{bc}{d}\ge c\ge 1\Rightarrow b'\ge 1\Rightarrow \frac{n'}{bd'}\lt \frac{c}{d}\). Final problem, here we know that \(bd\) overflows. Thanks to the reduction at the end of the Peasant multiplication/division \(d'\) may be smaller than \(d\) and avoid the overflow. If not there is also the possibility that \(\frac{n'}{b}\) is reducable, if so this also may save us from the overflow. And if all of that fails, then it means that \(\frac{n'}{bd'}\) can't be represented exactly within the accuracy of 1/UINT64_MAX. We have to approximate it, which can be safely calculated as follow with an error less than 1/UINT64_MAX: \(\frac{n'}{bd'}\simeq \frac{n'/\lfloor b/\lfloor UINT64\_MAX/d'\rfloor\rfloor}{UINT64\_MAX}\)

The code becomes:

// Reduce a CapyRatio // Input: // that: the CapyRatio to reduce // Output: // Return a new CapyRatio equals to the reduced CapyRatio // Exception: // May raise CapyExc_NumericalOverflow CapyRatio CapyRatioReduce(CapyRatio const that) { // Trying to reduce NaN gives NaN if (CapyRatioIsNaN(that)) return capyRatioNaN; // Declare the variable to memorise the result CapyRatio res = that; // If the numerator is greater than the denominator if (res.num >= res.den) { // If the reduction triggers an overflow of the base uint64_t inc = res.num / res.den; uint64_t threshold = 0; if (res.base < 0) threshold = (uint64_t)INT64_MAX + (uint64_t)(-res.base); else threshold = (uint64_t)INT64_MAX - (uint64_t)(res.base); if (inc > threshold) { // Raise an exception, or return NaN if the exception is not // caught raiseExc(CapyExc_NumericalOverflow); return capyRatioNaN; } // Update the components to keep the numerator smaller than the // denominator res.base += inc; res.num -= inc * res.den; } // If the numerator is not null if (res.num != 0) { // Divide the numerator and denominator by their gcd uint64_t gcd = CapyGcd(res.num, res.den); res.num /= gcd; res.den /= gcd; // Else, the numerator is null, ensure the denominator is equal to 1 // by convention } else res.den = 1; // Return the result return res; } // Add two CapyRatios // Inputs: // x: the first CapyRatio (must be in reduced form) // y: the second CapyRatio (must be in reduced form) // Output: // Return a new CapyRatio (in reduced form) equal to x+y // Exception: // May raise CapyExc_NumericalOverflow CapyRatio CapyRatioAdd( CapyRatio const x, CapyRatio const y) { // Trying to add NaN gives NaN if (CapyRatioIsNaN(x) || CapyRatioIsNaN(y)) return capyRatioNaN; // Declare the variable to memorise the result CapyRatio res = capyRatioZero; // Declare a variable to manage a residual part when adding fraction CapyRatio residual = capyRatioZero; // Variable to memorise the eventual increment of the base int64_t inc = 0; // If the multiplication of denominators doesn't overflow if (x.den <= (UINT64_MAX / y.den)) { // Calculate the result denominator res.den = x.den * y.den; // Calculate the two halves of the numerator (they can't overflow) uint64_t firstHalf = x.num * y.den; uint64_t sndHalf = y.num * x.den; // If the sum of the two halves doesn't overflow, calculate the // result numerator at once if (firstHalf <= (UINT64_MAX - sndHalf)) res.num = firstHalf + sndHalf; // Else, the sum of the two halves overflows else { // If the base cannot be incremented if (res.base == INT64_MAX) { // Raise an exception or return NaN in case the exception isn't // caught raiseExc(CapyExc_NumericalOverflow); return capyRatioNaN; } // Bring back the numerator in range by increasing the base // and correcting the numerator accordingly inc += 1; res.num = firstHalf - (res.den - sndHalf); } // Else, the multiplication of denominators overflows } else { // Variables to commonalise code uint64_t a, b, c, d; // Set the commonalised variables according to the smallest // component of the denominator if (x.den > y.den) { a = x.num; b = x.den; c = y.num; d = y.den; res.den = x.den; } else { a = y.num; b = y.den; c = x.num; d = x.den; res.den = y.den; } // Calculate bc/d first using the modified Peasant multiplication // algorithm. Here the fractional part is below the available // precision level, it is then ignored. CapyPeasantMulDivRes pmd = CapyPeasantMulDiv(b, c, d); if (CapyRatioIsNaN(pmd.frac)) { // Raise an exception or return NaN in case the exception isn't // caught raiseExc(CapyExc_NumericalOverflow); return capyRatioNaN; } // The base is normally in pmd.base, but for the trivial case it // may be in pmd.frac.base. Adding both cover all cases at once. res.num = pmd.base + pmd.frac.base; // Avoid doing extra useless calculation if the residual is null if (pmd.frac.num > 0) { // The numerator of the residual may be reducable by the result // denominator. If it is, do it right now as it can help to avoid // overflow in the next step. uint64_t gcd = CapyGcd(pmd.frac.num, res.den); if (gcd > 1) { pmd.frac.num /= gcd; residual.den = res.den / gcd; // Else it wasn't reducable, leave it as it is } else residual.den = res.den; // If the residual of the Peasant multiplication can be multiplied // by the denominator uint64_t r = UINT64_MAX / pmd.frac.den; if (residual.den <= r) { // Memorise the residual to add it later residual.num = pmd.frac.num; residual.den *= pmd.frac.den; // Else, the residual of the Peasant multiplication can't be // multiplied by the denominator } else { // Approximate the residual to the highest possible precision residual.num = pmd.frac.num / (res.den / r); residual.den = UINT64_MAX; } // Reduce the residual residual = CapyRatioReduce(residual); } // If the remaining addition does not overflow, do it if (res.num <= UINT64_MAX - a) res.num += a; // Else, the remaining addition overflows, increment the base // and perform the corrected addition in the numerator else { inc += 1; res.num = a - (res.den - res.num); } } // Check if the reduction would increase the base, if it does do // it right now to be able to check for overflow when performing // the addition of the two bases while (res.num >= res.den) { inc += 1; res.num -= res.den; } // Check for overflow on the addition of the bases if ((x.base >= 0 && y.base >= 0 && x.base > (INT64_MAX - y.base - inc)) || (x.base < 0 && y.base < 0 && (x.base + inc) < (INT64_MIN - y.base))) { // Raise an exception or return NaN in case the exception isn't caught raiseExc(CapyExc_NumericalOverflow); return capyRatioNaN; } // Add the bases and the eventual increment in the proper order to // avoid overflow if (x.base < y.base) res.base = (x.base + inc) + y.base; else res.base = (y.base + inc) + x.base; // Reduce the result res = CapyRatioReduce(res); // If the residual is not null, add it to the result if (residual.num != 0) res = CapyRatioAdd(res, residual); // Return the result return res; }

Negation.

The negation of a rational is easy: \(-\left(b+\frac{n}{d}\right)=-b-\frac{n}{d}=(-b-1)+\frac{d-n}{d}\). The only two things we have to take care of is the eventual overflow of \(-b-1\), and that INT64_MAX=-INT64_MIN-1. The code is as follow:

// Get the negative of a CapyRatio // Input: // x: the CapyRatio (must be in reduced form) // Output: // Return a new CapyRatio (in reduced form) equal to -x // Exception: // May raise CapyExc_NumericalOverflow CapyRatio CapyRatioNeg(CapyRatio const x) { // Negative of NaN gives NaN if (CapyRatioIsNaN(x)) return capyRatioNaN; // Check for overflow if (x.base == INT64_MIN && x.num == 0) { // Raise an exception or return NaN in case the exception isn't caught raiseExc(CapyExc_NumericalOverflow); return capyRatioNaN; } // Calculate the negative CapyRatio res; if (x.base > 0) res.base = -x.base - 1; else res.base = -(x.base + 1); res.num = x.den - x.num; res.den = x.den; // Reduce the result res = CapyRatioReduce(res); // Return the result return res; }

Substraction.

For the substraction, I simply reuse the two previous functions: \(x-y=x+(-y)\). The code is as follow:

// Substract two CapyRatios // Inputs: // x: the first CapyRatio (must be in reduced form) // y: the second CapyRatio (must be in reduced form) // Output: // Return a new CapyRatio (in reduced form) equal to x-y // Exception: // May raise CapyExc_NumericalOverflow CapyRatio CapyRatioSub( CapyRatio const x, CapyRatio const y) { return CapyRatioAdd(x, CapyRatioNeg(y)); }

Comparison.

To compare two rationals we can first compare their base, if they are different the result is trivial. If they are equals we need to compare the numerators, which can be done by multiplying them in such a way they become fractions with same numerator or denumerator, which in turn can easily be compared. Unfortunately, with this multiplication method comes again complication due to overflow. Comparing the difference of the two rationals to 0 instead, as we have already solved any complication in the substraction function (or rather the underlying addition and negation). One more catch here: the difference of the rationals can overflow too, it is the difference of the fractions only that must be checked. The code becomes:

// Compare two CapyRatios // Inputs: // x: the first CapyRatio (must be in reduced form) // y: the second CapyRatio (must be in reduced form) // Output: // Return -1 if x<y, else 0 if x==y, else 1 if x>y. Return 2 if an // exception occured and wasn't caught. // Exception: // May raise CapyExc_NumericalOverflow int8_t CapyRatioCmp( CapyRatio const x, CapyRatio const y) { if (CapyRatioIsNaN(x) || CapyRatioIsNaN(y)) { // Raise an exception, or return 2 in case the exception isn't caught raiseExc(CapyExc_NumericalOverflow); return 2; } // Compare the base to eliminate trivial cases if (x.base < y.base) return -1; else if (x.base > y.base) return 1; else if (x.num == y.num && x.den == y.den) return 0; // Else, the bases are equals else { // Calculate the difference of the fractional parts CapyRatio diff = CapyRatioSub( (CapyRatio){.base=0, .num=x.num, .den=x.den}, (CapyRatio){.base=0, .num=y.num, .den=y.den}); // Return the result according to the sign of the difference if (diff.base == 0 && diff.num == 0) return 0; else if (diff.base < 0) return -1; else return 1; } }

Absolute value.

The absolute value is trivial: if the rational is positive return it, else return its negative.

// Get the absolute value of a CapyRatio // Input: // x: the CapyRatio (must be in reduced form) // Output: // Return a new CapyRatio (in reduced form) equal to abs(x) // Exception: // May raise CapyExc_NumericalOverflow CapyRatio CapyRatioAbs(CapyRatio const x) { // Absolute value of NaN gives NaN if (CapyRatioIsNaN(x)) return capyRatioNaN; // Return the absolute value if (x.base < 0) return CapyRatioNeg(x); else return x; }

Multiplication.

The multiplication of \(x\) by \(y\) is equal to \(x_by_b+\frac{x_by_n}{y_d}+\frac{y_bx_n}{x_d}+\frac{x_ny_n}{x_dy_d}\). Calculating the 3 fractions separately instead of grouping them into one simplifies things, so that's what I'll do: calculate each component (\(x_by_b\), \(\frac{x_by_n}{y_d}\), \(\frac{y_bx_n}{x_d}\) and \(\frac{x_ny_n}{x_dy_d}\)) of the multiplication separately to get four rationals (plus a residual and three eventual corrective coefficients, see below), then add them all.

Same as for the addition of ratios, there are lot of overflow chances here, but it gets even worst due to the mix of signed and unsigned values. First, detecting the overflow on the multiplication of the bases \(x_by_b\) is a bit tricky due to the sign, the fractional part and the fact that \(INT64\_MIN\ne -INT64\_MAX\). If \(x_b=0\) or \(y_b=0\), the result is trivial. If \(x_b\gt 0\) and \(y_b\gt 0\), everything's working fine: there is overflow iff \(x_b\gt \lfloor\frac{INT64\_MAX}{y_b}\rfloor\). Other cases are processed as follow.

If \(x_b\gt 0\) and \(y_b\lt 0\), without taking account of the fractional part the multiplication of the bases only may indicates an incorrect overflow. Consider for example \((6+\frac{0}{1})*(-2+\frac{1}{2})\) with an overflow at -10. \(6*-2=-12\) does overflow, but the correct result with the fractional part is \(6*-1.5=-9\), which doesn't overflow. So what I'll do is rewrite \(x_by_b\) as \(x_b(y_b+1)-x_b\). \(x_b\gt 0\) guarantees \(-x_b\) exists, \(y_b\lt 0\) guarantees \(y_b+1\) exists. Then I check for overflow based on \(y_b+1\lt\lfloor\frac{INT64\_MIN}{x_b}\rfloor\). If the inequality is true, the multiplication is guaranteed to overflow, if not we don't know yet but we can move on to calculate it and the eventual overflow will be detected when adding the extra term \(-x_b\) (more on that later). The same goes for \(x_b\lt 0\) and \(y_b\gt 0\): I rewrite it as \(y_b(x_b+1)-y_b\) and use \(x_b+1\lt\lfloor\frac{INT64\_MIN}{y_b}\rfloor\).

If \(x_b\lt 0\) and \(y_b\lt 0\), the previous trick doesn't work any more: there is no term we can safely multiply by -1. One more step will save the day, if \(a\lt 0\), \(-a\) is not guaranteed to exist but \(-(a+1)\) is ! So we can use \(x_by_b=(x_b+1)(y_b+1)-(x_b+1)-(y_b+1)+1\) and check for \(-(x_b+1)\gt\lfloor\frac{INT64\_MAX}{-(y_b+1)}\rfloor\) (no overflow if \(y_b+1=0\)).

The second and third components can be calculated in the same manner, so I'll commonalise the code. We already know how to calculate them: it's the modified Peasant multiplication used in CapyRatioAdd. However they may be negative, while the modified Peasant multiplication handles only positive values. Easy, I'll ignore the sign, calculate, and if needed correct the result as follow: \(-(a+\frac{b}{c})=-(a+1)+\frac{c-b}{c}\).

The fourth component is resolved exactly as \(\frac{x_ny_d+y_nx_d}{x_dy_d}\) in the addition, with one half of the numerator equal to zero, which simplify things a bit.

As you can see we have already solved almost all the problems. There is only one new problem to tackle: how to add/substract the components without intermediate overflow. Given that adding a pair of positive and negative values cannot overflow, we can add the component by pairs of opposite sign values until there is only one value left or all the remaining values have same signedness (in which case we add them up).

The code becomes:

// Multiply two CapyRatios // Inputs: // x: the first CapyRatio (must be in reduced form) // y: the second CapyRatio (must be in reduced form) // Output: // Return a new CapyRatio (in reduced form) equal to x*y // Exception: // May raise CapyExc_NumericalOverflow CapyRatio CapyRatioMul( CapyRatio const x, CapyRatio const y) { // Trying to multiply NaN gives NaN if (CapyRatioIsNaN(x) || CapyRatioIsNaN(y)) return capyRatioNaN; // Declare the variable to memorise the result CapyRatio res = capyRatioZero; // Intermediate CapyRatio-s to memorise each component of the // multiplication #define nbComps 8 CapyRatio comps[nbComps] = {capyRatioZero, capyRatioZero, capyRatioZero, capyRatioZero, capyRatioZero, capyRatioZero, capyRatioZero, capyRatioZero}; // If the base product is not null if (x.base != 0 && y.base != 0) { // Variables to memorise the corrected bases int64_t a = x.base; int64_t b = y.base; // Check for overflow from the multplication of the bases bool overflow = false; if (a > 0 && b > 0) { if (a > INT64_MAX / b) overflow = true; } else if (a < 0 && b < 0) { a += 1; b += 1; comps[5].base = -a; comps[6].base = -b; comps[7].base = 1; if (b != 0 && -a > INT64_MAX / -b) overflow = true; } else if (a < 0 && b > 0) { a += 1; comps[6].base = -b; if (a < INT64_MIN / b) overflow = true; } else if (a > 0 && b < 0) { b += 1; comps[5].base = -a; if (b < INT64_MIN / a) overflow = true; } if (overflow == true) { // Raise an exception or return NaN in case the exception isn't // caught raiseExc(CapyExc_NumericalOverflow); return capyRatioNaN; } // Calculate the first component with the eventually corrected bases comps[0].base = a * b; } // Variables to commonalise code in the calculation of the second // and third components struct {uint64_t b; uint64_t n; uint64_t d; int s;} v[2]; v[0].b = labs(y.base); v[0].n = x.num; v[0].d = x.den; v[0].s = (y.base < 0 ? -1 : 1); v[1].b = labs(x.base); v[1].n = y.num; v[1].d = y.den; v[1].s = (x.base < 0 ? -1 : 1); // Loop to commonalise code of the calculation of the second and third // components loop (i, 2) { // If the numerator is not null if (v[i].b != 0 && v[i].n != 0) { // Calculate the intermediate fraction using the modified Peasant // algorithm. CapyPeasantMulDivRes pmd = CapyPeasantMulDiv(v[i].b, v[i].n, v[i].d); // Memorise the result. pmd.base + pmd.frac.base is guaranteed // to not overflow: it's the multiplication of an integer by a // fraction less than 1. comps[1 + i] = pmd.frac; comps[1 + i].base += pmd.base; // Correct the sign if the base was negative. if (v[i].s == -1) { comps[1 + i].base = -(comps[1 + i].base + 1); comps[1 + i].num = comps[1 + i].den - comps[1 + i].num; comps[1 + i] = CapyRatioReduce(comps[1 + i]); } } } // If the denominator of the third fraction doesn't overflow if (x.den <= UINT64_MAX / y.den) { // Calculate the denominator of the third fraction uint64_t den = x.den * y.den; // Now that the denominator is calculated we can reuse the modified // Peasant multiplication to calculate the third fraction CapyPeasantMulDivRes pmd = CapyPeasantMulDiv(x.num, y.num, den); // Memorise the result. pmd.base + pmd.frac.base is guaranteed // to be null: it's the multiplication of two fractions less // than 1. comps[3] = pmd.frac; // Else, the denominator of the third fraction overflows } else { // Variables to commonalise code uint64_t a, b; // Set the commonalised variables according to the smallest // component of the denominator if (x.den > y.den) { a = y.den; b = x.den; } else { a = x.den; b = y.den; } // Calculate x_n*y_n/a first using the modified Peasant // multiplication algorithm. Here the fractional part is below the // available precision level, it is then ignored. CapyPeasantMulDivRes pmd = CapyPeasantMulDiv(x.num, y.num, a); if (CapyRatioIsNaN(pmd.frac)) { // Raise an exception or return NaN in case the exception isn't // caught. raiseExc(CapyExc_NumericalOverflow); return capyRatioNaN; } // The base is normally in pmd.base, but for the trivial case it // may be in pmd.frac.base. Adding both cover all cases at once. comps[3].num = pmd.frac.base + pmd.base; comps[3].den = b; comps[3] = CapyRatioReduce(comps[3]); // Avoid doing extra useless calculation if the residual is null if (pmd.frac.num > 0) { // Initialise the residual comps[4] = pmd.frac; comps[4].base = 0; // The numerator of the residual may be reducable by the result // denominator. If it is, do it right now as it can help to avoid // overflow in the next step. uint64_t gcd = CapyGcd(comps[4].num, b); if (gcd > 1) { comps[4].num /= gcd; b /= gcd; } // If the residual of the Peasant multiplication can be multiplied // by the denominator uint64_t r = UINT64_MAX / b; if (comps[4].den <= r) { // Multiply by the denominator comps[4].den *= b; // Else, the residual of the Peasant multiplication can't be // multiplied by the denominator } else { // Approximate the residual to the highest possible precision comps[4].num /= comps[4].den / r; comps[4].den = UINT64_MAX; } // Reduce the residual comps[4] = CapyRatioReduce(comps[4]); } } // Variables to memorise the components gathered by signedness CapyRatio posComps[nbComps]; CapyRatio negComps[nbComps]; int nbPos = 0; int nbNeg = 0; // Gather the components by signedness loop (i, nbComps) if (comps[i].base >= 0) posComps[nbPos++] = comps[i]; else negComps[nbNeg++] = comps[i]; // Add the components by pair of positive/negative values int iPos = 0; int iNeg = 0; while (iPos < nbPos && iNeg < nbNeg) { CapyRatio p = CapyRatioAdd(posComps[iPos], negComps[iNeg]); if (p.base >= 0) { posComps[iPos] = p; ++iNeg; } else { negComps[iNeg] = p; ++iPos; } } // Here there is no more pair, just sum up the remaining values CapyRatio* remainVals; int iRemain; int nbRemain; if (iPos < nbPos) { remainVals = posComps; iRemain = iPos; nbRemain = nbPos; } else { remainVals = negComps; iRemain = iNeg; nbRemain = nbNeg; } while (iRemain < nbRemain) res = CapyRatioAdd(res, remainVals[iRemain++]); // Return the result return res; }

Inverse.

The inverse of \(x\) is equal to \(\frac{x_d}{x_bx_d+x_n}\). Lets first get rid of the sign problem, if \(x_b\lt 0\) I'll calculate the inverse of \(-x\) and then return the negative of the result. Now the only problem is the possible overflow in the calcualtion of \(x_bx_d+x_n\), i.e. \(x_bx_d+x_n\gt UINT64\_MAX\), equivalent to \(\frac{x_bx_d+x_n}{UINT64\_MAX}\gt 1\). This can be easily calculated thanks to the Peasant modified multiplication plus an addition. If there is no overflow, the result is calculated directly, else we need to approximate the results, which can be done with a dichotomic search. There is a trick though in the update of the bounds: we need a condition to compare the median to the searched value. This search value being the inverse of the input, we can multiply the median with the input, if it's greater than 1 the median is too large, else it's too small. The convergence is reached when multiplying the median by the input gives exactly 1, or when it isn't possible anymore to calculate the median without overflowing.

And that's it for the inverse. The code becomes:

// Get the inverse of a CapyRatio // Input: // x: the CapyRatio (must be in reduced form) // Output: // Return a new CapyRatio (in reduced form) equal to 1/that // Exception: // May raise CapyExc_NumericalOverflow CapyRatio CapyRatioInv(CapyRatio const x) { // Eliminate the trivial case of x equal to 0 if (x.base == 0 && x.num == 0) return capyRatioNaN; // Variable to memorise the input eventually corrected to be positive CapyRatio y; if (x.base < 0) y = CapyRatioNeg(x); else y = x; // Trying to inverse NaN gives NaN if (CapyRatioIsNaN(y)) return capyRatioNaN; // Check for overflow CapyRatio threshold = {.base=0, .num=1, .den=INT64_MAX}; if (CapyRatioCmp(y, threshold) < 0) return capyRatioNaN; // Variable to memorise the result CapyRatio res = capyRatioZero; // Variable to check for overflow of the denominator CapyPeasantMulDivRes pmd = CapyPeasantMulDiv(y.base, y.den, UINT64_MAX); CapyRatio z = {.base=0, .num=y.num, .den=UINT64_MAX}; pmd.frac = CapyRatioAdd(pmd.frac, z); // If the result can be calculated directly if (pmd.base == 0 && pmd.frac.base == 0) { // Calculate the inverse res.num = y.den; res.den = (uint64_t)(y.base) * y.den + y.num; // Else the result can't be calculated directly, approximate it by // dichotomic search. } else { // Bounds of the dichotomic search. CapyRatio low = {.base = 0, .num = 0, .den = 1}; CapyRatio high = {.base = 0, .num = 1, .den = 1}; // Loop until convergence. while (low.den <= UINT64_MAX / 2 && high.den <= UINT64_MAX && (low.num != high.num || low.den != high.den)) { // Calculate the median CapyRatio a = low; a.den *= 2; CapyRatio b = high; b.den *= 2; CapyRatio ab = CapyRatioAdd(a, b); // If the median could not be calculated, end the search if (CapyRatioIsNaN(ab)) low = high; // Else the median could be calculated else { // Update the result res = ab; // Update the bounds. We can't calculate the value we are // approximating, but we have a comparison condition on the // median by multiplying it with the input. If the result is // greater than 1 the median is greater than the approximated // value, if it's lower than 1 it is lower than the approximated // value. CapyRatio mul = CapyRatioMul(y, res); if (res.num == low.num && res.den == low.den) low = high; else if (mul.base == 1 && mul.num == 0) low = high; else if (mul.base == 0) low = res; else high = res; } } } // Reduce the result res = CapyRatioReduce(res); // If the input was negative, correct the result if (x.base < 0) res = CapyRatioNeg(res); // Return the result return res; }

Edited on 2023/12/12 to refactor the test (pmd.base + pmd.frac.base == 0) into (pmd.base == 0 && pmd.frac.base == 0).

Division.

For the division I simply reuse the two previous functions: \(\frac{a}{b}=a*(b)^{-1}\).

// Divide two CapyRatios // Inputs: // x: the first CapyRatio (must be in reduced form) // y: the second CapyRatio (must be in reduced form) // Output: // Return a new CapyRatio (in reduced form) equal to x/y // Exception: // May raise CapyExc_NumericalOverflow CapyRatio CapyRatioDiv( CapyRatio const x, CapyRatio const y) { // Reuse the inversion and multiplication return CapyRatioMul(x, CapyRatioInv(y)); }

Unit tests.

Given the complexity of the algorithms introduced here, the probability for a bug is not low ! To avoid it as much as possible I've checked my code carefully with the following unit tests. I've considered the set of 28 test values made of the bases {INT64_MIN, INT64_MIN+1, -1, 0, 1, INT64_MAX-1, INT64_MAX} combined with the fractional parts {0/1, 1/UINT64_MAX, 1/2, (UINT64_MAX-1)/UINT64_MAX}. Then I've checked the addition and multiplication of each pair, as well as the negation and the inversion of each one, against truth tables (calculated manually and using Wolfram Alpha). For the negation I checked that \(x+(-x)=0\), and for the inversion \((x^{-1})^{-1}=x)\) and \(x*(x^{-1})=1\). Reduction and conversion from/to double were also tested on a different, smaller, set of values. Unit tests' result were as expected, within an inaccuracy lower than 8/UINT64_MAX for the inverse of the inverse, 6/UINT64_MAX for \(x*(x^{-1})=1\), and 2/UINT64_MAX for the other tests.

And that's all folks for this time ! There is a lot more about rationals to have fun with, and I'll probably come back to it one day. I'll also probably add a small example of pratical use in the near future, in a separate article which I'll link here.

Edit 2022/02/01: Check how CapyRatio can be used to solve the n-body problem in this article.