Count Disjoint Arrays
Source: codechef
Difficulty: 2709

Count Disjoint Arrays
Given an integer$M$, count the number of non-empty arrays$A$satisfying the following conditions$1 \le A_i \le M$for all$1 \le i \le |A|$$A_i$$\&$$A_j = 0$for all$1 \le i < j \le |A|$, where$\&$denotes Bitwise AND operator.
Since the answer may be large, find it modulo$998244353$.
Input Format
The first line of input will contain a single integer$T$, denoting the number of test cases.
The first and only line of each test case contains a single integer$M$.
Output Format
For each test case, output on a new line the total number of arrays satisfying the conditions modulo$998244353$.
Constraints$1 \le T \le 10^4$$1 \le M \le 10^{18}$Sample 1:
Explanation:
Test case$1$: There is only$1$array:$[1]$.
Test case$2$: There are$4$arrays as enlisted below:$[1]$$[2]$$[1, 2]$$[2, 1]$

Sample Input:
9
1
2
3
4
6
8
13
16
1000000000000000000

Sample Output:
1
4
5
18
24
100
145
690
771115550

C++ solution:
#include <bits/stdc++.h>
using namespace std;
#define all(v) (v).begin(),(v).end()
#define pb(a) push_back(a)
#define rep(i, n) for(int i=0;i<n;i++)
#define foa(e, v) for(auto&& e : v)
using ll = long long;
const ll MOD7 = 1000000007, MOD998 = 998244353, INF = (1LL << 60);
#define dout(a) cout<<fixed<<setprecision(10)<<a<<endl;

template<int MOD> struct Modint {
    long long val;
    constexpr Modint(long long v = 0) noexcept : val(v % MOD) { if (val < 0) val += MOD; }
    constexpr int mod() const { return MOD; }
    constexpr long long value() const { return val; }
    constexpr Modint operator - () const noexcept { return val ? MOD - val : 0; }
    constexpr Modint operator + (const Modint& r) const noexcept { return Modint(*this) += r; }
    constexpr Modint operator - (const Modint& r) const noexcept { return Modint(*this) -= r; }
    constexpr Modint operator * (const Modint& r) const noexcept { return Modint(*this) *= r; }
    constexpr Modint operator / (const Modint& r) const noexcept { return Modint(*this) /= r; }
    constexpr Modint& operator += (const Modint& r) noexcept {
        val += r.val;
        if (val >= MOD) val -= MOD;
        return *this;
    }
    constexpr Modint& operator -= (const Modint& r) noexcept {
        val -= r.val;
        if (val < 0) val += MOD;
        return *this;
    }
    constexpr Modint& operator *= (const Modint& r) noexcept {
        val = val * r.val % MOD;
        return *this;
    }
    constexpr Modint& operator /= (const Modint& r) noexcept {
        long long a = r.val, b = MOD, u = 1, v = 0;
        while (b) {
            long long t = a / b;
            a -= t * b, swap(a, b);
            u -= t * v, swap(u, v);
        }
        val = val * u % MOD;
        if (val < 0) val += MOD;
        return *this;
    }
    constexpr bool operator == (const Modint& r) const noexcept { return this->val == r.val; }
    constexpr bool operator != (const Modint& r) const noexcept { return this->val != r.val; }
    friend constexpr istream& operator >> (istream& is, Modint<MOD>& x) noexcept {
        is >> x.val;
        x.val %= MOD;
        if (x.val < 0) x.val += MOD;
        return is;
    }
    friend constexpr ostream& operator << (ostream& os, const Modint<MOD>& x) noexcept {
        return os << x.val;
    }
    constexpr Modint<MOD> pow(long long n) noexcept {
        if (n == 0) return 1;
        if (n < 0) return this->pow(-n).inv();
        Modint<MOD> ret = pow(n >> 1);
        ret *= ret;
        if (n & 1) ret *= *this;
        return ret;
    }
    constexpr Modint<MOD> inv() const noexcept {
        long long a = this->val, b = MOD, u = 1, v = 0;
        while (b) {
            long long t = a / b;
            a -= t * b, swap(a, b);
            u -= t * v, swap(u, v);
        }
        return Modint<MOD>(u);
    }
};

const int MOD = MOD998;
using mint = Modint<MOD>;

struct Combination {
    long long C_MOD;
    vector<long long> fac, finv, inv;
    Combination(long long n, long long mod) noexcept : C_MOD(mod) {
        n = max(n, 2LL);
        fac.resize(n, 0);
        finv.resize(n, 0);
        inv.resize(n, 0);
        fac[0] = fac[1] = finv[0] = finv[1] = inv[1] = 1;
        for(int i = 2; i < n; i ++) {
            fac[i] = fac[i - 1] * i % mod;
            inv[i] = mod - inv[mod % i] * (mod / i) % mod;
            finv[i] = finv[i - 1] * inv[i] % mod;
        }
    }
    long long com(long long n, long long k) {
        if(n < k || n < 0 || k < 0) return 0;
        return fac[n] * (finv[k] * finv[n - k] % C_MOD) % C_MOD;
    }
};

mint memo[100][100];
bool thr[100][100];

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    Combination comb(10000, MOD998);
    memset(thr, 0, sizeof(thr));
    int t; cin >> t;

    auto calc = [&](auto calc, int n, int k) -> mint {
        if(n < 0 or k < 0) return mint(0);
        if(thr[n][k]) return memo[n][k];
        // n bit not used
        // k want use
        if(k == 0) {
            thr[n][k] = 1;
            return memo[n][k] = 1;
        } 
        mint ret = 0;
        for(int i = 1; i <= n; i ++) {
            ret += calc(calc, n - i, k - 1) * mint(comb.com(n, i));
        }
        thr[n][k] = 1;
        return memo[n][k] = ret;
    };
    bool ex[100][100];
    memset(ex, 0, sizeof(ex));
    mint dp[100][100];
    auto C = [&](int cnt, int i) -> mint {
        if(ex[cnt][i]) return dp[cnt][i];
        mint ans = 0;
        for(int j = 0; j <= i; j ++) {
            for(int k = 0; k <= 60; k ++) {
                ans += calc(calc, cnt - j, k) * mint(comb.com(i, j)) * mint(k + 1);
            }
        }
        ex[cnt][i] = 1;
        return dp[cnt][i] = ans;
    };
    while(t --) {
        ll m; cin >> m;
        int id = 0;
        for(int i = 0; i < 60; i ++) {
            if(m >> i & 1) {
                id = i;
            }
        }
        mint ans = 0;
        for(int i = 1; i <= 60; i ++) {
            ans += calc(calc, id, i);
        }
        int cnt = id;
        for(ll i = id - 1; i >= 0; i --) {
            if(m >> i & 1) {
                ans += C(cnt, i);
                cnt --;
            }
        }
        for(int j = 0; j <= 60; j ++) {
            ans += calc(calc, cnt, j) * mint(j + 1);
        }
        cout << ans << endl;
    }
    return 0;
}

Go solution:
package main

import (
    "bufio"
    "os"
    "strconv"
    "strings"
)

const mod = 998244353
const maxn = 60

var fact, invfact [maxn + 1]int64
var binom [maxn + 1][maxn + 1]int64
var stirling [maxn + 1][maxn + 1]int64
var obel, full, Uarr [maxn + 1]int64

func main() {
    precompute()
    in := bufio.NewReader(os.Stdin)
    out := bufio.NewWriter(os.Stdout)
    defer out.Flush()
    line, _ := in.ReadString('\n')
    T, _ := strconv.Atoi(strings.TrimSpace(line))
    for t := 0; t < T; t++ {
        l, _ := in.ReadString('\n')
        M, _ := strconv.ParseInt(strings.TrimSpace(l), 10, 64)
        L := bits(M)
        R := M - (1 << (L - 1))
        n := L - 1
        part1 := full[n] % mod
        cnt := countSub(n, R)
        part2 := int64(0)
        for s := 0; s <= n; s++ {
            part2 = (part2 + cnt[s]*Uarr[n-s]) % mod
        }
        ans := (part1 + part2) % mod
        out.WriteString(strconv.FormatInt(ans, 10) + "\n")
    }
}

func bits(x int64) int {
    b := 0
    for x > 0 {
        b++
        x >>= 1
    }
    return b
}

func precompute() {
    fact[0] = 1
    for i := 1; i <= maxn; i++ {
        fact[i] = (fact[i-1] * int64(i)) % mod
    }
    invfact[maxn] = modInv(fact[maxn])
    for i := maxn; i >= 0; i-- {
        if i > 0 {
            invfact[i-1] = (invfact[i] * int64(i)) % mod
        }
    }
    for i := 0; i <= maxn; i++ {
        for j := 0; j <= i; j++ {
            binom[i][j] = (((fact[i] * invfact[j]) % mod) * invfact[i-j]) % mod
        }
    }
    stirling[0][0] = 1
    for n := 1; n <= maxn; n++ {
        for k := 1; k <= n; k++ {
            stirling[n][k] = (stirling[n-1][k-1] + int64(k)*stirling[n-1][k]) % mod
        }
    }
    for n := 1; n <= maxn; n++ {
        sum := int64(0)
        for k := 1; k <= n; k++ {
            sum = (sum + (fact[k] * stirling[n][k]) % mod) % mod
        }
        obel[n] = sum % mod
    }
    full[0] = 0
    for n := 1; n <= maxn; n++ {
        s := int64(0)
        for r := 1; r <= n; r++ {
            s = (s + binom[n][r]*obel[r]) % mod
        }
        full[n] = s % mod
    }
    Uarr[0] = 1
    for n := 1; n <= maxn; n++ {
        s := int64(0)
        for r := 1; r <= n; r++ {
            inner := int64(0)
            for k := 1; k <= r; k++ {
                inner = (inner + (int64(k+1)*fact[k]%mod*stirling[r][k])%mod) % mod
            }
            s = (s + binom[n][r]*inner) % mod
        }
        Uarr[n] = (1 + s) % mod
    }
}

func modInv(a int64) int64 {
    return modPow(a, mod-2)
}

func modPow(a, b int64) int64 {
    res := int64(1)
    a %= mod
    for b > 0 {
        if b&1 == 1 {
            res = (res * a) % mod
        }
        a = (a * a) % mod
        b >>= 1
    }
    return res
}

func countSub(n int, R int64) []int64 {
    cnt := make([]int64, n+1)
    ones := 0
    for i := 0; i < n; i++ {
        bit := (R >> (n - 1 - i)) & 1
        if bit == 1 {
            for j := 0; j <= n-1-i; j++ {
                cnt[ones+j] = (cnt[ones+j] + binom[n-1-i][j]) % mod
            }
            ones++
        }
    }
    cnt[ones] = (cnt[ones] + 1) % mod
    return cnt
}

Java solution:
import java.util.*;
import java.lang.*;
import java.io.*;

class Codechef {
    static final int MOD = 998244353;
    static final int BIT_MAX = 60;
    static long[] factorial = new long[BIT_MAX + 1];
    static long[] invFactorial = new long[BIT_MAX + 1];
    static long[][] binomial = new long[BIT_MAX + 1][BIT_MAX + 1];
    static long[][] stirling = new long[BIT_MAX + 1][BIT_MAX + 1];
    static long[] orderedBell = new long[BIT_MAX + 1];
    static long[] completeSet = new long[BIT_MAX + 1];
    static long[] auxArray = new long[BIT_MAX + 1];

    public static void main(String[] args) throws java.lang.Exception {
        Scanner scan = new Scanner(System.in);
        initializeTables();
        
        int testCases = scan.nextInt();
        StringBuilder result = new StringBuilder();
        
        int caseIdx = 0;
        while (caseIdx < testCases) {
            long targetVal = scan.nextLong();
            int bitLen = Long.SIZE - Long.numberOfLeadingZeros(targetVal);
            long offsetVal = targetVal - (1L << (bitLen - 1));
            int width = bitLen - 1;
            
            long firstPart = completeSet[width] % MOD;
            long[] validSets = calculateValidSets(width, offsetVal);
            long secondPart = 0;
            
            int idx = 0;
            while (idx <= width) {
                secondPart = (secondPart + validSets[idx] * auxArray[width - idx]) % MOD;
                idx++;
            }
            
            long answer = (firstPart + secondPart) % MOD;
            result.append(answer).append("\n");
            caseIdx++;
        }
        
        System.out.print(result);
        scan.close();
    }

    static long modPow(long base, long exponent) {
        long res = 1;
        base %= MOD;
        
        while (exponent > 0) {
            if ((exponent & 1) == 1) {
                res = multiply(res, base);
            }
            base = multiply(base, base);
            exponent >>= 1;
        }
        
        return res;
    }

    static long modInverse(long number) {
        return modPow(number, MOD - 2);
    }

    static long multiply(long a, long b) {
        return (a * b) % MOD;
    }

    static void initializeTables() {
        factorial[0] = 1;
        
        int idx = 1;
        while (idx <= BIT_MAX) {
            factorial[idx] = multiply(factorial[idx - 1], idx);
            idx++;
        }
        
        invFactorial[BIT_MAX] = modInverse(factorial[BIT_MAX]);
        
        idx = BIT_MAX;
        while (idx > 0) {
            invFactorial[idx - 1] = multiply(invFactorial[idx], idx);
            idx--;
        }
        
        int row = 0;
        while (row <= BIT_MAX) {
            int col = 0;
            while (col <= row) {
                binomial[row][col] = multiply(multiply(factorial[row], invFactorial[col]), invFactorial[row - col]);
                col++;
            }
            row++;
        }
        
        stirling[0][0] = 1;
        
        row = 1;
        while (row <= BIT_MAX) {
            int col = 1;
            while (col <= row) {
                stirling[row][col] = (stirling[row - 1][col - 1] + multiply(col, stirling[row - 1][col])) % MOD;
                col++;
            }
            row++;
        }
        
        idx = 1;
        while (idx <= BIT_MAX) {
            long sum = 0;
            int k = 1;
            while (k <= idx) {
                sum = (sum + multiply(factorial[k], stirling[idx][k])) % MOD;
                k++;
            }
            orderedBell[idx] = sum;
            idx++;
        }
        
        completeSet[0] = 0;
        
        idx = 1;
        while (idx <= BIT_MAX) {
            long sum = 0;
            int r = 1;
            while (r <= idx) {
                sum = (sum + multiply(binomial[idx][r], orderedBell[r])) % MOD;
                r++;
            }
            completeSet[idx] = sum;
            idx++;
        }
        
        auxArray[0] = 1;
        
        idx = 1;
        while (idx <= BIT_MAX) {
            long outerSum = 0;
            int r = 1;
            while (r <= idx) {
                long innerSum = 0;
                int k = 1;
                while (k <= r) {
                    innerSum = (innerSum + multiply(multiply(k + 1, factorial[k]), stirling[r][k])) % MOD;
                    k++;
                }
                outerSum = (outerSum + multiply(binomial[idx][r], innerSum)) % MOD;
                r++;
            }
            auxArray[idx] = (1 + outerSum) % MOD;
            idx++;
        }
    }

    static long[] calculateValidSets(int bits, long target) {
        long[] dp = new long[bits + 1];
        int oneCount = 0;
        
        int pos = 0;
        while (pos < bits) {
            long currentBit = (target >> (bits - 1 - pos)) & 1;
            
            if (currentBit == 1) {
                int j = 0;
                while (j < bits - pos) {
                    dp[oneCount + j] = (dp[oneCount + j] + binomial[bits - 1 - pos][j]) % MOD;
                    j++;
                }
                oneCount++;
            }
            pos++;
        }
        
        dp[oneCount] = (dp[oneCount] + 1) % MOD;
        return dp;
    }
}

JavaScript solution:
const MOD = 998244353n;
const MAXN = 60;

// Precompute binomial coefficients C[n][k]
const C = Array.from({ length: MAXN + 1 }, () => Array(MAXN + 1).fill(0n));
for (let n = 0; n <= MAXN; n++) {
    C[n][0] = 1n;
    for (let k = 1; k <= n; k++) {
        C[n][k] = (C[n - 1][k - 1] + C[n - 1][k]) % MOD;
    }
}

// dp[n][k] = Stirling numbers of second kind S(n,k) multiplied by k!
const dp = Array.from({ length: MAXN + 1 }, () => Array(MAXN + 1).fill(0n));
dp[0][0] = 1n;
for (let n = 1; n <= MAXN; n++) {
    for (let k = 1; k <= n; k++) {
        dp[n][k] = (BigInt(k) * dp[n - 1][k] + BigInt(k) * dp[n - 1][k - 1]) % MOD;
    }
}

// a[n] = total ordered partitions = sum_{k} dp[n][k]
const a = Array(MAXN + 1).fill(0n);
a[0] = 1n;
for (let n = 1; n <= MAXN; n++) {
    let sum = 0n;
    for (let k = 1; k <= n; k++) {
        sum = (sum + dp[n][k]) % MOD;
    }
    a[n] = sum;
}

// c[n] = sum_{k=0..n} dp[n][k] * (k+1)
const c = Array(MAXN + 1).fill(0n);
for (let n = 0; n <= MAXN; n++) {
    let s = 0n;
    for (let k = 0; k <= n; k++) {
        s = (s + dp[n][k] * BigInt(k + 1)) % MOD;
    }
    c[n] = s;
}

// S2[t] = sum_{u=0..t} C[t][u] * c[u]
const S2 = Array(MAXN + 1).fill(0n);
for (let t = 0; t <= MAXN; t++) {
    let s = 0n;
    for (let u = 0; u <= t; u++) {
        s = (s + C[t][u] * c[u]) % MOD;
    }
    S2[t] = s;
}

// F0[W] = sum_{t=1..W} C[W][t]*a[t]
const F0 = Array(MAXN + 1).fill(0n);
for (let W = 1; W <= MAXN; W++) {
    let s = 0n;
    for (let t = 1; t <= W; t++) {
        s = (s + C[W][t] * a[t]) % MOD;
    }
    F0[W] = s;
}
F0[0] = 0n;

// Read input and process test cases
const readline = require('readline');
const rl = readline.createInterface({
    input: process.stdin,
    output: process.stdout,
    terminal: false
});

const input = [];
rl.on('line', (line) => {
    input.push(line.trim());
}).on('close', () => {
    const T = parseInt(input[0], 10);
    const output = [];
    
    for (let i = 1; i <= T; i++) {
        const M = BigInt(input[i]);
        if (M === 0n) {
            output.push("0");
            continue;
        }
        
        // Calculate number of bits
        const W = M.toString(2).length - 1;
        const M_low = M - (1n << BigInt(W));
        
        // Count numbers with specific bit patterns
        const cnt = Array(W + 1).fill(0n);
        let ones = 0;
        for (let i = W - 1; i >= 0; i--) {
            if ((M_low >> BigInt(i)) & 1n) {
                for (let k = 0; k <= i; k++) {
                    cnt[ones + k] = (cnt[ones + k] + C[i][k]) % MOD;
                }
                ones++;
            }
        }
        cnt[ones] = (cnt[ones] + 1n) % MOD;
        
        // Calculate S1
        let S1 = 0n;
        for (let m = 0; m <= W; m++) {
            S1 = (S1 + cnt[m] * S2[W - m]) % MOD;
        }
        
        const ans = (F0[W] + S1) % MOD;
        output.push(ans.toString());
    }
    
    console.log(output.join('\n'));
});

Python solution:
# cook your dish here
MOD = 998244353
MAXN = 60

# 1) Precompute binomial C[n][k]
C = [[0]*(MAXN+1) for _ in range(MAXN+1)]
for n in range(MAXN+1):
    C[n][0] = 1
    for k in range(1, n+1):
        C[n][k] = (C[n-1][k-1] + C[n-1][k]) % MOD

# 2) dp[n][k] = S(n,k) * k!  (ordered partitions into k blocks)
dp = [[0]*(MAXN+1) for _ in range(MAXN+1)]
dp[0][0] = 1
for n in range(1, MAXN+1):
    for k in range(1, n+1):
        # S(n,k) = k*S(n-1,k) + S(n-1,k-1)
        # multiply by k! gives
        dp[n][k] = (k * dp[n-1][k] + k * dp[n-1][k-1]) % MOD

# 3) a[n] = total ordered partitions = sum_{k} dp[n][k]
a   = [0]*(MAXN+1)
a[0] = 1
for n in range(1, MAXN+1):
    a[n] = sum(dp[n][k] for k in range(1, n+1)) % MOD

# 4) c[n] = sum_{k=0..n} dp[n][k] * (k+1)
c = [0]*(MAXN+1)
for n in range(0, MAXN+1):
    s = 0
    for k in range(0, n+1):
        s = (s + dp[n][k] * (k+1)) % MOD
    c[n] = s

# 5) S2[t] = sum_{u=0..t} C[t][u] * c[u]
S2 = [0]*(MAXN+1)
for t in range(MAXN+1):
    s = 0
    for u in range(t+1):
        s = (s + C[t][u] * c[u]) % MOD
    S2[t] = s

# 6) F0[W] = sum_{t=1..W} C[W][t]*a[t]   (the f_full(W) for W‐bits)
F0 = [0]*(MAXN+1)
for W in range(1, MAXN+1):
    s = 0
    for t in range(1, W+1):
        s = (s + C[W][t] * a[t]) % MOD
    F0[W] = s
F0[0] = 0  # no nonempty partition on 0 bits

import sys
input = sys.stdin.readline

T = int(input())
out = []
for _ in range(T):
    M = int(input())
    if M == 0:
        out.append("0")
        continue
    # W = floor(log2 M)
    W = M.bit_length() - 1
    # remainder after stripping top bit
    M_low = M - (1 << W)

    # build cnt[0..W]
    cnt = [0]*(W+1)
    ones = 0
    for i in range(W-1, -1, -1):
        if (M_low >> i) & 1:
            # we choose to put 0 at bit i, then any k of the lower i bits
            for k in range(i+1):
                cnt[ones + k] = (cnt[ones + k] + C[i][k]) % MOD
            ones += 1
    # include the exact M_low itself
    cnt[ones] = (cnt[ones] + 1) % MOD

    # S1 = sum_{m=0..W} cnt[m]*S2[W-m]
    S1 = 0
    for m in range(W+1):
        S1 = (S1 + cnt[m] * S2[W-m]) % MOD

    ans = (F0[W] + S1) % MOD
    out.append(str(ans))

print("\n".join(out))

Ruby solution:
MOD = 998244353
MAXN = 60

# Precompute binomial coefficients C[n][k]
C = Array.new(MAXN + 1) { Array.new(MAXN + 1, 0) }
(0..MAXN).each do |n|
  C[n][0] = 1
  (1..n).each do |k|
    C[n][k] = (C[n-1][k-1] + C[n-1][k]) % MOD
  end
end

# dp[n][k] = Stirling numbers of second kind S(n,k) multiplied by k!
dp = Array.new(MAXN + 1) { Array.new(MAXN + 1, 0) }
dp[0][0] = 1
(1..MAXN).each do |n|
  (1..n).each do |k|
    dp[n][k] = (k * dp[n-1][k] + k * dp[n-1][k-1]) % MOD
  end
end

# a[n] = total ordered partitions = sum_{k} dp[n][k]
a = Array.new(MAXN + 1, 0)
a[0] = 1
(1..MAXN).each do |n|
  sum = 0
  (1..n).each do |k|
    sum = (sum + dp[n][k]) % MOD
  end
  a[n] = sum
end

# c[n] = sum_{k=0..n} dp[n][k] * (k+1)
c = Array.new(MAXN + 1, 0)
(0..MAXN).each do |n|
  sum = 0
  (0..n).each do |k|
    sum = (sum + dp[n][k] * (k + 1)) % MOD
  end
  c[n] = sum
end

# S2[t] = sum_{u=0..t} C[t][u] * c[u]
S2 = Array.new(MAXN + 1, 0)
(0..MAXN).each do |t|
  sum = 0
  (0..t).each do |u|
    sum = (sum + C[t][u] * c[u]) % MOD
  end
  S2[t] = sum
end

# F0[w] = sum_{t=1..w} C[w][t]*a[t]
F0 = Array.new(MAXN + 1, 0)
(1..MAXN).each do |w|
  sum = 0
  (1..w).each do |t|
    sum = (sum + C[w][t] * a[t]) % MOD
  end
  F0[w] = sum
end
F0[0] = 0

# Read input and process
input = $stdin.read.split("\n")
test_cases = input[0].to_i
output = []

(1..test_cases).each do |i|
  m = input[i].to_i
  if m == 0
    output << "0"
    next
  end
  
  bits = m.bit_length - 1
  remainder = m - (1 << bits)

  # Build count array
  count = Array.new(bits + 1, 0)
  ones_count = 0
  (bits - 1).downto(0) do |i|
    if (remainder >> i) & 1 == 1
      (0..i).each do |k|
        count[ones_count + k] = (count[ones_count + k] + C[i][k]) % MOD
      end
      ones_count += 1
    end
  end
  count[ones_count] = (count[ones_count] + 1) % MOD

  # Calculate the sum
  current_sum = 0
  (0..bits).each do |m|
    current_sum = (current_sum + count[m] * S2[bits - m]) % MOD
  end

  result = (F0[bits] + current_sum) % MOD
  output << result.to_s
end

puts output.join("\n")