$(1+x^M)^{1/M}$ need to be able to calculate any order derivatives vs. $x$ and $M$ for $x\ge 0$ and $M\gt 2$

cannot delete my own question, so I try to overwrite it instead...

• Do you have an idea what the maximum orders of derivatives needed will be? Jan 22 '16 at 3:13
• This is to be used for a compiler. Thus the max. order of derivatives depends on the input source code written by the user.
– user7098
Jan 22 '16 at 17:23

One approach is to use the Faà di Bruno formula: $$\big((1+x^M)^{1/M}\big)^{(n)} = \sum \frac{n!}{m_1!1!^{m_1}\cdots m_n! n!^{m_n}} (1/M)^{\underline{m_1+\cdots+m_n}}(1+x^M)^{1/M-m_1-\cdots-m_n} \prod_{j=1}^{n} \big(M^{\underline{\phantom{,}j}} x^{M-j}\big)^{m_j}, \qquad M\notin\mathbb{Z}$$ Here $a^{\underline{b}}$ is the falling factorial, and the sum is over all $n$-tuples of integers $m_1,\ldots,m_n\geq0$ such that $\sum_k km_k = n$. This formula is likely to be unstable for very large values of $n$.

A decent alternative is to evaluate the contour integral $$\frac{n!}{2\pi \mathrm{i}}\int_\gamma (1+z^M)^{1/M}\frac{dz}{(z-x)^{n+1}},$$ where the contour $\gamma$ is a circle of some radius $R$ surrounding $x$, $z = x+Re^{\mathrm{i}t}$, with $R$ chosen such that all the branch points of $(1+x^M)^{1/M}$ in the complex plane are outside the circle. The choice of $R$ affects the stability of this formula. This is a straightforward Fourier integral, so one can directly use an existing library to compute it, and one can also evaluate many derivatives at once efficiently using FFT.

• I don't see how the first method avoids calculating the derivative vs. M at x==0. The problem is, that I always end up with a log(0.0)...
– user7098
Jan 22 '16 at 18:53
• @GrapschKnutsch In the first formula, when $M$ is an integer, $x^{M-j}$ should be replaced with $x^{M-j}[M\geq j]$ because it's really $(x^M)^{(j)}$ (see the full formula). At $x=0$, one just has $0^0=1, 0^k=0$. Not sure where you get $\log$. When $M$ is not an integer, the function is not infinitely differentiable at $x=0$, anyway, so the formula has to break down. The second (Cauchy integral) formula is usually slower, but more robust, and doesn't need such tinkering. Jan 22 '16 at 19:07
• "the function is not infinitely differentiable at x=0": Does this mean, that this function should not be part of a transistor model or it should not be calculated a x=0 when being part of such a model?
– user7098
Jan 22 '16 at 19:28
• @GrapschKnutsch You lost me there, I have no idea what you mean. What model? Can you perhaps create a separate question about what you actually want to do? Jan 22 '16 at 19:37
• This function is part of some famous transistor model (a mathematical description of electrical currents as a function of voltages) used in simulation of electrical circuits. Can I assume, that this function is not suitable for this purpose, as it is not physical.
– user7098
Jan 22 '16 at 19:48