Hello! On Sun, Dec 16, 2007 at 03:15:20PM +0100, Jan Bos wrote:
Please find attached an arbitrary precision implementation of AiryAi, AiryBi and their derivatives.
First of all, thanks a lot for your contribution! I've got some questions and remarks. First, I doubt the numerical integration is a good method for arbitrary-precision calculation. Although the paper claims the method can be used to compute Airy functions "to any desired accuracy... by modifying the step size", the actual analysis of *roundoff* errors is absent (I'll re-read the paper once more, may be I've missed it). I'd expect smaller step sizes to increase roundoff errors in uncontrollable manner. Something like (asymptotic) series expansion would be a better method. That said, it's just my gut feeling, I'm not really an expert. [skipped]
static ex AiryAi_evalf(const ex& z)
Few minor coding style related comments: - please avoid using caps in function names - FunctionAndVariableNamesLikeThisOneReallyHurtEyesAndAreNotReallyWellcome
if (is_a<numeric>(z)) { numeric zn = ex_to<numeric>(z); // Ai(z) -> conjugate(Ai(conjugate(z))) for Im(z) < 0 if (imag(zn).is_negative()) { numeric ret = ex_to<numeric>( AiryAi(conjugate(zn)).evalf()); return conjugate(ret); }
// for ph(z) in <2Pi/3,Pi] use // [1] 10.4.7 Ai(z) = -( zp*Ai(zp*z) + zm*Ai(zm*z) ), // zm = exp(-2/3*Pi*I), zp = exp(2/3*Pi*I) so // zm*z and conjugate(zp*z) are in [0,2Pi/3] if (!AiryInPrincipleSector(zn)) { numeric zm=ex_to<numeric>(AiryZm.evalf()); numeric zp=ex_to<numeric>(AiryZp.evalf()); return (-(zm*AiryAi(zm*zn)+zp*AiryAi(zp*zn))).evalf(); }
Using GiNaC::numeric for numeric calculations incurs significant overhead (yes I know the irony). In fact, GiNaC::numeric is not even intended for numeric intensive code. I suggest to use proper CLN types (i.e. cl_N for complex floating point numbers) instead (that is, if the numerical integration method is really suitable for arbitrary-precision calculation).
// ph(z) in [0,2Pi/3] perform contour integration AiryContour* pcntr=createAiryContour(zn);
Pointers are evil. What is the reason for explicit allocation/deallocation?
template <typename FUNC> numeric ModifiedTrapezoidalRule (const FUNC& f, numeric h, numeric d)
GiNaC already has numerical integration routine. Why re-invent the wheel? @ developers: BTW, the `integral' class sounds like re-inventing the wheel too. There are a *lot* of good numerical libraries out of there, and seamless integration with existing C/C++ code is one of the main reasons why GiNaC exists.
numeric answ=f(d);
// small number used to prevent div by 0 numeric TINY=ex_to<numeric>(pow(ex(10),ex(-(Digits+10))));
Please don't paper over the errors. Either handle singular points properly, or throw an exception if you can't (or don't want) to do it.
// approximation for smallest number so that 1+tol > 1 numeric tol=ex_to<numeric>(3*pow(ex(10),ex(-(Digits-2))));
Not really.
numeric j0=1;
You'd better use appropriate integer type for counters. Incrementing it is a single CPU instruction. Incrementing GiNaC::numeric is *much* more expensive. Best regards, Alexei -- All science is either physics or stamp collecting.