Десмос , 62 59 байт
Это порт
Решение Кристиана Сивера с некоторыми изменениями для обеспечения числовой точности..
Попробуйте это на Десмосе!
-3 байта благодаря
Эйден Чоу
Последний столбец таблицы показывает относительную ошибку, которая достигает 5,71×10^-13 в худшем случае при c=0,2. Это едва соответствует требованию относительной ошибки 10^-12.
0.2
can be changed to from functools import*
#let's start with the comments. I'll exclude them from total count
#from functools I need only reduce
#I've been coding in python2 for 5 years so I'm used to built-in reduce much
#and it's very handy here
g=range
#just keep it in mind while reading) ranGe -> g
#and we're getting rid of classes in favor of functions, that gives
#us ~100 bytes more. It's how it was:
#class P:
#it's a Polynomial, self.c are coefficients and self.d is a divisor. So all polynomial computations are exact
#we need polynomials to implement https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss%E2%80%93Legendre_quadrature
# c=[0]
# d=1
#don't even know if these initial values are necessary
# def __init__(self,c,d=1):
# self.c=list(c)
# self.d=d
#n'th derivative
# def dn(self,n=1):
# c=self.c[:]
# for _ in g(n):
# c=[c[i]*i for i in g(1,len(c))]
# d=reduce(gcd,c,self.d)
#gcd is imported from math along with exp,log and gamma
# return P([i//d for i in c],self.d//d)
# def m(self,x):
#the value of poly at point x. __call__ may look better
# return sum((x**i)*self.c[i] for i in range(len(self.c)-1,-1,-1))/self.d
#the reverse order was for floats and really needed there for precision,
#but python3 have exact Fractions, so maybe it's to be golfed too
#Now let's convert this class to functions
#derivative=lambda c:reduce(lambda c,_:[c[i]*i for i in g(1,len(c))],g(n),c)
#the_gcd=lambda c,d:reduce(gcd,c,d)
#the_result=lambda c,d,m:([i//m for i in c],d//m)
#where c is now the derivative
#and combined:
v=lambda c,d,n=1:(lambda c:
(lambda c,d,m:([i//m for i in c],d//m))(c,d,reduce(gcd,c,d))
)\
(reduce(lambda c,_:[c[i]*i for i in g(1,len(c))],g(n),c))
#the computing of a plynomial's value:
m=lambda c,d,x:sum((x**i)*c[i] for i in g(len(c)-1,-1,-1))/d
d={}
#https://en.wikipedia.org/wiki/Combination written recursively with Pascal's triangle, with the dict d for speed
def c(n,k):
if (n,k) not in d:
d[(n,k)]=1 if k==0 or k==n or n==0 else c(n-1,k)+c(n-1,k-1)
return d[(n,k)]
#p=30
#it's the floating point precision for Decimal
#it was needed for computing the polynomial roots, you can see the raw version
#maybe I'll get rid of it if it's used in the only place
from math import*
from fractions import Fraction as F
n=61
l=v(sum(([c(n,i)*(-1)**(n-i),0] for i in g(n+1)),[]),
2**n*reduce(lambda x,y:x*y,g(1,n+1),1),n)
#Legendre polynomial of order n, pretty straightforward, -3bytes for space
#you can't see these bytes in the tio version, there is 1 line
d=v(*l)
#and the derivative for weights
#see we can use c and d now 'cause we don't need them more
c=2**100
#maybe the hardest-to-read part, the roots of l
#taking too long to compute (2-3 mins) so hard-stored
#b(b'GKSmRyJoJz...') is the 390-bytes byte array,
#representing the <0 roots (the function is odd, so -x and x are both roots and we have to store only half)
#30 13-byte(100-bit, for 10^{-30} precision) integers are stored in lsbf order of bytes
#and c is common multiplier for all computed roots' denominators
from base64 import b64decode as b
#extract numbers from base-64,
r=(lambda b:[sum(2**(n*8)*i for n,i in enumerate(b[j:j+13])) for j in g(0,391,13)])(b(b'GKSmRyJoJzZNd978D7hBJ/SAXF6k3Q+D7w9eWlMN1hwU7v+nhNcPP3k+PPY3ALm75/G0DwxkhnBW7vYeTLvhhw8XTAWijBbPeSIxclAPxUQBLJIlU0r8S8gODz6W/ntGb0p0WeUPww6bS5oF+R0NzuqOe20OpkOWAx5DBcFecEQODqjUMBgWzmCaCiGqpQ2xekr44w56Xl998jMNN0/S8UF6TE04eGm5DL0i6Wm5wA/nHOhgNgxOR3kgsJNHFpZQMKsLdK2AMa1xRv65pzQYCxSYft+TX736FBjQfQqY8CuipXP9Bxm/adwJEKYLfykKeTs/aG00CVuQ6lnO2j04CUVLhggLyUIxC/eIqRKid9IHdSe6aSQRS8hkmWoZB2rhfvNKb9bHP8KfWwb379JeWbyiuY/elZkFsIKVLmMYZwdEhs7TBAvB47snQ/UMwtDNCgRb9soEqlNzn6z8GT8DTfSlRAVDCG07FjtxAheAQ01VIxgWXpy6oQEELmiRzDIqr+ckI9EA'))
#and make fractions
#0 is in r for magically mistaken 391 for 390 )
#maybe +[0]+r[::-1] is 1 byte shorter (or longer?)
#but at least more readable )
#roots computation code is in the raw below
r=[F(i,c) for i in [-j for j in r]+r[-2::-1]]
#the weights
w=[2/((1-x**2)*(m(*d,x))**2) for x in r]
from decimal import*
D=Decimal
setcontext(Context(prec=30, rounding=ROUND_HALF_DOWN))
#that was float=
#it's only converting Fraction to Decimal. no built-in
f=lambda r:(lambda n,d:D(n)/D(d))(*str(r).split('/'))
#number of sub-intervals of (0,1)
d=2**8
s=input().strip()
c=D(s)
#the inverse of c, not the math.exp(1) ))
e=D(1)/c
#pretty straightforward, I is the integrator,
I=lambda h:sum(sum(f(v)*h(f(((x+1)/2+i)/d))/(2*d) for v,x in zip(w,r))
for i in range(d))
#the rest follows (at math:)
print(D(gamma(1/float(s)+1))+(c.ln())-I(lambda t:(D('0.5')+5*((-9*((-((t).ln()))**(e))).exp())).ln())-D('0.5772156649015328606065120900824')*(1-e)-1)
Что касается требований к скорости, на моем компьютере это занимает менее 1500 мс.