以下は テスターした問題 解いたときの思考回路のメモ書きです．

とりあえず愚直解を書いておきましょう．

ll T, N, M, X, Y;
ll A[1d6], B[];

ll baka1(){
  ll i, j, res = 0;
  rep(i,N) A[i] = (X*i + Y) % M;
  rep(i,N) rep(j,i+1,N) if(A[i] > A[j]) res++;
  return res;
}

丁寧に愚直解と比較しながら進めていきたいです．

さて，項数 $N$ の数列なので何かしら一気に計算できないと駄目ですね．
とりあえず，思いつくところでは，$j$ を固定したとき，$i < j$ かつ $A[i] > A[j]$ なる $i$ の個数を求めようとすると，$(Xi+Y)\ \text{mod}\ M$ がある値以下になる個数が欲しいです．
これは floor_sum でできるってのが maspy さんが何回か言及してた気がしますね（リンク1，リンク2）．

// 他に縦横2分割にして分割統治とか，$O(\sqrt{\max(N,M)})$ みたいなのとかないかなとか考えそうにもなりましたが，あまり考えなかった．

これだけできたところでしょうがないのですが，とりあえずこの方針で進めてみましょう．
$j$ 固定で一気に計算してもしょうがないので，$A[i] > A[j]$ を式で書いてみます．
$A[i] > A[j]$ なら 1，そうでないなら 0 になる式を作りましょう．

maspy さんの教えの通りにやります．
$A[i]$ に $1$ ～ $M - A[j] + 1$ 足したときに $M$ を上回って，$M$ で割ったときに1増えるようにすることを考えると良いです．
$B[i] = Xi+Y$ として（modを取らないもの）以下のコードが書けます．

ll baka2(){
  ll i, j, res = 0;
  rep(i,N) B[i] = X*i + Y;
  rep(i,N) rep(j,i+1,N) res += fDiv(B[i] + M - 1 - (B[j] % M), M) - fDiv(B[i], M);
  return res;
}

一応，ここまで愚直解と一致です．
fDiv(a,b) は floor(a/b) です．今は非負なので / でも問題ないですが，今後のために．

で，とにかく，B[j] % M がやばそうです．
てか，$M$ で割った余りとか商とか考えているので，この辺は分離できますね．

ll baka3(){
  ll i, j, res = 0;
  rep(i,N) B[i] = X*i + Y;
  rep(i,N) rep(j,i+1,N) res += fDiv(B[i] + M - 1 - B[j], M) - fDiv(B[i], M) + fDiv(B[j], M);
  return res;
}

分離しました．
fDiv(B[i], M) は何回出てくるかカウントして一気にできますね．
後，$B[i] - B[j]$ は $X(i-j)$ です．

ll baka4(){
  ll i, j, res = 0;
  rep(i,N) B[i] = X*i + Y;
  rep(i,N) res += fDiv(B[i], M) * (2 * i + 1 - N);
  rep(i,N) rep(j,i+1,N) res += fDiv(X * (i-j) + M - 1, M);
  return res;
}

ということで，こうなります．
$i-j$ のところも何回出てくるかカウントして一気に計算したり，整理したりすると，結局こうなります．

ll baka5(){
  ll i, res = 0;
  rep(i,N) res += fDiv(X * i + Y, M) * (2 * i + 1 - N);
  rep(i,N) res += fDiv(-X * i + M - 1, M) * (N - i);
  return res;
}

このループを一気に計算したいですが，後ろの 2 * i + 1 - N とか N - i とかが定数なら floor_sum ですが i が入っています．
floor_sum は何だったか考えてみると，kyopro_friendsさんの絵（リンク）を思い出すと，格子点の数みたいな感じでしたね．
別に点に重み i を載せても重みの和を計算とかでもできそうな気もしますね．
x軸とy軸を反転させるので，x座標，y座標ともに計算できそうな重みを載せて拡張して modified_floor_sum なるものを作ってみましょう．

最初，格子点 (X,Y) に $w(X,Y) = (x_1 X + x_0)(y_1 Y + y_0)$ とか何故か非線形な重みを載せようとして失敗しましたが，普通に線形で $w(X,Y) = xX + yY + c$ の重みを載せましょう（$x,y,c$ が重みを決めるパラメータ）．

// ordinary_floor_sum(ll n, ll m, ll a, ll b)
// = sum[X=0..n-1] floor( (a*X+b) / m )
// = num of lattice points [0 <= X < n] && [0 < Y <= (Ax+B)/M]

// modified_floor_sum(ll n, ll m, ll a, ll b, ll x, ll y, ll c)
// = sum of lattice points [0 < X <= n] && [0 < Y <= (Ax+B)/M] : weight of (X,Y) = x*X + y*Y + c
// = sum[X=1..n] [let Y=floor((a*X+b)/m)] Y * (Y+1) / 2 * y + Y * (c + x*X)

ll modified_floor_sum(ll n, ll m, ll a, ll b, ll x, ll y, ll c){
  ll res = 0;
  ll h, d, w, aa, bb;

  if(a < 0){
    return modified_floor_sum(n, m, -a, b + (n+1)*a, -x, y, c+x*(n+1));
  }

  h = fDiv(b, m);
  if(h){
    res += (h * n * ( (h+1)*y + (n+1)*x )) / 2 + h * n * c;
    b -= m * h;
    c += y * h;
  }

  d = fDiv(a, m);
  res += n * (n+1) / 2 * d * ( 6*c + y*(2*d*n+d+3) + (4*n+2) * x ) / 6;
  a -= d * m;
  x += d * y;

  if(a == 0) return res;

  w = fDiv(a*n+b, m);
  aa = -m;
  bb = b + a * (n+1);
  res += modified_floor_sum(w, a, aa, bb, y, -x, c+x*(n+1));

  return res;
}

ということで，modified_floor_sum の愚直解と比較しながら，kyopro_friendsさんの絵を見ながらぐだぐだと書きました（後でコードは整理とかしましょう）．
方針は同じです．三角形領域の和を求めるのがちょっとめんどいです．長方形とか三角形を削ったときにパラメータがどう変形するか考えます．

愚直解は

ll baka_modified_floor_sum(ll n, ll m, ll a, ll b, ll x, ll y, ll c){
  ll res = 0, X, Y;
  rep(X,1,n+1){
    Y = fDiv(a*X+b, m);
    res += Y * (Y+1) / 2 * y + Y * (c + x*X);
  }
  return res;
}

です．
長方形削ったら残りは愚直解に投げるのと，最初から愚直解に投げるのとで一致しますよねーとか言いながら書きました．

ちなみに，愚直解に baka とつけてるのは診断人師匠に教わりました．
格好良い名前で作ると，愚直解をサブミットしがちなので．

さて，元の問題に戻って，1-indexed に変えたので（失敗だった気もしてきた），解もそっちに揃えて

ll baka6(){
  ll i, res = 0;
  rep(i,1,N+1) res += fDiv(X * i + Y - X, M) * (2 * i - 1 - N);
  rep(i,1,N+1) res += fDiv(-X * i + M - 1 + X, M) * (N - i + 1);
  return res;
}

となり，modified_floor_sum が使える形になったので

ll baka7(){
  ll i, res = 0;
  res += modified_floor_sum(N, M, X, Y-X, 2, 0, -1-N);
  res += modified_floor_sum(N, M, -X, M-1+X, -1, 0, N+1);
  return res;
}

として，終了です．
ではなく，オーバーフローしそうな気がするので，128ビット整数で計算した場合と比較してみます．

はい，落ちました．ダメです．
64ビットだと $N,M$ が $10^6$ 程度は大丈夫ですが $2 \times 10^6$ でダメです．
途中で 3 乗ぐらいの値が出てくるのかな？
はい，出てきますね（2乗個の点があって重みが1乗ぐらい）．
なので128ビット整数でやりましょう（か中国剰余定理？）．

ところで，これ計算できるの面白いですね．

Current time: 2024年05月06日03時22分03秒
Last modified: 2023年06月23日23時24分25秒 (by laycrs)
Tags: no_tags
トップページに戻る

Logged in as: unknown user (not login)

ログイン：