« back to all changes in this revision
Viewing changes to bc-lib/complex.bc
- browse files at revision 6
- compare with another revision
- download diff
- download tarball
- view history from revision 6
- Committer: Jean-Michel Juzan
- Date: 2009-04-01 12:03:01 UTC
- Revision ID: jean-michel@juzan.org-20090401120301-j6g6xj1z0561zdg3
copie de la branche principale
added
removed
bc-lib/complex.bc
1
### Complex.BC - Rudimentary complex number handling for GNU BC
2
3
## Not to be regarded as suitable for any purpose
4
## Not guaranteed to return correct answers
5
6
# Read a digit from the decimal representation of a number
7
# read_digit(3210.987, -1) returns 9 as that is the digit in
8
# the 10^-1 position.
9
define read_digit(number, place) {
10
auto os,d;os = scale
11
number*=10^-place
12
scale=0
13
d = (number%10)/1
14
scale=os
15
return(d)
16
}
17
18
# MakeComplex: Turn two numbers into one by interleaving the digits
19
define makecomplex(r,i) {
20
auto c, rd, id, ri, ii, j, sign
21
22
r /= 1
23
i /= 1
24
ri = length(r)-scale
25
ii = length(i)-scale
26
27
sign = 1
28
if (r<0){r*=-1;sign+=2}
29
if (i<0){i*=-1;sign+=1}
30
31
if(ri<ii)ri=ii
32
scale *= 2
33
34
c = 0/1
35
36
for(j=-scale;j<ri;j++) {
37
rd = read_digit(r, j)
38
id = read_digit(i, j)
39
c += (rd*10+id)*100^j
40
}
41
42
c += 100^j*sign
43
44
scale /= 2
45
return(c)
46
47
}
48
49
# workhorse for real & imag
50
define realimag_(c,f) {
51
auto x, sign, cs, ci, j, os
52
os = scale
53
cs = scale(c);ci = length(c)-cs-1
54
55
sign = read_digit(c,ci)
56
57
scale = cs/2
58
59
x = 0
60
for(j=-scale;j*2<ci;j++) {
61
x += 10^j * read_digit(c,j*2+f)
62
}
63
64
if ( (sign == 2+f) || (sign == 4) )x*=-1
65
66
scale = os
67
return(x)
68
}
69
70
# Extract the real part of a generated complex number
71
define real(c) {
72
return realimag_(c,1)
73
}
74
75
# Extract the imaginary part of a generated complex number
76
define imag(c) {
77
return realimag_(c,0)
78
}
79
80
# Print a generated complex number
81
# (normal print functions don't work on complex numbers)
82
define printc(c) {
83
auto r,i
84
r = real(c)/1
85
i = imag(c)/1
86
print r
87
if(i<0){print " -";i*=-1}else{print " +"}
88
print " i*"
89
return(i)
90
}
91
92
## Basic math - ordinary bc operators don't work on these special
93
## interleaved complex numbers
94
95
# Add two complex numbers and return another complex number
96
define cadd(c1,c2) {
97
return makecomplex(real(c1)+real(c2),imag(c1)+imag(c2))
98
}
99
100
# Subtract a complex number from another, returning complex
101
define csub(c1,c2) {
102
return makecomplex(real(c1)-real(c2),imag(c1)-imag(c2))
103
}
104
105
# Multiply two complex, return complex
106
define cmul(c1,c2) {
107
auto r1,r2,i1,i2;
108
r1 = real(c1)
109
i1 = imag(c1)
110
r2 = real(c2)
111
i2 = imag(c2)
112
return makecomplex(r1*r2-i1*i2, r2*i1+r1*i2)
113
}
114
115
# Calculates the absolute value of a complex number
116
# NB: returns a standard bc number and not a new fangled 'complex'
117
define cabs(c) {
118
return sqrt(real(c)^2+imag(c)^2)
119
}
120
121
# Returns the complex conjugate of a complex number
122
# ...is faster than the hashed out return line
123
define cconj(c) {
124
#return makecomplex(real(c), -imag(c))
125
return c - 2*10^(length(c)-scale(c)-1)
126
}
127
128
# Divide one complex by another, returning a third
129
define cdiv(c1,c2) {
130
auto r1,r2,i1,i2,aa;
131
r1 = real(c1)
132
r2 = real(c2)
133
i1 = imag(c1)
134
i2 = imag(c2)
135
aa = r2^2+i2^2
136
return makecomplex((r1*r2+i1*i2)/aa, (r2*i1-r1*i2)/aa)
137
}
138
139
# Negate a complex number (i.e. multiply by -1)
140
# ...faster than cmul(c, makecomplex(-1,0))
141
# ...or the hashed out return line
142
define cneg(c) {
143
auto sign, t;
144
# return makecomplex(-real(c), -imag(c))
145
t = length(c) - scale(c) - 1
146
sign = read_digit(c,t)
147
c += (5-2*sign)*(10^t)
148
return(c)
149
}
150
151
## Basic Math II - fast[er] squaring, square roots and integer power
152
153
# Square a complex number
154
define csquare(c) {
155
auto r,i
156
r = real(c)
157
i = imag(c)
158
return makecomplex(r^2-i^2, 2*r*i)
159
}
160
161
# Find the primary square root of a complex number
162
define csqrt(c) {
163
auto r,i, sr, si
164
r = real(c)
165
i = imag(c)
166
if(i==0){if(r>=0){
167
return(makecomplex(sqrt(r),0))
168
} else {
169
return(makecomplex(0,sqrt(-r)))
170
}}
171
sr = sqrt((sqrt(r^2+i^2)+r)/2)
172
si = i/sr/2
173
return makecomplex(sr,si)
174
}
175
176
# subfunctions for use by cintpow
177
define mod(n,m) {auto s;s=scale;scale=0;n%=m;scale=s;return(n)}
178
define int(n) {auto s;s=scale;scale=0;n/=1;scale=s;return(n)}
179
180
# Raise a complex number [c] to an integer power [n]
181
# NB: n is a standard bc number and not a complex
182
define cintpow(c, n) {
183
auto x
184
x = makecomplex(1,0)
185
n = int(n)
186
if(n<0)return( cdiv(x,cintpow(c,-n)) )
187
if(n==0)return( x )
188
if(n==1)return( c )
189
if(n==2)return( csquare(c) )
190
while(n) {
191
if(mod(n,2)) {
192
x = cmul(x, c) ; n -= 1
193
#print "DEBUG: cintpow - odd_ (",n,")\n"
194
} else {
195
c = csquare(c) ; n = int(n/2)
196
#print "DEBUG: cintpow - even (",n,")\n"
197
}
198
}
199
return x
200
}
Loggerhead is a web-based interface for Breezy
Version: 2.0.1