1
2 /*********************************************************/
3 /* TAUCS */
4 /* Author: Sivan Toledo */
5 /* */
6 /* Simple complex arithmetic routines. */
7 /* They are called if the compiler does not support */
8 /* complex. GCC supports complex, and so do all C99 */
9 /* compilers. */
10 /* */
11 /*********************************************************/
12
13 #include <math.h>
14 #include "taucs.h"
15
16 #ifdef TAUCS_CORE_DOUBLE
taucs_get_nan()17 double taucs_get_nan()
18 {
19 double zero = 0.0;
20 double inf = 1.0 / zero;
21 double nan = inf - inf;
22 return nan;
23 }
24 #endif
25
26 #ifdef TAUCS_CORE_DOUBLE
27 taucs_double taucs_dzero_const = 0.0;
28 taucs_double taucs_done_const = 1.0;
29 taucs_double taucs_dminusone_const = -1.0;
30 #endif
31
32 #ifdef TAUCS_CORE_SINGLE
33 taucs_single taucs_szero_const = 0.0f;
34 taucs_single taucs_sone_const = 1.0f;
35 taucs_single taucs_sminusone_const = -1.0f;
36 #endif
37
38 /*#if defined(__GNUC__) && !defined(TAUCS_CONFIG_GENERIC_COMPLEX)*/
39 #ifdef TAUCS_C99_COMPLEX
40
41 #ifdef TAUCS_CORE_DCOMPLEX
42 taucs_dcomplex taucs_zzero_const = 0.0+0.0*_Complex_I;
43 taucs_dcomplex taucs_zone_const = 1.0+0.0*_Complex_I;
44 taucs_dcomplex taucs_zminusone_const = -1.0+0.0*_Complex_I;
45 #endif
46
47 #ifdef TAUCS_CORE_SCOMPLEX
48 taucs_scomplex taucs_czero_const = 0.0f+0.0f*_Complex_I;
49 taucs_scomplex taucs_cone_const = 1.0f+0.0f*_Complex_I;
50 taucs_scomplex taucs_cminusone_const = -1.0f+0.0f*_Complex_I;
51 #endif
52
53 #else /* TAUCS_C99_COMPLEX */
54
55 #ifdef TAUCS_CORE_DCOMPLEX
56 taucs_dcomplex taucs_zzero_const = { 0.0 , 0.0 };
57 taucs_dcomplex taucs_zone_const = { 1.0 , 0.0 };
58 taucs_dcomplex taucs_zminusone_const = {-1.0 , 0.0 };
59 #endif
60
61 #ifdef TAUCS_CORE_SCOMPLEX
62 taucs_scomplex taucs_czero_const = { 0.0f, 0.0f};
63 taucs_scomplex taucs_cone_const = { 1.0f, 0.0f};
64 taucs_scomplex taucs_cminusone_const = {-1.0f, 0.0f};
65 #endif
66
67 #ifdef TAUCS_CORE_COMPLEX
68
69 taucs_datatype
taucs_dtl(complex_create_fn)70 taucs_dtl(complex_create_fn)(taucs_real_datatype r, taucs_real_datatype i)
71 {
72 taucs_datatype c;
73 taucs_re(c) = r;
74 taucs_im(c) = i;
75 return c;
76 }
77
78 taucs_datatype
taucs_dtl(add_fn)79 taucs_dtl(add_fn)(taucs_datatype a, taucs_datatype b)
80 {
81 taucs_datatype c;
82 taucs_re(c) = taucs_re(a) + taucs_re(b);
83 taucs_im(c) = taucs_im(a) + taucs_im(b);
84 return c;
85 }
86
87 taucs_datatype
taucs_dtl(sub_fn)88 taucs_dtl(sub_fn)(taucs_datatype a, taucs_datatype b)
89 {
90 taucs_datatype c;
91 taucs_re(c) = taucs_re(a) - taucs_re(b);
92 taucs_im(c) = taucs_im(a) - taucs_im(b);
93 return c;
94 }
95
96 taucs_datatype
taucs_dtl(mul_fn)97 taucs_dtl(mul_fn)(taucs_datatype a, taucs_datatype b)
98 {
99 taucs_datatype c;
100 taucs_re(c) = taucs_re(a) * taucs_re(b) - taucs_im(a) * taucs_im(b);
101 taucs_im(c) = taucs_re(a) * taucs_im(b) + taucs_im(a) * taucs_re(b);
102 return c;
103 }
104
105 taucs_datatype
taucs_dtl(div_fn)106 taucs_dtl(div_fn)(taucs_datatype a, taucs_datatype b)
107 {
108 taucs_datatype c;
109 /*double r,den; omer*/
110 taucs_real_datatype r,den;
111
112 if (fabs(taucs_re(b)) >= fabs(taucs_im(b))) {
113 r = taucs_im(b) / taucs_re(b);
114 den = taucs_re(b) + r * taucs_im(b);
115 taucs_re(c) = (taucs_re(a) + r * taucs_im(a))/den;
116 taucs_im(c) = (taucs_im(a) - r * taucs_re(a))/den;
117 } else {
118 r = taucs_re(b) / taucs_im(b);
119 den = taucs_im(b) + r * taucs_re(b);
120 taucs_re(c) = (r * taucs_re(a) + taucs_im(a))/den;
121 taucs_im(c) = (r * taucs_im(a) - taucs_re(a))/den;
122 }
123 return c;
124 }
125
126 taucs_datatype
taucs_dtl(conj_fn)127 taucs_dtl(conj_fn)(taucs_datatype a)
128 {
129 taucs_datatype c;
130 taucs_re(c) = taucs_re(a);
131 taucs_im(c) = - taucs_im(a);
132 return c;
133 }
134
135 taucs_datatype
taucs_dtl(neg_fn)136 taucs_dtl(neg_fn)(taucs_datatype a)
137 {
138 taucs_datatype c;
139 taucs_re(c) = - taucs_re(a);
140 taucs_im(c) = - taucs_im(a);
141 return c;
142 }
143
144 double
taucs_dtl(abs_fn)145 taucs_dtl(abs_fn)(taucs_datatype a)
146 {
147 double x,y,temp;
148
149 #if 1
150 x = fabs(taucs_re(a));
151 y = fabs(taucs_im(a));
152
153 if (x==0.0) return y;
154 if (y==0.0) return x;
155
156 if (x > y) {
157 temp = y/x;
158 return ( x*sqrt(1.0+temp*temp) );
159 } else {
160 temp = x/y;
161 return ( y*sqrt(1.0+temp*temp) );
162 }
163 #else
164 return hypot(taucs_re(a), taucs_im(a));
165 #endif
166 }
167
168 taucs_datatype
taucs_dtl(sqrt_fn)169 taucs_dtl(sqrt_fn)(taucs_datatype a)
170 {
171 taucs_datatype c;
172 double x,y,t;/*,w; omer*/
173 taucs_real_datatype w;
174
175 if (taucs_re(a) == 0.0 && taucs_im(a) == 0.0) {
176 taucs_re(c) = 0.0;
177 taucs_im(c) = 0.0;
178 } else {
179 x = fabs((double) taucs_re(a));
180 y = fabs((double) taucs_im(a));
181 if (x >= y) {
182 t = y/x;
183 w = (taucs_real_datatype )(sqrt(x) * sqrt(0.5 * (1.0 + sqrt(1.0 + t * t))));
184 } else {
185 t = x/y;
186 w = (taucs_real_datatype )(sqrt(y) * sqrt(0.5 * (t + sqrt(1.0 + t * t))));
187 }
188
189 if (taucs_re(a) > 0.0) {
190 taucs_re(c) = w;
191 /*taucs_im(c) = taucs_im(a) / (2.0 * w); omer*/
192 taucs_im(c) = (taucs_real_datatype)(taucs_im(a) / (2.0 * w));
193 } else {
194 x = (taucs_im(a) >= 0.0) ? w : -w;
195 taucs_im(c) = (taucs_real_datatype )x;
196 /*taucs_re(c) = taucs_im(a) / (2.0 * x); omer*/
197 taucs_re(c) = (taucs_real_datatype)(taucs_im(a) / (2.0 * x));
198 }
199 }
200
201 return c;
202 }
203
204 #endif /* TAUCS_C99_COMPLEX */
205
206 #endif /* TAUCS_CORE_COMPLEX */
207
208
209
210
211
212
213