Precision issue in python

[mukesh tiwari]

Hello everyone. I think it is related to the precision with double
arithmetic so i posted here.I am trying with this problem (https://
www.spoj.pl/problems/CALCULAT) and the problem say that "Note : for
all test cases whose N>=100, its K<=15." I know precision of doubles
in c is 16 digits. Could some one please help me with this precision
issue.I used stirling (http://en.wikipedia.org/wiki/
Stirling's_approximation) to calculate the first k digits of N.
Thank you.

__author__="Administrator"
__date__ ="$Feb 19, 2010 3:16:47 PM$"
import math
if __name__ == "__main__":
n=int(raw_input())
while(n>0):
n-=1
raw_input()
l=raw_input();
m=l.split(" ")
N,K,L=int(m[0]),int(m[1]),int(m[2])
fwd,bkd=1,1
s_1,s_2="",""
if(N<=200):
for i in range(1,N+1):
fwd=fwd*i;
s=str(fwd)
s_1=s[:K]
else:
d_1=(N*math.log(N)-N+0.5*math.log(2*math.pi*N)+(1/12/N)-
(1/360/pow(N,3))+(1/1260/pow(N,5))-(1/1680/pow(N,7))+(1/1188/pow(N,9))-
(691/2730/12/11/pow(N,11))+(7/6/14/13/pow(N,13)))*math.log10(math.e)
d_2=d_1-int(d_1)+K-1
fwd=pow(10,d_2)
#print fwd
fwd=int(fwd)
s_1=str(fwd)

if(N<=500):
for i in range(1,N+1):
bkd=bkd*i
bkd=bkd%pow(10,L)
if(bkd==0):
s_2="0"*L
else:
s_2=str(bkd)
else:
s_2="0"*L

print s_1+" "+s_2

 

[Mark Dickinson]

Hello everyone. I think it is  related to the precision with double
arithmetic so i posted here.I am trying with this problem (https://www.spoj.pl/problems/CALCULAT) and the problem say that "Note : for
all test cases whose N>=100, its K<=15." I know precision of doubles
in c is 16 digits. Could some one please help me with this precision
issue.I used stirling (http://en.wikipedia.org/wiki/
Stirling's_approximation) to calculate the first k digits of N.
Thank you.

If I understand you correctly, you're trying to compute the first k
digits in the decimal expansion of N!, with bounds of k <= 15 and 100
<= N < 10**8. Is that right?

Python's floats simply don't give you enough precision for this:
you'd need to find the fractional part of log10(N!) with >= 15 digits
of precision. But the integral part needs ~ log10(log10(N!)) digits,
so you end up needing to compute log10(N!) with at least 15 +
log10(log10(N!)) ~ 15 + log10(N) + log10(log10(N)) digits (plus a few
extra digits for safety); so that's at least 25 digits needed for N
close to 10**8.

The decimal module would get you the results you need (are you allowed
imports?). Or you could roll your own log implementation based on
integer arithmetic.

 

[mukesh tiwari]

If I understand you correctly, you're trying to compute the first k
digits in the decimal expansion of N!, with bounds of k <= 15 and 100
<= N < 10**8.  Is that right?

Python's floats simply don't give you enough precision for this:
you'd need to find the fractional part of log10(N!) with >= 15 digits
of precision.  But the integral part needs ~ log10(log10(N!)) digits,
so you end up needing to compute log10(N!) with at least 15 +
log10(log10(N!)) ~ 15 + log10(N) + log10(log10(N)) digits (plus a few
extra digits for safety);  so that's at least 25 digits needed for N
close to 10**8.

The decimal module would get you the results you need (are you allowed
imports?).  Or you could roll your own log implementation based on
integer arithmetic.

Yes i am trying to first k digits of N!.I will try your method.
Thank you

 

[mukesh tiwari]

Yes i am trying to first k digits of N!.I will try your method.
Thank you

I am out of luck.I put d_2=d_1-int(d_1)+25 instead of d_2=d_1-
int(d_1)+K-1 but i am getting WA.
"The decimal module would get you the results you need (are you
allowed
imports?). Or you could roll your own log implementation based on
integer arithmetic. "
I don't know if is possible to import this decimal module but kindly
tell me.Also a bit about log implementation
Thank you

 

[Mark Dickinson]

I don't know if is possible to import this decimal module but kindly
tell me.Also a bit about log implementation

The decimal module is part of the standard library; I don't know what
the rules are for SPOJ, but you're already importing the math module,
so at least *some* imports are obviously permitted. Here's a quick
demonstration of how you might use the decimal module: you can
probably find ways to tweak it for speed and accuracy.

from decimal import Decimal as D
from decimal import getcontext

getcontext().prec = 100 # working precision = 100 digits

pi =
D('3.14159265358979323846264338327950288419716939937510582097494459'

'2307816406286208998628034825342117067982148086513282306647093844')
half = D('0.5')
log = D.ln

def logfac(x):
"""Approximation to log(x!), using first few terms of Stirling's
series."""

x = D(x)
return log(2*pi)/2 + (x + half)*log(x) - x + \
1/(12*x) - 1/(360*x**3) + 1/(1260*x**5)

def fac_first_digits(n, k):
"""Get first k decimal digits of n!."""

log10_nfac = logfac(n)/log(D(10))
frac = log10_nfac - int(log10_nfac)
return int(10**(frac - 1 + k))


With the above code I get, for example (with Python 2.6):

fac_first_digits(12345, 15) 344364246918678
from math import factorial
int(str(factorial(12345))[:15])

344364246918678

For small n, you'll need more terms of Stirling's series than I've
used above. And there's always a small chance that the intermediate
errors involved in the computation might change those first 15
digits; with a working precision of 100 and lots of terms of
Stirling's series it's a *very* small chance, but it's there. If you
want something 100% watertight, you'd need to compute strict upper
and lower bounds for the true result, and increase precision until
those bounds match sufficiently far that you can be sure of the first
k digits being correct.

 

[Mark Dickinson]

A quick solution I came out with, no stirling numbers and had tried to avoid
large integer multiplication as much as possible.

import math

for i in range(int(raw_input())):
    n, k, l = [int(i) for i in raw_input().split()]
    e = sum(math.log10(i) for i in range(1, n+1))
frac_e = e - math.floor(e)

This isn't going to give you enough accuracy when n gets large (and
the original problem seems to allow n to be as large as 10**8), for
the reasons I described earlier---that is, Python floats only give you
around 16 decimal digits of precision; your 'e' value will already
have used up some of those 16 digits before the point, leaving you
fewer than 16 digits of precision after the point, so the absolute
error in frac_e will likely be larger than 1e-15. That's not good
enough for getting the first 15 digits of 10**frac_e accurately. For
large n, you'll also be accumulating significant error in the sum.

a = str(10**frac_e * 10**(k - 1)).split('.')[0]

The str(...).split('.') here doesn't do a good job of extracting the
integer part when its argument is >= 1e12, since Python produces a
result in scientific notation. I think you're going to get strange
results when k >= 13.