アルゴリズムと数学 101

https://atcoder.jp/contests/math-and-algorithm/tasks/typical90_o

選んだボールの1、選ばれなかったボールを0と表すと、例えばN=5で1と4が選ばれたら、10010と表現できます。
出力例を見ると、下のほうは差が等差数列になっていて簡単そうですね。一方、kが小さいとDPが使えます。例えば、k=2なら、左から考えて、直前に1があった状態と次に1を置ける状態の2つがあって、

$\begin{pmatrix} 0 \\ 1 \\ \end{pmatrix}$

という状態に、

$\begin{pmatrix} 0 & 1 \\ 1 & 1 \\ \end{pmatrix}^N$

という行列を作用すれば答えが出るので、kが小さければ速いです。ただ、これだと個々の計算量が $O(k^3\log{N})$ なのであまり大きいkには使えません。
もう少し詳しく見てみましょう。k=5だとして、(a, b, c, d, e)を一歩進めると、(e, a, b, c, d+e)となって、計算量O(k)で処理できます。しかし、e, a, b, cは移動しているだけなので、これを固定して考えると、(a, b, c, d, e) -> (a, b, c, d+e, e) -> (a, b, c+d+e, d+e, e)となって、変化する場所が変わっていくだけなので、この処理はO(1)となって、一つのkの全体でO(N)で処理できます。ただ、全体で $O(N^2)$ になってしまうので、kがある程度小さくないといけません。

kが大きい方は、定番の考え方で、例えばk=5なら、10000を塊と考えて、0と10000の組み合わせを考えれば、選ばれたボールの距離は必ず5以上になります。ただし、これだと右端のほうで1がある場合が考えられないので、一番右の1の位置を固定して、その左をさきほどの組み合わせとします。10000のような塊がいくつあるかを固定して、まず0個で、一番右の1の位置を変えると、

$_{N-1}C_0 + \cdots _0C_0 = _NC_1$
となります。塊が1つある場合は、
$_{N-k}C_1 + \cdots _1C_1 = _{N-k+1}C_2$
となります。一般に塊がnある場合は、
$_{N-n(k-1)-1}C_n + \cdots _nC_n = _{N-n(k-1)}C_{n+1}$
となります。
kが大きければ小さなnだけ考えればよいです。また、kが大きいとき、nの最大が同じになるkが多いので、上の組み合わせをkの多項式と考えて、多項式の和を作っておけば一つのkの計算量が $O(N/k)$ なるので、 $k \ge N^{1/2}$ 全体で計算量が $O(N^{3/2})$ となります。kが小さい方も $O(N^{3/2})$ となるので、全体でも同じ表記になります。

// Don't be too close
#![allow(non_snake_case)]


//////////////////// Constants ////////////////////

const D: i64 = 10i64.pow(9) + 7;


//////////////////// library ////////////////////

fn read<T: std::str::FromStr>() -> T {
    let mut line = String::new();
    std::io::stdin().read_line(&mut line).ok();
    line.trim().parse().ok().unwrap()
}

// ax = by + 1 (a, b > 0)
fn linear_diophantine(a: i64, b: i64) -> Option<(i64, i64)> {
    if a == 1 {
        return Some((1, 0))
    }
    
    let q = b / a;
    let r = b % a;
    if r == 0 {
        return None
    }
    let (x1, y1) = linear_diophantine(r, a)?;
    Some((-q * x1 - y1, -x1))
}

fn inverse(a: i64, d: i64) -> i64 {
    let (x, _y) = linear_diophantine(a, d).unwrap();
    if x >= 0 {
        x % d
    }
    else {
        x % d + d
    }
}


//////////////////// Polynomial ////////////////////

use std::ops::{Add, Mul};

struct Polynomial {
    a: Vec<i64>
}

impl Add for Polynomial {
    type Output = Self;
    
    fn add(self, other: Self) -> Self {
        if self.len() < other.len() {
            other + self
        }
        else {
            let mut a = self.a.to_vec();
            for (i, c) in other.a.iter().enumerate() {
                a[i] += c
            }
            Self { a }
        }
    }
}

impl Mul for Polynomial {
    type Output = Self;
    
    fn mul(self, other: Self) -> Self {
        let L1 = self.len();
        let L2 = other.len();
        let mut a: Vec<i64> = vec![0; L1+L2-1];
        for i in 0..L1 {
            for j in 0..L2 {
                a[i+j] = (a[i+j] + self.a[i] * other.a[j]) % D
            }
        }
        Self { a }
    }
}

impl Polynomial {
    fn len(&self) -> usize {
        self.a.len()
    }
    
    fn apply(&self, x: i64) -> i64 {
        let v = self.a.iter().rev().fold(0, |y, &c| (y*x+c)%D);
        (v + D) % D
    }
}


//////////////////// process ////////////////////

fn F_small_each(k: usize, N: usize) -> i64 {
    let mut a: Vec<i64> = vec![0; k];
    a[k-1] = 1;
    for i in 0..N {
        let p = ((N+1)*k - i - 2) % k;
        a[p] = (a[p] + a[(p+1)%k]) % D
    }
    (a.into_iter().sum::<i64>() - 1) % D
}

fn F_small(M: usize, N: usize) {
    for k in 1..M+1 {
        println!("{}", F_small_each(k, N))
    }
}

fn F_large(M1: usize, N: usize) {
    let mut p = Polynomial { a: vec![N as i64] };
    let mut fact_invs = vec![1; M1+1];
    for i in 2..M1+1 {
        fact_invs[i] = fact_invs[i-1] * inverse(i as i64, D) % D
    }
    
    let mut results: Vec<i64> = vec![N as i64; N+1];
    for n in 1..M1 {
        let mut p1 = Polynomial { a: vec![1] };
        for i in 0..n+1 {
            p1 = p1 * Polynomial { a: vec![(N+n-i) as i64, -(n as i64)] };
        }
        p1 = p1 * Polynomial { a: vec![fact_invs[n+1]] };
        p = p + p1;
        for k in (N-1)/(n+1)+1..(N-1)/n+1 {
            results[k] = p.apply(k as i64)
        }
    }
    
    for k in (N-1)/M1+1..N+1 {
        println!("{}", results[k])
    }
}

fn F(N: usize) {
    let M1 = (1..).filter(|k| k*k >= N).next().unwrap();
    let M = (N-1)/M1;
    F_small(M, N);
    F_large(M1, N)
}

fn main() {
    let N = read();
    F(N)
}