1 /*
2 * Copyright (c) 2007 - 2015 Joseph Gaeddert
3 *
4 * Permission is hereby granted, free of charge, to any person obtaining a copy
5 * of this software and associated documentation files (the "Software"), to deal
6 * in the Software without restriction, including without limitation the rights
7 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
8 * copies of the Software, and to permit persons to whom the Software is
9 * furnished to do so, subject to the following conditions:
10 *
11 * The above copyright notice and this permission notice shall be included in
12 * all copies or substantial portions of the Software.
13 *
14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
19 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
20 * THE SOFTWARE.
21 */
22
23 //
24 // Complex math functions (trig, log, exp, etc.)
25 //
26
27 #include <stdio.h>
28 #include <stdlib.h>
29 #include <math.h>
30
31 #include "liquid.internal.h"
32
33 // complex square root
liquid_csqrtf(float complex _z)34 float complex liquid_csqrtf(float complex _z)
35 {
36 float r = cabsf(_z); // magnitude of _z
37 float a = crealf(_z); // real component of _z
38
39 float re = sqrtf(0.5f*(r+a)); // real component of return value
40 float im = sqrtf(0.5f*(r-a)); // imag component of return value
41
42 // return value, retaining sign of imaginary component
43 return cimagf(_z) > 0 ? re + _Complex_I*im :
44 re - _Complex_I*im;
45 }
46
47 // complex exponent
liquid_cexpf(float complex _z)48 float complex liquid_cexpf(float complex _z)
49 {
50 float r = expf( crealf(_z) );
51 float re = cosf( cimagf(_z) );
52 float im = sinf( cimagf(_z) );
53
54 return r * ( re + _Complex_I*im );
55 }
56
57 // complex logarithm
liquid_clogf(float complex _z)58 float complex liquid_clogf(float complex _z)
59 {
60 return logf(cabsf(_z)) + _Complex_I*cargf(_z);
61 }
62
63 // complex arcsin
liquid_casinf(float complex _z)64 float complex liquid_casinf(float complex _z)
65 {
66 return 0.5f*M_PI - liquid_cacosf(_z);
67 }
68
69 // complex arccos
liquid_cacosf(float complex _z)70 float complex liquid_cacosf(float complex _z)
71 {
72 // return based on quadrant
73 int sign_i = crealf(_z) > 0;
74 int sign_q = cimagf(_z) > 0;
75
76 if (sign_i == sign_q) {
77 return -_Complex_I*liquid_clogf(_z + liquid_csqrtf(_z*_z - 1.0f));
78 } else {
79 return -_Complex_I*liquid_clogf(_z - liquid_csqrtf(_z*_z - 1.0f));
80 }
81
82 // should never get to this state
83 return 0.0f;
84 }
85
86 // complex arctan
liquid_catanf(float complex _z)87 float complex liquid_catanf(float complex _z)
88 {
89 float complex t0 = 1.0f - _Complex_I*_z;
90 float complex t1 = 1.0f + _Complex_I*_z;
91
92 return 0.5f*_Complex_I*liquid_clogf( t0 / t1 );
93 }
94
95 // approximation to cargf() but faster
liquid_cargf_approx(float complex _x)96 float liquid_cargf_approx(float complex _x)
97 {
98 float theta;
99 float xi = crealf(_x);
100 float xq = cimagf(_x);
101
102 if (xi == 0.0f) {
103 if (xq == 0.0f)
104 return 0.0f;
105 return xq > 0.0f ? M_PI_2 : -M_PI_2;
106 } else {
107 theta = xq / fabsf(xi);
108 }
109
110 if (theta > M_PI_2)
111 theta = M_PI_2;
112 else if (theta < -M_PI_2)
113 theta = -M_PI_2;
114 return theta;
115 }
116
117