#4850 realpart/imagpart/rectform of erf and other special functions give useless or unsimplifiable results for complex arguments

[David Scherfgen]

(co-authored by Claude)

For functions that satisfy f(conjugate(z)) = conjugate(f(z)) — erf, erfc, erfi, erf_generalized, fresnel_s, fresnel_c, and likewise gamma, zeta, beta, the Airy functions, expintegral_ei and others — rectform, realpart and imagpart often fail as soon as the argument is complex but not literally of the form a+b*%i:

(%i1) rectform(erf((-1)^(1/4)*x));
(%o1) 'realpart(erf((-1)^(1/4)*x))+%i*'imagpart(erf((-1)^(1/4)*x))

(%i2) float(realpart(erf((-1)^(1/4))));
(%o2) 'realpart(erf(1.0*(-1)^(1/4)))

The second example is a plain numerical constant that cannot be brought to a number. Even when a result is produced, it can be unsatisfactory in two ways. An "imaginary part" may still contain %i:

(%i3) imagpart(erf(%i*t));
(%o3) -(%i*erf(%i*t))

(erfi(t) would be the natural answer). And a successful split can mix two different representations of the same value, which then never cancel against each other in later computations, because Maxima does not identify the various forms of e^(%i*%pi/4)-type constants:

(%i4) rectform(erf(sqrt(-%i)*x));
(%o4) (erf((%i/sqrt(2)+1/sqrt(2))*x)+erf(sqrt(-%i)*x))/2-(erf((%i/sqrt(2)+1/sqrt(2))*x)-erf(sqrt(-%i)*x))/2

(%i5) ratsimp(erf(%e^(-%i*%pi/4)*x) - erf(sqrt(-%i)*x));
(%o5) -erf(((sqrt(2)*%i-sqrt(2))*x)/2)-erf(sqrt(-%i)*x)

%e^(-%i*%pi/4) and sqrt(-%i) are the same number, so %o5 is exactly zero, but no simplifier gets there. The same failures occur for other functions:

(%i6) rectform(gamma((-1)^(1/4)*x));
(%o6) 'realpart(gamma((-1)^(1/4)*x))+%i*'imagpart(gamma((-1)^(1/4)*x))

and likewise for zeta, airy_ai, expintegral_ei, ...

This is also the reason integrate produces noun-form-polluted antiderivatives:

(%i7) integrate(sin(x^2), x);
(%o7) (sqrt(%pi)*((sqrt(2)*%i+sqrt(2))*'realpart(erf((-1)^(1/4)*x))+(sqrt(2)*%i-sqrt(2))*'imagpart(erf((-1)^(1/4)*x))+(sqrt(2)-sqrt(2)*%i)*erf(sqrt(-%i)*x)))/8

which does not differentiate back to sin(x^2).

Cause: The mirror-symmetry branch of risplit (rpart.lisp) builds the second kernel by calling conjugate(). That either fails outright (conjugate((-1)^(1/4)) does not simplify) or returns the conjugated coefficient in a different representation than the one already present.

I have a patch: Rebuild both kernels from risplit's own rectangular split (a . b) of the argument, so realpart(f(a+%i*b)) = (f(a+%i*b)+f(a-%i*b))/2 and imagpart(f(a+%i*b)) = (f(a+%i*b)-f(a-%i*b))/(2*%i). Both parts are real by the symmetry, and every representation of the same argument lands on the same pair of kernels. The change is in the mirror-symmetry branch of risplit itself, so it covers every function carrying the commutes-with-conjugate property at once. The erf family additionally gets small risplit-function handlers (the same hook bessel.lisp uses) so that purely imaginary arguments split cleanly — with them, (%o3) comes out as erfi(t). (The erf simplifier does contain an erf(%i*b) -> %i*erfi(b) rule behind the undocumented, default-false flag erf_%iargs, analogous to %iargs for sin and the Fresnel functions; but enabling that flag globally makes other parts of the test suite fail, so the handlers provide this only where realpart/imagpart/rectform need it.) With the patch, %i7 returns a clean combination of erf((1±%i)*x/sqrt(2)) whose derivative trigrat's back to sin(x^2), and the test suite passes.