| Bytes | Lang | Time | Link |
|---|---|---|---|
| 050 | APL Dyalog APL | 250724T115336Z | panadest |
| 098 | C gcc | 240724T225036Z | engineer |
| 026 | Charcoal | 240720T070451Z | Neil |
| 022 | Pyth | 240725T230147Z | CursorCo |
| 141 | Python | 240722T035101Z | Bubbler |
| 096 | K ngn/k | 240722T100340Z | Bubbler |
| 058 | JavaScript Node.js | 240722T084413Z | l4m2 |
APL (Dyalog APL), 50 bytes
{⊃⌽⊃(1∘↓⍤⊢,⊢/⍤⊢-2∘↑⍤⊢+.÷2×0 1+⊣)/(⍺+⍳⍺×⍵-1),⊂1⍺/1}
-9 bytes thanks to @Adám.
BQN, 50 bytes
1e9{⊢´(1⥊˜𝕗+1)(⊢«˜⊢´∘⊢-2÷˜2⊸↑+´∘÷+⟜0‿1)´𝕗+↕𝕗-˜𝕗×𝕩}
Explanation
Both programs implement Heun's method, basically this recurrence:
\$p_{i +1} = p_i - \frac{1}{2} \left( \frac{p_{i -N}}{i} + \frac {p_{i+1-N }}{i+1} \right)\$
The accuracy in the REPL link has been reduced to 1e4 to prevent the browser from crashing.
C (gcc), 112 110 109 98 bytes
-2 thanks to @l4m2
-1 thanks to @ceilingcat
-11 thanks to @att
V=2e4,i;double k(double u){double L[V*8],l=1;for(i=V;i<u*V;L[i]=l-=(i<2*V?:L[i-V])/i++);return l;}
Splits the domain \$[1,8]\$ across \$(8-1)(2*10^4)=140000\$ slices and computes each slice's y-value successively based off of the previous y-values.
Charcoal, 26 bytes
≧×⁸φＦ×Ｎφ⊞υ∨‹ιφ⁻⌊υ∕§υ±φιＩ⊟υ
Try it online! Link is to verbose version of code. Note that the version shown is only within 0.1% for x up to about 6; for more accuracy the ⁸ could be changed to χ or even φ but it would make the resulting program inordinately slow. Explanation:
≧×⁸φ
Divide the interval into strips of width 1/8000 (or 1/10000 or 1/1000000).
Ｆ×Ｎφ
Calculate the number of strips needed to reach an input of x.
⊞υ∨‹ιφ⁻⌊υ∕§υ±φι
For strips of 0<x<1 use a height of 1, otherwise estimate the slope of the graph and use that to estimate the height of the next strip.
Ｉ⊟υ
Output the height of the last strip.
A port of @Bubbler's Python (backport of @att's K) answer weighs in at a massive 72 63 bytes:
ＮθＦ⁹⁹⊞υ¬ιＦθ«≔⟦⁰⟧ηＦυ⊞η∕⁺κ§η±¹⁺²ιＵＭη∕κ∨¬λλ≔Ｅη⎇λκ⁻§υ⁰Σηυ»Ｉ↨⁻¹﹪θ¹⮌υ
Try it online! Link is to verbose version of code. Explanation:
Ｎθ
Input x.
Ｆ⁹⁹⊞υ¬ι
Start with the polynomial for 0<=x<=1.
Ｆθ«
Repeat ⌊x⌋ times.
≔⟦⁰⟧η
Start a new polynomial with a dummy constant term.
Ｆυ⊞η∕⁺κ§η±¹⁺²ι
For each term of the previous polynomial, subtract the latest term of the new polynomial, divide it by the endpoint of the interval, and append it to the new polynomial.
ＵＭη∕κ∨¬λλ
Divide each term (except the dummy term) of the new polynomial by its index.
≔Ｅη⎇λκ⁻§υ⁰Σηυ
Calculate the constant term to replace the dummy term of the new polynomial and replace the old polynomial.
»Ｉ↨⁻¹﹪θ¹⮌υ
Evaluate the polynomial at ⌊x⌋+1-x.
Pyth, 22 bytes
L|g1b-y-b^T_6cytb*b^T6
Try it online! (link is to a faster version of the code which has way too much error, but will actually run online)
Defines a recursive function y which computes the Dickman function of its input.
Explanation
This answer is in essence, the same as all the other answers which slice the domain up into intervals and compute the function discretely, except that here it is done recursively. This method can be thought of as a piecewise function \$f(x) = 1\$ for \$x \leq 1\$ and \$f(x) = f(x-d)-\frac{d\cdot f(x-1)}{x}\$ for \$1 < x\$ where \$d\$ is a constant which controls the accuracy. The smaller the \$d\$ the less error, but the computation becomes very expensive very quickly. Using \$d=10^{-6}\$ should be enough to meet the requirements, but the link above is to \$d=10^{-3}\$ as that's as much as will finish online. If the function was memoized, it would actually be about the same as computing with slices, but without it it's about as bad as recursive fibonacci.
L # define y(b)
g1b # if 1 >= b
| # return 1
# otherwise return
y-b^T_6 # y(b-10^(-6))
- # minus
ytb # y(b-1)
c # divided by
*b^T6 # b * 10^6
Python, 141 bytes
g=lambda x,p=[1]+[0]*99,I=2,z=0,i=0:x<1and sum(a*(1-x)**~-(i:=i+1)for a in p)or g(x-1,[p[0]-sum(A:=[(z:=(z+w)/I)/(i:=i+1)for w in p])]+A,I+1)
Ungolfed
def N(p,I):
A=p[:];l=len(A)
for i in range(1,l):A[i]+=A[i-1]/I
for i in range(l):A[i]/=I*-~i
return[p[0]-sum(A)]+A
def g(x):
f=[1]+[0]*99
for i in range(int(x)):f=N(f,i+2)
return sum(a*(1-x%1)**i for i,a in enumerate(f))
A backport of att's short K solution. Simply using the interval [0, 1] degrades the error massively, so I wondered why att's code worked and what att meant by interval [1, 0] written backwards.
First, rewrite the math below using \$-1 \le x \le 0\$ instead:
- \$f_0 (x) = 0\$
- \$(x + 2) f_1' (x) + f_0 (x) = 0\$, \$f_1(-1) = f_0(0)\$
- \$(x + i + 1) f_i' (x) + f_{i-1} (x) = 0\$, \$f_i(-1) = f_{i-1}(0)\$
Then substitute \$-x\$ for \$x\$:
- \$-(-x + i + 1) f_i' (-x) + f_{i-1} (-x) = 0\$ (the leading
-comes from the chain rule for \$f_i'\$)
Then redefine \$f_i(x)\$ to \$f_i(-x)\$:
- \$-(i + 1 - x) f_i'(x) + f_{i-1}(x) = 0\$, \$f_i(1) = f_{i-1}(0)\$
Using the same logic, we get
$$f_i(x) = C + \int \frac{f_{i-1}(x)}{i + 1 - x} dx$$
Now the Maclaurin series expansion of \$\frac{1}{k-x}\$ is [k^-1, k^-2, k^-3, ...] (notice the lack of alternating signs). Also, for the boundary condition, \$f(1)\$ is the sum of coefficients and \$f(0)\$ is the first term, so they can be simplified.
Python, 254 bytes
def N(p,I):
A=p[:];l=len(A)
for i in range(1,l):A[i]-=A[i-1]/I
for i in range(l):A[i]/=I*~i
return[E(p,.5)+E(A,-.5)/2]+A
E=lambda f,x:sum(a*x**i for i,a in enumerate(f))
def g(x):
f=[1]+[0]*20
for i in range(int(x)):f=N(f,i+1.5)
return E(f,x%1-.5)
g(x) computes the Dickman function \$ \rho(x) \$ with relative error under 5e-7 for up to x <= 8.
Instead of a single continuous function \$ \rho(x) \$, let's define a sequence of functions \$f_0(x), f_1(x), f_2(x), \cdots\$ which are defined on the interval \$-\frac{1}{2} \le x \le \frac{1}{2}\$:
- \$f_0(x) = 1\$.
- \$(x + \frac{3}{2}) f_1'(x) + f_0(x) = 0\$.
- \$(x + i + \frac{1}{2}) f_i'(x) + f_{i-1}(x) = 0\$ in general.
Now we will approximate these functions using truncated Maclaurin series. Let's write f(x) = [a0, a1, a2, ...] as a convenient notation for \$f(x) = a_0 x^0 + a_1 x^1 + a_2 x^2 + \cdots\$.
Solving \$(x + i + \frac{1}{2}) f_i'(x) + f_{i-1}(x) = 0\$ for \$f_i\$ gives
$$f_i(x) = C + \int \frac{f_{i-1}(x)}{x + i + \frac{1}{2}} dx$$
so we will multiply the Maclaurin series of \$f_{i-1}(x)\$ with that of \$\frac{1}{x + i + \frac{1}{2}}\$, integrate it, and find the constant term using the boundary condition \$f_{i-1}(\frac{1}{2}) = f_{i}(-\frac{1}{2})\$.
The Maclaurin series expansion of \$\frac{1}{x + k}\$ is [k^-1, -(k^-2), k^-3, -(k^-4), ...]. Factor out k^-1 (we will multiply it back later), and multiplying with [a0, a1, a2, a3, ...] gives
[
a0,
-k^-1 * a0 + a1,
k^-2 * a0 - k^-1 * a1 + a2,
-k^-3 * a0 + k^-2 * a1 - k^-1 * a2 + a3,
...
]
which can be converted to a simple scan "multiply the previous cumulative result with -1/k and add current term", i.e. if the resulting series is [b0, b1, b2, b3, ...] then it is equal to
[a0, -1/k * b0 + a1, -1/k * b1 + a2, -1/k * b2 + a3, ...]
Integration is also easy: since integration of x^n is x^(n+1)/(n+1), we can divide each term with its 1-based index and prepend a constant term.
The constant term is simply \$f_{i-1}(\frac{1}{2}) - f_{i}(-\frac{1}{2})\$ where \$f_i\$ is evaluated without the constant term.
E(f,x) evaluates the polynomial f at x. N(p,i) takes the series for \$f_i\$ and the value of \$i+1.5\$ to compute the series for \$f_{i+1}\$. g(x) finds the appropriate function for the interval that x belongs to, and then evaluates the function at x offset by an appropriate amount.
I had to read a few papers cited on the Wikipedia page to get the idea, but computing with the idea turned out to be not that complex, yet surprisingly efficient and accurate. I intentionally didn't golf the program hard to keep the intents relatively clear, in the hopes that others can port this without too much problem.
K (ngn/k), 96 bytes
{E[-0.5+1!x;(1,&20){(E[0.5;x]+E[-0.5;w]%2),w:({z-y%x}[y]\x)%y*|!-#x}/1.5+!_x]}
E:{{z+y*x}[x]/|y}
Port of my own Python answer. It could have been a bit shorter if base conversion allowed non-integer bases, and those 0.5s hurt too.
JavaScript (Node.js), 58 bytes
f=x=>f[x*1e6|0]
for(i=0;i<9e6;)f[1+i]=f[i]-f[i-1e6]/i++||1
Delta matters lol