[ARC167] B - Product of Divisors

[Q] https://atcoder.jp/contests/arc167/tasks/arc167_b

2で割り切れない場合をなおすのに1時間以上かけて、脳がとろけてしまった。

考察
1. B=1のときを考える
2. Aを素因数分解する。
素因数の総積が分かったとして、それがA何回分になるかを考える。素直に全部掛け算するとあっという間にオーバーフローしてしまうので、小さくできないか考える。
3. Aの素因数のうち"適当に"素因数aを1つ選んで考える。aの指数をbとする。Aの約数の総積のうち、素因数aの指数cがわかれば、c/bが答えになるはず。

A = 12 = 2 * 2 * 3 を考える
A = {1,2,3,4,6,12}
  = 2^(0,1,2) * 3^(0,1) のバリエーションで6通り。

Aのある素因数a = 2を選ぶ。aの指数は2(b)。
Aの約数の総積 1*2*3*4*6*12 について、2の指数cはいくつになる？

c = {2の指数の総和 * 3のバリエーション}の取り方なので、
(0+1+2) * 2通り = 6(c)
c/bが答えになるので、6/2 = 3回割り切れる。

Q. 本当に素因数を適当に選んでいいの？
A. いいと思う。とりあえず手計算してみて、一般化を考えてみる。

A = 2*2*2*3*3*5を考える。

Aの約数の総積について、素因数2^X, 3^Y, 5^Z がいくつ含まれているかを調べる
  2      3     5 // 素因数
------------------
  3      2     1 // 指数
  6      3     1 // 指数の総和

2^X = {2の指数の総和(6) * 3の指数(2+1) * 5の指数(1+1)}
    = 6 * 3 * 2
    = 36
3^Y = {3の指数の総和(3) * 2の指数(3+1) * 5の指数(1+1)}
    = 3 * 4 * 2
    = 24
5^Z = {5の指数の総和(1) * 2の指数(3+1) * 3の指数(2+1)}
    = 1 * 4 * 3
    = 12

ある素因数aについて
a=2のとき、X / (2の指数) = 36 / 3 = 12
a=3のとき、Y / (3の指数) = 24 / 2 = 12
a=5のとき、X / (5の指数) = 12 / 1 = 12

どれを選んでも12回割り切れることがわかる。


////////////// 一般化 /////////////////////////
A = 2^x * 3^y * 5^zを考える
  2         3          5 // 素因数
------------------------
  x         y          z // 指数
 (x+1)*x/2  (y+1)*y/2  (z+1)*z/2   // 指数の総和

2^X = (x+1)*x/2*(y+1)*(z+1)
3^Y = (y+1)*y/2*(x+1)*(z+1)
5^Z = (z+1)*z/2*(x+1)*(y+1)

求めたい回答は X/x, Y/y, Z/zのどれか。本当に同じか整理すると、
X/x = (x+1) * (y+1) * (z+1) / 2
Y/y = (x+1) * (y+1) * (z+1) / 2
Z/z = (x+1) * (y+1) * (z+1) / 2

いずれも同じ解答になることが分かった。

4, B>1を考える。特に難しくない。
Aの約数総積で求まった指数をB倍すればいい。

A=12 B=2を考える

A = 12^2
  = 144
  = 2^4 * 3^2
  = {(2^2) * 3^1)}^2

2の指数2を2倍して4
3の指数1を2倍して2になるだけ。

2の指数部 = 2の指数の総和 * 3の指数のバリエーション
         = (0+1+2+3+4) / 2 * {0,1,2}
         = 10 / 2 * 3
         = 30

答えは 30/2 = 15回割り切れる。

////////////// 一般化 /////////////////////////
A = 2^x * 3^y * 5^zを考える
  2         3          5 // 素因数
------------------------
  x*B       y*B        z*B // 指数
 (x*B+1)*x*B/2  // 指数の総和

2^X = (x*B+1)*x*B/2*(y*B+1)*(z*B+1)

求めたい回答は X/x
X/x = (x*B+1) * (y*B+1) * (z*B+1) * B / 2

5, サンプル通るが、まだ罠があるようだ。21ケース落ちてる。

6, 壊れるサンプルを探す。A=4 B=1で変だ。どうやら2で割れないときに壊れてしまうようだ。

さっきの一般式をながめる。
X/xが解答。

X/x = (x*B+1) * (y*B+1) * (z*B+1) * B / 2
      ~~~~~~~   ~~~~~~~   ~~~~~~~   ~

壊れないときは次のルールがあるっぽい。
1) 指数部のバリエーション{(x*B+1) * (y*B+1)*…}に、偶数の取り方が1つでも含まれていれば問題ない。
2) Bが偶数なら問題ない。
3) 壊れた場合、指数部をxとするとx/2あまる。

指数部のバリエーションに偶数通りのものがなく、Bが偶数でないときに壊れて、その余りは指数部/2。

A=4 B=1のときNG
4=2^2
2の指数部: (0+2)*3/2 = 3
3/2 = 1...1 になってしまう。あまりは2/2

A=16 B=1のときNG
16=2^4
2の指数部: (0+4)*5/2 = 10
10/4 = 2...2 になってしまう。あまりは4/2

A=64 B=1のときNG
64=2^6
2の指数部: (0+6)*7/2 = 21
21/6 = 3...3 になってしまう。あまりは6/2

指数部をxとすると、あまりはいつもx/2になる。


A=12 B=1のときOK
12=2*2*3
2の指数部: (0+2)*3/2*2 = 6
6/2 = 3でセーフ // 指数部のバリエーションが偶数の箇所があれば問題ない

A=4 B=2のときOK
A=16なので、約数総積の2の指数部は10
Aの2の指数部は2なので
10/2 = 5回割り切れる。

////////////// 一般化 /////////////////////////
A = 2^x * 3^y * 5^zを考える
  2         3          5 // 素因数
------------------------
  x*B       y*B        z*B // 指数
 (x*B+1)*x*B/2  // 指数の総和

2^X = (x*B+1)*x*B/2*(y*B+1)*(z*B+1)

求めたい回答は X/x
X/x = (x*B+1) * (y*B+1) * (z*B+1) * B / 2

7, 2で割り切れず壊れる場合は、Aの約数総積の指数部cから、もとのAの指数bのうち、cからb/2を引いてから、c/bを答えれば大丈夫。あまりはb/2になる。
実装

int main() {
  MOD = 998244353;
  cincout();
  ll A, B;
  cin >> A >> B;

  // Aを素因数分解。指数の配列を作る
  vector<ll> exponents;
  {
    ll a = A;
    for (ll i = 2; i * i <= a; ++i) {
      ll ex = 0;
      while (a % i == 0) {
        a /= i;
        ++ex;
      }
      if (ex > 0) exponents.emplace_back(ex);
    }
    if (a != 1) {
      exponents.emplace_back(1);
    }
  }

  // exponents[0]で割った値を答えとする。
  mint all = 1; // exponents[1]~の総積
  bool is_div2 = false;  // 2で割り切れるかどうか
  for(ll i = 1; i < (ll)exponents.size(); ++i) {
    all *= (mint)exponents[i] * B + 1;
    if (exponents[i] % 2 == 1 && B % 2 == 1) is_div2 = true; // exponents[i] * B + 1 が偶数 B<=1e18なので判定わけ
  }
  if (B % 2 == 0) is_div2 = true;
  mint b = (mint)exponents[0] * B + 1; // 項数
  mint sum = (b - 1) * b / 2;  // 指数の総和 2^(0 + 1 +...+ b-1) 初項0, 末項b-1, 項数b
  mint ans = sum * all; // 2^X 個数ある
  // 2で割り切れない場合がある
  if (exponents[0] % 2 == 0 && !is_div2) {
    ans -= exponents[0] / 2; // 割り切れないときの余りはexponents[0]/2
  }
  ans /= exponents[0];
  cout << ans << endl;
}

https://atcoder.jp/contests/arc167/submissions/46644631