Re^7: Introspection into floats/NV (math behind periodic numbers)

[LanX]

> I think there is a generic formula to calculate these using prime factors.

I looked into it, there is a trivial and less trivial part to it.

(I'm writing it down for therapeutically reasons, because wanna stop reinvent Euler's wheels again)

for shortened p/q* and with q* =  q * 2^i for i in N₀, and p,q in N.

(i.a.W. ignore prime factors of 2)

The period is the smallest k for which there is a factor l such that q*l= 2^k -1 (aka ord₂(q) said order of 2 modulo q)

The periodic fraction will be p*l , because 1/ 2^k-1 is a periodic fraction of the length k and sequence '0...01'. (that's similar to 1/9, 1/99, ... in decimal, same math, we did this in school)

Example 1/5 = 3/15 =  3/2**4-1 = 0.[0011] because 1/15 = 0.0001

  DB<15> $, =","; say unpack "A1A11A*",unpack "B*", pack "F>", 1/5
0,01111111100,1001100110011001100110011001100110011001100110011010

  DB<16> $, =","; say unpack "A1A11A*",unpack "B*", pack "F>", 1/15
0,01111111011,0001000100010001000100010001000100010001000100010001

NB: The exponents are different to normalize the mantissa.

That was the trivial part. ;-)

The rest dives into number theory

for a concrete formula° to deduce k from q one needs q's prime factors and Eulers phi function (it eludes me why this is called Totient in English)

for many - not all - prime numbers q holds k = q-1 (e.g 3, 5, but not 7)

In the case of 1/121 (NB 121=11**2) we have a period of 110 = 10 ** 11 digits, easily cracking the 53 bits of a double float.

This is also exemplifies why rational numbers can't be fixed by adding the period length to floats, storing 1 and 121 is far more precise. OTOH floats are capable of approximating many non-rational numbers like pi.

DISCLAIMER: not a mathematical text, I wrote it down as fast as possible, here are typos, dragons and typodragons. ;-)

°) There are detailed Math Exchange on the topic see https://math.stackexchange.com/questions/1025578/is-there-a-better-way-of-finding-the-order-of-a-number-modulo-n#1025593

Update

here a brutal but O(n) efficient way to calculate it even for very long periods:

sub period_length { 
    my $q = shift ; 

    until ( $q % 2 ) { $q /= 2 }; 
    my $r=1; 
    # say "q:$q"; 
    for (1..$q) { 
        $r = $r * 2 % $q; 
        return $_ if $r == 1; 
        # say "r:$r" 
    } 
}