[Numpy-discussion] Re: indexing problem

[Tim Hochberg]

Ryan Krauss wrote:

>At the risk of sounding silly, can you explain to me in simple terms
>why s**2 is less accurate than s*s.  I can sort of intuitively
>appreciate that that would be true, but but might like just a little
>more detail.
>  
>
I don't know that it has to be *less* accurate, although it's unlikely 
to be more accurate since s*s should be nearly as accurate as you get 
with floating point. Multiplying two complex numbers in numpy is done in 
the most straightforward way imaginable:

    result.real    = z1.real*z2.real - z1.imag*z2.imag
    result.image = z1.real*z2.imag + z1.imag*z2.real

The individual results lose very little precision and the overall result 
will be nearly exact to within the limits of floating point.

On the other hand, s**2 is being calculated by a completely different 
route. Something that will look like:

    result = pow(s, 2.0)

Pow is some general function that computes the value of s to any power. 
As such it's a lot more complicated than the above simple expression. I 
don't think that there's any reason in principle that pow(s,2) couldn't 
be as accurate as s*s, but there is a tradeoff between accuracy, speed 
and simplicity of implementation.

That being said, it may be worthwhile having a look at complex pow and 
see if there's anything suspicious that might make the error larger than 
it needs to be.

If all of that sounds a little bit like "I really know", there's some of 
that in there too.

Regards,

-tim


>Thanks,
>
>Ryan
>
>On 2/13/06, Tim Hochberg <tim.hochberg at cox.net> wrote:
>  
>
>>>>Ryan Krauss wrote:
>>>>
>>>>
>>>>
>>>>        
>>>>
>>>>>This may only be a problem for ridiculously large numbers.  I actually
>>>>>meant to be dealing with these values:
>>>>>
>>>>>In [75]: d
>>>>>Out[75]:
>>>>>array([   246.74011003,    986.96044011,   2220.66099025,   3947.84176044,
>>>>>       6168.50275068,   8882.64396098,  12090.26539133,  15791.36704174,
>>>>>      19985.94891221,  24674.01100272])
>>>>>
>>>>>In [76]: s=d[-1]*1.0j
>>>>>
>>>>>In [77]: s
>>>>>Out[77]: 24674.011002723393j
>>>>>
>>>>>In [78]: type(s)
>>>>>Out[78]: <type 'complex128scalar'>
>>>>>
>>>>>In [79]: s**2
>>>>>Out[79]: (-608806818.96251547+7.4554869875188623e-08j)
>>>>>
>>>>>So perhaps the previous difference of 26 orders of magnitude really
>>>>>did mean that the imaginary part was negligibly small, that just got
>>>>>obscured by the fact that the real part was order 1e+135.
>>>>>
>>>>>On 2/13/06, Ryan Krauss <ryanlists at gmail.com> wrote:
>>>>>
>>>>>
>>>>>          
>>>>>
>>I got myself all tied up in a knot over this because I couldn't figure
>>out how multiplying two purely complex numbers was going to result in
>>something with a complex portion. Since I couldn't find the complex
>>routines my imagination went wild: perhaps, I thought, numpy uses the
>>complex multiplication routine that uses 3 multiplies instead of the
>>more straightforward one that uses 4 multiples, etc, etc. None of these
>>panned out, and of course they all evaporated when I got pointed to the
>>code that implements this which is pure vanilla.  All the time I was
>>overlooking the obvious:
>>
>>Ryan is using s**2, not s*s.
>>
>>So the obvious answer, is that he's just seeing normal error in the
>>function that is implementing pow.
>>
>>If this is inacuracy is problem, I'd just replace s**2 with s*s. It will
>>probably be both faster and more accurate anyway
>>
>>Foolishly,
>>
>>-tim
>>
>>
>>
>>
>>
>>
>>
>>-------------------------------------------------------
>>This SF.net email is sponsored by: Splunk Inc. Do you grep through log files
>>for problems?  Stop!  Download the new AJAX search engine that makes
>>searching your log files as easy as surfing the  web.  DOWNLOAD SPLUNK!
>>http://sel.as-us.falkag.net/sel?cmd=lnk&kid=103432&bid=230486&dat=121642
>>_______________________________________________
>>Numpy-discussion mailing list
>>Numpy-discussion at lists.sourceforge.net
>>https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>>
>>    
>>
>
>
>-------------------------------------------------------
>This SF.net email is sponsored by: Splunk Inc. Do you grep through log files
>for problems?  Stop!  Download the new AJAX search engine that makes
>searching your log files as easy as surfing the  web.  DOWNLOAD SPLUNK!
>http://sel.as-us.falkag.net/sel?cmd=k&kid3432&bid#0486&dat1642
>_______________________________________________
>Numpy-discussion mailing list
>Numpy-discussion at lists.sourceforge.net
>https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>
>
>  
>