3M multiplication for t_COMPLEX

[Karim Belabas <Karim.Belabas@math.u-psud.fr>]

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         
"obviously" faster: 3 mul + 5 add instead of 4 mul + 2 add, assuming            
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]