-
Notifications
You must be signed in to change notification settings - Fork 10
/
complex.h
157 lines (134 loc) · 3.66 KB
/
complex.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
/* complex.h define type complex and arithmetic operations */
#ifndef COMPLEX_H_
#define COMPLEX_H_
typedef struct {double re; double im;} complex;
complex cmplx(double a, double b); /* create a complex number */
double real (complex a); /* get real part */
double imag (complex a); /* get imaginary part */
complex cxadd(complex a, complex b); /* add two complex numbers */
complex cxadd4(complex a, complex b, complex c, complex d);
complex cxsub(complex a, complex b); /* subtract b from a, complex */
complex cxmul(complex a, complex b); /* multiply two complex numbers */
complex cmuld(complex a, double b); /* multiply complex times double */
complex cxdiv(complex a, complex b); /* divide complex by complex */
complex cdivd(complex a, double b); /* divide complex by double */
double cxabs(complex a); /* return magnitude of complex munber */
double cxarg(complex a); /* return argument of complex number */
complex cxneg(complex a); /* negate a complex number */
complex cxsqr(complex a); /* complex square */
complex cxsqrt(complex a); /* complex square root */
complex cxsin(complex a); /* complex sin */
complex cxcos(complex a); /* complex cos */
complex cmplx(double a, double b) /* make a complex number */
{
complex c;
c.re = a;
c.im = b;
return c;
}
double real(complex a) /* get real part */
{
return a.re;
}
double imag(complex a) /* get imaginary part */
{
return a.im;
}
complex cxadd(complex a, complex b) /* add two complex numbers */
{
complex c;
c.re = a.re+b.re;
c.im = a.im+b.im;
return c;
}
complex cxadd4(complex a, complex b, complex c, complex d)
{
complex e;
e.re = a.re+b.re+c.re+d.re;
e.im = a.im+b.im+c.im+d.im;
return e;
}
complex cxsub(complex a, complex b) /* subtract b from a, complex */
{
complex c;
c.re = a.re-b.re;
c.im = a.im-b.im;
return c;
}
complex cxmul(complex a, complex b) /* multiply two complex numbers */
{
complex c;
c.re = a.re*b.re - a.im*b.im;
c.im = a.re*b.im + a.im*b.re;
return c;
}
complex cmuld(complex a, double b) /* multiply complex times double */
{
complex c;
c.re = a.re*b;
c.im = a.im*b;
return c;
}
complex cxdiv(complex a, complex b) /* divide complex by complex */
{
complex c;
double r;
r = b.re*b.re+b.im*b.im;
c.re = (a.re*b.re+a.im*b.im)/r;
c.im = (a.im*b.re-a.re*b.im)/r;
return c;
}
complex cdivd(complex a, double b) /* divide complex by double */
{
complex c;
c.re = a.re/b;
c.im = a.im/b;
return c;
}
double cxabs(complex a) /* return magnitude of complex munber */
{
return sqrt(a.re*a.re+a.im*a.im);
}
double cxarg(complex a) /* return argument of complex munber */
{
return atan2(a.im,a.re);
}
complex cxneg(complex a) /* negate a complex */
{
complex c;
c.re = -a.re;
c.im = -a.im;
return c;
}
complex cxsqr(complex a)
{
complex c;
c.re = a.re*a.re - a.im*a.im;
c.im = 2 * a.re * a.im;
return c;
}
complex cxsqrt(complex a)
{
complex c;
c.re = sqrt((cxabs(a)+real(a))/2);
if(imag(a)<0)
c.im = -sqrt((cxabs(a)-real(a))/2);
else
c.im = sqrt((cxabs(a)-real(a))/2);
return c;
}
complex cxsin(complex a)
{
complex c;
c.re = sin(a.re)*cosh(a.im);
c.im = cos(a.re)*sinh(a.im);
return c;
}
complex cxcos(complex a)
{
complex c;
c.re = cos(a.re)*cosh(a.im);
c.im = -sin(a.re)*sinh(a.im);
return c;
}
#endif /* COMPLEX_H_ */