|Bill Allombert on Thu, 14 Apr 2005 15:09:04 +0200|
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
|Re: 3M multiplication for t_COMPLEX|
On Wed, Apr 13, 2005 at 08:09:49PM +0200, Karim Belabas wrote: > 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. That might be completly irrelevant, but a standard trick here is to pick the 3M formula that does not cause cancellation amongst the 4 possible formulas below: Im = (xr + xi) * (yr + yi) - (xr * yr + xi * yi), Re accordingly Im = (xr - xi) * (yi - yr) + (xr * yr + xi * yi), Re accordingly Re = (xr + xi) * (yr - yi) + (xr * yi - xi * yr), Im accordingly Re = (xr - xi) * (yr + yi) + (xi * yr - xr * yi), Im accordingly In the case at hand, the 2 last formulas will give the good result. What do mpc here ? Cheers, Bill.