REAL Numbers: Exact Decimal Arithmetic

REAL Numbers: Exact Decimal Arithmetic

Most programmers assume numbers behave exactly as written. In practice, this is not always true.

Most mainstream programming languages use IEEE binary floating-point for their standard floating-point numeric types.

For example:

In most mainstream languages, this expression does not produce exactly 0. Instead, it produces a very small leftover value because binary floating-point cannot represent many ordinary decimal fractions, such as 0.1, exactly.

Sheerpower's REAL data type is designed around decimal arithmetic rather than binary floating-point. Ordinary decimal values are represented exactly, and within REAL's normal fixed-decimal range of up to 54 digits before and 16 digits after the decimal point, addition and subtraction are exact. This prevents the small accumulated errors commonly associated with binary floating-point and allows results to be compared directly using =, without epsilon tolerances.

Rounding: Half-to-Even, and round()

Some operations, particularly multiplication and division, can produce more digits than REAL can retain. When rounding is necessary, Sheerpower uses half-to-even rounding, producing predictable results while preserving up to 64 significant decimal digits.

For example:

Output:

sqrt(2) has an infinite decimal expansion, so REAL keeps it correctly rounded to sixteen places: 1.4142135623730950. The printed 1.414213562373095 shows fifteen of them because the sixteenth is a zero, and print leaves off trailing zeros. If that value is later used in another calculation, any additional rounding that becomes necessary is performed consistently according to REAL's decimal arithmetic rules.

Half-to-even sends an exact tie to the even digit, so ties go up as often as down and rounding adds no bias to a long column of results. The round() function is different on purpose: it rounds a tie away from zero (half-up), the rule business code expects, so round(2.5, 0) is 3 where the arithmetic rule would give 2:

Output:

For business calculations such as currency, results can still be rounded explicitly to the required working scale — for example, two decimal places — when that scale is part of the application's requirements.

A whole-number power is computed on all of its digits and rounded only once, at the end, so even a small base raised to a large power keeps every significant digit it can, and a negative power is a single divide. (x ^ 2 is the same one multiply as x * x.)

Output:

How Many Digits a Value Keeps

Very small results keep 64 digits. The sixteen places are an absolute grid: near 1 they hold sixteen significant digits, but a value like 0.0000000000000012 would keep only two of them. So a result smaller than 0.00000001 (10^-8) that the sixteen places cannot hold exactly is kept in scientific form instead, with 64 significant digits, rounded once. A small value the places do hold exactly stays on them, and from 10^-8 up a result uses the sixteen places as always:

Output:

How many digits a rounded result keeps. Above 0.00000001 the sixteen places count from the decimal point, so the size of a result decides how many significant digits it keeps. A third, at several sizes:

Output:

That matters when a small ratio will be multiplied back up. A 5% annual rate per minute lies between 0.00000001 and 0.0000001, so it keeps nine digits (the ninth is a zero), and multiplied back by the 525,600 minutes in a year it misses 0.05 in the eleventh place. The same rate per second is below 0.00000001, keeps 64 digits, and comes back as 0.05. divr() (Exact Fractions, below) keeps the exact fraction at any size:

Output:

The rule: when a small ratio will be multiplied back up, make it with divr().

How many of those digits are right. Every result is the correctly rounded answer for the values it was given. When one of those values was itself rounded earlier, the extra digits carry that rounding forward: they are exact for the rounded value, not for the number you had in mind. Here x holds 1/3 to sixteen places, so x / 700000000 agrees with 1/2100000000 for sixteen digits and then follows the rounded x. isrounded() says a value was rounded somewhere along the way; when every digit matters, keep the exact fraction with divr() (below) or compute from the original numbers:

Output:

Overflow and Underflow

Overflow does not silently wrap. If a result is larger than REAL's supported numeric range, Sheerpower detects it and raises a runtime exception. A result too small for the range — below 1e-1000000000 — becomes 0, the way IEEE floating point and Python's decimal type handle underflow: it is closer to 0 than to any value a REAL can hold.

Output:

Functions Computed Through a Double

The statistics built on a square root — stats$pcorr(), stats$spearman() and stats$linreg_r2() — work from exact sums and round only once, so every digit they print is right. A very small result that comes from a C double (exp(), log(), the trigonometric functions) keeps the double's fifteen digits, not 64. Near 1 the logarithms do better than the double alone could: 1.000000001 is not exact in binary, and a language that hands its double to the logarithm — NumPy, for one — prints 1.000000082240371e-9, right to seven digits. Sheerpower still holds the exact digits of x, so log(), log2() and log10() add back the error of the conversion: log(1.000000001) is 9.999999995e-10, every digit right.

Functions computed through a double. exp(), log(), log2(), log10(), sin(), cos(), tan(), asin(), acos(), atn(), sinh(), cosh() and tanh() work through a C double, so their results carry the double's fifteen significant digits and are marked *Rounded* unless the double is exact (cos(0) is a plain 1). sec(), csc() and cot() divide one of those results, so they print sixteen places but are no more accurate than the double. A double's range is not theirs, though: past it exp(), a non-whole power and the logarithms split the value into its power of ten and the rest, and only the rest goes through the double, so exp(1000) is 1.97007111401705e+434 and 10 ^ 400.5 is 3.16227766016838e+400, still to fifteen digits. sinh() and cosh() do stop there: sinh(1000) raises "Numeric value out of range for this function". sqrt() is not on the list: it works in whole-number arithmetic and is correctly rounded.

Output:

Exact Fractions: division that stays exact

Divide two whole numbers and the answer often will not fit in the sixteen digits a decimal real provides: 1/3 becomes 0.3333333333333333, and three of those add up to 0.9999999999999999 -- not 1. Many languages hand you the rounded number and leave you to notice. Sheerpower offers the exact fraction when you ask for it by name.

divr(a, b) divides and keeps the EXACT FRACTION when the quotient does not terminate, and the fraction stays exact through the add, subtract, multiply, divide and whole-number powers that follow, while its denominator is at most 4,294,967,295 (2³² - 1). The / operator is always the sixteen-place decimal divide, so a program says exactly where it wants exactness:

Output:

What joins a fraction exactly. Whole numbers of any size (past 10⁵⁴, in scientific form, too), other fractions, and short decimals: a decimal counts as the exact fraction it is written as (0.25 is 1/4, 0.05 is 1/20) whenever that fraction's denominator is at most 4,294,967,295 — every decimal of nine places or fewer. That holds for the arguments of divr() and for anything a fraction meets:

Output:

A quotient that terminates is simply a decimal: divr(1, 2) is 0.5 and divr(6, 3) is 2 (typeof$() says *Narrow*). Because 0.5 is a short decimal, x + divr(1, 2) above is still exactly 5/6.

The same holds for every result, not only for divr(): when arithmetic on a fraction lands on a value that terminates, the value is stored as that decimal — still exact, and fraction$() still shows it as a fraction. A monthly rate of 1/240 times 7 is still a fraction; times 12 it is 12/240 = 1/20, the decimal 0.05; and divr(1, 3) + divr(1, 6) is the decimal 0.5:

Output:

Where a fraction becomes a decimal. A LONG decimal — one whose fraction would need a bigger denominator, such as the sixteen places of a rounded result like 0.3333333333333333 — is taken as the decimal it is: beside it, in arithmetic and in comparison, a fraction is the sixteen-place decimal it rounds to. So a value printed and read back compares equal to the fraction that printed it. A result whose denominator would pass 4,294,967,295 is rounded the same way — once, from its exact value, however large the power that made it — and every print, str$() and f$() shows the decimal:

Output:

How long a fraction lasts. A sum of fractions needs a denominator every term divides, so with prime denominators they multiply and the sum reaches the limit quickly. The sum of 1/p over the primes from 3 to 29 is still exact; adding 1/31 needs a denominator of 100,280,245,065, and the sum becomes its sixteen-place decimal:

Output:

A fraction lives only in memory. Written to a file, or anywhere as text, it is its sixteen-place decimal, and that is what reads back: equal to the fraction, but times 3 it no longer gives 1:

Output:

Very small results. A fraction smaller than 0.00000001 follows the rule for every very small result: its decimal is what the / divide of the same two numbers gives, 64 significant digits in scientific form (*Sci*), marked *Rounded* unless those digits are the whole value. A tiny fraction does not turn into 0 (only a result below 1e-1000000000, the end of REAL's range, does):

Output:

typeof$() says *Rational*, fraction$() shows the fraction, numerator() / denominator() take it apart, decimal() rounds it, and divr() maps over an array (divr(v, 3)).

Those functions take a value in scientific form too. A whole one is its own fraction, written as str$() writes it. A rounded value smaller than 0.00000001 has 64 significant digits, so its denominator has 72 digits or more: fraction$() writes it in e-notation, the way str$() writes so large a number, and fraction() of a value below about 10^-10 is 0, the closest fraction any allowed denominator can make:

Output:

The Modes and the *Rounded* Mark

typeof$() names the mode a value is in: *Narrow* (18 integer digits and 16 places), *Wide* (54 integer digits), *Sci* (64 significant digits with an exponent), *Rational* (an exact fraction), *Complex*, *NaN* and *Inf*. The modes are automatic: a value moves between them by its size and its origin, never by a declaration.

One more marker can follow the mode: *Rounded*, when Sheerpower had to round in making the value — a multiply or divide that needed more than sixteen places, a result past 64 significant digits, a square root that is not exact, a function computed through a double whose answer is inexact (sin(1), but not cos(0)), a fraction turned into its decimal, or text with more digits than fit. Arithmetic on a rounded value carries the mark on. The functions that round because you asked — round(), truncate(), int(), floor(), ceil() — never add it; they keep their input's, because a value rounded earlier can sit on the wrong side of the boundary: (1 / 3) * 7.5 is 2.4999999999999998, not 2.5, so round() gives 2 where 7.5 / 3 rounds to 3. Zero never carries the mark: (1 / 3) * 0 is a plain 0. The mark says how a value was made, not what it is: a rounded 0.3333333333333333 equals one typed in, and lookups and unique() treat the two as one value. It lives only in memory — text, files and table fields do not carry it.

Arrays, statistics and complex values follow the same rule. sort() keeps each element's mark; a statistic — stats$max(), stats$median(), stats$percentile() and the rest — is marked when its data holds a rounded value; matmul() marks a product that needed more than sixteen places. Complex arithmetic keeps thirty digits in each part and marks a result it had to round to fit them, and a complex function computed through doubles is marked when the double is inexact.

Output:

Output:

isrounded(x) asks the question directly: it returns 1 when x carries the mark and 0 otherwise, and like isnan() it works on a whole array — stats$sum(isrounded(v)) counts the rounded values and filter(v, not isrounded(v)) keeps the exact ones. The mark also decides how a very large value prints: an exact 64-digit integer, such as 1e63 or the 304th Fibonacci number, prints every digit, while one that had to be rounded to fit prints in scientific notation.

Output:

Not-a-Number and Infinity

A REAL can hold two special values besides ordinary numbers: NaN ("not-a-number") and infinity. They are written _nan, _inf, and -_inf.

They are values you put in. A constant, val("NaN"), a CSV cell that reads Inf, a float field read from a table — each of those can bring one into a program. What Sheerpower does not do is manufacture them from ordinary numbers: 1 / 0 still raises DIVBY0, sqrt(-1) still raises a range error, and assigning either value to an INTEGER raises NUMOVER. So a NaN or an infinity in your data is always something that was put there, never a calculation that went wrong behind your back.

Once one is present, it propagates through arithmetic the way IEEE 754 specifies:

• any expression containing a NaN is a NaN
• _inf + 1 is _inf
• _inf - _inf is _nan
• 1 / _inf is 0

div0inf(a, b) is the division that produces these values on purpose: _inf or -_inf for a nonzero number over zero, _nan for 0 / 0, and the ordinary exact quotient otherwise. It is the companion of div0(a, b), which returns 0 for division by zero instead. Because you have to call it by name, every place an infinity can enter a calculation is visible in the source.

A NaN is never equal to anything, not even to another NaN. Two variables that both hold _nan compare unequal, so an if x = y test quietly passes over exactly the values you were looking for. Sheerpower will not let you write the literal form at all: if x = _nan is a compile-time error, because the answer is always "no" and the comparison is therefore always a mistake. Use isnan(var), which returns 1 when the value is a NaN and 0 otherwise.

Infinity is different. It compares normally, so if x = _inf is legal and does what it says. But it tests for one sign only — isinf(var) returns 1 for both _inf and -_inf, which is usually the question you meant to ask.

Complex Numbers

Sheerpower REALs can also represent complex numbers. Complex-number syntax, arithmetic, functions, precision, and FFT support are covered in the separate Processing Complex Numbers tutorial Processing Complex Numbers.

Inside a REAL: the Internal Representation

Sheerpower stores REAL numbers using multiple 64-bit integers. In its normal fixed-point representation, a REAL consists of:

• IP — The integer portion of the number (digits to the left of the decimal point), supporting approximately 54 decimal digits.
• FP — The fractional portion of the number (digits to the right of the decimal point), stored as an integer scaled by 10¹⁶. This preserves up to 16 decimal digits exactly.

Together, these provide up to 70 significant decimal digits.

When still larger values are required, REAL automatically uses scientific notation while retaining up to 64 significant digits. The exponent can have a magnitude of up to 1,000,000,000, allowing values on the order of 10^1000000000 as well as correspondingly tiny values.

That is why the two limits differ: just under 10⁵⁴ a fixed-decimal value holds 70 digits (54 + 16), while past it the scientific form holds 64, so at 10⁵⁴ only nine decimal places are left. 10 ^ 54 + 0.123456789123 prints as 1000000000000000000000000000000000000000000000000000000.123456789 and is marked *Rounded*.

The sign is kept once, as the top bit of the highest word; IP and FP always hold the magnitude, so -13.599 and 13.599 differ in that one bit. The ip() and fp() functions give the parts with the sign: ip(-13.599) is -13 and fp(-13.599) is -0.599.

This separation allows Sheerpower to maintain complete decimal accuracy capped at 54 digits before and 16 digits after the decimal point, predictably.

In standard floating-point math, "magnitude mismatch" errors occur when numbers of vastly different scales are mixed. For example, adding a micro-transaction to a national balance sheet might result in the smaller value simply disappearing due to the limitations of binary representation—technically known as loss of significance (absorption). By storing the integer and fractional portions separately, Sheerpower REALs prevent magnitude mismatch. You can safely add one quintillion and 0.000000000000123 without losing a single digit of precision.

Visualizing the Structure

The illustration below shows how a REAL value is stored:

Why This Is Important

1. Reliable financial systems — every penny is preserved
2. Accounting that stays accurate—$0.01 stays $0.01 forever, no drift over millions of transactions
3. Clean and predictable results for statistical and data-processing workloads
4. Comparisons you can rely on — sums and differences compare exactly with =; a rounded product or quotient compares as the sixteen-place decimal it is, so 1 / 3 * 3 = 1 is false and divr(1, 3) * 3 = 1 is true
5. Totals that do not depend on order — exact addition gives the same sum however the values are grouped, so a threaded stats$sum() or stats$mean() returns the same digits however the array is split
6. Data interchange that keeps every digit — str$() writes every digit a value holds (a fraction, its sixteen-place decimal) and val() reads it back equal, across JSON, CSV, databases and APIs. The other side must read decimal text exactly too: many JSON parsers read numbers into doubles

Exact REAL Math Across Extreme Magnitudes

In IEEE doubles the same loop ends at exactly 105,000,000,000,000: at that size neighboring doubles are 0.015625 apart, so each 0.00001 added is rounded away and the whole $100 is lost.

Summary