]>
Commit | Line | Data |
---|---|---|
6de9cd9a DN |
1 | /* Complex exponential functions |
2 | Copyright 2002, 2004 Free Software Foundation, Inc. | |
3 | Contributed by Paul Brook <paul@nowt.org> | |
4 | ||
57dea9f6 | 5 | This file is part of the GNU Fortran 95 runtime library (libgfortran). |
6de9cd9a DN |
6 | |
7 | Libgfortran is free software; you can redistribute it and/or | |
57dea9f6 | 8 | modify it under the terms of the GNU General Public |
6de9cd9a | 9 | License as published by the Free Software Foundation; either |
57dea9f6 TM |
10 | version 2 of the License, or (at your option) any later version. |
11 | ||
12 | In addition to the permissions in the GNU General Public License, the | |
13 | Free Software Foundation gives you unlimited permission to link the | |
14 | compiled version of this file into combinations with other programs, | |
15 | and to distribute those combinations without any restriction coming | |
16 | from the use of this file. (The General Public License restrictions | |
17 | do apply in other respects; for example, they cover modification of | |
18 | the file, and distribution when not linked into a combine | |
19 | executable.) | |
6de9cd9a DN |
20 | |
21 | Libgfortran is distributed in the hope that it will be useful, | |
22 | but WITHOUT ANY WARRANTY; without even the implied warranty of | |
23 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
57dea9f6 | 24 | GNU General Public License for more details. |
6de9cd9a | 25 | |
57dea9f6 TM |
26 | You should have received a copy of the GNU General Public |
27 | License along with libgfortran; see the file COPYING. If not, | |
6de9cd9a DN |
28 | write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, |
29 | Boston, MA 02111-1307, USA. */ | |
30 | #include <math.h> | |
31 | #include "libgfortran.h" | |
32 | ||
33 | ||
34 | /* z = a + ib */ | |
35 | /* Absolute value. */ | |
36 | GFC_REAL_4 | |
37 | cabsf (GFC_COMPLEX_4 z) | |
38 | { | |
39 | return hypotf (REALPART (z), IMAGPART (z)); | |
40 | } | |
41 | ||
1e38f159 PB |
42 | /* Complex argument. The angle made with the +ve real axis. |
43 | Range -pi-pi. */ | |
6de9cd9a DN |
44 | GFC_REAL_4 |
45 | cargf (GFC_COMPLEX_4 z) | |
46 | { | |
47 | GFC_REAL_4 arg; | |
48 | ||
1e38f159 | 49 | return atan2f (IMAGPART (z), REALPART (z)); |
6de9cd9a DN |
50 | } |
51 | ||
52 | /* exp(z) = exp(a)*(cos(b) + isin(b)) */ | |
53 | GFC_COMPLEX_4 | |
54 | cexpf (GFC_COMPLEX_4 z) | |
55 | { | |
56 | GFC_REAL_4 a; | |
57 | GFC_REAL_4 b; | |
58 | GFC_COMPLEX_4 v; | |
59 | ||
60 | a = REALPART (z); | |
61 | b = IMAGPART (z); | |
62 | COMPLEX_ASSIGN (v, cosf (b), sinf (b)); | |
63 | return expf (a) * v; | |
64 | } | |
65 | ||
66 | /* log(z) = log (cabs(z)) + i*carg(z) */ | |
67 | GFC_COMPLEX_4 | |
68 | clogf (GFC_COMPLEX_4 z) | |
69 | { | |
70 | GFC_COMPLEX_4 v; | |
71 | ||
72 | COMPLEX_ASSIGN (v, logf (cabsf (z)), cargf (z)); | |
73 | return v; | |
74 | } | |
75 | ||
76 | /* log10(z) = log10 (cabs(z)) + i*carg(z) */ | |
77 | GFC_COMPLEX_4 | |
78 | clog10f (GFC_COMPLEX_4 z) | |
79 | { | |
80 | GFC_COMPLEX_4 v; | |
81 | ||
82 | COMPLEX_ASSIGN (v, log10f (cabsf (z)), cargf (z)); | |
83 | return v; | |
84 | } | |
85 | ||
86 | /* pow(base, power) = cexp (power * clog (base)) */ | |
87 | GFC_COMPLEX_4 | |
88 | cpowf (GFC_COMPLEX_4 base, GFC_COMPLEX_4 power) | |
89 | { | |
90 | return cexpf (power * clogf (base)); | |
91 | } | |
92 | ||
93 | /* sqrt(z). Algorithm pulled from glibc. */ | |
94 | GFC_COMPLEX_4 | |
95 | csqrtf (GFC_COMPLEX_4 z) | |
96 | { | |
97 | GFC_REAL_4 re; | |
98 | GFC_REAL_4 im; | |
99 | GFC_COMPLEX_4 v; | |
100 | ||
101 | re = REALPART (z); | |
102 | im = IMAGPART (z); | |
103 | if (im == 0.0) | |
104 | { | |
105 | if (re < 0.0) | |
106 | { | |
107 | COMPLEX_ASSIGN (v, 0.0, copysignf (sqrtf (-re), im)); | |
108 | } | |
109 | else | |
110 | { | |
111 | COMPLEX_ASSIGN (v, fabsf (sqrt (re)), | |
112 | copysignf (0.0, im)); | |
113 | } | |
114 | } | |
115 | else if (re == 0.0) | |
116 | { | |
117 | GFC_REAL_4 r; | |
118 | ||
119 | r = sqrtf (0.5 * fabs (im)); | |
120 | ||
121 | COMPLEX_ASSIGN (v, copysignf (r, im), r); | |
122 | } | |
123 | else | |
124 | { | |
125 | GFC_REAL_4 d, r, s; | |
126 | ||
127 | d = hypotf (re, im); | |
128 | /* Use the identity 2 Re res Im res = Im x | |
129 | to avoid cancellation error in d +/- Re x. */ | |
130 | if (re > 0) | |
131 | { | |
132 | r = sqrtf (0.5 * d + 0.5 * re); | |
133 | s = (0.5 * im) / r; | |
134 | } | |
135 | else | |
136 | { | |
137 | s = sqrtf (0.5 * d - 0.5 * re); | |
138 | r = fabsf ((0.5 * im) / s); | |
139 | } | |
140 | ||
141 | COMPLEX_ASSIGN (v, r, copysignf (s, im)); | |
142 | } | |
143 | return v; | |
144 | } | |
145 |