[GiNaC-list] Re: Airy functions
Alexei Sheplyakov
varg at theor.jinr.ru
Fri Dec 21 15:44:20 CET 2007
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.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 827 bytes
Desc: Digital signature
Url : http://www.cebix.net/pipermail/ginac-list/attachments/20071221/c338af37/attachment.pgp
More information about the GiNaC-list
mailing list