So I can generate many tuples like this:
(601550405185810455248373798733610900689885946410558295383908863020551447581889414152035914344864580636662293641050147614154610394724089543305418716041082523503171641011728703744273399267895810412812627682686305964507416778143771218949050158028407021152173879713433156038667304976240165476457605035649956901133077856035193743615197184,
191479441008508487760634222418439911957601682008868450843945373670464368694409556412664937174591858285324642229867265839916055393493144203677491629737464170928066273172154431360491037381070068776978301192672069310596051608957593418323738306558817176090035871566224143565145070495468977426925354101694409791889484040439128875732421875)
They are all tuples of two ints, they represent fractions with (nearly) infinite precision (bounded only by computer memory), the first number is the numerator, the second denominator.
If we divide them, we get the first 101 decimal places of π:
'3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798'
I did everything using integer math and without using any libraries. I didn't use any floating point because I know they are all of finite precision, in contrast to int's infinite precision in Python.
Python's float uses this, and it only has log10(2) * 52 = 15.653559774527022 decimal places, far less than what I wanted.
I wrote two correct functions to get n specified decimal places of any fraction:
from typing import Dict, List, Tuple
def decimalize(dividend: int, divisor: int, places: int) -> str:
div, mod = divmod(dividend, divisor)
result = [f"{div}."]
while mod and places:
div, mod = divmod(mod * 10, divisor)
result.append(str(div))
places -= 1
integral, first, *others = result
return integral + first + "".join(others).rstrip("0")
def pad_cycles(mod: int, places: int, pairs: Dict[int, str], result: List[str]) -> None:
if mod and places:
i = list(pairs).index(mod)
cycle = "".join(list(pairs.values())[i:])
div, mod = divmod(places, len(cycle))
result.append(cycle * div + cycle[:mod])
def decimalize1(dividend: int, divisor: int, places: int) -> str:
div, mod = divmod(dividend, divisor)
result = [f"{div}."]
pairs = {}
while mod and places and mod not in pairs:
div, mod1 = divmod(mod * 10, divisor)
pairs[mod] = (div := str(div))
result.append(div)
mod = mod1
places -= 1
pad_cycles(mod, places, pairs, result)
integral, first, *others = result
return integral + first + "".join(others).rstrip("0")
They work but both are inefficient, as they are literally doing long division.
They both adhere to the following rules:
They should return the fraction expanded to the desired width, unless the fraction terminates before reaching the width (we had an exact finite decimal representation before reaching the width);
If the last few digits of the decimal expansion cutoff are 0, in which case the result should not contain any trailing zeros, unless the very first decimal digit is 0. For example, it should be '0.0' instead of '0.'
Additionally, the second function breaks out of the loop once it has all digits of one cycle, and constructs the remaining digits by repeating the cycle, though in each iteration more work is done, as a result it is faster if the cycle is short, and slower if the cycle is longer. Worst case is if the cycle is longer than length, so it does all the extra work without terminating early:
In [364]: %timeit decimalize(1, 3, 100)
25.8 μs ± 371 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
In [365]: %timeit decimalize1(1, 3, 100)
2.07 μs ± 20.4 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
In [366]: %timeit decimalize(1, 137, 100)
26.8 μs ± 209 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
In [367]: %timeit decimalize1(1, 137, 100)
4.94 μs ± 38.1 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
In [368]: %timeit decimalize1(1, 1337, 100)
43.4 μs ± 280 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
In [369]: %timeit decimalize(1, 1337, 100)
28.6 μs ± 389 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
In [370]: %timeit decimalize1(1, 123456789, 100)
42.4 μs ± 309 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
In [371]: %timeit decimalize(1, 123456789, 100)
29.7 μs ± 494 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
In [372]: a = 60155040518581045524837379873361090068988594641055829538390886302055144758188941415203591434486458063666229364105014761415461039472408954330541871604108252350317164101172870374427339926789581041
...: 2812627682686305964507416778143771218949050158028407021152173879713433156038667304976240165476457605035649956901133077856035193743615197184
In [373]: b = 19147944100850848776063422241843991195760168200886845084394537367046436869440955641266493717459185828532464222986726583991605539349314420367749162973746417092806627317215443136049103738107006877
...: 6978301192672069310596051608957593418323738306558817176090035871566224143565145070495468977426925354101694409791889484040439128875732421875
In [374]: decimalize(a, b, 101)
Out[374]: '3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798'
In [375]: decimalize1(a, b, 101) == decimalize(a, b, 101)
Out[375]: True
In [376]: %timeit decimalize1(a, b, 101)
96.3 μs ± 295 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
In [377]: %timeit decimalize(a, b, 101)
64.4 μs ± 402 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
How can we do better than long division, so that the execution time is reduced drastically, while achieving the same result (using integer math to get infinite precision representation of fractions)? Preferably this should be done without using any libraries as I would like to know the algorithm.
Just to prove a point:
In [398]: decimalize1(1, 5, 99999999999999999999999999999999999999999)
Out[398]: '0.2'
In [399]: decimalize1(1, 5, 10**1000-1)
Out[399]: '0.2'
In my solution, no matter how many digits are asked, if there is an exact finite decimal representation shorter than asked, it will always be returned, the number of digits asked doesn't affect it.
I chose long division because it is the only thing I know of that works, I want infinite precision, and that method spits out the first N correct decimal places for any N.
I had read this Wikipedia article, but as you can see none of them can be used for my purposes. I want to know what algorithm can be used for infinite precision decimal.
To address some comment, code taken directly from decimalize1:
In [409]: places = 4096
...: div, mod = divmod(1, 1337)
...: result = [f"{div}."]
...: pairs = {}
...: while mod and places and mod not in pairs:
...: div, mod1 = divmod(mod * 10, 1337)
...: pairs[mod] = (div := str(div))
...: result.append(div)
...: mod = mod1
...: places -= 1
In [410]: len(result)
Out[410]: 571
In [411]: len(pairs)
Out[411]: 570
Python has a built-in decimal.Decimal type and you can set the total number of places of precision for operations. I don't know if you would call this "using a library" because the Decimal type is just as fundamental to Python as the int type except implementing its math on CPUs that do not have decimal math capabilities means that a computational algorithm has to be involved. But this is true of the int type also because Python supports essentially infinite precision and that requires multiple-precision algorithms to perform its computations. If you are interested in how multiple-precision arithmetic is done, there are probably many sources of information available to you. As an example I refer to you to Donald Knuth's classic volume on the subject, The Art of Computer Programming: Seminumerical Algorithms, Volume 2.
from decimal import Decimal, getcontext
# Total number of places (allow 1 place for the integer part):
getcontext().prec = 102
n = Decimal('601550405185810455248373798733610900689885946410558295383908863020551447581889414152035914344864580636662293641050147614154610394724089543305418716041082523503171641011728703744273399267895810412812627682686305964507416778143771218949050158028407021152173879713433156038667304976240165476457605035649956901133077856035193743615197184')
d = Decimal('191479441008508487760634222418439911957601682008868450843945373670464368694409556412664937174591858285324642229867265839916055393493144203677491629737464170928066273172154431360491037381070068776978301192672069310596051608957593418323738306558817176090035871566224143565145070495468977426925354101694409791889484040439128875732421875')
assert str(n / d) == '3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798'
print(n / d)
Prints:
3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798
Since getcontext().prec specifies the precision for all digits, then if you are looking for a specific number of digits for the fractional part, you have to first determine how many digits of precision the integer portion of the quotient requires:
from decimal import Decimal, getcontext
n = Decimal('601550405185810455248373798733610900689885946410558295383908863020551447581889414152035914344864580636662293641050147614154610394724089543305418716041082523503171641011728703744273399267895810412812627682686305964507416778143771218949050158028407021152173879713433156038667304976240165476457605035649956901133077856035193743615197184')
d = Decimal('191479441008508487760634222418439911957601682008868450843945373670464368694409556412664937174591858285324642229867265839916055393493144203677491629737464170928066273172154431360491037381070068776978301192672069310596051608957593418323738306558817176090035871566224143565145070495468977426925354101694409791889484040439128875732421875')
integer_places = len(str(int(n / d)))
# Total number of places (allow 101 places for the fractional part):
getcontext().prec = integer_places + 101
assert str(n / d) == '3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798'
print(n / d)