Fractions and decimals

"""Adding fractions: the right rule, a tempting wrong rule, and computer decimals.

Python's fractions.Fraction does exact arithmetic with fractions. We use it to
compare the rule a/b + c/d = (ad + cb)/(bd) with the tempting wrong rule
"add the tops and add the bottoms", (a + c)/(b + d), on every pair of
fractions between 0 and 1 with denominators up to 12; to check the
cross-multiplication test ad = cb against exact equality; and to see why a
computer's ordinary decimals (floats) can disagree with exact fractions.
"""

from fractions import Fraction

def add_rule(a, b, c, d):
    """Numerator and denominator of a/b + c/d by the common-denominator rule."""
    return a * d + c * b, b * d

def wrong_rule(a, b, c, d):
    """The tempting mistake: add the numerators and add the denominators."""
    return a + c, b + d

def equivalent(a, b, c, d):
    """Cross-multiplication test: a/b and c/d are the same number."""
    return a * d == c * b

def fractions_in_unit_interval(max_den):
    """All pairs (a, b) with 1 <= a <= b <= max_den: fractions in (0, 1]."""
    return [(a, b) for b in range(1, max_den + 1) for a in range(1, b + 1)]

def compare_rules(max_den):
    """Return (pairs, right rule correct, wrong rule correct, wrong in between)."""
    fracs = fractions_in_unit_interval(max_den)
    right = wrong = between = 0
    for a, b in fracs:
        for c, d in fracs:
            exact = Fraction(a, b) + Fraction(c, d)
            right += Fraction(*add_rule(a, b, c, d)) == exact
            mediant = Fraction(*wrong_rule(a, b, c, d))
            wrong += mediant == exact
            low, high = sorted([Fraction(a, b), Fraction(c, d)])
            between += low <= mediant <= high
    return len(fracs) ** 2, right, wrong, between

def cross_test_agrees(max_num, max_den):
    """Return (pairs, agreements, equal pairs) for the cross-multiplication test."""
    fracs = [(a, b) for a in range(max_num + 1) for b in range(1, max_den + 1)]
    exact = {f: Fraction(*f) for f in fracs}
    pairs = agree = equal = 0
    for a, b in fracs:
        for c, d in fracs:
            pairs += 1
            same = exact[(a, b)] == exact[(c, d)]
            agree += equivalent(a, b, c, d) == same
            equal += same
    return pairs, agree, equal

def main():
    num, den = add_rule(1, 2, 1, 3)
    wn, wd = wrong_rule(1, 2, 1, 3)
    print(f"1/2 + 1/3 by the common-denominator rule: {num}/{den}")
    print(f"1/2 + 1/3 by 'add tops, add bottoms': {wn}/{wd}"
          f" = {float(Fraction(wn, wd))}, but 1/2 = 0.5")
    print("  The wrong answer is smaller than one of the parts!")
    print()
    pairs, right, wrong, between = compare_rules(12)
    print(f"All {pairs:,} pairs of fractions in (0, 1], denominators up to 12:")
    print(f"  right rule equals the exact sum:   {right:,} of {pairs:,}")
    print(f"  wrong rule equals the exact sum:   {wrong:,} of {pairs:,}")
    print(f"  wrong rule lands between the two:  {between:,} of {pairs:,}")
    pairs, agree, equal = cross_test_agrees(24, 12)
    print(f"Cross-multiplication test versus exact equality, {pairs:,} pairs:")
    print(f"  they agree {agree:,} times ({equal:,} pairs are equal numbers)")
    print()
    print(f"Floats:    0.1 + 0.2 == 0.3 is {0.1 + 0.2 == 0.3};"
          f" 0.1 + 0.2 = {0.1 + 0.2!r}")
    exact = Fraction(1, 10) + Fraction(2, 10) == Fraction(3, 10)
    print(f"Fractions: 1/10 + 2/10 == 3/10 is {exact}")
    print()
    print("Epistemic status: the checks are complete for the fractions listed,")
    print("so they prove nothing about larger denominators; the lesson's argument")
    print("(same-sized pieces, then count) covers every fraction. The float")
    print("result is a real counterexample: binary floats cannot store 0.1 exactly.")

if __name__ == "__main__":
    main()

"""Which fractions have decimals that stop?

Long division turns a fraction into a decimal. The remainders can only be
0, 1, ..., b - 1, so either a remainder of 0 appears (the decimal stops) or
a remainder repeats (the digits repeat forever). We print the decimals of
1/2 to 1/20 and then test the rule: a fraction in lowest terms has a decimal
that stops exactly when its denominator is built from 2s and 5s only.
"""

from math import gcd

def decimal_expansion(a, b):
    """Long division of a/b (b > 0, a >= 0).

    Returns (whole part, digits before the repeat, repeating digits); the
    repeating list is empty when the decimal stops.
    """
    if b <= 0 or a < 0:
        raise ValueError("need a >= 0 and b > 0")
    whole, r = divmod(a, b)
    digits, seen = [], {}
    while r and r not in seen:
        seen[r] = len(digits)
        digit, r = divmod(10 * r, b)
        digits.append(digit)
    if r == 0:
        return whole, digits, []
    start = seen[r]
    return whole, digits[:start], digits[start:]

def show(a, b):
    """The decimal as text, with the repeating block in brackets: 0.1(6)."""
    whole, fixed, repeat = decimal_expansion(a, b)
    text = "".join(map(str, fixed))
    if repeat:
        text += "(" + "".join(map(str, repeat)) + ")"
    return f"{whole}.{text}" if text else str(whole)

def only_twos_and_fives(n):
    """True if n > 0 has no prime factors other than 2 and 5."""
    if n <= 0:
        raise ValueError("n must be positive")
    for p in (2, 5):
        while n % p == 0:
            n //= p
    return n == 1

def check_rule(max_den):
    """Return (fractions checked, rule correct, decimals that stop)."""
    checked = correct = stops = 0
    for b in range(2, max_den + 1):
        for a in range(1, b):
            checked += 1
            terminates = decimal_expansion(a, b)[2] == []
            stops += terminates
            correct += terminates == only_twos_and_fives(b // gcd(a, b))
    return checked, correct, stops

def main():
    print("Long division of 1/b; brackets mark digits that repeat forever:")
    cells = [f"1/{b:<2} = {show(1, b)}" for b in range(2, 21)]
    for i in range(0, len(cells), 2):
        print("  " + "".join(f"{c:<36}" for c in cells[i:i + 2]).rstrip())
    print()
    checked, correct, stops = check_rule(100)
    print(f"All {checked:,} fractions a/b with 1 <= a < b <= 100:")
    print(f"  {stops:,} have decimals that stop; the rest repeat forever.")
    print("  Rule: 'stops exactly when the lowest-terms denominator uses only")
    print(f"  2s and 5s' is right in {correct:,} of {checked:,} cases.")
    print(f"  Careful with lowest terms: 3/12 = {show(3, 12)} stops, because"
          f" 3/12 = 1/4.")
    print()
    print("Epistemic status: long division here is exact (whole-number")
    print("remainders, no rounding), so each expansion is correct; the rule is")
    print("checked completely for these 4,950 fractions only. Why it holds for")
    print("every fraction: a stopping decimal is m/10^k, and the only prime")
    print("factors of 10^k are 2 and 5.")

if __name__ == "__main__":
    main()