Dear Alexei, Thank you for your interest and careful examination of the code. I made an effort to fit the functions in the GiNaC framework and got series expansions etc. more or less for free. I'm looking forward to seeing the result when it meets the GiNaC standards. To answer your questions a fundamental property of the algorithm is needed: Depending on the argument a so called "steepest descent" contour is determined along which the integrand is monotonically decreasing. This ensures that a) roundoff does not occur, b) the range over which the integrand has to be evaluated to obtain a certain accuracy is minimized and c) the dominant factor (possibly very larger or small) can be placed outside the integral. Furthermore all derivatives at +/- infinity are zero. For functions with this property the trapezoidal rule has exponential convergence, that is, halving the stepsize more than doubles the number of significant digits of the result. For an integrator with non-constant step size this property would be lost.
AiryContour* pcntr=createAiryContour(zn);
Pointers are evil. What is the reason for explicit
allocation/deallocation?
Again depending on the argument (zn) but more involved, 3 different types of contours are needed of which one is the truly steepest descent contour and the other 2 are approximations thereof. AiryContour is the base class of them. I hope this answers your questions Best regards, Jan