dp[i]:=gcdがiになるような部分列の総数とすると以下の関係式が成り立つ。

\begin{equation}
dp[i]=\text{gcdがiの倍数の部分列の総数}-dp[2i]-dp[3i]-dp[4i]-\cdots
\end{equation}

gcdがiの倍数になる部分列の総数を得るために、\( \sum_{k=1}^{n} k\binom{n}{k} \) を計算する必要が出てくるが、\((1+x)^n\) を微分すれば分かる。

#include <iostream>
#include <algorithm>
#include <vector>

using namespace std;

constexpr int MOD = 1e9 + 7;

struct modint {
  int n;
  modint(int n = 0) : n(n) {}
};

modint operator+(modint a, modint b) { return modint((a.n += b.n) >= MOD ? a.n - MOD : a.n); }
modint operator-(modint a, modint b) { return modint((a.n -= b.n) < 0 ? a.n + MOD : a.n); }
modint operator*(modint a, modint b) { return modint(1LL * a.n * b.n % MOD); }
modint &operator+=(modint &a, modint b) { return a = a + b; }
modint &operator-=(modint &a, modint b) { return a = a - b; }
modint &operator*=(modint &a, modint b) { return a = a * b; }

template<typename T>
modint modpow(modint a, T b) {
  modint ret = 1;
  for (; b > 0; b >>= 1, a *= a) if (b & 1) ret *= a;
  return ret;
}

int main() {
  constexpr int N = 1.01e6;
  int n;
  cin >> n;
  vector<modint> dp(N);
  for (int i = 0; i < n; i++) {
    int a;
    scanf("%d", &a);
    dp[a] += 1;
  }
  for (int i = 1; i < N; i++) {
    for (int j = i * 2; j < N; j += i) dp[i] += dp[j];
    dp[i] *= modpow(2, dp[i].n - 1);
  }
  modint ans = 0;
  for (int i = N - 1; i >= 2; i--) {
    for (int j = i * 2; j < N; j += i) dp[i] -= dp[j];
    ans += i * dp[i];
  }
  cout << ans.n << endl;
}

2017-08-04

Educational Codeforces Round 26: G. Functions On The Segments

http://codeforces.com/contest/837/problem/G

解法

二次元の総和クエリに落とし込む。高次元クエリを扱うテクニックとして永続的データ構造を用いるという方法がある。永続 segtree により解くことができるが、wavelet matrix を用いて解くこともできる。

計算量は、構築O(n log x) 質問O(log x) である。

#include <iostream>
#include <algorithm>
#include <vector>

using namespace std;

// wavelet-matrix
constexpr int N=75000*2;
constexpr int H=18;

struct foo {
  int key;
  long long a, b;
};

foo dat[N], mat[H][N+1];

void build(int l, int r, int d) {
  if(d==H || r-l==0) return;
  for (int i=l; i<r; i++) {
    mat[d][i+1].key=mat[d][i].key+(dat[i].key>>H-1-d&1);
  }
  int m=stable_partition(dat+l, dat+r, [&](foo x) { return ~x.key>>(H-1-d)&1; })-dat;
  for (int i=l; i<r; i++) {
    mat[d][i+1].a=mat[d][i].a+dat[i].a;
    mat[d][i+1].b=mat[d][i].b+dat[i].b;
  }
  build(l, m, d+1);
  build(m, r, d+1);
}

int count1(int d, int l, int r) { 
  return mat[d][r].key-mat[d][l].key;
}

int count0(int d, int l, int r) {
  return (r-l)-count1(d, l, r);
}

long long sum(int l, int r, long long x, int k) {
  int ll=0, rr=N;
  long long ret=0;
  for (int i=0; i<H; i++) {
    int mm=count0(i, ll, rr)+ll;
    int l0=count0(i, ll, l)+ll;
    int r0=count0(i, ll, r)+ll;
    if (~k>>H-1-i&1){
      l=l0;
      r=r0;
      rr=mm;
    } else {
      ret+=mat[i][r0].a*x+mat[i][r0].b;
      ret-=mat[i][l0].a*x+mat[i][l0].b;
      l=count1(i, ll, l)+mm;
      r=count1(i, ll, r)+mm;
      ll=mm;
    }
  }
  return ret;
}

int main() {
  int n;
  cin>>n;

  vector<long long> base(n+1);
  for (int i=0; i<n; i++) {
    int x1, x2, y1, A, B, y2;
    scanf("%d %d %d %d %d %d", &x1, &x2, &y1, &A, &B, &y2);
    base[i+1]=base[i]+y1;
    dat[i*2+0]={x1, A, -y1+B};
    dat[i*2+1]={x2, -A, -B+y2};
  }
  build(0, N, 0);

  long long last=0;
  int q;
  cin>>q;
  while (q--) {
    int l, r, x;
    scanf("%d %d %d", &l, &r, &x);
    l--;
    int xi=(x+last)%int(1e9);
    int t=min(xi, 200005);
    last=sum(l*2, r*2, xi, t)+base[r]-base[l];
    printf("%lld\n", last);
  }
}

2017-08-04

Educational Codeforces Round 26: E. Vasya's Function

http://codeforces.com/contest/837/problem/E

解法

b-gcd(a,b) は gcd(a,b) の倍数である。つまり次の操作でも gcd は gcd(a,b) の倍数になるはずである。何だかんだでf(a,b)=f(a/g,b/g)が成り立つ。よってf(a,b)の計算は

a,bをgcd(a,b)で割る
gcd(a,b)≠1 になるまでbを引き続ける

の繰り返しでできる。前者が起きる回数は O(log a) であるから、もし後者を高速にシミュレーションできれば解くことができる。

後者はあらかじめ a の素因数リストを計算しておけば分かる。

#include <iostream>
#include <algorithm>
#include <vector>

using namespace std;

long long gcd(long long a, long long b) { 
  while (b!=0) {
    long long tmp=a;
    a=b;
    b=tmp%b;
  }
  return a;
}

vector<long long> primes(long long n) {
  vector<long long> ret;
  for (long long i=2; i*i<=n; i++) {
    if (n%i==0){ 
      ret.push_back(i);
      while (n%i==0) n/=i;
    }
  }
  if (n!=1) ret.push_back(n);
  return ret;
}

int main() {
  long long a, b;
  cin>>a>>b;  

  auto ps=primes(a);
  long long ans=0;
  while (b!=0) {
    long long g=gcd(a, b);
    a/=g;
    b/=g;
    long long next=0;
    for (long long p:ps) if (a%p==0) {
      next=max(next, (b-1)/p*p);
    }
    ans+=b-next;
    b=next;
  }
  cout<<ans<<endl;
}

2017-07-29

yukicoder 550: 夏休みの思い出（１）

https://yukicoder.me/problems/no/550

解法

x³+ax²+bx+c=0 を解く前に x³+ax²+bx+c≡0 (mod 33331) で解く。これは全探索で解くことができ、解の個数が 3 個以下であることが保証されている。

得られた解を α,β,γ とすると、元の方程式の解は α+33331k,β+33331k,γ+33331k の形で表せる。解の範囲が分かっているため、こちらも全探索できる。

#include <iostream>
#include <set>

int main() {
    long long a, b, c;
    std::cin >> a >> b >> c;

    constexpr int MOD = 33331;
    std::set<long long> ans;
    for (__int128 i = 0; i < MOD; i++) {
        if ((i*i*i + a*i*i + b*i + c) % MOD == 0) {
            for (__int128 j = i - MOD*MOD; j <= MOD*MOD; j += MOD) {
                if (j*j*j + a*j*j + b*j + c == 0) ans.insert(j);
            }
        }
    }

    for (long long x : ans) {
        std::cout << x << ' ';
    }
}

python3 (153 bytes)

a,b,c=map(int,input().split())
P=33331
f=lambda x:(x+a)*x*x+x*b+c
print(*sorted(j for i in range(P) if f(i)%P==0 for j in range(i-P*P,P*P,P) if f(j)==0))

解の範囲が狭いことが分かっていれば、より高次の方程式に対しても同じ方法が適用できる。

2017-07-25

Codeforces #425 (Div.2) E. Vasya and Shifts

http://codeforces.com/contest/832/problem/E

問題概要

文字列 s₁,s₂,...,s_n がそれぞれ 4 つずつ――全部で4n個――与えられる
使われている文字はabcdeの5種類で、文字列長はすべて m
4n個の文字列の中からいくつか選び、その総和が b になるようなパターン数を求めるクエリが飛んでくる
文字列の和は'a'=0,...,'e'=4として (s+t)[i]=(s[i]+t[i])%5で定義する
s_i を複数回使うとき、使う順序は区別しない。

解法

文字列で考える必要なまったくなくて、単に GF(5) 上の m 次元ベクトルと考えればいい。この問題は

\begin{equation}
\mathbf{b} = c_1 \mathbf{s}_1 + c_2 \mathbf{s}_2 + \cdots +c_n \mathbf{s}_n
\end{equation}

を満たすような \( (c_1,c_2,\ldots,c_n) \) が何通りあるか、という問題に還元できる。これは連立方程式であるから、線形代数の理論を用いると良い。係数行列を A とする。解が存在するならば \(5^{n-\mathrm{rank}A}\) 通り存在する。

したがって、連立方程式が解を持つかどうかを判定できれば良い。\(b\) が与えられるたびに掃き出し法を行うのでは遅いため、\( (A\,\mathbf{b_1}\,\mathbf{b_2}\,\cdots\,\mathbf{b_q}) \)という拡大係数行列を構築しまとめて掃き出し法を行う。

掃き出し法はまあ良いと思うんだけど、解の個数があまり自明ではないので証明のアイディアをつらつらと書いておく。

連立方程式\(Ac=b\) が解 \(Ac'=b\) を持ってるなら、\(A(c-c')=0\) と変形できる
上のように変形できるのであれば、解の個数は連立方程式\(Ax=0\) の解の個数と等しいはず
\(Ax=0\)となるような\(x\)全体を行列の\(A\)のカーネルと呼ぶ
解の個数は \( |\mathrm{Ker}A| \) と表せる
\( \mathrm{Ker} A \) は線形空間である
次元定理により \( \mathrm{rank}A+\mathrm{dim}\mathrm{Ker} A = \mathrm{dim}V \)――ここで A は V→W の線形写像
\(m\) 次元の線形空間は \(m\) 個の基底が存在し、すべての元は \( \lambda_1 u_1 + \lambda_2 u_2 + \cdots + \lambda_m u_m \) の形で一意に定まる
\(\mathbb{F}_5\)上の\(m\)次元の線形空間の要素数は \(5^m\) である

#include <iostream>
#include <algorithm>
#include <vector>
#include <string>

constexpr int MOD = 5;

struct modint {
    int n;
    modint(int n = 0) : n(n) {}
};

std::string input() {
    static char buf[1000];
    scanf("%s", buf);
    return std::string(buf);
}

modint inv[5];

modint operator+(modint a, modint b) { return modint((a.n += b.n) >= MOD ? a.n - MOD : a.n); }
modint operator-(modint a, modint b) { return modint((a.n -= b.n) < 0 ? a.n + MOD : a.n); }
modint operator*(modint a, modint b) { return modint(1LL * a.n * b.n % MOD); }
modint operator/(modint a, modint b) { return a * inv[b.n]; }
modint &operator+=(modint &a, modint b) { return a = a + b; }
modint &operator-=(modint &a, modint b) { return a = a - b; }
modint &operator*=(modint &a, modint b) { return a = a * b; }
modint &operator/=(modint &a, modint b) { return a = a / b; }

int main() {
    inv[1] = 1;
    for (int i = 2; i < MOD; i++) {
        inv[i] = inv[MOD % i] * (MOD - MOD / i);
    }

    int n, m;
    std::cin >> n >> m;

    static modint a[500][800];

    for (int i = 0; i < n; i++) {
        std::string s = input();
        for (int j = 0; j < m; j++) {
            a[j][i] = s[j] - 'a';
        }
    }
    int q;
    std::cin >> q;
    for (int i = 0; i < q; i++) {
        std::string s = input();
        for (int j = 0; j < m; j++) {
            a[j][i + n] = s[j] - 'a';
        }
    }

    const int h = m;
    const int w = n + q;
    int r = 0;
    for (int j = 0; j < n && r < h; j++) {
        int p = r;
        for (int i = r; i < h; i++) {
            if (a[i][j].n != 0) {
                p = i;
                break;
            }
        }
        std::swap(a[p], a[r]);
        if (a[r][j].n == 0) continue;
        for (int jj = j + 1; jj < w; jj++) {
            a[r][jj] /= a[r][j];
        }
        a[r][j] = 1;
        for (int ii = r + 1; ii < h; ii++) {
            for (int jj = j + 1; jj < w; jj++) {
                a[ii][jj].n += 25 - a[ii][j].n * a[r][jj].n;
                a[ii][jj].n %= MOD;
            }
            a[ii][j] = 0;
        }
        r++;
    }

    constexpr long long MOD2 = 1e9 + 7;
    long long pw = 1;
    for (int i = 0; i < n - r; i++) {
        pw *= 5;
        pw %= MOD2;
    }

    for (int j = 0; j < q; j++) {
        long long ans = pw;
        for (int i = r; i < m; i++) {
            if (a[i][j + n].n != 0) {
                ans = 0;
            }
        }
        printf("%d\n", (int)ans);
    }
}

187ms。

pekempey's blog

Codeforces Round #428 (Div. 2): D. Winter is here