09-23-2025, 06:53 AM
here's my translation of pi1.c from https://bellard.org/pi/pi1.c
it prints 9 digits of Pi starting at position n
it prints 9 digits of Pi starting at position n
Code: (Select All)
_Title "pi1_c"
$Console:Only
_Dest _Console
Option _Explicit
Dim As Double t
Dim As Long i, n
t = Timer(.0001)
n = -8
pi1 1
pi1 10
pi1 32
pi1 100 - 9
pi1 1000 - 9
Print
Print " Pi = 3."
For i = 1 To 111
n = n + 9
pi1 n
Next i
t = Timer(.0001) - t
Print "runtime in seconds = "; t
'Print "press return to end"
'Sleep
Function inv_mod& (x As Long, y As Long)
Dim As Long q, u, v, a, c, t
u = x
v = y
c = 1
a = 0
Do
q = v \ u
t = c
c = a - (q * c)
a = t
t = u
u = v - (q * u)
v = t
Loop While u <> 0
a = a Mod y
If a < 0 Then
a = y + a
End If
inv_mod = a
End Function
Function inv_mod2& (u As Long, v As Long)
Dim As Long u1, u3, v1, v3, t1, t3
u1 = 1
u3 = u
v1 = v
v3 = v
If (u And 1) <> 0 Then
t1 = 0
t3 = -v
GoTo Y4
Else
t1 = 1
t3 = u
End If
Do
Do
If (t1 And 1) = 0 Then
t1 = _ShR(t1, 1)
t3 = _ShR(t3, 1)
Else
t1 = _ShR((t1 + v), 1)
t3 = _ShR(t3, 1)
End If
Y4:
Loop While (t3 And 1) = 0
If t3 >= 0 Then
u1 = t1
u3 = t3
Else
v1 = v - t1
v3 = -t3
End If
t1 = u1 - v1
t3 = u3 - v3
If t1 < 0 Then
t1 = t1 + v
End If
Loop While t3 <> 0
inv_mod2 = u1
End Function
Function pow_mod& (a_in As Long, b_in As Long, m_in As Long)
Dim As _Integer64 a, b, m, r, aa
a = a_in
b = b_in
m = m_in
r = 1
aa = a
While 1
If b And 1 Then
r = (r * aa) Mod m
End If
b = _ShR(b, 1)
If b = 0 Then Exit While
aa = (aa * aa) Mod m
Wend
pow_mod = r
End Function
Function is_prime& (n As Long)
Dim As Long r, i
If (n Mod 2) = 0 Then
is_prime = 0
Exit Function
End If
r = Int(Sqr(CDbl(n)))
i = 3
While i <= r
If (n Mod i) = 0 Then
is_prime = 0
Exit Function
End If
i = i + 2
Wend
is_prime = 1
End Function
Function next_prime& (n As Long)
Do
n = n + 1
Loop While is_prime(n) = 0
next_prime = n
End Function
Sub fmod (f_mod As Double, num As Double, denom As Double)
Dim As Double q, fm 'make copy in case the destination and source are the same
fm = num / denom
q = Fix(fm)
fm = fm - q
fm = fm * denom
f_mod = fm
End Sub
Sub pi1 (n As Long)
Dim As _Integer64 av, a, vmax, N2, num, den, k
Dim As _Integer64 kq1, kq2, kq3, kq4, t, v
Dim As _Integer64 s, i, t1
Dim As Double sum
Dim As String sn
N2 = Int(((n + 20#) * Log(10#)) / Log(13.5#))
sum = 0
a = 2
While a <= (3 * N2)
vmax = Int(Log(3# * N2) / Log(CDbl(a)))
If a = 2 Then
vmax = vmax + (N2 - n)
If (vmax <= 0) Then
GoTo Continue_While
End If
End If
av = 1
i = 0
While i < vmax
av = av * a
i = i + 1
Wend
s = 0
den = 1
kq1 = 0
kq2 = -1
kq3 = -3
kq4 = -2
If a = 2 Then
num = 1
v = -n
Else
num = pow_mod(2, n, av)
v = 0
End If
k = 1
While k <= N2
t = 2 * k
kq1 = kq1 + 2
If kq1 >= a Then
Do
kq1 = kq1 - a
Loop While kq1 >= a
If kq1 = 0 Then
Do
t = t \ a
v = v - 1
Loop While (t Mod a) = 0
End If
End If
num = (num * t) Mod av
t = (2 * k) - 1
kq2 = kq2 + 2
If kq2 >= a Then
Do
kq2 = kq2 - a
Loop While kq2 >= a
If kq2 = 0 Then
Do
t = t \ a
v = v - 1
Loop While (t Mod a) = 0
End If
End If
num = (num * t) Mod av
t = 3 * ((3 * k) - 1)
kq3 = kq3 + 9
If kq3 >= a Then
Do
kq3 = kq3 - a
Loop While kq3 >= a
If kq3 = 0 Then
Do
t = t \ a
v = v + 1
Loop While (t Mod a) = 0
End If
End If
den = (den * t) Mod av
t = (3 * k) - 2
kq4 = kq4 + 3
If kq4 >= a Then
Do
kq4 = kq4 - a
Loop While kq4 >= a
If kq4 = 0 Then
Do
t = t \ a
v = v + 1
Loop While (t Mod a) = 0
End If
End If
If a <> 2 Then
t = t * 2
Else
v = v + 1
End If
den = (den * t) Mod av
If v > 0 Then
If a <> 2 Then
t = inv_mod2(den, av)
Else
t = inv_mod(den, av)
End If
t = (t * num) Mod av
i = v
While i < vmax
t = (t * a) Mod av
i = i + 1
Wend
t1 = (25 * k) - 3
t = (t * t1) Mod av
s = s + t
If s >= av Then
s = s - av
End If
End If
k = k + 1
Wend
t = pow_mod(5, n - 1, av)
s = (s * t) Mod av
fmod sum, sum + (CDbl(s) / CDbl(av)), 1#
Continue_While:
a = next_prime(a)
Wend
sn = _Trim$(Str$(Int(sum * 1000000000)))
While Len(sn) < 9
sn = "0" + sn
Wend
Print Using "Decimal digits of pi at position ########"; n;: Print ": "; sn
End Sub
