in reply to Re^6: Introspection into floats/NV
in thread Introspection into floats/NV
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 N0, 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 ord2(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. ;-)
Cheers Rolf
(addicted to the Perl Programming Language :)
see Wikisyntax for the Monastery
°) 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
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" } }
|
---|