| Title: | Mathematics at DEC |
| Moderator: | RUSURE::EDP |
| Created: | Mon Feb 03 1986 |
| Last Modified: | Fri Jun 06 1997 |
| Last Successful Update: | Fri Jun 06 1997 |
| Number of topics: | 2083 |
| Total number of notes: | 14613 |
I have a program which does COMPLEX calculations. It gives me the
following result:
(612323399.573678,1.000000000000003E+025)
^ ^
| |
| |
Real part Imaginary part
The result should be:
(0.0,1000000000000003E+025)
What am I doing wrong? Is there a bug in the Fortran compiler?
If I replace the x**0.5 expression with CDSQRT(x) then I get the correct
result. What is the reason for this? Why are both results so far appart?
-----------------------------------------------------------------
program test
complex*16 x, y
real*8 a /-1E50/
real*8 b /0.D0/
x = DCMPLX (a,b)
y = x**0.5
type *,y
call exit
end
-----------------------------------------------------------------
FORTRAN V5.0 (same result on V4.8)
Any help appreciated.
M.
| T.R | Title | User | Personal Name | Date | Lines |
|---|---|---|---|---|---|
| 1033.1 | HPSTEK::XIA | Mon Feb 27 1989 18:31 | 6 | ||
re -1,
The relative error in the real part is only 10E-14. This is pretty
good considering you are only using double precision.
Eugene
| |||||
| 1033.2 | CTCADM::ROTH | If you plant ice you'll harvest wind | Tue Feb 28 1989 06:45 | 11 | |
Mathematically, X**0.5 and SQRT(X) are the same, but in complex G floating
FORTRAN will call the complex exponentiation library routine for the
former and will call the complex square root routine for the latter.
Since these both involve approximations there can be a slight relative
error.
In real E and D floating, the compiler will special case exponentiation
by 0.5 to a call to SQRT, so the results are identical in that case.
- Jim
| |||||
| 1033.3 | Strange, but fairly accurate! | VINO::EKLUND | Dave Eklund | Tue Feb 28 1989 11:35 | 21 |
I did not understand .1 at first, but it is basically correct!
While I have not looked at the routines involved, a good guess is
that you have found your way into a very general purpose exponentiation
routine, and have obtained a result that is accurate to within a
few bits (despite the way it looks). Without getting into a long
numerical analysis answer, it may be easier to understand if you
think about the answer as not an X and Y coordinate answer, but
one expressed as magnitude and angle (in the complex plane). When
viewed that way, the angle is VERY small (as it should be), but
not zero (as intuition tells us it must be). One other small comment
should be made - your -1E50 is not exact; it cannot be exactly
represented within the real*8 world. This is not part of why your
answer is the way it is, just a comment. While I doubt that you
have reason to cause the exponentiation routine to change, it is
possible to special case such arguments with a zero imaginary part,
to get a more "expected" result. However, the answer you got is
numerically almost as good - it just doesn't fit well with what
we "know" is the "real" answer!
Dave E
| |||||