10
\$\begingroup\$

WARNING: This challenge may need 128 bit floats.1

The task is to perform numerical integration. Consider the following three functions.

\$ f(x) = cx^{c - 1}e^{-x^c} \$

\$ g_1(x) = 0.5e^{-x} \$

\$ g_2(x) = 5 e^{-10 x} \$

We will have that \$c \geq 0.2\$. Your code should be correct for any value between 0.2 and 1.0.

The task is to perform the integral:

\$ \int_{0}^{\infty} f(x) \ln\left(\frac{f(x)}{g_1(x)+g_2(x)}\right) \mathrm{d}x. \$

Your output should be correct to within plus or minus \$10^{-10}\$ percent of the correct answer.

Input

A single real number \$c\$.

Output shown to 15 significant places

c = 0.2. Output: 119.2233798550179
c = 0.3. Outout: 8.077346771987397
c = 0.4. Output: 2.07833944013377
c = 0.5. Output: 0.8034827042913205
c = 0.6. Output: 0.3961739639940003
c = 0.7. Output: 0.2585689391629249
c = 0.8. Output: 0.2287758419066705
c = 0.9. Output: 0.2480070283065089
c = 1.0. Output: 0.2908566108890576

Timing

Your code must complete on TIO before it times out.

As is usual, this is a competition per language. In other words, Java coders should not be put off by the single character solution in a language they can't believe exists. You may use any language/library you choose.

Questions and notes

  • Does this integral actually exist/converge?

Yes. It is the Kullback Leibler divergence between two probability distributions.

  • Show me any code that can compute these integrals.

See this Wolfram Alpha example.

  • Why such high precision that we need 128 bit floats?

It is quite easy to solve the integral to low precision using simple ideas. However, to get high precision seems to need careful thought both about the integral itself but also the quadrature method that will be used. I felt this made the question more challenging and interesting.

1. For lots of languages this is easy. E.g Fortran has real128, gcc/clang have _Float128 for C/C++, Python has Decimal (and many other options), Dyalog APL has 128 bit (decimal) floats, C# has the Decimal type, Haskell has rounded, PHP has decimal, MATLAB has the Multiprecision Computing Toolbox, Julia has BigFloat, Java has BigDecimal, JavaScript has big.js, Ruby has BigDecimal, R has Rmpfr, Perl has Math::BigFloat, SOGL has 200 digit decimal precision, etc. If your favourite golfing language doesn't, maybe consider adding it?

Update

Alexey Burdin has shown that

\$ \int_{0}^{\infty} f(x) \ln\left(\frac{f(x)}{g_1(x)+g_2(x)}\right) \mathrm{d}x = \int\limits_0^1 \left(\ln c+\left(\ln \left(-\ln t\right)\right)\cdot\frac{c-1}{c}+ \ln t- \ln\left(0.5+5e^{-9x(t)}\right) +(-\ln{t})^{1/c}\right)dt\$

This potentially greatly simplifies the problem.

\$\endgroup\$
  • 6
    \$\begingroup\$ Because of the very high precision requirements here, I think the answer to "How compactly can your language perform numerical integration?" is largely based on "Does your language have a high-precision real-number representation?" \$\endgroup\$ – xnor Nov 2 at 19:29
  • \$\begingroup\$ @xnor It's a little more difficult than that :) But I think most (all?) main stream languages have at least libraries for 128 bit floats. \$\endgroup\$ – Anush Nov 2 at 19:40
  • 5
    \$\begingroup\$ Many answers on this site are not written in mainstream languages. \$\endgroup\$ – pppery Nov 2 at 19:59
  • 6
    \$\begingroup\$ @pppery that is true. But for every challenge there are a set of languages most suited and a set that are ill suited. Also the challenge is per language so even if you have to implement 128 bit multiplication/division yourself in your language, that's ok. \$\endgroup\$ – Anush Nov 2 at 20:08
4
\$\begingroup\$

Python 3, 150 133 122 bytes

(thanks to Delta and Alexey Burdin for further golfing this)

SymPy's mpmath module is a great way to handle this problem. All math functions (exp, log, euler, gamma and quad for integration) are run in high precision.

lambda c:gamma(1/c+1)+l(c)-euler*(1-1/c)-1-quad(lambda t:l(.5+5*exp(-9*((-l(t))**(1/c)))),[0,1])
from mpmath import*
l=log

Try it online!

\$\endgroup\$
  • \$\begingroup\$ This is great. Thank you. \$\endgroup\$ – Anush Nov 7 at 7:21
  • 2
    \$\begingroup\$ Try it online! you can save some bytes if you use a lamda instead of the input() \$\endgroup\$ – Delta Nov 7 at 13:48
  • \$\begingroup\$ @Delta thanks! I thought we needed to output a fixed number of digits, thus I hang on to the nprint, but indeed that is not required. That saves 17 bytes. Nice one. \$\endgroup\$ – agtoever Nov 7 at 21:00
  • 1
    \$\begingroup\$ 128 bytes \$\endgroup\$ – Alexey Burdin Nov 8 at 23:51
  • \$\begingroup\$ @AlexeyBurdin thanks! I further golfed it to 122 bytes. \$\endgroup\$ – agtoever Nov 9 at 12:04
4
+50
\$\begingroup\$

python3, 3982 2183 1866 1821 1677, ~14s

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)

Try it online! -- takes input as a valid Decimal, like 0.2
Math:
0. We do google numerical integration.
1. Let's see provocative \$f(x)=c\cdot x^{c-1}\cdot e^{-x^c}= -\left(e^{-x^c}\right)'\$
2. So, we have \$-\int\limits_0^\infty \ln\left(\frac{c\cdot x^{c-1}\cdot e^{-x^c}}{0.5e^{-x}+5e^{-10x}}\right)\left(e^{-x^c}\right)'dx\$,
we change the variable to \$t=e^{-x^c}\$, \$x^c=-\ln(t)\$, \$x=\left(-\ln t\right)^\frac1c\$,
\$t(0)=1\$, \$t(+\infty)=0\$ and the integral is
\$\int\limits_0^1 \ln\left(\frac{c\cdot x(t)^{c-1}\cdot t}{\left(0.5+5e^{-9x(t)}\right)e^{-x(t)}}\right)dt=\$
\$\int\limits_0^1 \left(\ln c+\left(\ln x(t)\right)\cdot(c-1)+ \ln t- \ln\left(0.5+5e^{-9x(t)}\right) +x(t)\right)dt=\$
\$\int\limits_0^1 \left(\ln c+\left(\ln \left(-\ln t\right)\right)\cdot\frac{c-1}{c}+ \ln t- \ln\left(0.5+5e^{-9x(t)}\right) +x(t)\right)dt\$
3. Next we wolframalpha these things:
\$\int\limits_0^1\ln(-\ln t)dt\$, \$\int\limits_0^1 (-\ln t)^a dt\$, \$\int\limits_0^1 \ln t\,dt\$ and do the rest numerically (the latter thing maybe can be done with trapeziums or rectangles or Simpson -- the idea for short submissions).
The raw thing:

from copy import copy
from functools import reduce
class Poly:
    c=[0]
    d=1
    def __init__(self,c,d=1):
        self.c=list(c)
        self.d=d
    def dn(self,n=1):
        c=copy(self.c)
        for _ in range(n):
            c=[c[i]*i for i in range(1,len(c))]
        d=reduce(gcd,c,self.d)
        return Poly([i//d for i in c],self.d//d)
    def m(self,x):
        return sum((x**i)*self.c[i] for i in range(len(self.c)-1,-1,-1))/self.d
def gcd(x,y):
    return abs(gcd(y,x%y)) if y else x
cdict={}
def c(n,k):
    if (n,k) not in cdict:
        cdict[(n,k)]=1 if k==0 or k==n or n==0 else c(n-1,k)+c(n-1,k-1)
    return cdict[(n,k)]
def root(x,f,f_):
    x_old=0
    while x!=x_old:# and abs(f(x))>1e-15:
        print(x,f(x),f_(x))
        x_old=x
        x=x-f(x)/f_(x)
    print('===')
    return x
precision=30
def b(b,e,f,eps=1e-17):
    mold,m=0,1
    while e-b>eps and mold!=m:
        mold=m
        m=(e+b)/2
        if (l.m(m)>0)!=(l.m(b)>0):
            e=m
        else:
            b=m
        #print(*map(float,[b,e,f(b),f(e)]))
    print('.',sep='',end='')
    return b
##class Frac:
##    n=0
##    d=1
##    def __init__(self,n,d):
##        self.n=n
##        self.d=d
##    def
lcm=lambda x,y:x*y//gcd(x,y)
import math
from fractions import Fraction
degree=61
n=degree
l=Poly(sum(([c(n,i)*(-1)**(n-i),0] for i in range(n+1)),[]),
       2**n*reduce(lambda x,y:x*y,range(1,n+1),1)).dn(n)
ld=l.dn()
#print(l.coeffs,l.d)
#r=[root(math.cos(math.pi*(4*i-1)/(4*n+2)),l.m,ld.m) for i in range(1,n+1)]
##d=2**8#4096
##while True:
##    i_=[(b,e) for b,e in [(Fraction(i,d)-1,Fraction(i+1,d)-1)
##                          for i in range(2*d)] if (l.m(b)>0)!=(l.m(e)>0)]
##    print(d,len(i_))
##    if len(i_)==n:
##        break
##    d*=2
##r=[b(*i,l.m,Fraction(1,10**precision)) for i in i_]
##print('\n',sep='',end='')
cm=1267650600228229401496703205376
#cm=reduce(lambda x,i:lcm(int((str(i).split('/')+[1])[1]),x),r,1);cm
#[str(i*cm) for i in r]
r=[Fraction(int(i),cm) for i in ['-1266681605106811426160584860697', '-1262547799267707740863710314936', '-1255122086390021131342967954015', '-1244422184933863703004183427391', '-1230475806835976018691325584397', '-1213319288075032955237836344344', '-1192997371846808208190052844742', '-1169563069074221666262667138623', '-1143077514053466175428649175963', '-1113609802945881973204995949479', '-1081236812728978305124643427497', '-1046043000289730543348383120049', '-1008120181921745998647915597624', '-967567293702390287305984778942', '-924490133333634877286205966159', '-879001084101412705885157895540', '-831218821664589525108846925845', '-781268004433877062952267018393', '-729278948345971539276192261648', '-675387286880006687633990914140', '-619733617202494872284783888651', '-562463133363442075271048079222', '-503725247500289525435137450347', '-443673200037730123766618583031', '-382463659900222496924931293872', '-320256315780122313563780137228', '-257213459527711400732339664476', '-193499562749975057795207525453', '-129280847722706546595631169560', '-64724853735360852478581550597', '0', '64724853735360852478581550596', '129280847722706546595631169559', '193499562749975057795207525452', '257213459527711400732339664475', '320256315780122313563780137227', '382463659900222496924931293871', '443673200037730123766618583030', '503725247500289525435137450346', '562463133363442075271048079221', '619733617202494872284783888650', '675387286880006687633990914139', '729278948345971539276192261647', '781268004433877062952267018392', '831218821664589525108846925844', '879001084101412705885157895539', '924490133333634877286205966158', '967567293702390287305984778941', '1008120181921745998647915597623', '1046043000289730543348383120048', '1081236812728978305124643427496', '1113609802945881973204995949478', '1143077514053466175428649175962', '1169563069074221666262667138622', '1192997371846808208190052844741', '1213319288075032955237836344343', '1230475806835976018691325584396', '1244422184933863703004183427390', '1255122086390021131342967954014', '1262547799267707740863710314935', '1266681605106811426160584860696']]
#r=[Fraction(-72002513042367095, 72057594037927936), Fraction(-35883766692782835, 36028797018963968), Fraction(-35672715239059617, 36028797018963968), Fraction(-141474422995755015, 144115188075855872), Fraction(-69944451686068079, 72057594037927936), Fraction(-137938433007892963, 144115188075855872), Fraction(-135628098615462501, 144115188075855872), Fraction(-66481963419518119, 72057594037927936), Fraction(-32488216960856169, 36028797018963968), Fraction(-126602776952709289, 144115188075855872), Fraction(-61461197033678349, 72057594037927936), Fraction(-59460660411885473, 72057594037927936), Fraction(-7163124720376473, 9007199254740992), Fraction(-109999665903886841, 144115188075855872), Fraction(-105102359763536115, 144115188075855872), Fraction(-99930853605361871, 144115188075855872), Fraction(-23624659822433771, 36028797018963968), Fraction(-88819888837164975, 144115188075855872), Fraction(-41454708727199633, 72057594037927936), Fraction(-19195661220692463, 36028797018963968), Fraction(-70455555169530143, 144115188075855872), Fraction(-63944654967081315, 144115188075855872), Fraction(-28633465234423143, 72057594037927936), Fraction(-6304975386764939, 18014398509481984), Fraction(-10870271009372985, 36028797018963968), Fraction(-36408927801417385, 144115188075855872), Fraction(-29241784833142379, 144115188075855872), Fraction(-21998353389560073, 144115188075855872), Fraction(-14697530755564293, 144115188075855872), Fraction(-7358363943167303, 144115188075855872), Fraction(0, 1), Fraction(3679181971583651, 72057594037927936), Fraction(3674382688891073, 36028797018963968), Fraction(2749794173695009, 18014398509481984), Fraction(14620892416571189, 72057594037927936), Fraction(4551115975177173, 18014398509481984), Fraction(43481084037491939, 144115188075855872), Fraction(50439803094119511, 144115188075855872), Fraction(57266930468846285, 144115188075855872), Fraction(31972327483540657, 72057594037927936), Fraction(35227777584765071, 72057594037927936), Fraction(76782644882769851, 144115188075855872), Fraction(82909417454399265, 144115188075855872), Fraction(44409944418582487, 72057594037927936), Fraction(94498639289735083, 144115188075855872), Fraction(49965426802680935, 72057594037927936), Fraction(52551179881768057, 72057594037927936), Fraction(13749958237985855, 18014398509481984), Fraction(114609995526023567, 144115188075855872), Fraction(118921320823770945, 144115188075855872), Fraction(122922394067356697, 144115188075855872), Fraction(15825347119088661, 18014398509481984), Fraction(129952867843424675, 144115188075855872), Fraction(132963926839036237, 144115188075855872), Fraction(33907024653865625, 36028797018963968), Fraction(68969216503946481, 72057594037927936), Fraction(139888903372136157, 144115188075855872), Fraction(70737211497877507, 72057594037927936), Fraction(142690860956238467, 144115188075855872), Fraction(143535066771131339, 144115188075855872), Fraction(144005026084734189, 144115188075855872)]
#r=[Fraction(-144104924687379029, 144115188075855872), Fraction(-144061113633417193, 144115188075855872), Fraction(-71991152681636629, 72057594037927936), Fraction(-143868501524294767, 144115188075855872), Fraction(-35929931908622757, 36028797018963968), Fraction(-143536019344824829, 144115188075855872), Fraction(-71658710569760757, 72057594037927936), Fraction(-143063986084641675, 144115188075855872), Fraction(-142775775751233995, 144115188075855872), Fraction(-35613215044775229, 36028797018963968), Fraction(-142095317851501427, 144115188075855872), Fraction(-141703235672459401, 144115188075855872), Fraction(-141276708943911307, 144115188075855872), Fraction(-70407920670826367, 72057594037927936), Fraction(-140320744889653401, 144115188075855872), Fraction(-34947884983137447, 36028797018963968), Fraction(-139228355106226913, 144115188075855872), Fraction(-138631327306447751, 144115188075855872), Fraction(-69000300827755959, 72057594037927936), Fraction(-68668165733468573, 72057594037927936), Fraction(-136638678208163739, 144115188075855872), Fraction(-135907811461286037, 144115188075855872), Fraction(-67571954440908669, 72057594037927936), Fraction(-134347156155496359, 144115188075855872), Fraction(-16689718369143059, 18014398509481984), Fraction(-132655882883583863, 144115188075855872), Fraction(-131761773444627501, 144115188075855872), Fraction(-65417817986076227, 72057594037927936), Fraction(-64938847793634331, 72057594037927936), Fraction(-32222046285398959, 36028797018963968), Fraction(-127867345160661637, 144115188075855872), Fraction(-31703855946358703, 36028797018963968), Fraction(-62866338356003711, 72057594037927936), Fraction(-15577420891180097, 18014398509481984), Fraction(-123475765655790143, 144115188075855872), Fraction(-61151075136161875, 72057594037927936), Fraction(-121098806255952059, 144115188075855872), Fraction(-119866026109883705, 144115188075855872), Fraction(-118604109492524957, 144115188075855872), Fraction(-117313363144639965, 144115188075855872), Fraction(-115994100814789491, 144115188075855872), Fraction(-28661660795766563, 36028797018963968), Fraction(-28317829445786347, 36028797018963968), Fraction(-111868458922669023, 144115188075855872), Fraction(-110438407601984251, 144115188075855872), Fraction(-27245377857813575, 36028797018963968), Fraction(-26874531136490757, 36028797018963968), Fraction(-105988607520833283, 144115188075855872), Fraction(-26113331820545011, 36028797018963968), Fraction(-102892657018719659, 144115188075855872), Fraction(-50653488045428427, 72057594037927936), Fraction(-49848334969235775, 72057594037927936), Fraction(-98062129987227927, 144115188075855872), Fraction(-96403753553428489, 144115188075855872), Fraction(-94721943747436259, 144115188075855872), Fraction(-23254277343922147, 36028797018963968), Fraction(-45644832420663197, 72057594037927936), Fraction(-89540030043462969, 144115188075855872), Fraction(-10971078784389609, 18014398509481984), Fraction(-85975896119833679, 144115188075855872), Fraction(-42081131673510875, 72057594037927936), Fraction(-20582043201506859, 36028797018963968), Fraction(-40237035159487741, 72057594037927936), Fraction(-39300203286200335, 72057594037927936), Fraction(-19176909251924261, 36028797018963968), Fraction(-74796221710411347, 144115188075855872), Fraction(-72866625298407583, 144115188075855872), Fraction(-35459658404464939, 72057594037927936), Fraction(-68954769584591107, 144115188075855872), Fraction(-33486730579157493, 72057594037927936), Fraction(-32487936568629803, 72057594037927936), Fraction(-31481245542875313, 72057594037927936), Fraction(-60933804407252565, 144115188075855872), Fraction(-58890306225406907, 144115188075855872), Fraction(-56832493264165915, 144115188075855872), Fraction(-54760865727051299, 144115188075855872), Fraction(-26337963587783547, 72057594037927936), Fraction(-25289092203398143, 72057594037927936), Fraction(-48468147330210951, 144115188075855872), Fraction(-46346328843725859, 144115188075855872), Fraction(-11053311177256411, 36028797018963968), Fraction(-21034706713097939, 72057594037927936), Fraction(-19957678053844251, 72057594037927936), Fraction(-4718949543956531, 18014398509481984), Fraction(-35578660114658865, 144115188075855872), Fraction(-33397075583856055, 144115188075855872), Fraction(-15603686524289111, 72057594037927936), Fraction(-29010084771446223, 144115188075855872), Fraction(-1675359053686717, 9007199254740992), Fraction(-24594889131807827, 144115188075855872), Fraction(-22378054994346827, 144115188075855872), Fraction(-20155781304247947, 144115188075855872), Fraction(-4482152060343901, 36028797018963968), Fraction(-15697077176510773, 144115188075855872), Fraction(-6730865269878557, 72057594037927936), Fraction(-5611555844344803, 72057594037927936), Fraction(-4490882388138869, 72057594037927936), Fraction(-6738234618615355, 144115188075855872), Fraction(-2246533281244659, 72057594037927936), Fraction(-2246806352819167, 144115188075855872), Fraction(0, 1), Fraction(1123403176409583, 72057594037927936), Fraction(4493066562489317, 144115188075855872), Fraction(3369117309307677, 72057594037927936), Fraction(8981764776277737, 144115188075855872), Fraction(11223111688689605, 144115188075855872), Fraction(13461730539757113, 144115188075855872), Fraction(3924269294127693, 36028797018963968), Fraction(17928608241375603, 144115188075855872), Fraction(10077890652123973, 72057594037927936), Fraction(11189027497173413, 72057594037927936), Fraction(12297444565903913, 72057594037927936), Fraction(26805744858987471, 144115188075855872), Fraction(14505042385723111, 72057594037927936), Fraction(31207373048578221, 144115188075855872), Fraction(16698537791928027, 72057594037927936), Fraction(2223666257166179, 9007199254740992), Fraction(37751596351652247, 144115188075855872), Fraction(39915356107688501, 144115188075855872), Fraction(42069413426195877, 144115188075855872), Fraction(44213244709025643, 144115188075855872), Fraction(23173164421862929, 72057594037927936), Fraction(24234073665105475, 72057594037927936), Fraction(50578184406796285, 144115188075855872), Fraction(52675927175567093, 144115188075855872), Fraction(27380432863525649, 72057594037927936), Fraction(28416246632082957, 72057594037927936), Fraction(29445153112703453, 72057594037927936), Fraction(15233451101813141, 36028797018963968), Fraction(62962491085750625, 144115188075855872), Fraction(64975873137259605, 144115188075855872), Fraction(66973461158314985, 144115188075855872), Fraction(34477384792295553, 72057594037927936), Fraction(70919316808929877, 144115188075855872), Fraction(36433312649203791, 72057594037927936), Fraction(37398110855205673, 72057594037927936), Fraction(76707637007697043, 144115188075855872), Fraction(78600406572400669, 144115188075855872), Fraction(80474070318975481, 144115188075855872), Fraction(82328172806027435, 144115188075855872), Fraction(84162263347021749, 144115188075855872), Fraction(42987948059916839, 72057594037927936), Fraction(87768630275116871, 144115188075855872), Fraction(11192503755432871, 18014398509481984), Fraction(91289664841326393, 144115188075855872), Fraction(93017109375688587, 144115188075855872), Fraction(47360971873718129, 72057594037927936), Fraction(12050469194178561, 18014398509481984), Fraction(49031064993613963, 72057594037927936), Fraction(99696669938471549, 144115188075855872), Fraction(101306976090856853, 144115188075855872), Fraction(51446328509359829, 72057594037927936), Fraction(104453327282180043, 144115188075855872), Fraction(52994303760416641, 72057594037927936), Fraction(107498124545963027, 144115188075855872), Fraction(108981511431254299, 144115188075855872), Fraction(55219203800992125, 72057594037927936), Fraction(55934229461334511, 72057594037927936), Fraction(113271317783145387, 144115188075855872), Fraction(114646643183066251, 144115188075855872), Fraction(57997050407394745, 72057594037927936), Fraction(29328340786159991, 36028797018963968), Fraction(29651027373131239, 36028797018963968), Fraction(14983253263735463, 18014398509481984), Fraction(60549403127976029, 72057594037927936), Fraction(122302150272323749, 144115188075855872), Fraction(61737882827895071, 72057594037927936), Fraction(124619367129440775, 144115188075855872), Fraction(125732676712007421, 144115188075855872), Fraction(126815423785434811, 144115188075855872), Fraction(31966836290165409, 36028797018963968), Fraction(128888185141595835, 144115188075855872), Fraction(129877695587268661, 144115188075855872), Fraction(130835635972152453, 144115188075855872), Fraction(32940443361156875, 36028797018963968), Fraction(66327941441791931, 72057594037927936), Fraction(133517746953144471, 144115188075855872), Fraction(67173578077748179, 72057594037927936), Fraction(135143908881817337, 144115188075855872), Fraction(33976952865321509, 36028797018963968), Fraction(68319339104081869, 72057594037927936), Fraction(137336331466937145, 144115188075855872), Fraction(138000601655511917, 144115188075855872), Fraction(69315663653223875, 72057594037927936), Fraction(4350886097069591, 4503599627370496), Fraction(139791539932549787, 144115188075855872), Fraction(17540093111206675, 18014398509481984), Fraction(140815841341652733, 144115188075855872), Fraction(70638354471955653, 72057594037927936), Fraction(17712904459057425, 18014398509481984), Fraction(71047658925750713, 72057594037927936), Fraction(142452860179100915, 144115188075855872), Fraction(71387887875616997, 72057594037927936), Fraction(71531993042320837, 72057594037927936), Fraction(143317421139521513, 144115188075855872), Fraction(35884004836206207, 36028797018963968), Fraction(143719727634491027, 144115188075855872), Fraction(71934250762147383, 72057594037927936), Fraction(143982305363273257, 144115188075855872), Fraction(18007639204177149, 18014398509481984), Fraction(36026231171844757, 36028797018963968)]
w=[2/((1-x**2)*(ld.m(x))**2) for x in r]
from decimal import Decimal,Context,setcontext,ROUND_HALF_DOWN
setcontext(Context(prec=30, rounding=ROUND_HALF_DOWN))
float=lambda r:(lambda x:Decimal(x[0])/Decimal(x[1]))(str(r).split('/'))
d=2**8
s='{120.000000000000000000000000000,9.26052826812554737559992621571,3.32335097044784255118406403126,2.00000000000000000000000000000,1.50457548825155601882809780906,1.26582350605728327736200759570,1.13300309631934634747833911121,1.05218372089129332722399903882,1.00000000000000000000000000000}'
l=list(map(Decimal,s[1:-1].split(',')))
gamma=Decimal('0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495')
for (c,c_1,real),l_ in zip([
    (Decimal('0.1')*i,Decimal('10')/i,10/i) for i in range(2,11)],l):
    #let t=e^{-x^c},\ln t=-x^c,x=(-\ln t)^{\frac1c}
    #t(0)=1,t(\infty)=0
    #dt=-e^{-x^c}\cdot x^{c-1}\cdot c
    #f(x)=e^{-x^c}\cdot x^{c-1}\cdot c=t\cdot (-\ln t)^{\frac{c-1}c}\cdot c
    xt=lambda t:(-((t).ln()))**(c_1)
##    f=lambda x:c*(x**(c-1))#*(math.exp(-(x**c)))
    fln=lambda t:((-((t).ln())).ln())#+(t.ln())
##    g=lambda x:Decimal('0.5')+5*((-9*x).exp())#*math.exp(-x)
    g=lambda t:Decimal('0.5')+5*((-9*xt(t)).exp())#*math.exp(-x)
##    h=lambda x:f(x)*((-(x**c)).exp())*\
##       ((f(x)).ln()-(x**c)-((g(x)).ln()-x))
##    #0,inf->0,1: x=t/(t-1), xt-x=t, (x-1)t=x, t=x/(x-1)
##    h1=lambda t:h(t/(1-t))/((1-t)**2) # 0,1
##    h2=lambda v:h1((v+1)/2)/2
    h1=lambda t:g(t).ln()#f(t).ln()-((g(t).ln())-xt(t))
    int_=lambda h1:sum(
          sum(float(w_)*h1(float(((x+1)/2+i)/d))/(2*d) for w_,x in zip(w,r))
          for i in range(d))
    print(c,Decimal(math.gamma(real+1))+(c.ln())-int_(h1)-gamma*(1-c_1)-1)
h1_=r'''0.2 0.476044892154085268532810576285
0.3 0.326045243249124607083205564208
0.4 0.194544295792217104304557844994
0.5 0.0805857800502670858642542661562
0.6 -0.0176136562407463311404194282348
0.7 -0.102042234943715186336125809639
0.8 -0.174612380676135222634071174702
0.9 -0.237048749195042896134197564780
1.0 -0.290856610889057599214369459143'''
x=[i.split(': ') for i in r'''c = 0.2. Output: 119.2233798550179
c = 0.3. Outout: 8.077346771987397
c = 0.4. Output: 2.07833944013377
c = 0.5. Output: 0.8034827042913205
c = 0.6. Output: 0.3961739639940003
c = 0.7. Output: 0.2585689391629249
c = 0.8. Output: 0.2287758419066705
c = 0.9. Output: 0.2480070283065089
c = 1.0. Output: 0.2908566108890576'''.split('\n')]
y=[i.split(' ') for i in r'''0.2 119.223379855017945799292478451
0.3 8.07734677198740392689540304712
0.4 2.07833944013376923348715605315
0.5 0.80348270429132046532502570246
0.6 0.39617396399400051391328713871
0.7 0.25856893916292307961043114604
0.8 0.22877584190665494099938058972
0.9 0.24800702830645810748432393277
1.0 0.29085661088905759921436945914'''.split('\n')]
#[(lambda l:(l and min(l)) or 0)([n for n,(d1,d2) in enumerate(zip(i[1],j[1])) if d1!=d2]) for i,j in zip(x,y)]
#[0, 14, 15, 17, 17, 16, 15, 14, 17]
\$\endgroup\$
3
\$\begingroup\$

Wolfram Language (Mathematica), 105 bytes ~7.5 seconds

Mathematica is good at automatically determining the best method needed to perform the integration. Setting the working precision high gets us to results that match the stated answer. The time stated is to perform all the example integrals.

NIntegrate[((f=#*x^(#-1))*Log[f/(E^x^#*(5/E^(10*x)+0.5/E^x))])/E^x^#,{x,0,-Log@0},WorkingPrecision->128]&

Try it online!

\$\endgroup\$
  • \$\begingroup\$ At least**!** ))) \$\endgroup\$ – Alexey Burdin Nov 5 at 15:27
  • 1
    \$\begingroup\$ At least someone posted wolfram mathematica solution. The language is perfect for numerical integration) or for all math) \$\endgroup\$ – Alexey Burdin Nov 5 at 15:51
  • 1
    \$\begingroup\$ 62 bytes. In Mathematica, precision is specified in decimal digits. MachinePrecision offers 15.9 decimal digits (53 bits) of precision by default on TIO, which is sufficient for the stated requirement. \$\endgroup\$ – attinat Nov 6 at 8:12
  • 1
    \$\begingroup\$ @Anush You only required accuracy to within +/-10^-10. 119.22337985501592 is sufficiently close. \$\endgroup\$ – attinat Nov 6 at 22:43
  • 1
    \$\begingroup\$ @attinat, I think you've got some good edits, but I'm not convinced that dropping the WorkingPrecision option works when c=0.8. As the problem states, 10^-10% implies a relative error of 10^-12 or less, but the relative difference between your answer and the stated correct answer is 7.5*10^-11 which is just a tad bit too much error. I figured that was a corner case left in by the questioner to force the use of extra precision. \$\endgroup\$ – Kelly Lowder Nov 7 at 22:48
3
\$\begingroup\$

Pari/GP, 88 87 bytes, < 1 sec

Great work from @AlexeyBurdin!

I also saw the "provocative" way in which \$f\$ was a derivative of an easier function, but only had the idea to use this for partial integration, without success.

For golfing reasons I recombined all the logarithms, except the one giving \$-1\$ because it makes no difference for the code length and subtracting one seemed slightly more efficient than multiplying the integrand by \$t\$. I also took the logarithms of \$1/t\$ instead of \$t\$ to be able to get rid of the otherwise needed parentheses for the minus sign. The 2 in the last argument to intnum reduces the integration step size. That is needed to achieve the required precision.

l=log;i(c)=intnum(t=0,1,l(c*l(1/t)^(1-1/c)/(.5+5*exp(-9*l(1/t)^c^-1)))+l(1/t)^c^-1,2)-1

Try it online!

\$\endgroup\$
  • 1
    \$\begingroup\$ This is pretty awesome. Thank you. \$\endgroup\$ – Anush Nov 6 at 16:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.