Plot Riemann Zeta Function ζ(1/2 +ti) for PicoCalc / MMBasic

zeta.bas

CLS
Const LOWER_THRESHOLD = 1E-6
Const UPPER_BOUND = 1E4
Const INFINITY = 50
Const X_CENTER = 100
Const Y_CENTER = 150
Const X_SCALE = 50
Const Y_SCALE = 50

'------------------------------------------------
' Complex
'------------------------------------------------
Sub ComplexAdd(ar, ai, br, bi, BYREF rr, BYREF ri)
    rr = ar + br
    ri = ai + bi
End Sub

Sub ComplexDiv(ar, ai, br, bi, BYREF rr, BYREF ri)
    Local denom As FLOAT
    denom = br*br + bi*bi
    rr = (ar*br + ai*bi) / denom
    ri = (ai*br - ar*bi) / denom
End Sub

Sub ComplexPowNeg(j As FLOAT, sr As FLOAT, si As FLOAT, BYREF rr, BYREF ri)
    Local lr As FLOAT, a As FLOAT, b As FLOAT, rmag As FLOAT
    lr = Log(j)
    a = -sr * lr
    b = -si * lr
    rmag = Exp(a)
    rr = rmag * Cos(b)
    ri = rmag * Sin(b)
End Sub

'------------------------------------------------
' nCk
'------------------------------------------------
Function Binom(n, k)
    Local res As FLOAT, i As INTEGER
    If k > n - k Then k = n - k
    If k = 0 Then Binom = 1: Exit Function
    res = 1
    For i = 1 To k
        res = res * (n - k + i) / i
    Next i
    Binom = res
End Function

'------------------------------------------------
' Zeta
' s = sr + i*si
'------------------------------------------------
Sub Zeta(sr, si, BYREF zr, BYREF zi)
    Local m As INTEGER, j As INTEGER
    Local inner_r As FLOAT, inner_i As FLOAT
    Local pr As FLOAT, P_i As FLOAT
    Local outer_r As FLOAT, outer_i As FLOAT
    Local prev_r As FLOAT, prev_i As FLOAT
    Local denom_r As FLOAT, denom_i As FLOAT

    outer_r = 0: outer_i = 0
    prev_r = 1E9: prev_i = 1E9

    Local pow2_r As FLOAT, pow2_i As FLOAT
    Local a As FLOAT, b As FLOAT, mag As FLOAT
    a = (1 - sr) * Log(2)
    b = -si * Log(2)
    mag = Exp(a)
    pow2_r = mag * Cos(b)
    pow2_i = mag * Sin(b)
    denom_r = 1 - pow2_r
    denom_i = -pow2_i
    Local c1 As FLOAT, c2 As FLOAT

    For m = 1 To INFINITY
        inner_r = 0: inner_i = 0
        For j = 1 To m
            If ((j - 1) Mod 2) = 0 Then c1 = 1 Else c1 = -1
            c2 = Binom(m - 1, j - 1)
            ComplexPowNeg(j, sr, si, pr, P_i)
            pr = pr * c1 * c2
            P_i = P_i * c1 * c2
            ComplexAdd(inner_r, inner_i, pr, P_i, inner_r, inner_i)
        Next j
        inner_r = inner_r / (2^m)
        inner_i = inner_i / (2^m)
        ComplexDiv(inner_r, inner_i, denom_r, denom_i, pr, P_i)
        ComplexAdd(outer_r, outer_i, pr, P_i, outer_r, outer_i)

        If Abs(prev_r - pr) + Abs(prev_i - P_i) < LOWER_THRESHOLD Then Exit For
        If Abs(outer_r) + Abs(outer_i) > UPPER_BOUND Then Exit For
        prev_r = pr: prev_i = P_i
    Next m

    zr = outer_r
    zi = outer_i
End Sub

'------------------------------------------------
Dim t As INTEGER
Dim s_r As FLOAT, s_i As FLOAT
Dim z_r As FLOAT, z_i As FLOAT
Dim x As INTEGER, y As INTEGER
Dim px As integer, py As integer

For t = 0 To 500
    s_r = 0.5
    s_i = t / 10.0
    Zeta(s_r, s_i, z_r, z_i)
    x = X_CENTER + z_r * X_SCALE
    y = Y_CENTER - z_i * Y_SCALE
    If t > 0 Then
      If t < 256 Then
        c = RGB(255-t,t,0)
      Else
        c = RGB(0,511-t,t-256)
      End If
      Line px, py, x, y, 1, c
    End If
    px = x
    py = y
Next t

Save IMAGE "zeta"