Karim Belabas on Wed, 13 Apr 2005 20:16:32 +0200

 3M multiplication for t_COMPLEX

```Hi pari-dev,

I have been sent the following example:

\p36
x = 3.32307372544194849414053517840951697E-19 - 1.12813451938986614095048572644641324E-38*I;
y = 3.05042759366319004212246229138501820E36 + 103557517877613079.124097118142914269*I;

xr = real(x); xi = imag(x);
yr = real(y); yi = imag(y);

? imag(x * y)
%1 = -3.388131790 E-21

? xr*yi + xi*yr
%2 = 1.563040145 E-37

The problem comes from the 3M formula :

Im(y * x) := (xr + xi) * (yr + yi) - (xr * yr + xi * yi)
Re(y * x) := (xr * yr) - (xi * yi)  \\ no new multiplication here

Which has the side effect of combining the imaginary and real accuracy.
Note that

z = x*y;      \\ 3M product
pr = xr*yr - xi*yi;
pi = xr*yi + xi*yr;
p = pr + I*pi; \\ standard product

abs(z / p - 1)

--> 3.801531906 E-39

so the relative error as a complex number is OK. Not so for the real and
imaginary parts individually.

I am not sure what can be done here. Obviously 3M is less stable when real
and imaginary parts of the result have different magnitudes. It is also
add << mul, that's a 33% speedup.

Actually, the speed situation is not clear-cut: disregarding the accuracy
problem, a general complex multiplication with 3M is, on my machine
[ an aging P4 1.6MHz ]

-- a little slower (!) with the GMP kernel until ~1000 decimal digits,
where it becomes progressively faster up to the expected 33% around 10000
digits.
-- about 33% faster with the native PARI kernel, from the start

I am not sure what is the correct solution here (besides precisely deocumenting the problem...).

Any thoughts ?

Best wishes,

Karim.
--
Karim Belabas                     Tel: (+33) (0)1 69 15 57 48
Dep. de Mathematiques, Bat. 425   Fax: (+33) (0)1 69 15 60 19
Universite Paris-Sud              http://www.math.u-psud.fr/~belabas/
F-91405 Orsay (France)            http://pari.math.u-bordeaux.fr/  [PARI/GP]

```