Yukicoder No.456 Millions of Submits!

,

http://yukicoder.me/problems/no/456

面白かった。

この整列するテクは予想以上に効いたので覚えておきたい。 ただし二分探索だと、区間がパラメタなので適用が難しい。

solution

Newton法 + 精度検査 + バッチ処理。計算量は分からないが、平均イテレーション回数は$2.34$回であった。

Newton法。 $f(x) = y$を求めるには$x’ \gets x - \frac{f(x) - y}{f’(x)}$を繰り返す。 今回$f(x) = x^a(\log x)^b$で、$f’(x) = ax^{a-1} \cdot (\log x)^b + x^a \cdot b\frac{1}{x}(\log x)^{b-1}$。 二次収束であり速いが、二分探索より汎用性は劣る。

固定の回数でなく精度を検査して収束してそうなら打ち切りにする。 そうしないと高い収束速度が生かせない。

しかしまだ定数倍間に合わない。 $0 \le a, b \le 10$と小さいので、$(a_i, b_i)$が同じ$t_i$をまとめて整列させ、端から順にNewton法を行い、その際初期値に隣のものの収束先を利用する。 手元の乱択ケースで$\frac{1}{10}$倍ぐらいに高速化されて通る。

implementation

#include <cstdio>
#include <vector>
#include <algorithm>
#include <array>
#include <cmath>
#include <cassert>
#define repeat(i,n) for (int i = 0; (i) < int(n); ++(i))
#define repeat_reverse(i,n) for (int i = (n)-1; (i) >= 0; --(i))
#define whole(f,x,...) ([&](decltype((x)) whole) { return (f)(begin(whole), end(whole), ## __VA_ARGS__); })(x)
using namespace std;
const double eps = 1e-10;
double newton(int a, int b, double t, double n) {
    assert (a >= 1 and b >= 1);
    auto f = [&](double n) {
        return pow(n, a) * pow(log(n), b);
    };
    auto f1 = [&](double n) {
        return a * pow(n, a-1) * pow(log(n), b) + pow(n, a) * b * 1/n * pow(log(n), b-1);
    };
    for (double delta = INFINITY; delta > eps; ) {
        delta = (f(n) - t) / f1(n);
        n -= delta;
    }
    return n;
}
int main() {
    // input
    int m; scanf("%d", &m);
    vector<int> as(m), bs(m);
    vector<double> ts(m);
    repeat (i,m) scanf("%d%d%lf", &as[i], &bs[i], &ts[i]);
    // solve
    array<array<vector<int>, 10+1>, 10+1> xs = {};
    repeat (i,m) xs[as[i]][bs[i]].push_back(i);
    vector<double> n(m);
    repeat (a,10+1) {
        repeat (b,10+1) {
            if (a == 0) {
                for (int x : xs[a][b]) {
                    n[x] = exp(pow(ts[x], 1.0/b));
                }
            } else if (b == 0) {
                for (int x : xs[a][b]) {
                    n[x] = pow(ts[x], 1.0/a);
                }
            } else {
                whole(sort, xs[a][b], [&](int x, int y) {
                    return ts[x] < ts[y];
                });
                double prev = 22027; // exp(10)
                repeat_reverse (i, int(xs[a][b].size())) {
                    int x = xs[a][b][i];
                    double t = ts[x];
                    prev = n[x] = newton(a, b, t, prev);
                }
            }
        }
    }
    // output
    repeat (i,m) {
        printf("%.12lf\n", n[i]);
    }
    return 0;
}