On Thu, May 14, 2009 at 10:12 PM, Carl Witty <[email protected]> wrote:
> On Thu, May 14, 2009 at 4:48 PM, Lisandro Dalcin <[email protected]> wrote:
>> However, interestingly enough, my C++ <complex> header (from GCC
>> 4.3.2) seems to implements the naive way (although warned in the
>> comment)... Am I missing something here??
>
> GSL (the GNU Scientific Library) also uses the naive algorithm.
> Perhaps everybody else has decided that speed is more important than
> accuracy, or perhaps they haven't thought of the issues at all?
>
I bet that the second one is the case ;-)
>
> How much difference in accuracy are we talking about?
>
Basically, we have two issues: overflow and cancelation.
* Overflow will occur if the real or imaginary part of the denominator
has exponent > 1/2 max_exponent; for example real|imag parts are above
1e20 for float or 1e155 for double
* Cancellation can occur on subtraction when computing the imaginary component.
See yourself the effects of all this in single precision:
import numpy
def cdiv(z1, z2):
z1 = numpy.array(z1, dtype='F')
z2 = numpy.array(z2, dtype='F')
a, b = z1.real, z1.imag
c, d = z2.real, z2.imag
nrm = c*c + d*d
z = numpy.array(0, dtype='F')
z.real = (a*c + b*d) / nrm
z.imag = (b*c - a*d) / nrm
return z
def test(z1, z2):
print ('z1 = %s' % z1)
print ('z2 = %s' % z2)
print ('z1/z2 (naive) = %s' % cdiv(z1,z2))
print ('z1/z2 (smart) = %s' % (z1/z2))
print("")
if __name__ == "__main__":
z1 = 1 + 2j
z2 = 3 + 4j
test(z1, z2)
z1 = 0 + 1e20j
z2 = 1e20 + 0j
test(z1, z2)
z1 = 1 + (1+1e-7)*1j
z2 = 1 + 1j
test(z1, z2)
I get this output (Win32, Py2.6, NumPy 1.3.0rc1)
z1 = (1+2j)
z2 = (3+4j)
z1/z2 (naive) = (0.439999997616+0.0799999982119j)
z1/z2 (smart) = (0.44+0.08j)
z1 = 1e+20j
z2 = (1e+20+0j)
z1/z2 (naive) = nan*j
z1/z2 (smart) = 1j
z1 = (1+1.0000001j)
z2 = (1+1j)
z1/z2 (naive) = (1+5.96046447754e-08j)
z1/z2 (smart) = (1.00000005+5.00000000292e-08j)
In the second bit, overflow generated a NaN. In the third bit,
cancellation produced an imaginary component 20% larger.
> (BTW, the Sage arbitrary-precision complex and complex interval
> arithmetic also use the naive algorithm; if it is significantly bad,
> then Sage should be fixed too.)
>
> On the other hand, I think that gcc's "complex double" does not use
> the naive algorithm. (I haven't traced through the assembly language
> to see exactly what it does, but it does have conditional branches, if
> I'm reading the x86 assembly correctly.)
>
> Carl
> _______________________________________________
> Cython-dev mailing list
> [email protected]
> http://codespeak.net/mailman/listinfo/cython-dev
>
--
Lisandro Dalcín
---------------
Centro Internacional de Métodos Computacionales en Ingeniería (CIMEC)
Instituto de Desarrollo Tecnológico para la Industria Química (INTEC)
Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET)
PTLC - Güemes 3450, (3000) Santa Fe, Argentina
Tel/Fax: +54-(0)342-451.1594
_______________________________________________
Cython-dev mailing list
[email protected]
http://codespeak.net/mailman/listinfo/cython-dev
