Karim Belabas on Wed, 13 Apr 2005 20:16:32 +0200 |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
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 "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]