0

How do I go about calculating exp(i * x) where |x| >> 2 pi?

I know that it can be done with an arbitrary precision library, but is there an algorithm to do it stably (without roundoff errors) that doesn't require a library?

Andrew Spott
  • 3,457
  • 8
  • 33
  • 59

1 Answers1

4

exp(i*x) is cos(x) + i*sin(x). A good math library will calculate cos(x) and sin(x) correctly even when x is large. E.g., you should get correct results with the OS X or iOS standard math library, within a few ULP.

In C, this should work:

#include <complex.h>
…
    double complex Y = cexp(I * x);
Eric Postpischil
  • 195,579
  • 13
  • 168
  • 312
  • how large a x would a good library handle? I thought that once ulp(x) > pi, most library would not bother... – aka.nice Dec 19 '12 at 20:13
  • 1
    @aka.nice: A good library handles all values. OS X and iOS handle all values. It is not proper for the library to decide “Well, your input is big, you must have gotten some error computing it, so I will not try very hard.” All the library knows is that it was passed a floating-point object that very definitely (it is specified in a very good standard) exactly represents a certain value, and its job is to return the cosine or sine of that value. It may well be that the user was careful in producing the input, and the library should avoid introducing new error. – Eric Postpischil Dec 19 '12 at 20:17
  • I like this principle which was yet behind Sun math library. I'm going to try it for fun. Is there any step needed with gcc to make sure that proper library is used rather than broken x86 hardwired cos/sin? – aka.nice Dec 19 '12 at 20:23
  • @aka.nice: No, the default compilation for OS X or iOS should give you the right stuff. Example: `sin((2-0x1p-52)*0x1p1023)` (the largest finite double) returns 0.00496195478918406117363470997361218906007707118988037109375, which is within .72 ULP of the exact sine (OS X 10.8.2, MacPro4,1). So there is a reduction of over a thousand bits in there. And it is still fast. – Eric Postpischil Dec 19 '12 at 20:28
  • I confirm it works very well on OSX with no specific gcc flag – aka.nice Dec 19 '12 at 20:28
  • I also want to brag :-) For Eric's test case, CUDA 5.0 returns a result of 4.9619547891840612e-3 (0x3f7452fc98b34e96), with an error < 0.712 ulps. As Eric says, a good math library will handle large inputs to trigonometric functions correctly. Because the granularity of the results increases with the magnitude of the argument, it is usually not a good idea to let the arguments to trig functions grow very large, even with a good math library. – njuffa Dec 20 '12 at 03:19