Asterism Stream (Codeforces 1864H)
Tags: Math, matrices

Bogocubic is playing a game with amenotiomoi. First, Bogocubic fixed an integer $$$n$$$, and then he gave amenotiomoi an integer $$$x$$$ which is initially equal to $$$1$$$.In one move amenotiomoi performs one of the following operations with the same probability:  increase $$$x$$$ by $$$1$$$;  multiply $$$x$$$ by $$$2$$$. Bogocubic wants to find the expected number of moves amenotiomoi has to do to make $$$x$$$ greater than or equal to $$$n$$$. Help him find this number modulo $$$998\,244\,353$$$.Formally, let $$$M = 998\,244\,353$$$. It can be shown that the answer can be expressed as an irreducible fraction $$$\frac{p}{q}$$$, where $$$p$$$ and $$$q$$$ are integers and $$$q \not \equiv 0 \pmod{M}$$$. Output the integer equal to $$$p \cdot q^{-1} \bmod M$$$. In other words, output such an integer $$$y$$$ that $$$0 \le y < M$$$ and $$$y \cdot q \equiv p \pmod{M}$$$.

Input: Each test contains multiple test cases. The first line contains the number of test cases $$$t$$$ ($$$1 \le t \le 100$$$). The description of the test cases follows.The only line of each test case contains one integer $$$n$$$ ($$$1 \le n \le 10^{18}$$$).

Output: For each test case, output a single integer — the expected number of moves modulo $$$998\,244\,353$$$.

Note: In the first test case, $$$n\le x$$$ without any operations, so the answer is $$$0$$$.In the second test case, for $$$n = 4$$$, here is the list of all possible sequences of operations and their probabilities:  $$$1\stackrel{+1}{\longrightarrow}2\stackrel{+1}{\longrightarrow}3\stackrel{+1}{\longrightarrow}4$$$, the probability is $$$\frac{1}{8}$$$;  $$$1\stackrel{\times 2}{\longrightarrow}2\stackrel{+1}{\longrightarrow}3\stackrel{+1}{\longrightarrow}4$$$, the probability is $$$\frac{1}{8}$$$;  $$$1\stackrel{+1}{\longrightarrow}2\stackrel{+1}{\longrightarrow}3\stackrel{\times 2}{\longrightarrow}6$$$, the probability is $$$\frac{1}{8}$$$;  $$$1\stackrel{\times 2}{\longrightarrow}2\stackrel{+1}{\longrightarrow}3\stackrel{\times 2}{\longrightarrow}6$$$, the probability is $$$\frac{1}{8}$$$;  $$$1\stackrel{+1}{\longrightarrow}2\stackrel{\times 2}{\longrightarrow}4$$$, the probability is $$$\frac{1}{4}$$$;  $$$1\stackrel{\times 2}{\longrightarrow}2\stackrel{\times 2}{\longrightarrow}4$$$, the probability is $$$\frac{1}{4}$$$. So the expected number of moves is $$$4 \cdot \left(3 \cdot \frac{1}{8}\right) + 2 \cdot \left(2 \cdot \frac{1}{4} \right) =\frac{5}{2} \equiv 499122179 \pmod{998244353}$$$. In the third test case, for $$$n = 8$$$, the expected number of moves is $$$\frac{137}{32}\equiv 717488133\pmod{998244353}$$$. In the fourth test case, for $$$n = 15$$$, the expected number of moves is $$$\frac{24977}{4096} \equiv 900515847 \pmod{998244353}$$$. 

Sample Inputs:
7
1
4
8
15
998244353
296574916252563317
494288321850420024

Sample Outputs:
0
499122179
717488133
900515847
93715054
44488799
520723508

Solution (Python):
from functools import cache
mod = 998244353

@cache
def pol(d, e=0, mod=mod):
    if d == 1:
        return [1, -1]
    b = pow(2, -2 ** (d - 2 + e), mod)
    c = [*pol(d - 1, e=e, mod=mod), 0]
    for i in range(len(c) - 1)[::-1]:
        c[i + 1] -= b * c[i]
        c[i + 1] %= mod
    return c

def solve(n, mod=mod):
    bounds = [n]
    while bounds[-1] > 1:
        bounds.append(bounds[-1] + 1 >> 1)
    s = set()

    @cache
    def e(v):
        i = 0
        while bounds[i] > v:
            i += 1
        if i == 0:
            return 0
        w = bounds[i - 1]
        d = i + 1
        if w - v <= d:
            return (1 + (e(v + 1) + e(v * 2)) * pow(2, -1, mod)) % mod
        seq = lin(w, d)
        n = w - 1 - v
        if n >= len(seq):
            seq = [*seq]
            ee = 0
            while n >= len(seq):
                c = pol(d, ee, mod=mod)
                while len(seq) < 2 * d:
                    seq.append(sum((seq[-i - 1] * -c[i + 1] for i in range(d))) % mod)
                seq = seq[n & 1::2]
                n >>= 1
                ee += 1
        return seq[n]

    @cache
    def lin(w, d):
        c = pol(d, mod=mod)
        bas = [e(w - 1 - i) for i in range(d)]
        for i in range(d + 2):
            bas.append(sum((bas[-i - 1] * -c[i + 1] for i in range(d))) % mod)
        return bas
    return e(1) % mod

def main():
    for cas in range(int(input())):
        print(solve(int(input())) % mod)
if __name__ == '__main__':
    main()

Solution (Python):
from fractions import Fraction as F
from heapq import *
from collections import defaultdict
MOD = 998244353

def inv(x):
    return pow(x, MOD - 2, MOD)
h2 = inv(2)

def modmul(a, b):
    return a * b % MOD
t = int(input())
for _ in range(t):
    n = int(input()) - 1
    if n == 0:
        print(0)
        continue
    p2 = 1
    while p2 <= n:
        p2 *= 2
    state = defaultdict(int)
    state[n, h2] = 1
    q = [(-n, h2)]
    out = 0
    uuv = 0
    seen = set()

    def upd(a, b, c):
        if (a, c) not in state:
            heappush(q, (-a, c))
        state[a, c] += b
    while q:
        uuv += 1
        (i, r) = heappop(q)
        i = -i
        v = state[i, r]
        ct = i
        nex = pow(r, ct, MOD)
        out += modmul(v, modmul(1 - nex, inv(1 - r)))
        if i % 2 == 0:
            upd(i // 2, modmul(v, h2), h2)
            v *= r
            v %= MOD
            i -= 1
        if i == 1:
            continue
        rr = 2 * r * r % MOD
        den = inv(2 - 2 * rr)
        vv = modmul(v, den)
        vv = modmul(vv, 1 + r)
        upd(i // 2, vv, h2)
        upd(i // 2, modmul(vv, -rr), modmul(h2, rr))
    print(out % MOD)

Solution (Python):
from fractions import Fraction as F
from heapq import *
from collections import defaultdict
MOD = 998244353

def inv(x):
    return pow(x, MOD - 2, MOD)
h2 = inv(2)

def modmul(a, b):
    return a * b % MOD
t = int(input())
for _ in range(t):
    n = int(input()) - 1
    if n == 0:
        print(0)
        continue
    p2 = 1
    while p2 <= n:
        p2 *= 2
    state = defaultdict(int)
    state[n, h2] = 1
    q = [(-n, h2)]
    out = 0
    uuv = 0
    seen = set()

    def upd(a, b, c):
        if (a, c) not in state:
            heappush(q, (-a, c))
        state[a, c] += b
    while q:
        uuv += 1
        (i, r) = heappop(q)
        i = -i
        v = state[i, r]
        ct = i
        nex = pow(r, ct, MOD)
        out += modmul(v, modmul(1 - nex, inv(1 - r)))
        if i % 2 == 0:
            upd(i // 2, modmul(v, h2), h2)
            v *= r
            v %= MOD
            i -= 1
        if i == 1:
            continue
        rr = 2 * r * r % MOD
        den = inv(2 - 2 * rr)
        vv = modmul(v, den)
        vv = modmul(vv, 1 + r)
        upd(i // 2, vv, h2)
        upd(i // 2, modmul(vv, -rr), modmul(h2, rr))
    print(out % MOD)