[GiNaC-list] Re: Airy functions

[Alexei Sheplyakov]

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