102 lines
2.2 KiB
Text
102 lines
2.2 KiB
Text
#define PI 3.141592653589793238462643383279502884197169399375105821
|
|
#define MAXITER 12
|
|
|
|
'---------------------------------------
|
|
' complex numbers and their arithmetic
|
|
'---------------------------------------
|
|
|
|
type complex
|
|
r as double
|
|
i as double
|
|
end type
|
|
|
|
function conj( a as complex ) as complex
|
|
dim as complex c
|
|
c.r = a.r
|
|
c.i = -a.i
|
|
return c
|
|
end function
|
|
|
|
operator + ( a as complex, b as complex ) as complex
|
|
dim as complex c
|
|
c.r = a.r + b.r
|
|
c.i = a.i + b.i
|
|
return c
|
|
end operator
|
|
|
|
operator - ( a as complex, b as complex ) as complex
|
|
dim as complex c
|
|
c.r = a.r - b.r
|
|
c.i = a.i - b.i
|
|
return c
|
|
end operator
|
|
|
|
operator * ( a as complex, b as complex ) as complex
|
|
dim as complex c
|
|
c.r = a.r*b.r - a.i*b.i
|
|
c.i = a.i*b.r + a.r*b.i
|
|
return c
|
|
end operator
|
|
|
|
operator / ( a as complex, b as complex ) as complex
|
|
dim as double bcb = (b*conj(b)).r
|
|
dim as complex acb = a*conj(b), c
|
|
c.r = acb.r/bcb
|
|
c.i = acb.i/bcb
|
|
return c
|
|
end operator
|
|
|
|
sub printc( a as complex )
|
|
if a.i>=0 then
|
|
print using "############.############### + ############.############### i"; a.r; a.i
|
|
else
|
|
print using "############.############### - ############.############### i"; a.r; -a.i
|
|
end if
|
|
end sub
|
|
|
|
function intc( n as integer ) as complex
|
|
dim as complex c
|
|
c.r = n
|
|
c.i = 0.0
|
|
return c
|
|
end function
|
|
|
|
function absc( a as complex ) as double
|
|
return sqr( (a*conj(a)).r )
|
|
end function
|
|
|
|
'-----------------------
|
|
' the algorithm
|
|
' Uses a rapidly converging continued
|
|
' fraction expansion for e^z and recursive
|
|
' expressions for its convergents
|
|
'-----------------------
|
|
|
|
dim as complex pii, pii2, curr, A2, A1, A0, B2, B1, B0
|
|
dim as complex ONE, TWO
|
|
dim as integer i, k = 2
|
|
pii.r = 0.0
|
|
pii.i = PI
|
|
pii2 = pii*pii
|
|
|
|
B0 = intc(2)
|
|
A0 = intc(2)
|
|
B1 = (intc(2) - pii)
|
|
A1 = B0*B1 + intc(2)*pii
|
|
printc( A0/B0)
|
|
print " Absolute error = ", absc(A0/B0)
|
|
printc( A1/B1)
|
|
print " Absolute error = ", absc(A1/B1)
|
|
|
|
for i = 1 to MAXITER
|
|
k = k + 4
|
|
A2 = intc(k)*A1 + pii2*A0
|
|
B2 = intc(k)*B1 + pii2*B0
|
|
curr = A2/B2
|
|
A0 = A1
|
|
A1 = A2
|
|
B0 = B1
|
|
B1 = B2
|
|
printc( curr )
|
|
print " Absolute error = ", absc(curr)
|
|
next i
|