Re: [Numpy-discussion] Re: indexing problem

[Tim H. <tim...@co...>]

Tim Hochberg wrote:

> 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.


To add a little more detail to this, the formula that numpy uses to 
compute a**b is:

    exp(b*log(a))

where:

    log(x) = log(|x|) + arctan2(x.imag, x,real)
    exp(x) = exp(x.real) * (cos(x.imag) + 1j*sin(x.imag))

With these definitions in hand, it should be apparent what's happening 
when *a* = |a|j and *b* = 2. First, let's compute 2*log(a) for *a* = 
24674.011002723393j:

   2*log(a) = (20.227011565110185+3.1415926535897931j)
  
Now it's clear what's happening: ideally the sine of the imaginary part 
of the above number should be zero. However:

  sin(3.1415926535897931) = 1.2246063538223773e-016

And this in turn leads to the error we see here.

-tim


>
> Regards,
>
> -tim
>
>
>> Thanks,
>>
>> Ryan
>>
>> On 2/13/06, Tim Hochberg <tim...@co...> 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 <rya...@gm...> 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
>>> Num...@li...
>>> 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
>> Num...@li...
>> 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=lnk&kid=103432&bid=230486&dat=121642
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>
>