On May 14, 2009, at 6:12 PM, Carl Witty 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?
>
> How much difference in accuracy are we talking about?
Not much, but a little:
C = ComplexField(100)
total = 100000
naive_avg = 0
naive_worst = 0
naive_better = 0
other_avg = 0
other_worst = 0
other_better = 0
for i in range(total):
a = C.random_element()
b = C.random_element()
naive_error = abs(C(CDF(a)/CDF(b)) - a/b)/abs(a/b)
other_error = abs(C(complex(a)/complex(b)) - a/b)/abs(a/b)
#print naive_error / other_error
naive_better += naive_error < other_error
naive_avg += naive_error / total
naive_worst = max(naive_worst, naive_error)
other_better += naive_error > other_error
other_avg += other_error / total
other_worst = max(other_worst, other_error)
print "naive\n better", naive_better, "\n avg", naive_avg,
"\n worst", naive_worst
print "other\n better", other_better, "\n avg", other_avg,
"\n worst", other_worst
naive
better 26187
avg 1.4940916064705992895601085724e-16
worst 5.7659574333851360909896621025e-16
other
better 63116
avg 9.5414951065097745299683547276e-17
worst 3.9584519821557590591785765975e-16
> (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.)
I think the biggest issue is overflow/underflow, which we
(essentially) don't have to worry about in the arbitrary-precision case.
> 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.)
Interesting. In terms of branching on abs, we can do a bitwise and
(but I'm not sure how to write good generic code to do this in both
the floating point and integer case).
- Robert
_______________________________________________
Cython-dev mailing list
[email protected]
http://codespeak.net/mailman/listinfo/cython-dev
