cmplxl.c 7.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461
  1. /* cmplxl.c
  2. *
  3. * Complex number arithmetic
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * typedef struct {
  10. * long double r; real part
  11. * long double i; imaginary part
  12. * }cmplxl;
  13. *
  14. * cmplxl *a, *b, *c;
  15. *
  16. * caddl( a, b, c ); c = b + a
  17. * csubl( a, b, c ); c = b - a
  18. * cmull( a, b, c ); c = b * a
  19. * cdivl( a, b, c ); c = b / a
  20. * cnegl( c ); c = -c
  21. * cmovl( b, c ); c = b
  22. *
  23. *
  24. *
  25. * DESCRIPTION:
  26. *
  27. * Addition:
  28. * c.r = b.r + a.r
  29. * c.i = b.i + a.i
  30. *
  31. * Subtraction:
  32. * c.r = b.r - a.r
  33. * c.i = b.i - a.i
  34. *
  35. * Multiplication:
  36. * c.r = b.r * a.r - b.i * a.i
  37. * c.i = b.r * a.i + b.i * a.r
  38. *
  39. * Division:
  40. * d = a.r * a.r + a.i * a.i
  41. * c.r = (b.r * a.r + b.i * a.i)/d
  42. * c.i = (b.i * a.r - b.r * a.i)/d
  43. * ACCURACY:
  44. *
  45. * In DEC arithmetic, the test (1/z) * z = 1 had peak relative
  46. * error 3.1e-17, rms 1.2e-17. The test (y/z) * (z/y) = 1 had
  47. * peak relative error 8.3e-17, rms 2.1e-17.
  48. *
  49. * Tests in the rectangle {-10,+10}:
  50. * Relative error:
  51. * arithmetic function # trials peak rms
  52. * DEC cadd 10000 1.4e-17 3.4e-18
  53. * IEEE cadd 100000 1.1e-16 2.7e-17
  54. * DEC csub 10000 1.4e-17 4.5e-18
  55. * IEEE csub 100000 1.1e-16 3.4e-17
  56. * DEC cmul 3000 2.3e-17 8.7e-18
  57. * IEEE cmul 100000 2.1e-16 6.9e-17
  58. * DEC cdiv 18000 4.9e-17 1.3e-17
  59. * IEEE cdiv 100000 3.7e-16 1.1e-16
  60. */
  61. /* cmplx.c
  62. * complex number arithmetic
  63. */
  64. /*
  65. Cephes Math Library Release 2.3: March, 1995
  66. Copyright 1984, 1995 by Stephen L. Moshier
  67. */
  68. #include <math.h>
  69. /*
  70. typedef struct
  71. {
  72. long double r;
  73. long double i;
  74. }cmplxl;
  75. */
  76. #ifdef ANSIPROT
  77. extern long double fabsl ( long double );
  78. extern long double cabsl ( cmplxl * );
  79. extern long double sqrtl ( long double );
  80. extern long double atan2l ( long double, long double );
  81. extern long double cosl ( long double );
  82. extern long double sinl ( long double );
  83. extern long double frexpl ( long double, int * );
  84. extern long double ldexpl ( long double, int );
  85. extern int isnanl ( long double );
  86. void cdivl ( cmplxl *, cmplxl *, cmplxl * );
  87. void caddl ( cmplxl *, cmplxl *, cmplxl * );
  88. #else
  89. long double fabsl(), cabsl(), sqrtl(), atan2l(), cosl(), sinl();
  90. long double frexpl(), ldexpl();
  91. int isnanl();
  92. void cdivl(), caddl();
  93. #endif
  94. extern double MAXNUML, MACHEPL, PIL, PIO2L, INFINITYL, NANL;
  95. cmplx czerol = {0.0L, 0.0L};
  96. cmplx conel = {1.0L, 0.0L};
  97. /* c = b + a */
  98. void caddl( a, b, c )
  99. register cmplxl *a, *b;
  100. cmplxl *c;
  101. {
  102. c->r = b->r + a->r;
  103. c->i = b->i + a->i;
  104. }
  105. /* c = b - a */
  106. void csubl( a, b, c )
  107. register cmplxl *a, *b;
  108. cmplxl *c;
  109. {
  110. c->r = b->r - a->r;
  111. c->i = b->i - a->i;
  112. }
  113. /* c = b * a */
  114. void cmull( a, b, c )
  115. register cmplxl *a, *b;
  116. cmplxl *c;
  117. {
  118. long double y;
  119. y = b->r * a->r - b->i * a->i;
  120. c->i = b->r * a->i + b->i * a->r;
  121. c->r = y;
  122. }
  123. /* c = b / a */
  124. void cdivl( a, b, c )
  125. register cmplxl *a, *b;
  126. cmplxl *c;
  127. {
  128. long double y, p, q, w;
  129. y = a->r * a->r + a->i * a->i;
  130. p = b->r * a->r + b->i * a->i;
  131. q = b->i * a->r - b->r * a->i;
  132. if( y < 1.0L )
  133. {
  134. w = MAXNUML * y;
  135. if( (fabsl(p) > w) || (fabsl(q) > w) || (y == 0.0L) )
  136. {
  137. c->r = INFINITYL;
  138. c->i = INFINITYL;
  139. mtherr( "cdivl", OVERFLOW );
  140. return;
  141. }
  142. }
  143. c->r = p/y;
  144. c->i = q/y;
  145. }
  146. /* b = a
  147. Caution, a `short' is assumed to be 16 bits wide. */
  148. void cmovl( a, b )
  149. void *a, *b;
  150. {
  151. register short *pa, *pb;
  152. int i;
  153. pa = (short *) a;
  154. pb = (short *) b;
  155. i = 16;
  156. do
  157. *pb++ = *pa++;
  158. while( --i );
  159. }
  160. void cnegl( a )
  161. register cmplxl *a;
  162. {
  163. a->r = -a->r;
  164. a->i = -a->i;
  165. }
  166. /* cabsl()
  167. *
  168. * Complex absolute value
  169. *
  170. *
  171. *
  172. * SYNOPSIS:
  173. *
  174. * long double cabsl();
  175. * cmplxl z;
  176. * long double a;
  177. *
  178. * a = cabs( &z );
  179. *
  180. *
  181. *
  182. * DESCRIPTION:
  183. *
  184. *
  185. * If z = x + iy
  186. *
  187. * then
  188. *
  189. * a = sqrt( x**2 + y**2 ).
  190. *
  191. * Overflow and underflow are avoided by testing the magnitudes
  192. * of x and y before squaring. If either is outside half of
  193. * the floating point full scale range, both are rescaled.
  194. *
  195. *
  196. * ACCURACY:
  197. *
  198. * Relative error:
  199. * arithmetic domain # trials peak rms
  200. * DEC -30,+30 30000 3.2e-17 9.2e-18
  201. * IEEE -10,+10 100000 2.7e-16 6.9e-17
  202. */
  203. /*
  204. Cephes Math Library Release 2.1: January, 1989
  205. Copyright 1984, 1987, 1989 by Stephen L. Moshier
  206. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  207. */
  208. /*
  209. typedef struct
  210. {
  211. long double r;
  212. long double i;
  213. }cmplxl;
  214. */
  215. #ifdef UNK
  216. #define PRECL 32
  217. #define MAXEXPL 16384
  218. #define MINEXPL -16384
  219. #endif
  220. #ifdef IBMPC
  221. #define PRECL 32
  222. #define MAXEXPL 16384
  223. #define MINEXPL -16384
  224. #endif
  225. #ifdef MIEEE
  226. #define PRECL 32
  227. #define MAXEXPL 16384
  228. #define MINEXPL -16384
  229. #endif
  230. long double cabsl( z )
  231. register cmplxl *z;
  232. {
  233. long double x, y, b, re, im;
  234. int ex, ey, e;
  235. #ifdef INFINITIES
  236. /* Note, cabs(INFINITY,NAN) = INFINITY. */
  237. if( z->r == INFINITYL || z->i == INFINITYL
  238. || z->r == -INFINITYL || z->i == -INFINITYL )
  239. return( INFINITYL );
  240. #endif
  241. #ifdef NANS
  242. if( isnanl(z->r) )
  243. return(z->r);
  244. if( isnanl(z->i) )
  245. return(z->i);
  246. #endif
  247. re = fabsl( z->r );
  248. im = fabsl( z->i );
  249. if( re == 0.0 )
  250. return( im );
  251. if( im == 0.0 )
  252. return( re );
  253. /* Get the exponents of the numbers */
  254. x = frexpl( re, &ex );
  255. y = frexpl( im, &ey );
  256. /* Check if one number is tiny compared to the other */
  257. e = ex - ey;
  258. if( e > PRECL )
  259. return( re );
  260. if( e < -PRECL )
  261. return( im );
  262. /* Find approximate exponent e of the geometric mean. */
  263. e = (ex + ey) >> 1;
  264. /* Rescale so mean is about 1 */
  265. x = ldexpl( re, -e );
  266. y = ldexpl( im, -e );
  267. /* Hypotenuse of the right triangle */
  268. b = sqrtl( x * x + y * y );
  269. /* Compute the exponent of the answer. */
  270. y = frexpl( b, &ey );
  271. ey = e + ey;
  272. /* Check it for overflow and underflow. */
  273. if( ey > MAXEXPL )
  274. {
  275. mtherr( "cabsl", OVERFLOW );
  276. return( INFINITYL );
  277. }
  278. if( ey < MINEXPL )
  279. return(0.0L);
  280. /* Undo the scaling */
  281. b = ldexpl( b, e );
  282. return( b );
  283. }
  284. /* csqrtl()
  285. *
  286. * Complex square root
  287. *
  288. *
  289. *
  290. * SYNOPSIS:
  291. *
  292. * void csqrtl();
  293. * cmplxl z, w;
  294. *
  295. * csqrtl( &z, &w );
  296. *
  297. *
  298. *
  299. * DESCRIPTION:
  300. *
  301. *
  302. * If z = x + iy, r = |z|, then
  303. *
  304. * 1/2
  305. * Im w = [ (r - x)/2 ] ,
  306. *
  307. * Re w = y / 2 Im w.
  308. *
  309. *
  310. * Note that -w is also a square root of z. The root chosen
  311. * is always in the upper half plane.
  312. *
  313. * Because of the potential for cancellation error in r - x,
  314. * the result is sharpened by doing a Heron iteration
  315. * (see sqrt.c) in complex arithmetic.
  316. *
  317. *
  318. *
  319. * ACCURACY:
  320. *
  321. * Relative error:
  322. * arithmetic domain # trials peak rms
  323. * DEC -10,+10 25000 3.2e-17 9.6e-18
  324. * IEEE -10,+10 100000 3.2e-16 7.7e-17
  325. *
  326. * 2
  327. * Also tested by csqrt( z ) = z, and tested by arguments
  328. * close to the real axis.
  329. */
  330. void csqrtl( z, w )
  331. cmplxl *z, *w;
  332. {
  333. cmplxl q, s;
  334. long double x, y, r, t;
  335. x = z->r;
  336. y = z->i;
  337. if( y == 0.0L )
  338. {
  339. if( x < 0.0L )
  340. {
  341. w->r = 0.0L;
  342. w->i = sqrtl(-x);
  343. return;
  344. }
  345. else
  346. {
  347. w->r = sqrtl(x);
  348. w->i = 0.0L;
  349. return;
  350. }
  351. }
  352. if( x == 0.0L )
  353. {
  354. r = fabsl(y);
  355. r = sqrtl(0.5L*r);
  356. if( y > 0.0L )
  357. w->r = r;
  358. else
  359. w->r = -r;
  360. w->i = r;
  361. return;
  362. }
  363. /* Approximate sqrt(x^2+y^2) - x = y^2/2x - y^4/24x^3 + ... .
  364. * The relative error in the first term is approximately y^2/12x^2 .
  365. */
  366. if( (fabsl(y) < 2.e-4L * fabsl(x))
  367. && (x > 0) )
  368. {
  369. t = 0.25L*y*(y/x);
  370. }
  371. else
  372. {
  373. r = cabsl(z);
  374. t = 0.5L*(r - x);
  375. }
  376. r = sqrtl(t);
  377. q.i = r;
  378. q.r = y/(2.0L*r);
  379. /* Heron iteration in complex arithmetic */
  380. cdivl( &q, z, &s );
  381. caddl( &q, &s, w );
  382. w->r *= 0.5L;
  383. w->i *= 0.5L;
  384. cdivl( &q, z, &s );
  385. caddl( &q, &s, w );
  386. w->r *= 0.5L;
  387. w->i *= 0.5L;
  388. }
  389. long double hypotl( x, y )
  390. long double x, y;
  391. {
  392. cmplxl z;
  393. z.r = x;
  394. z.i = y;
  395. return( cabsl(&z) );
  396. }