This site is developed to XHTML and CSS2 W3C standards.
If you see this paragraph, your browser does not support those standards and you
need to upgrade. Visit WaSP
for a variety of options.
Paste #96
Posted by: Pi
Posted on: 2026-11-03 08:20:40
Age: 4 days ago
Views: 6
import sys
import time
import numpy as np
from math import isqrt, log2, log10, ceil
SH = 96
def primes(m):
s = np.ones(m + 1, dtype=bool)
s[:2] = False
for i in range(2, isqrt(m) + 1):
if s[i]:
s[i * i::i] = False
return np.nonzero(s)[0].astype(np.int64)
def setup_euler(n, g):
N = ceil((n + g + 5) * log2(10)) + 5
pr = primes(2 * N + 1)
pr = pr[(pr != 2) & (pr != 5)]
big = pr * pr > 2 * N + pr
return N, pr[~big], pr[big]
def adic_sum_euler(pr, n, N):
w = len(pr)
if w == 0:
return 0
Ea = np.array([max(1, sum(1 for j in range(1, 40) if int(a) ** j <= 2 * N + 1)) for a in pr], dtype=np.int64)
qs = [int(a) ** int(e) for a, e in zip(pr, Ea)]
assert max(qs) ** 2 < 2 ** 63
q = np.array(qs, dtype=np.int64)
Em = int(Ea.max())
pt = np.ones((w, Em + 1), dtype=np.int64)
for t in range(1, Em + 1):
pt[:, t] = pt[:, t - 1] * pr % q
c0 = np.array([2 * pow(10, n, qq) % qq for qq in qs], dtype=np.int64)
num = np.ones(w, dtype=np.int64)
den = num.copy()
v = np.zeros(w, dtype=np.int64)
A = np.zeros(w, dtype=np.int64)
B = num.copy()
ar = np.arange(w)
for k in range(1, N + 1):
m1 = (k % pr == 0)
m2 = ((2 * k + 1) % pr == 0)
if m1.any() or m2.any():
kv = np.full(w, k, dtype=np.int64)
sv = np.full(w, 2 * k + 1, dtype=np.int64)
for i in np.nonzero(m1)[0]:
a = int(pr[i])
x = k
while x % a == 0:
x //= a
v[i] += 1
kv[i] = x
for i in np.nonzero(m2)[0]:
a = int(pr[i])
x = 2 * k + 1
while x % a == 0:
x //= a
v[i] -= 1
sv[i] = x
num = num * (kv % q) % q
den = den * (sv % q) % q
else:
num = num * k % q
den = den * (2 * k + 1) % q
neg = v < 0
if neg.any():
e = np.where(neg, -v, 1)
T = num * (c0 * pt[ar, Ea - e] % q) % q
A = np.where(neg, (A * den + B * T) % q, A)
B = np.where(neg, B * den % q, B)
tot = 0
for i in range(w):
qq = qs[i]
r = int(A[i]) * pow(int(B[i]), -1, qq) % qq
tot += (r << SH) // qq
return tot
def block_weights(a, M):
num = 1
den = 1
fs = []
for m in range(M + 1):
fs.append(num * pow(den, -1, a) % a)
num = num * (m + 1) % a
den = den * (2 * m + 3) % a
return sum(fs[:M]) % a, fs[M]
def euler_blocks(pr, n, N):
if len(pr) == 0:
return 0
M = (N + 1) // pr
J = N - M * pr
h = (pr - 1) // 2
act = J >= h
L = np.where(act, pr - 1 - J, 0)
order = np.argsort(-L, kind='stable')
pa = pr[order]
La = L[order]
aa_act = act[order]
Lmax = int(La[0])
A = np.zeros(len(pa), dtype=np.int64)
P = np.ones_like(A)
Q = P.copy()
cnt = len(La) - np.searchsorted(La[::-1], np.arange(Lmax), side='right')
for i in range(Lmax):
c = int(cnt[i])
m = pa[:c]
A[:c] = (A[:c] + P[:c]) * (i + 1) % m
P[:c] = P[:c] * (2 * i + 1) % m
Q[:c] = Q[:c] * (i + 1) % m
tot = 0
for idx in range(len(pa)):
a = int(pa[idx])
sg = -1 if ((a - 1) // 2) & 1 else 1
F, fm = block_weights(a, (N + 1) // a)
t = sg * F
if aa_act[idx]:
Pv = int(A[idx]) * pow(int(Q[idx]), -1, a) % a if La[idx] > 0 else 0
t += fm * (sg - Pv)
c = 2 * pow(10, n, a) * t % a
tot += (c << SH) // a
return tot
def digits_euler(n, g=14):
N, small, big = setup_euler(n, g)
tot = adic_sum_euler(small, n, N) + euler_blocks(big, n, N)
return finish(tot, g)
R = log10(13.5)
def setup_g(n, g):
N = ceil((n + g + 6) / R) + 3
pr = primes(3 * N + 2)
pr = pr[(pr != 2) & (pr != 5)]
big = pr * pr > 3 * N + 2 * pr
return N, pr[~big], pr[big]
def adic_sum_g(pr, n, N):
w = len(pr)
if w == 0:
return 0
Ea = np.array([max(1, sum(1 for j in range(1, 40) if int(a) ** j <= 3 * N + 2)) for a in pr], dtype=np.int64)
qs = [int(a) ** int(e) for a, e in zip(pr, Ea)]
assert max(qs) ** 2 < 2 ** 63
q = np.array(qs, dtype=np.int64)
Em = int(Ea.max())
pt = np.ones((w, Em + 1), dtype=np.int64)
for t in range(1, Em + 1):
pt[:, t] = pt[:, t - 1] * pr % q
c0 = np.array([pow(10, n, qq) for qq in qs], dtype=np.int64)
num = np.ones(w, dtype=np.int64)
den = num.copy()
v = np.zeros(w, dtype=np.int64)
A = np.zeros(w, dtype=np.int64)
B = num.copy()
ar = np.arange(w)
def strip(x, mask, sign):
out = np.full(w, x, dtype=np.int64)
for i in np.nonzero(mask)[0]:
a = int(pr[i])
y = x
while y % a == 0:
y //= a
v[i] += sign
out[i] = y
return out
for k in range(1, N + 1):
fs = [k, 2 * k - 1]
ds = [3, 3 * k - 1, 3 * k - 2]
v0 = v.copy()
for x in fs:
mk = (x % pr == 0)
if mk.any():
num = num * (strip(x, mk, 1) % q) % q
else:
num = num * x % q
for x in ds:
mk = (x % pr == 0)
if mk.any():
den = den * (strip(x, mk, -1) % q) % q
else:
den = den * x % q
y = 50 * k - 6
mk = (y % pr == 0)
vt = v.copy()
if mk.any():
sv = np.full(w, y, dtype=np.int64)
ext = np.zeros(w, dtype=np.int64)
for i in np.nonzero(mk)[0]:
a = int(pr[i]); z = y
while z % a == 0:
z //= a; ext[i] += 1
sv[i] = z
vt = v + ext
yv = sv % q
else:
yv = y % q
neg = vt < 0
if neg.any():
e = np.where(neg, -vt, 1)
T = num * yv % q * (c0 * pt[ar, Ea - e] % q) % q
A = np.where(neg, (A * den + B * T) % q, A)
B = np.where(neg, B * den % q, B)
tot = 0
for i in range(w):
qq = qs[i]
r = int(A[i]) * pow(int(B[i]), -1, qq) % qq
tot += (r << SH) // qq
return tot
def suffix_sums(pa, lens, kind):
n_ = len(pa)
out = np.zeros(n_, dtype=np.int64)
if n_ == 0 or lens.max() <= 0:
return out
order = np.argsort(-lens, kind='stable')
a = pa[order]
L = lens[order]
if kind == 1:
h = (a - 1) // 2
top = h.copy()
P = np.array([(-int(x) * pow(2, -int(x), int(y))) % int(y) for x, y in zip(h, a)], dtype=np.int64)
else:
top = a - 1
P = np.full(n_, 2, dtype=np.int64)
Q = np.ones(n_, dtype=np.int64)
A = np.zeros(n_, dtype=np.int64)
Lm = int(L[0])
cnt = n_ - np.searchsorted(L[::-1], np.arange(Lm), side='right')
for s in range(Lm):
c = int(cnt[s])
m = a[:c]
j = top[:c] - s
A[:c] = (A[:c] + (50 * j - 6) % m * P[:c]) % m
if kind == 1:
d = 3 * j - m
nu = (d % m) * ((d - 1) % m) % m * ((d - 2) % m) % m
de = j * j % m * (2 * j - 1) % m
else:
d = 3 * j - 2 * m
e2 = 2 * j - m
nu = 2 * (d % m) % m * ((d - 1) % m) % m * ((d - 2) % m) % m
de = j * (e2 % m) % m * ((e2 - 1) % m) % m
P[:c] = P[:c] * nu % m
Q[:c] = Q[:c] * de % m
A[:c] = A[:c] * de % m
res = np.zeros(n_, dtype=np.int64)
for i in range(n_):
mm = int(a[i])
res[order[i]] = int(A[i]) * pow(int(Q[i]), -1, mm) % mm if L[i] > 0 else 0
return res
THR = 150000
def gosper_blocks(pr, n, N):
if len(pr) == 0:
return 0
M = (N + 1) // pr
J = N - M * pr
h = (pr - 1) // 2
lo1 = (pr + 2) // 3
lo2 = (2 * pr + 2) // 3
len1 = np.where((J >= lo1) & (J < h), h - J, 0)
len2 = np.where(J >= lo2, pr - 1 - J, 0)
S1 = suffix_sums(pr, np.where(len1 > THR, 0, len1), 1)
S2 = suffix_sums(pr, np.where(len2 > THR, 0, len2), 2)
for i in np.nonzero(len1 > THR)[0]:
S1[i] = tail_bgs(int(pr[i]), 1, int(len1[i]))
for i in np.nonzero(len2 > THR)[0]:
S2[i] = tail_bgs(int(pr[i]), 2, int(len2[i]))
tot = 0
for idx in range(len(pr)):
a = int(pr[idx]); Mi = int(M[idx]); Ji = int(J[idx]); hi = int(h[idx])
sg = -1 if hi & 1 else 1
i2 = pow(2, -1, a)
w1n, w1d, w2n, w2d, p2 = 1, 1, 1, 2, 1
F = 0
for m in range(Mi + 1):
u1 = w1n * pow(w1d, -1, a) % a * p2 % a
u2 = w2n * pow(w2d, -1, a) % a * p2 % a
if m < Mi:
F += u1 + u2
w1n = w1n * (m + 1) * (2 * m + 1) * (2 * m + 2) % a
w1d = w1d * (3 * m + 2) * (3 * m + 3) * (3 * m + 4) % a
w2n = w2n * (m + 1) * (2 * m + 2) * (2 * m + 3) % a
w2d = w2d * (3 * m + 3) * (3 * m + 4) * (3 * m + 5) % a
p2 = p2 * i2 % a
H1 = 0 if Ji < lo1[idx] else (-2 * sg if Ji >= hi else (-2 * sg - int(S1[idx])))
H2 = 0 if Ji < lo2[idx] else (-2 * sg - int(S2[idx]))
c = pow(10, n, a) * (2 * sg * F - u1 * H1 - u2 * H2) % a
tot += (c << SH) // a
return tot
def conv_mod(a, b, p):
n = len(a) + len(b) - 1
sz = 1 << (n - 1).bit_length()
a0 = a & 32767; a1 = a >> 15
b0 = b & 32767; b1 = b >> 15
fa0 = np.fft.rfft(a0.astype(np.float64), sz); fa1 = np.fft.rfft(a1.astype(np.float64), sz)
fb0 = np.fft.rfft(b0.astype(np.float64), sz); fb1 = np.fft.rfft(b1.astype(np.float64), sz)
r00 = np.rint(np.fft.irfft(fa0 * fb0, sz)[:n]).astype(np.int64) % p
r01 = np.rint(np.fft.irfft(fa0 * fb1 + fa1 * fb0, sz)[:n]).astype(np.int64) % p
r11 = np.rint(np.fft.irfft(fa1 * fb1, sz)[:n]).astype(np.int64) % p
return (r00 + (r01 << 15) % p + ((r11 << 15) % p << 15)) % p
def prefix_prod(x, p):
x = x.copy()
s = 1
while s < len(x):
x[s:] = x[s:] * x[:-s] % p
s <<= 1
return x
def batch_inv(x, p):
pre = prefix_prod(x, p)
inv_all = pow(int(pre[-1]), -1, p)
suf = prefix_prod(x[::-1].copy(), p)[::-1]
out = np.empty_like(x)
out[0] = inv_all * (suf[1] if len(x) > 1 else 1) % p
out[1:-1] = pre[:-2] * suf[2:] % p * inv_all % p if len(x) > 2 else out[1:-1]
if len(x) > 1:
out[-1] = pre[-2] * inv_all % p
return out
class Ctx:
def __init__(self, K, p):
self.K = K
self.p = p
f = [1] * (K + 1)
for i in range(1, K + 1):
f[i] = f[i - 1] * i % p
inv = [1] * (K + 1)
inv[K] = pow(f[K], -1, p)
for i in range(K, 0, -1):
inv[i - 1] = inv[i] * i % p
self.delta = np.array([(inv[k] * inv[K - k]) % p * (1 if (K - k) % 2 == 0 else p - 1) % p for k in range(K + 1)], dtype=np.int64)
def shift(vals, a, ctx):
K, p = ctx.K, ctx.p
t = np.arange(-K, K + 1, dtype=np.int64)
base = (a + t) % p
if (base == 0).any():
return None
e = batch_inv(base, p)
c = vals * ctx.delta % p
full = conv_mod(c, e, p)[K:2 * K + 1]
i = np.arange(K + 1)
pre0 = 1
for l in range(K + 1):
pre0 = pre0 * ((a - l) % p) % p
ratio_num = (a + i[:-1] + 1) % p
ratio_den_inv = e[i[:-1]]
steps = ratio_num * ratio_den_inv % p
prefix = np.concatenate([[pre0], steps])
pre = prefix_prod(prefix, p)
return pre * full % p
def make_steps(a, kind):
if kind == 1:
def f(j):
d = 3 * j - a
num = (d % a) * ((d - 1) % a) % a * ((d - 2) % a) % a
den = (j % a) * (j % a) % a * ((2 * j - 1) % a) % a
return den, (50 * j - 6) % a * den % a, num
h = (a - 1) // 2
P0 = (-h * pow(2, -h, a)) % a
top = h
else:
def f(j):
d = 3 * j - 2 * a
e2 = 2 * j - a
num = 2 * (d % a) % a * ((d - 1) % a) % a * ((d - 2) % a) % a
den = (j % a) * (e2 % a) % a * ((e2 - 1) % a) % a
return den, (50 * j - 6) % a * den % a, num
P0 = 2
top = a - 1
return f, P0, top
def comp(t1, t2, p):
D1, B1, N1 = t1
D2, B2, N2 = t2
return D1 * D2 % p, (D2 * B1 + B2 * N1) % p, N1 * N2 % p
def tail_bgs(a, kind, L):
f, P0, top = make_steps(a, kind)
v = max(2, isqrt(L // 3) + 1)
while True:
res = _bgs(a, kind, L, v, f, P0, top)
if res is not None:
return res
v += 1
def _bgs(a, kind, L, v, f, P0, top):
p = a
bits = bin(v)[3:]
d = 1
K = 4
x = np.arange(K + 1, dtype=np.int64)
G = tuple(np.asarray(t, dtype=np.int64) for t in f(top - v * x))
for bit in bits:
Kd = 3 * d + 1
K2 = 6 * d + 1
ctx = Ctx(Kd, p)
delta = d * pow(v, -1, p) % p
comps_first = []
comps_second = []
for c in range(3):
vals = G[c] % p
ext = shift(vals, Kd + 1, ctx)
if ext is None:
return None
first = np.concatenate([vals, ext])[:K2 + 1]
sh = []
for m in range(3):
r = shift(vals, (delta + m * (Kd + 1)) % p, ctx)
if r is None:
return None
sh.append(r)
second = np.concatenate(sh)[:K2 + 1]
comps_first.append(first)
comps_second.append(second)
G = comp(tuple(comps_first), tuple(comps_second), p)
d *= 2
if bit == '1':
Kc = 3 * d + 1
ctx = Ctx(Kc, p)
ext = []
for c in range(3):
e = shift(G[c] % p, Kc + 1, ctx)
if e is None:
return None
ext.append(np.concatenate([G[c] % p, e])[:Kc + 4])
xs = np.arange(Kc + 4, dtype=np.int64)
st = f((top - v * xs - d) % p)
G = comp(tuple(ext), tuple(np.asarray(t, dtype=np.int64) for t in st), p)
d += 1
w = L // v
tot = (1, 0, 1)
for xblk in range(w):
tot = comp(tot, (int(G[0][xblk]), int(G[1][xblk]), int(G[2][xblk])), p)
Dt, Bt, Nt = tot
A, P, Q = Bt * P0 % p, Nt * P0 % p, Dt
for s in range(v * w, L):
D, B, N = f(top - s)
A = (D * A + B * P) % p
P = P * N % p
Q = Q * D % p
return A * pow(Q, -1, p) % p
def finish(tot, g):
fr = tot & ((1 << SH) - 1)
return str(fr * 10 ** g >> SH).zfill(g)
def digits(n, g=14, mode='gosper'):
if mode == 'euler':
return digits_euler(n, g)
assert n >= 400
N, small, big = setup_g(n, g)
tot = adic_sum_g(small, n, N)
tot += gosper_blocks(big, n, N) if mode == 'gosper' else adic_sum_g(big, n, N)
return finish(tot, g)
if __name__ == "__main__":
n = int(sys.argv[1])
g = int(sys.argv[2]) if len(sys.argv) > 2 else 14
mode = sys.argv[3] if len(sys.argv) > 3 else 'gosper'
if len(sys.argv) > 4:
THR = int(sys.argv[4])
t = time.time()
d = digits(n, g, mode)
print(d[:g - 4], round(time.time() - t, 2), "s")
Download raw |
Create new paste
