'******************************************************************************
'Subject: Graphically exploring the Euler gamma- and Riemann zeta-functions.
'Author : Sjoerd.J.Schaper - vspickelen
'Date : 09-21-2008
'Code : QBasic, PDS, FreeBasic (-lang fblite)
'******************************************************************************
'This program is copyright (c) 2008 by the author. It is made available as is,
'and no warranty - about the program, its performance, or its conformity to any
'specification - is given or implied. It may be used, modified, and distributed
'freely, as long as the original author is credited as part of the final work.
'******************************************************************************
DEFINT A-Z
TYPE Cmpx
r AS DOUBLE
i AS DOUBLE
END TYPE
DECLARE SUB cadd (a AS Cmpx, b AS Cmpx)
DECLARE SUB csub (a AS Cmpx, b AS Cmpx)
DECLARE SUB cmul (a AS Cmpx, b AS Cmpx)
DECLARE SUB cdiv (a AS Cmpx, b AS Cmpx)
'complex arithmetic and functions
DECLARE SUB cexp (a AS Cmpx)
DECLARE SUB clog (a AS Cmpx)
DECLARE SUB c2sn (a AS Cmpx) ' 2*sin a
DECLARE SUB cpow (c AS Cmpx, z AS Cmpx)
'a + bi:= a ^ -(x + yi)
DECLARE SUB scadd (c AS Cmpx, z AS Cmpx, sc AS DOUBLE, k)
'c:= c + (-1)^k-1 * z * scalar
DECLARE FUNCTION seca% (a AS Cmpx, b AS Cmpx, d AS SINGLE)
'secant approximate root in [a, b]
DECLARE FUNCTION norm# (a AS Cmpx)
'return norm(a)
DECLARE FUNCTION cs# (sn AS DOUBLE, a AS DOUBLE)
'return cos(a), sn:= sin(a)
DECLARE FUNCTION magn! (a AS DOUBLE, b AS SINGLE)
'return log_b(a) (mod 1)
DECLARE FUNCTION rond$ (a AS DOUBLE, k)
'round(a > 0) to k decimal digits
DECLARE FUNCTION mkc$ (a AS Cmpx)
'convert argument to string
DECLARE SUB prntf (s AS Cmpx, t AS Cmpx)
DECLARE SUB bino (n AS INTEGER)
'binomial coefficients (n, k)
DECLARE SUB lgamma (t AS Cmpx, s AS Cmpx)
'let t = log(gamma(s)), s.r > 0
DECLARE SUB zetap (t AS Cmpx, s AS Cmpx)
'let t = zeta(s), s.r > -1, s <> 1
DECLARE SUB gamma (t AS Cmpx, s AS Cmpx)
'meromorphic gamma-function
DECLARE SUB zeta (t AS Cmpx, s AS Cmpx)
'meromorphic Riemann-zeta
DECLARE SUB gfxin (a AS SINGLE, b AS SINGLE, c AS SINGLE, d AS SINGLE)
DECLARE SUB plot (x AS DOUBLE, y AS DOUBLE, t)
DECLARE SUB eval (a AS DOUBLE, b AS DOUBLE, sw)
DECLARE SUB real (a, b, h AS SINGLE, p AS LONG, sw)
DECLARE SUB tcir (c AS SINGLE, r AS SINGLE, h AS SINGLE, p AS LONG, sw)
DECLARE SUB zlin (a, b, h AS SINGLE, p AS LONG)
DECLARE SUB zmod (a, b, h AS SINGLE, p AS LONG)
DECLARE SUB zflo (c AS SINGLE, d AS SINGLE, h AS SINGLE, a, b)
DECLARE SUB tpot (c AS SINGLE, d AS SINGLE, h AS SINGLE, sw)
'comments in CASE-block below
DATA 12,360,1260,1680,1188,522
CONST Lx = 1023, Gc = 6 ' gamma-coefficients
CONST Eps = 1D-18, Inf = 10000000#
CONST MODE = 12, dW = 640, dH = 480 ' screenmode, width, height
'CONST MODE = 19, dW = 800, dH = 600 ' FreeBasic only
DIM SHARED pi AS DOUBLE: pi = 4 * ATN(1)
DIM SHARED tpi AS DOUBLE: tpi = 2 * pi ' global constants
DIM SHARED L2pi AS DOUBLE: L2pi = LOG(tpi)
DIM SHARED L10 AS DOUBLE: L10 = 1 / LOG(10)
DIM SHARED bc(Lx) AS DOUBLE, m AS INTEGER
DIM c AS SINGLE, d AS SINGLE, h AS SINGLE
DIM f AS STRING, g(5) AS STRING, p AS LONG
g(0) = " zeta": g(1) = " gamma": g(2) = ", [Esc] to break"
g(3) = " eval, real, circ, tlin, modu, flow, port"
g(4) = " window height h, plotting scale p"
g(5) = " even(sw) to evaluate zeta, odd(sw) for gamma"
DO
SCREEN 0: CLS
LOCATE 2, 1: PRINT g(3)
DO: LOCATE 3, 1: INPUT " function"; f
LOOP WHILE INSTR(f, "'")
IF f = "0" OR f = "" THEN EXIT DO
PRINT : PRINT " argument s, value t"
SELECT CASE (INSTR(g(3), f) - 2) \ 6
CASE 1
PRINT " plot(re(s), re(t)) with re(s) in [a, b], im = 0"
PRINT g(4): PRINT g(5)
INPUT " a, b, h, p, sw"; a, b, h, p, sw
PRINT : PRINT " real valued"; g(sw AND 1)
real a, b, h, p, sw
CASE 2
PRINT " plot(t) with s on circle(c, 0), radius r"
PRINT g(4): PRINT g(5)
INPUT " c, r, h, p, sw"; c, d, h, p, sw
PRINT : PRINT g(sw AND 1); " transformation of a circle"
tcir c, d, h, p, sw
CASE 3
PRINT " plot(t) with s on line re = 1/2, im(s) in [a, b]"
PRINT g(4): INPUT " a, b, h, p"; a, b, h, p
PRINT : PRINT " transformation of the critical line"
zlin a, b, h, p
CASE 4
PRINT " plot(im(s), abs(t)) with re(s) = 1/2, im(s) in [a, b] <4200"
PRINT g(4): PRINT " zeroes are printed in file 'zetazero.log'"
PRINT " increase p if zeroes don't separate properly"
INPUT " a, b, h, p"; a, b, h, p
PRINT : PRINT " emulating Hardy's function Z(t)"
zmod a, b, h, p
CASE 5
PRINT " plot solutions to ds/dt = zeta(s)"
PRINT " center(c, d), height h, meshgrid (2a+1)x(2b+1)"
INPUT " c, d, h, a, b"; c, d, h, a, b
PRINT : PRINT " phase portrait of the zeta flow"; g(2)
zflo c, d, h, a, b
CASE 6
PRINT " dual color portrait, center(c, d), height h": PRINT g(5)
PRINT " sw (mod 4) <= 1: left arg(t) mod 2pi, right abs(t) mod 1,"
PRINT " else signed magnitude (mod 1): left im(t), right re(t) reflected"
PRINT " set sw>= 4 to map the secant transformation associated with zeta"
INPUT " c, d, h, sw"; c, d, h, sw
IF (sw AND 5) = 4 THEN
PRINT : PRINT " zeta secant iteration"; g(2)
ELSE
PRINT : PRINT g(sw AND 1); " range portrait"; g(2)
END IF
tpot c, d, h, sw
CASE ELSE
PRINT g(5): INPUT " re(s), im(s), sw"; r#, i#, sw
PRINT : PRINT " t ="; g(sw AND 1); "(s)"
eval r#, i#, sw
END SELECT
SLEEP
LOOP
SYSTEM
SUB bino (n AS INTEGER)
bc(0) = 1: m = n
IF m > Lx THEN m = Lx
IF m < 16 THEN m = 16
FOR k = 0 TO m \ 2
bc(m - k) = bc(k)
bc(k + 1) = (m - k) / (k + 1) * bc(k)
NEXT k
END SUB
SUB c2sn (a AS Cmpx)
DIM x AS DOUBLE, y AS DOUBLE, sn AS DOUBLE
x = EXP(a.i): y = 1 / x
a.i = (x - y) * cs(sn, a.r)
a.r = (x + y) * sn
END SUB
SUB cadd (a AS Cmpx, b AS Cmpx)
a.r = a.r + b.r: a.i = a.i + b.i
END SUB
SUB cdiv (a AS Cmpx, b AS Cmpx)
DIM tmp AS DOUBLE, d AS DOUBLE, q AS DOUBLE
IF ABS(b.r) < ABS(b.i) THEN ' ACM algorithm 116
q = b.r / b.i
d = b.r * q + b.i
tmp = (a.r * q + a.i) / d
a.i = (a.i * q - a.r) / d
ELSE
q = b.i / b.r
d = b.r + b.i * q
tmp = (a.r + a.i * q) / d
a.i = (a.i - a.r * q) / d
END IF
a.r = tmp
END SUB
SUB cexp (a AS Cmpx)
DIM r AS DOUBLE, sn AS DOUBLE
r = EXP(a.r)
a.r = r * cs(sn, a.i)
a.i = r * sn
END SUB
SUB clog (a AS Cmpx)
DIM fi AS DOUBLE
fi = ATN(a.i / a.r)
IF a.r < 0 THEN fi = fi + SGN(a.i) * pi
a.r = LOG(norm(a)) * .5: a.i = fi
END SUB
SUB cmul (a AS Cmpx, b AS Cmpx)
DIM tmp AS DOUBLE
tmp = a.r * b.r - a.i * b.i
a.i = a.i * b.r + a.r * b.i
a.r = tmp
END SUB
SUB cpow (a AS Cmpx, b AS Cmpx)
DIM r AS DOUBLE, sn AS DOUBLE
a.r = -LOG(a.r): r = EXP(a.r * b.r)
a.r = r * cs(sn, a.r * b.i)
a.i = r * sn
END SUB
FUNCTION cs# (s AS DOUBLE, a AS DOUBLE)
a = a / tpi: a = (a - INT(a)) * tpi
s = COS(a): cs = s
s = SQR(1 - s * s)
IF a > pi THEN s = -s
END FUNCTION
SUB csub (a AS Cmpx, b AS Cmpx)
a.r = a.r - b.r: a.i = a.i - b.i
END SUB
SUB eval (a AS DOUBLE, b AS DOUBLE, sw)
DIM dx AS SINGLE, dy AS SINGLE
DIM s AS Cmpx, t AS Cmpx
s.r = a: s.i = b
IF sw AND 1 THEN
gamma t, s
ELSE
IF ABS(b) >= m THEN bino CINT(ABS(b))
zeta t, s
END IF
'
dx = ABS(t.r): dy = ABS(t.i)
IF dx < Inf THEN
IF ABS(s.r) > dx THEN dx = ABS(s.r)
IF ABS(s.i) > dy THEN dy = ABS(s.i)
dx = dx * 1.1: dy = dy * 1.1
IF dx * dH / dW > dy THEN
dy = dH / dW * dx
ELSE
dx = dW / dH * dy
END IF
gfxin -dx, -dy, dx, dy
LINE (0, 0)-(s.r, s.i), 12
CIRCLE (s.r, s.i), dx * .007, 12
LINE (0, 0)-(t.r, t.i), 14
CIRCLE (t.r, t.i), dx * .007, 14
ELSE
SCREEN MODE
END IF
prntf s, t
END SUB
SUB gamma (t AS Cmpx, z AS Cmpx)
DIM x AS Cmpx
IF z.r > 0 THEN
lgamma t, z: cexp t
ELSE
'
IF (z.r - INT(z.r)) < Eps THEN
IF ABS(z.i) < Eps THEN ' integer poles
t.r = Inf: t.i = z.i
IF CINT(z.r) AND 1 THEN t.r = -Inf
EXIT SUB
END IF
END IF
z.r = 1 - z.r: z.i = -z.i ' reflection formula 6.1.17
lgamma t, z: cexp t ' t = gamma(1-z)
z.r = 1 - z.r: z.i = -z.i
x.r = z.r * pi
x.i = z.i * pi: c2sn x ' x = 2sin(pi*z)
cmul x, t: t.r = tpi: t.i = 0
cdiv t, x ' 2pi / t*x
END IF
END SUB
SUB gfxin (a AS SINGLE, b AS SINGLE, c AS SINGLE, d AS SINGLE)
DIM mx AS SINGLE, my AS SINGLE, dx AS SINGLE, dy AS SINGLE
mx = (a + c) * .5: my = (b + d) * .5
dx = (c - a) * .51: IF dx = 0 THEN dx = 1.02 * dW / dH
dy = (d - b) * .51: IF dy = 0 THEN dy = 1.02
SCREEN MODE: WINDOW (mx - dx, my - dy)-(mx + dx, my + dy)
LINE (mx - dx, 0)-(mx + dx, 0)
LINE (0, my - dy)-(0, my + dy)
FOR t = -1 TO 0
k = 1: j = INT(mx - dx)
IF dx > 20 THEN k = 10: j = (j \ k - 1) * k
FOR i = j TO 1 + INT(mx + dx) STEP k
IF t THEN
LINE (i, -.02 * dy)-(i, .02 * dy)
ELSE
LINE (-.02 * dx, i)-(.02 * dx, i)
END IF
NEXT i: mx = my: SWAP dx, dy
NEXT t
END SUB
SUB lgamma (t AS Cmpx, z AS Cmpx)
'Abramowitz & Stegun, Handbook of Mathematical Functions, 6.1.40
DIM x AS Cmpx, y AS Cmpx
sw = z.r < 8
IF sw THEN ' compute lgamma(z + n)
n = 8 - INT(z.r): z.r = z.r + n
END IF
'
t = z: x = z
t.r = t.r - .5: clog x
cmul t, x: csub t, z ' (z -.5)*log(z) - z
t.r = t.r + .5 * L2pi ' + log(sqrt(2pi))
'
x = z: y.r = 1: y.i = 0
cdiv y, z: cmul y, y ' y = z^-2
FOR k = 1 TO Gc ' Stirling's series
cmul x, y: READ d
scadd t, x, 1 / d, k
NEXT k: RESTORE
'
IF sw THEN
z.r = z.r - n: x = z: y = z
FOR k = 2 TO n
y.r = y.r + 1: cmul x, y ' rising factorial 6.1.22
NEXT k: clog x
csub t, x ' remove (z)n
END IF
END SUB
FUNCTION magn! (a AS DOUBLE, b AS SINGLE)
s = SGN(a): a = b * LOG(ABS(a) + Eps)
a = a - INT(a): IF s < 0 THEN a = a - 1
magn = a
END FUNCTION
FUNCTION mkc$ (a AS Cmpx)
DIM x AS DOUBLE, g AS STRING, s AS STRING
s = " ": IF a.r < 0 THEN s = "-"
g = s + rond$(ABS(a.r), 12)
IF ABS(a.r) = Inf THEN g = s + "#INF"
x = ABS(a.i)
IF x > 0 THEN
s = "+": IF a.i < 0 THEN s = "-"
IF ABS(x - 1) < Eps THEN
g = g + s + "i"
ELSE
g = g + s + rond$(x, 12) + "*i"
END IF
END IF
mkc = g
END FUNCTION
FUNCTION norm# (a AS Cmpx)
norm = a.r * a.r + a.i * a.i
END FUNCTION
SUB plot (x AS DOUBLE, y AS DOUBLE, t)
IF t THEN
LINE -(x, y)
ELSE
PSET (x, y): t = -1
END IF
END SUB
SUB prntf (a AS Cmpx, b AS Cmpx)
DIM g AS STRING * 80
MID$(g, 1) = "s=" + mkc(a)
MID$(g, 41) = "t=" + mkc(b)
LOCATE 1, 1: PRINT g;
END SUB
SUB real (a, b, h AS SINGLE, s AS LONG, sw)
DIM x AS Cmpx, y AS Cmpx, d AS SINGLE, i AS LONG
gfxin CSNG(a), -h, CSNG(b), h
bino 16: IF s = 0 THEN s = 15
FOR i = a * s TO b * s
x.r = i / s: x.i = 0 ' real axis
IF sw AND 1 THEN
gamma y, x
ELSE
zeta y, x
END IF
prntf x, y
IF ABS(y.r - d) > h + h THEN t = 0
d = y.r: plot x.r, y.r, t
NEXT i
END SUB
FUNCTION rond$ (a AS DOUBLE, k)
DIM b AS DOUBLE, e AS DOUBLE
IF a > 0 THEN
e = 1 + INT(LOG(a) * L10)
b = 10 ^ (k - e)
b = INT(a * b + .5) / b
END IF
rond = LTRIM$(STR$(b))
END FUNCTION
SUB scadd (a AS Cmpx, b AS Cmpx, s AS DOUBLE, k)
IF k AND 1 THEN
a.r = a.r + s * b.r
a.i = a.i + s * b.i
ELSE
a.r = a.r - s * b.r
a.i = a.i - s * b.i
END IF
END SUB
FUNCTION seca% (x AS Cmpx, y AS Cmpx, d AS SINGLE)
DIM f AS Cmpx, g AS Cmpx, z AS Cmpx
DIM p AS Cmpx, n AS SINGLE
IF d = 0 THEN
sw = 0: tx = 16 ' zmod,
ELSE
sw = 1: tx = 4: p.r = 1 ' tpot call
x.r = y.r - d * .1: x.i = y.i
END IF
'
zeta g, x
FOR t = 1 TO tx
f = g: zeta g, y
csub f, g: csub x, y
cdiv x, f: cmul x, g ' g*(x-y)/(f-g)
n = norm(x)
IF n < Eps THEN EXIT FOR
IF n > 225 THEN sw = 1: EXIT FOR
IF sw THEN cmul p, x ' trace(dx)
z = y: csub y, x: x = z ' update y, x
NEXT t: x = p
IF sw = 0 THEN t = t > 16
seca = t
END FUNCTION
SUB tcir (c AS SINGLE, r AS SINGLE, h AS SINGLE, s AS LONG, sw)
DIM a AS DOUBLE, sn AS DOUBLE, d AS SINGLE
DIM x AS Cmpx, y AS Cmpx, i AS LONG
bino CINT(r): a = dW / dH * h
gfxin .5 - a, -h, .5 + a, h
IF s = 0 THEN s = 180
FOR i = 0 TO 2 * s
a = pi * i / s
x.r = c + r * cs(sn, a) ' circle (c, 0), r
x.i = r * sn
IF sw AND 1 THEN
gamma y, x
ELSE
zeta y, x
END IF
prntf x, y
IF ABS(y.i - d) > h + h THEN t = 0
d = y.i: plot y.r, y.i, t
NEXT i
END SUB
SUB tpot (a AS SINGLE, b AS SINGLE, dy AS SINGLE, sw)
DIM sx AS SINGLE, sy AS SINGLE, dx AS SINGLE, h(255) AS LONG
DIM s AS Cmpx, t AS Cmpx, f AS SINGLE, g AS SINGLE
h(0) = &H40406: h(1) = &HC000B ' sepia, violet
h(2) = &H80010: h(3) = &H17 ' mauve, carmine
h(4) = &H61E: h(5) = &HC22 ' vermilion, orange
h(6) = &H1322: h(7) = &H1A23 ' deep yellow, lemon yellow
h(8) = &HE0006: h(9) = &H110500 ' blue violet, ultramarine
h(10) = &H100A00: h(11) = &HE0F00 ' cobalt blue, blue green
h(12) = &HB1400: h(13) = &H180D ' turquoise, light green
h(14) = &H181D: h(15) = &H3D3E3F ' greenish yellow, grey
SCREEN MODE: PALETTE USING h(0)
IF dy = 0 THEN dy = 3
bino ((ABS(b) + dy) * .5)
dx = dW / dH * dy * .5: f = .5 / LOG(15) ' gain 15^2
sx = dW * .25: dx = dx / sx
sy = dH * .5: dy = dy / sy
'
FOR p = -sx TO sx
FOR q = -sy TO sy
s.r = a + p * dx
s.i = b + q * dy
IF sw AND 1 THEN
g = f: gamma t, s
ELSEIF sw AND 4 THEN
g = f * 4 / seca(t, s, dx)
ELSE
g = f: zeta t, s
END IF
'
IF sw AND 2 THEN ' signed magnitude mod 1
kl = INT(7 * (1 + magn(t.i, g)))
kr = INT(7 * (1 + magn(t.r, g)))
ELSE
clog t: t.r = t.r * g
kl = INT(7 * (1 + t.i / pi)) ' argument mod 2pi
kr = INT(14 * (t.r - INT(t.r))) ' log(modulus) mod 1
END IF
PSET (sx + p, sy - q), 1 + kl
PSET (3 * sx - p, sy - q), 1 + kr ' reflect
NEXT q
IF INKEY$ = CHR$(27) THEN EXIT FOR
NEXT p
'
p = a / dx: q = b / dy
LINE (3 * sx + p, 0)-(3 * sx + p, 2 * sy), 0
LINE (sx - p, 0)-(sx - p, 2 * sy), 0
LINE (2 * sx, 0)-(2 * sx, 2 * sy), 0
END SUB
SUB zeta (t AS Cmpx, s AS Cmpx)
DIM x AS Cmpx
IF s.r > -1 THEN
IF ABS(s.r - 1) < Eps THEN
IF ABS(s.i) < Eps THEN ' harmonic pole
t.r = Inf: t.i = s.i: EXIT SUB
END IF
END IF
zetap t, s
ELSE
'
s.r = 1 - s.r: s.i = -s.i ' functional equation 23.2.6
lgamma x, s ' x = log(gamma(1-s))
scadd x, s, L2pi, 0: cexp x ' - (1-s)*log(2pi)
zetap t, s: cmul t, x ' t = exp(x) * zeta(1-s)
s.r = 1 - s.r: s.i = -s.i
x.r = s.r * .5 * pi
x.i = s.i * .5 * pi: c2sn x ' * 2sin(s*pi/2)
cmul t, x
END IF
END SUB
SUB zetap (t AS Cmpx, s AS Cmpx)
'Borwein's Algorithm 3 based on poly x^n * (1 - x)^n, in:
'An efficient algorithm for the Riemann zeta function' (1995)
DIM p2 AS DOUBLE, c AS DOUBLE
DIM p AS Cmpx, z AS Cmpx
t.r = 0: t.i = 0
FOR k = 1 TO m ' split alternating sums
p.r = k: cpow p, s ' (-1)^k-1 / k^s
IF k AND 1 THEN
cadd t, p
ELSE
csub t, p
END IF
NEXT k
'
z.r = 0: z.i = 0
p2 = 2 ^ m: c = p2 - 1
FOR k = 1 + m TO m + m ' scale with combi sum
p.r = k: cpow p, s
scadd z, p, c, k
c = c - bc(k - m)
NEXT k
scadd t, z, 1 / p2, 1 ' t + z / 2^m
'
p.r = .5: s.r = 1 - s.r
cpow p, s: s.r = 1 - s.r
p.r = 1 - p.r: cdiv t, p ' eta to zeta
END SUB
SUB zflo (a AS SINGLE, b AS SINGLE, dy AS SINGLE, sx, sy)
DIM s AS Cmpx, x AS Cmpx, z AS Cmpx, dx AS SINGLE
DIM q AS DOUBLE, h AS DOUBLE, d AS SINGLE, r AS SINGLE
IF dy = 0 THEN dy = 3
bino ABS(b) + dy: dx = dW / dH * dy
gfxin a - dx, b - dy, a + dx, b + dy
r = dx * dx + dy * dy: d = SQR(r) * .01
sx = ABS(sx): dx = dx / (1 + sx)
sy = ABS(sy): dy = dy / (1 + sy)
'
FOR i = -sx TO sx
FOR j = -sy TO sy
FOR k = 0 TO 1
s.r = a + i * dx
s.i = b + j * dy
h = Inf: t = 0
FOR p = 1 TO 256 ' solve ds/dt = zeta(s)
IF ABS(s.r) = Inf THEN t = 0
plot s.r, s.i, t
'
zeta x, s: q = SQR(norm(x))
IF 1.8 * h * q < d THEN EXIT FOR ' critical point
z = s: h = d / q
scadd z, x, .5 * h, k ' midpoint method 25.5.7
zeta x, z: scadd s, x, h, k
'
z.r = s.r - a: z.i = s.i - b
IF norm(z) > r THEN EXIT FOR
NEXT p
IF INKEY$ = CHR$(27) THEN EXIT SUB
NEXT k
NEXT j
NEXT i
END SUB
SUB zlin (a, b, h AS SINGLE, s AS LONG)
DIM x AS Cmpx, y AS Cmpx, i AS LONG
n = ABS(b): IF ABS(a) > n THEN n = ABS(a)
bino n: x.r = dW / dH * h
gfxin .5 - x.r, -h, .5 + x.r, h
IF s = 0 THEN s = 15
FOR i = a * s TO b * s
x.r = .5: x.i = i / s ' critical line
zeta y, x: prntf x, y
plot y.r, y.i, t
NEXT i
END SUB
SUB zmod (a, b, h AS SINGLE, s AS LONG)
DIM x AS Cmpx, y AS Cmpx, g AS STRING, i AS LONG
DIM p AS SINGLE, q AS SINGLE, t0 AS SINGLE
OPEN "zetazero.log" FOR OUTPUT AS #1
n = ABS(b): IF ABS(a) > n THEN n = ABS(a)
bino n: ln = dH \ 16 - 1: j = 1: sg = -1
gfxin CSNG(a), -h, CSNG(b), h
t0 = TIMER: IF s = 0 THEN s = 10
'
FOR i = a * s TO b * s
x.r = .5: x.i = i / s
zeta y, x: prntf x, y
q = p: p = SQR(norm(y)) ' modulus
'
IF fl AND (p >= q) THEN ' zero crossed
IF ABS(x.i) > 11 THEN
y.r = .5: y.i = (i - 2) / s
g = "r" + LTRIM$(STR$(j))
IF seca(y, x, 0) THEN g = "fail"
g = g + mkc(x) + " "
PRINT #1, g: LOCATE ln, 1: PRINT g;
j = j + 1: sg = -sg: fl = 0
END IF
ELSEIF (p < q) AND NOT fl THEN
fl = -1
END IF
plot i / s, sg * p, t
NEXT i
g = "Timer " + STR$(CSNG(TIMER - t0)) + " s"
PRINT #1, g: LOCATE ln, 41: PRINT g;
CLOSE 1
END SUB