diff options
| author | schnetter <schnetter@17b73243-c579-4c4c-a9d2-2d5706c11dac> | 2005-05-15 11:48:38 +0000 |
|---|---|---|
| committer | schnetter <schnetter@17b73243-c579-4c4c-a9d2-2d5706c11dac> | 2005-05-15 11:48:38 +0000 |
| commit | aaa09178ea02c19bc06f9c61a8a49b5a600e38a0 (patch) | |
| tree | 056644335d2f204195d5b46ce9cfeb29e8cb35e1 /src/main/Complex.c | |
| parent | 659355e21e9f6874aec4997c58eccbc158e853cf (diff) | |
Prevent unnecessary overflow or underflow in complex division by
rescaling the arguments. This is a standard algorithm.
git-svn-id: http://svn.cactuscode.org/flesh/trunk@4050 17b73243-c579-4c4c-a9d2-2d5706c11dac
Diffstat (limited to 'src/main/Complex.c')
| -rw-r--r-- | src/main/Complex.c | 18 |
1 files changed, 14 insertions, 4 deletions
diff --git a/src/main/Complex.c b/src/main/Complex.c index 9415c2f2..e0f6b4d0 100644 --- a/src/main/Complex.c +++ b/src/main/Complex.c @@ -297,7 +297,11 @@ cctk_complex CCTK_Cmplx##Mul (cctk_complex a, cctk_complex b) \ @date Sat Dec 4 12:11:04 1999 @author Gabrielle Allen @desc - Divides two complex numbers + Divide two complex numbers. + Make sure that the intermediate values do not overflow. + See e.g. the C99 compliant implementation in libgcc + (which ships with gcc). + This implementation here does not handle nan and inf. @enddesc @var a @@ -319,13 +323,19 @@ cctk_complex CCTK_Cmplx##Mul (cctk_complex a, cctk_complex b) \ #define DEFINE_CCTK_CMPLX_DIV(CCTK_Cmplx, cctk_real, cctk_complex) \ cctk_complex CCTK_Cmplx##Div (cctk_complex a, cctk_complex b) \ { \ - cctk_real factor; \ + cctk_real afact, bfact, factor; \ cctk_complex result; \ \ \ + afact = fabs(a.Re) + fabs(a.Im); \ + a.Re /= afact; \ + a.Im /= afact; \ + bfact = fabs(b.Re) + fabs(b.Im); \ + b.Re /= bfact; \ + b.Im /= bfact; \ factor = b.Re*b.Re + b.Im*b.Im; \ - result.Re = (a.Re*b.Re + a.Im*b.Im) / factor; \ - result.Im = (a.Im*b.Re - a.Re*b.Im) / factor; \ + result.Re = afact / bfact * (a.Re*b.Re + a.Im*b.Im) / factor; \ + result.Im = afact / bfact * (a.Im*b.Re - a.Re*b.Im) / factor; \ \ return (result); \ } |