Re: [pure-lang-users] math.pure

[Eddie R. <er...@bm...>]

On Sat 23/08/08  8:27 AM , Albert Graef Dr....@t-... sent:
> Albert Graef wrote:
> Ok, I did some minor cosmetic surgery to the complex trigs and hyps, 
> actually most of my definitions (originally pilfered from Bronstein) 
> seemed to be pretty much in line with that document already. I also 
> fixed up my completely nonsensical complex sqrt definition, wonder how 
> that slipped in. ;-)

I pilfered that definition from a paper I found on Internet. It is actually the same definition as the one you
gave but they used the half angle identities cos (t/2) = +-sqrt (1+cos t) / 2 and sin (t/2) = +-sqrt (1-cos t) / 2
to derive it. As soon as I can find that link I again I will send it to you. I suspect I left out a 
sign function for the imaginary part.

However, why didn't you use Kahan's suggestion for sqrt?
sqrt z = exp (0.5 * ln z);
or even 
sqrt z = z ^ 0.5;

It's prettier to look at ;-)

> What still remains to be done is to add some type guards on the mixed 
> real/complex rules, so that we can later (when Eddie's GSL module 
> materializes) overload the common operators for scalar+vector/matrix 
> kind of stuff.

I'm going to act on this real soon.

> One thing that keeps me wondering is Kahan's suggestion for the acosh:
> 
> acosh z@(x+:y)	|
> acosh z@(r<:t)	= 2*ln (sqrt ((z+1)/2)+sqrt ((z-1)/2));
> 
> This is the definition I currently use (from Mathworld):
> 
> acosh z@(x+:y)	|
> acosh z@(r<:t)	= ln (z+sqrt (z-1)*sqrt (z+1));


> This has the same branch cut (-inf,1), and seems to yield the same 
> results AFAICT (up to rounding of course). Kahan's definition looks more
> complicated and requires more operations, but given that Kahan wrote it,
> I guess that there must be some reason he wrote it that way. ;-)
> 
> So can anyone think of a reason why one might prefer Kahan's definition?
> Numerical stability maybe?

Well, the Mathworld definition for actanh z didn't even give correct results comapred to the Common Lisp version.
 I'm just guessing here but as we well know the actual math definition is not always the best to use for floating 
point. I don't have the means to test it right now but I bet that sqrt ((z+1)/2)+sqrt ((z-1)/2 gives better results 
than z+sqrt (z-1)*sqrt (z+1) when either the real part or the imaginary part is either really small or really big.
Again, I'm just guessing.