Codegolf. Насколько Компактно Ваш Язык Может Выполнять Точное Числовое Интегрирование?

  • Автор темы Dmitry52
  • Обновлено
  • 21, Oct 2024
  • #1

ВНИМАНИЕ: для этой задачи могут потребоваться 128-битные числа с плавающей запятой.1

Задача состоит в том, чтобы выполнить численное интегрирование. Рассмотрим следующие три функции.

\$

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

\$

\$

g_1(x) = 0,5e^{-x}

\$

\$

g_2(x) = 5 е^{-10 х}

\$

У нас будет \$c \geq 0.2\$. Ваш код должен быть корректным для любого значения от 0,2 до 1,0.

 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
 

Задача – выполнить интеграл:

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

Ваш результат должен быть правильным с точностью до плюс-минус \$10^{-10}\$ процентов от правильного ответа. Вход Единственное вещественное число \$c\$.

Вывод показан в 15 значимых местах

  • Тайминг

Ваш код должен завершиться

  • ТИО

прежде чем истечет время. Как обычно это конкурс за

  • язык. Другими словами, программистов Java не должно отпугивать односимвольное решение на языке, в существование которого они не могут поверить. Вы можете использовать любой язык/библиотеку по вашему выбору.

Вопросы и примечания

Действительно ли этот интеграл существует/сходится?

Да. Это расхождение Кульбака-Лейблера между двумя распределениями вероятностей.

Покажите мне любой код, который может вычислить эти интегралы.

Посмотреть это

Вольфрам Альфа

пример.

Dmitry52


Рег
09 Jan, 2011

Тем
86

Постов
196

Баллов
636
  • 26, Oct 2024
  • #2

python3, 3982 2183 1866 1821 1677, ~14 с

 
 
 
 
 
 (-lnt) 

Попробуйте онлайн! -- принимает ввод как допустимое десятичное число, например ln(1/t)
Математика:
0. Мы занимаемся численным интегрированием Google.
1. Давайте посмотрим на провокационные \$f(x)=c\cdot x^{c-1}\cdot e^{-x^c}=
-\left(e^{-x^c}\right)'\$
2. Итак, имеем \$-\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\$,
мы меняем переменную на \$t=e^{-x^c}\$, \$x^c=-\ln(t)\$,
\$x=\left(-\ln t\right)^\frac1c\$,
\$t(0)=1\$, \$t(+\infty)=0\$ и интеграл равен
\$\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. Далее вольфрамируем вот эти вещи:\$\int\limits_0^1\ln(-\ln t)dt\$ \$\int\limits_0^1 (-\ln t)^a dt\$ ,
\$\int\limits_0^1 \ln t\,dt\$ и сделайте все остальное

a=ln(1/t)^{1/c} f(c)=∫_0^1(ln(a^c/a/e(.5+5e^{-9a})c)+a)dt ||answer||

численно(последнее, возможно, можно сделать с помощью трапеций, прямоугольников или Симпсона — идея для коротких материалов).

Сырая вещь:

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

Язык Wolfram (Математика)

, 105 байт ~7,5 секунд

Mathematica умеет автоматически определять лучший метод, необходимый для выполнения интегрирования. Установка высокой точности работы дает нам результаты, соответствующие заявленному ответу. Указанное время предназначено для выполнения всех примеров интегралов.Попробуйте онлайн!

 

Sawka2k13


Рег
09 Apr, 2014

Тем
71

Постов
202

Баллов
587
  • 26, Oct 2024
  • #3

Париж/GP

, 88 87 байт, < 1 сек. quad in the last argument to gamma Отличная работа от @AlexeyBurdin!

euler

Я также видел «провокационный» способ, которым \$f\$ была производной более простой функции, но у меня возникла идея использовать это только для частичного интегрирования, но безуспешно.

По соображениям игры в гольф я перекомбинировал все логарифмы, кроме того, который дает \$-1\$, потому что это не имеет значения для длины кода, а вычитание одного казалось немного более эффективным, чем умножение подынтегральной функции на \$t\$. Я также взял логарифмы \$1/t\$ вместо \$t\$, чтобы избавиться от необходимых в противном случае круглых скобок для знака минус.

уменьшает размер шага интегрирования. Это необходимо для достижения необходимой точности.Попробуйте онлайн!

 

Xaoc85


Рег
09 Sep, 2006

Тем
66

Постов
187

Баллов
517
  • 26, Oct 2024
  • #4

Питон 3 , 150 133 122 байт (спасибо Дельте и Алексею Бурдину за дальнейшую игру в гольф) log , exp , 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 , intnum SymPy's 2 модуль мпмат

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

это отличный способ справиться с этой проблемой. Все математические функции (

и

для интеграции) выполняются с высокой точностью.Попробуйте онлайн!

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]

 

Lenali


Рег
26 Jun, 2011

Тем
66

Постов
184

Баллов
524
  • 26, Oct 2024
  • #5
Десмос , 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 мс.

 

Sur54


Рег
20 May, 2006

Тем
67

Постов
195

Баллов
550
Тем
403,760
Комментарии
400,028
Опыт
2,418,908

Интересно

Lumtu.com © 2024