Problem extracting base B fractional digits of a number

[Georgi Guninski]

Given SR constant, I want to extract the first L base B digits
of the fractional part.

So far the best solution I found is:
floor((SR(expression)*B**L).numerical_approx(digits=L)).digits(B)

Can I do better?

Why the following approach fails numerically:

```
def base_B_digits(A, B,prot=False):
digits=[]
fractional_part = A - int(A)
while fractional_part != 0:
digit = int(fractional_part * B)
digits.append(digit)
fractional_part = fractional_part * B - digit
if prot: print(fractional_part)
return digits
```
tt=base_B_digits(SR(1/3).numerical_approx(digits=10),10,1)

The fractional part starts:

0.3333333333
0.3333333331
0.3333333314
0.3333333139
0.3333331393
...
0.7812500000
0.8125000000
0.1250000000
0.2500000000
0.5000000000
0.0000000000

I am not sure the program must terminate.

For A=1/4, the base 10 digits are computed correctly.

[Gareth Ma]

I assume you should use floor division to preserve accuracy. As a numerical example, say you want the first 100 decimal digits of 3 / 173. That's the same as floor(10^100 * 3 / 173).

In [105]: from decimal import Decimal, getcontext
     ...: getcontext().prec = 300
     ...: print(10**100 * 91 // 173)
     ...: print(int(str(Decimal(91) / Decimal(173)).split(".")[1][:100]))
5260115606936416184971098265895953757225433526011560693641618497109826589595375722543352601156069364
5260115606936416184971098265895953757225433526011560693641618497109826589595375722543352601156069364

Gareth

[Georgi Guninski]

On Sat, Feb 10, 2024 at 6:24 PM Georgi Guninski <ggun...@gmail.com> wrote:
>
>
> ```
> def base_B_digits(A, B,prot=False):
> digits=[]
> fractional_part = A - int(A)
> while fractional_part != 0:
> digit = int(fractional_part * B)
> digits.append(digit)
> fractional_part = fractional_part * B - digit
> if prot: print(fractional_part)
> return digits
> ```
>

> I am not sure the program must terminate.
>

Just for the record, this this digit extraction algorithm doesn't terminate
for A=0.5 and B=3 in sage and pari/gp.

[Gareth Ma]

Yes, because 1 / 2 is mathematically infinite in base 3, just like how A=1/3 and B=10 doesn't terminate.

[Georgi Guninski]

On Mon, Feb 12, 2024 at 9:35 AM Gareth Ma <grh...@gmail.com> wrote:
>
> Yes, because 1 / 2 is mathematically infinite in base 3, just like how A=1/3 and B=10 doesn't terminate.
>

Indeed, it is periodic and infinite.
When I use code finite precision I expect it to terminate.
It works with python `float` too.
btw, this code was written by chatgpt and I have seen it on
other places.