[Numpy-discussion] Re: indexing problem
Tim Hochberg
tim.hochberg at cox.net
Mon Feb 13 15:34:01 EST 2006
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
>
>
>
>
More information about the NumPy-Discussion
mailing list