#858188 (C++17) No.2447 行列累乗根

提出ソース
結果

問題	No.2447 行列累乗根
ユーザー	hotman78
提出日時	2023-04-09 21:05:23
言語	C++17 (gcc 12.3.0 + boost 1.83.0)
結果	AC
実行時間	322 ms / 2,000 ms
コード長	19,197 bytes
コンパイル時間	4,589 ms
コンパイル使用メモリ	252,456 KB
実行使用メモリ	5,248 KB
最終ジャッジ日時	2024-11-15 03:50:35
合計ジャッジ時間	12,097 ms
ジャッジサーバーID （参考情報）	judge3 / judge4
このコードへのチャレンジ
（要ログイン）
テストケース

テストケース表示
入力	結果	実行時間実行使用メモリ
testcase_00	AC	2 ms 5,248 KB
testcase_01	AC	2 ms 5,248 KB
testcase_02	AC	2 ms 5,248 KB
testcase_03	AC	2 ms 5,248 KB
testcase_04	AC	2 ms 5,248 KB
testcase_05	AC	2 ms 5,248 KB
testcase_06	AC	2 ms 5,248 KB
testcase_07	AC	2 ms 5,248 KB
testcase_08	AC	2 ms 5,248 KB
testcase_09	AC	2 ms 5,248 KB
testcase_10	AC	2 ms 5,248 KB
testcase_11	AC	2 ms 5,248 KB
testcase_12	AC	2 ms 5,248 KB
testcase_13	AC	2 ms 5,248 KB
testcase_14	AC	2 ms 5,248 KB
testcase_15	AC	2 ms 5,248 KB
testcase_16	AC	2 ms 5,248 KB
testcase_17	AC	2 ms 5,248 KB
testcase_18	AC	6 ms 5,248 KB
testcase_19	AC	33 ms 5,248 KB
testcase_20	AC	320 ms 5,248 KB
testcase_21	AC	307 ms 5,248 KB
testcase_22	AC	309 ms 5,248 KB
testcase_23	AC	311 ms 5,248 KB
testcase_24	AC	317 ms 5,248 KB
testcase_25	AC	322 ms 5,248 KB
testcase_26	AC	316 ms 5,248 KB
testcase_27	AC	321 ms 5,248 KB
権限があれば一括ダウンロードができます
ソースコード

raw source code
// author: hotman78
// date: 2023/04/09-21:05:15

// --- begin raw code -----------------
// #include"cpplib/util/template.hpp"
// #include<atcoder/string.hpp>
// 
// struct matrix{
//     using D=long double;
//     using M=matrix;
//     array<array<D,2>,2>v;
//     matrix(){
//         rep(i,2)v[i].fill(0);
//     }
//     M& operator+=(M x){
//         rep(i,2)rep(j,2)v[i][j]+=x.v[i][j];
//         return *this;
//     }
//     M& operator*=(M x){
//         M res;
//         rep(i,2)rep(j,2)rep(k,2)res.v[i][j]+=v[i][k]*x.v[k][j];
//         return *this=res;
//     }
//     M operator+(M x){
//         return M(*this)+=x;
//     }
//     M operator*(M x){
//         return M(*this)*=x;
//     }
// 
//     static M transpose(M x){
//         M res;
//         rep(i,2)rep(j,2)res.v[i][j]=x.v[j][i];
//         return res;
//     }
// 
//     M cbrt(){
//         if(abs(v[1][0])<=1e-11){
//             M res;
//             res.v={array<D,2>{std::cbrt(v[0][0]),0l},array<D,2>{0l,std::cbrt(v[1][1])}};
//             return res;
//         }
//         D a=1,b=-v[0][0]-v[1][1],c=v[0][0]*v[1][1]-v[1][0]*v[0][1];
//         D d=b*b-4*a*c;
//         D x1=(-b-sqrt(d))/(2*a);
//         D x2=abs(x1)<=1e-11?-b:c/x1;
//         // cerr<<x1<<" "<<x2<<endl;
//         M p;
//         bool fst=0;
//         {
//             bool ok=0;
//             if(!ok){
//                 p.v[0][0]=-v[0][1];
//                 p.v[0][1]=v[0][0]-x1;
//                 D tmp=p.v[0][0]*p.v[0][0]+p.v[0][1]*p.v[0][1];
//                 if(abs(tmp)>=1e-11){
//                     p.v[0][0]/=sqrt(tmp);
//                     p.v[0][1]/=sqrt(tmp);
//                     ok=1;
//                 }
//             }
//             if(!ok){
//                 p.v[0][0]=v[1][1]-x1;
//                 p.v[0][1]=-v[1][0];
//                 D tmp=p.v[0][0]*p.v[0][0]+p.v[0][1]*p.v[0][1];
//                 if(abs(tmp)>=1e-11){
//                     p.v[0][0]/=sqrt(tmp);
//                     p.v[0][1]/=sqrt(tmp);
//                     ok=1;
//                 }
//             }
//             if(!ok){
//                 p.v[0][0]=1;
//                 p.v[0][1]=0;
//             }
//         }
//         p.v[1][0]=p.v[0][1];
//         p.v[1][1]=-p.v[0][0];
//         // rep(i,2){
//         //     rep(j,2){
//         //         cout<<p.v[i][j]<<" \n"[j];
//         //     }
//         // }
//         M m;
//         m.v={array<D,2>{std::cbrt(x1),0l},array<D,2>{0l,std::cbrt(x2)}};
//         // rep(i,2){
//         //     rep(j,2){
//         //         cout<<m.v[i][j]<<" \n"[j];
//         //     }
//         // }
//         // M xx=transpose(p);
//         // rep(i,2){
//         //     rep(j,2){
//         //         cout<<xx.v[i][j]<<" \n"[j];
//         //     }
//         // }
//         return transpose(p)*m*(p);
//     }
// };
// 
// int main(){
//     lint t;
//     cin>>t;
//     while(t--){
//         matrix x;
//         rep(i,2)rep(j,2)cin>>x.v[i][j];
//         rep(i,2)rep(j,2)x.v[i][j]*=1000000;
//         x=x.cbrt();
//         rep(i,2){
//             rep(j,2){
//                 cout<<fixed<<setprecision(5)<<x.v[i][j]*0.01<<" \n"[j];
//             }
//         }
//     }
// }

// --- end raw code -----------------

#line 2 "cpplib/util/template.hpp"
#pragma GCC optimize("Ofast")
#pragma GCC optimize("unroll-loops")
#pragma GCC target("avx2")
#include <bits/stdc++.h>
using namespace std;
#line 1 "cpplib/util/ioutil.hpp"
// template <class Head,class... Args>
// std::ostream& output(std::ostream& out,const Head& head,const Args&... args){
//     out>>head;
//     return output(head,args...);
// }
// template <class Head>
// std::ostream& output(std::ostream& out,const Head& head){
//     out>>head;
//     return out;
// }

template <typename T, typename E>
std::ostream &operator<<(std::ostream &out, std::pair<T, E> v) {
    out << "(" << v.first << "," << v.second << ")";
    return out;
}

// template <class... Args>
// ostream& operator<<(ostream& out,std::tuple<Args...>v){
//     std::apply(output,v);
//     return out;
// }
#line 8 "cpplib/util/template.hpp"
struct __INIT__ {
    __INIT__() {
        cin.tie(0);
        ios::sync_with_stdio(false);
        cout << fixed << setprecision(15);
    }
} __INIT__;
typedef long long lint;
constexpr long long INF = 1LL << 60;
constexpr int IINF = 1 << 30;
constexpr double EPS = 1e-10;
#ifndef REACTIVE
#define endl '\n';
#endif
typedef vector<lint> vec;
typedef vector<vector<lint>> mat;
typedef vector<vector<vector<lint>>> mat3;
typedef vector<string> svec;
typedef vector<vector<string>> smat;
template <typename T> using V = vector<T>;
template <typename T> using VV = V<V<T>>;
template <typename T> inline void output(T t) {
    bool f = 0;
    for (auto i : t) {
        cout << (f ? " " : "") << i;
        f = 1;
    }
    cout << endl;
}
template <typename T> inline void output2(T t) {
    for (auto i : t)
        output(i);
}
template <typename T> inline void debug(T t) {
    bool f = 0;
    for (auto i : t) {
        cerr << (f ? " " : "") << i;
        f = 1;
    }
    cerr << endl;
}
template <typename T> inline void debug2(T t) {
    for (auto i : t)
        debug(i);
}
#define loop(n) for (long long _ = 0; _ < (long long)(n); ++_)
#define _overload4(_1, _2, _3, _4, name, ...) name
#define __rep(i, a) repi(i, 0, a, 1)
#define _rep(i, a, b) repi(i, a, b, 1)
#define repi(i, a, b, c)                                                       \
    for (long long i = (long long)(a); i < (long long)(b); i += c)
#define rep(...) _overload4(__VA_ARGS__, repi, _rep, __rep)(__VA_ARGS__)
#define _overload3_rev(_1, _2, _3, name, ...) name
#define _rep_rev(i, a) repi_rev(i, 0, a)
#define repi_rev(i, a, b)                                                      \
    for (long long i = (long long)(b)-1; i >= (long long)(a); --i)
#define rrep(...) _overload3_rev(__VA_ARGS__, repi_rev, _rep_rev)(__VA_ARGS__)

#define all(n) begin(n), end(n)
template <typename T, typename E> bool chmin(T &s, const E &t) {
    bool res = s > t;
    s = min<T>(s, t);
    return res;
}
template <typename T, typename E> bool chmax(T &s, const E &t) {
    bool res = s < t;
    s = max<T>(s, t);
    return res;
}
const vector<lint> dx = {1, 0, -1, 0, 1, 1, -1, -1};
const vector<lint> dy = {0, 1, 0, -1, 1, -1, 1, -1};
#define SUM(v) accumulate(all(v), 0LL)
#if __cplusplus >= 201703L
template <typename T, typename... Args>
auto make_vector(T x, int arg, Args... args) {
    if constexpr (sizeof...(args) == 0)
        return vector<T>(arg, x);
    else
        return vector(arg, make_vector<T>(x, args...));
}
#endif
#define extrep(v, ...) for (auto v : __MAKE_MAT__({__VA_ARGS__}))
#define bit(n, a) ((n >> a) & 1)
vector<vector<long long>> __MAKE_MAT__(vector<long long> v) {
    if (v.empty())
        return vector<vector<long long>>(1, vector<long long>());
    long long n = v.back();
    v.pop_back();
    vector<vector<long long>> ret;
    vector<vector<long long>> tmp = __MAKE_MAT__(v);
    for (auto e : tmp)
        for (long long i = 0; i < n; ++i) {
            ret.push_back(e);
            ret.back().push_back(i);
        }
    return ret;
}
using graph = vector<vector<int>>;
template <typename T> using graph_w = vector<vector<pair<int, T>>>;

#if __cplusplus >= 201703L
constexpr inline long long powll(long long a, long long b) {
    long long res = 1;
    while (b--)
        res *= a;
    return res;
}
#endif

template <typename T, typename E>
pair<T, E> &operator+=(pair<T, E> &s, const pair<T, E> &t) {
    s.first += t.first;
    s.second += t.second;
    return s;
}
template <typename T, typename E>
pair<T, E> &operator-=(pair<T, E> &s, const pair<T, E> &t) {
    s.first -= t.first;
    s.second -= t.second;
    return s;
}
template <typename T, typename E>
pair<T, E> operator+(const pair<T, E> &s, const pair<T, E> &t) {
    auto res = s;
    return res += t;
}
template <typename T, typename E>
pair<T, E> operator-(const pair<T, E> &s, const pair<T, E> &t) {
    auto res = s;
    return res -= t;
}
#define BEGIN_STACK_EXTEND(size)                                               \
    void *stack_extend_memory_ = malloc(size);                                 \
    void *stack_extend_origin_memory_;                                         \
    char *stack_extend_dummy_memory_ = (char *)alloca(                         \
        (1 + (int)(((long long)stack_extend_memory_) & 127)) * 16);            \
    *stack_extend_dummy_memory_ = 0;                                           \
    asm volatile("mov %%rsp, %%rbx\nmov %%rax, %%rsp"                          \
                 : "=b"(stack_extend_origin_memory_)                           \
                 : "a"((char *)stack_extend_memory_ + (size)-1024));
#define END_STACK_EXTEND                                                       \
    asm volatile("mov %%rax, %%rsp" ::"a"(stack_extend_origin_memory_));       \
    free(stack_extend_memory_);
#line 2 "main.cpp"

#include <algorithm>
#include <cassert>
#include <numeric>
#include <string>
#include <vector>

namespace atcoder {

namespace internal {

std::vector<int> sa_naive(const std::vector<int> &s) {
    int n = int(s.size());
    std::vector<int> sa(n);
    std::iota(sa.begin(), sa.end(), 0);
    std::sort(sa.begin(), sa.end(), [&](int l, int r) {
        if (l == r)
            return false;
        while (l < n && r < n) {
            if (s[l] != s[r])
                return s[l] < s[r];
            l++;
            r++;
        }
        return l == n;
    });
    return sa;
}

std::vector<int> sa_doubling(const std::vector<int> &s) {
    int n = int(s.size());
    std::vector<int> sa(n), rnk = s, tmp(n);
    std::iota(sa.begin(), sa.end(), 0);
    for (int k = 1; k < n; k *= 2) {
        auto cmp = [&](int x, int y) {
            if (rnk[x] != rnk[y])
                return rnk[x] < rnk[y];
            int rx = x + k < n ? rnk[x + k] : -1;
            int ry = y + k < n ? rnk[y + k] : -1;
            return rx < ry;
        };
        std::sort(sa.begin(), sa.end(), cmp);
        tmp[sa[0]] = 0;
        for (int i = 1; i < n; i++) {
            tmp[sa[i]] = tmp[sa[i - 1]] + (cmp(sa[i - 1], sa[i]) ? 1 : 0);
        }
        std::swap(tmp, rnk);
    }
    return sa;
}

template <int THRESHOLD_NAIVE = 10, int THRESHOLD_DOUBLING = 40>
std::vector<int> sa_is(const std::vector<int> &s, int upper) {
    int n = int(s.size());
    if (n == 0)
        return {};
    if (n == 1)
        return {0};
    if (n == 2) {
        if (s[0] < s[1]) {
            return {0, 1};
        } else {
            return {1, 0};
        }
    }
    if (n < THRESHOLD_NAIVE) {
        return sa_naive(s);
    }
    if (n < THRESHOLD_DOUBLING) {
        return sa_doubling(s);
    }

    std::vector<int> sa(n);
    std::vector<bool> ls(n);
    for (int i = n - 2; i >= 0; i--) {
        ls[i] = (s[i] == s[i + 1]) ? ls[i + 1] : (s[i] < s[i + 1]);
    }
    std::vector<int> sum_l(upper + 1), sum_s(upper + 1);
    for (int i = 0; i < n; i++) {
        if (!ls[i]) {
            sum_s[s[i]]++;
        } else {
            sum_l[s[i] + 1]++;
        }
    }
    for (int i = 0; i <= upper; i++) {
        sum_s[i] += sum_l[i];
        if (i < upper)
            sum_l[i + 1] += sum_s[i];
    }

    auto induce = [&](const std::vector<int> &lms) {
        std::fill(sa.begin(), sa.end(), -1);
        std::vector<int> buf(upper + 1);
        std::copy(sum_s.begin(), sum_s.end(), buf.begin());
        for (auto d : lms) {
            if (d == n)
                continue;
            sa[buf[s[d]]++] = d;
        }
        std::copy(sum_l.begin(), sum_l.end(), buf.begin());
        sa[buf[s[n - 1]]++] = n - 1;
        for (int i = 0; i < n; i++) {
            int v = sa[i];
            if (v >= 1 && !ls[v - 1]) {
                sa[buf[s[v - 1]]++] = v - 1;
            }
        }
        std::copy(sum_l.begin(), sum_l.end(), buf.begin());
        for (int i = n - 1; i >= 0; i--) {
            int v = sa[i];
            if (v >= 1 && ls[v - 1]) {
                sa[--buf[s[v - 1] + 1]] = v - 1;
            }
        }
    };

    std::vector<int> lms_map(n + 1, -1);
    int m = 0;
    for (int i = 1; i < n; i++) {
        if (!ls[i - 1] && ls[i]) {
            lms_map[i] = m++;
        }
    }
    std::vector<int> lms;
    lms.reserve(m);
    for (int i = 1; i < n; i++) {
        if (!ls[i - 1] && ls[i]) {
            lms.push_back(i);
        }
    }

    induce(lms);

    if (m) {
        std::vector<int> sorted_lms;
        sorted_lms.reserve(m);
        for (int v : sa) {
            if (lms_map[v] != -1)
                sorted_lms.push_back(v);
        }
        std::vector<int> rec_s(m);
        int rec_upper = 0;
        rec_s[lms_map[sorted_lms[0]]] = 0;
        for (int i = 1; i < m; i++) {
            int l = sorted_lms[i - 1], r = sorted_lms[i];
            int end_l = (lms_map[l] + 1 < m) ? lms[lms_map[l] + 1] : n;
            int end_r = (lms_map[r] + 1 < m) ? lms[lms_map[r] + 1] : n;
            bool same = true;
            if (end_l - l != end_r - r) {
                same = false;
            } else {
                while (l < end_l) {
                    if (s[l] != s[r]) {
                        break;
                    }
                    l++;
                    r++;
                }
                if (l == n || s[l] != s[r])
                    same = false;
            }
            if (!same)
                rec_upper++;
            rec_s[lms_map[sorted_lms[i]]] = rec_upper;
        }

        auto rec_sa =
            sa_is<THRESHOLD_NAIVE, THRESHOLD_DOUBLING>(rec_s, rec_upper);

        for (int i = 0; i < m; i++) {
            sorted_lms[i] = lms[rec_sa[i]];
        }
        induce(sorted_lms);
    }
    return sa;
}

} // namespace internal

std::vector<int> suffix_array(const std::vector<int> &s, int upper) {
    assert(0 <= upper);
    for (int d : s) {
        assert(0 <= d && d <= upper);
    }
    auto sa = internal::sa_is(s, upper);
    return sa;
}

template <class T> std::vector<int> suffix_array(const std::vector<T> &s) {
    int n = int(s.size());
    std::vector<int> idx(n);
    iota(idx.begin(), idx.end(), 0);
    sort(idx.begin(), idx.end(), [&](int l, int r) { return s[l] < s[r]; });
    std::vector<int> s2(n);
    int now = 0;
    for (int i = 0; i < n; i++) {
        if (i && s[idx[i - 1]] != s[idx[i]])
            now++;
        s2[idx[i]] = now;
    }
    return internal::sa_is(s2, now);
}

std::vector<int> suffix_array(const std::string &s) {
    int n = int(s.size());
    std::vector<int> s2(n);
    for (int i = 0; i < n; i++) {
        s2[i] = s[i];
    }
    return internal::sa_is(s2, 255);
}

template <class T>
std::vector<int> lcp_array(const std::vector<T> &s,
                           const std::vector<int> &sa) {
    int n = int(s.size());
    assert(n >= 1);
    std::vector<int> rnk(n);
    for (int i = 0; i < n; i++) {
        rnk[sa[i]] = i;
    }
    std::vector<int> lcp(n - 1);
    int h = 0;
    for (int i = 0; i < n; i++) {
        if (h > 0)
            h--;
        if (rnk[i] == 0)
            continue;
        int j = sa[rnk[i] - 1];
        for (; j + h < n && i + h < n; h++) {
            if (s[j + h] != s[i + h])
                break;
        }
        lcp[rnk[i] - 1] = h;
    }
    return lcp;
}

std::vector<int> lcp_array(const std::string &s, const std::vector<int> &sa) {
    int n = int(s.size());
    std::vector<int> s2(n);
    for (int i = 0; i < n; i++) {
        s2[i] = s[i];
    }
    return lcp_array(s2, sa);
}

template <class T> std::vector<int> z_algorithm(const std::vector<T> &s) {
    int n = int(s.size());
    if (n == 0)
        return {};
    std::vector<int> z(n);
    z[0] = 0;
    for (int i = 1, j = 0; i < n; i++) {
        int &k = z[i];
        k = (j + z[j] <= i) ? 0 : std::min(j + z[j] - i, z[i - j]);
        while (i + k < n && s[k] == s[i + k])
            k++;
        if (j + z[j] < i + z[i])
            j = i;
    }
    z[0] = n;
    return z;
}

std::vector<int> z_algorithm(const std::string &s) {
    int n = int(s.size());
    std::vector<int> s2(n);
    for (int i = 0; i < n; i++) {
        s2[i] = s[i];
    }
    return z_algorithm(s2);
}

} // namespace atcoder

struct matrix {
    using D = long double;
    using M = matrix;
    array<array<D, 2>, 2> v;
    matrix() { rep(i, 2) v[i].fill(0); }
    M &operator+=(M x) {
        rep(i, 2) rep(j, 2) v[i][j] += x.v[i][j];
        return *this;
    }
    M &operator*=(M x) {
        M res;
        rep(i, 2) rep(j, 2) rep(k, 2) res.v[i][j] += v[i][k] * x.v[k][j];
        return *this = res;
    }
    M operator+(M x) { return M(*this) += x; }
    M operator*(M x) { return M(*this) *= x; }

    static M transpose(M x) {
        M res;
        rep(i, 2) rep(j, 2) res.v[i][j] = x.v[j][i];
        return res;
    }

    M cbrt() {
        if (abs(v[1][0]) <= 1e-11) {
            M res;
            res.v = {array<D, 2>{std::cbrt(v[0][0]), 0l},
                     array<D, 2>{0l, std::cbrt(v[1][1])}};
            return res;
        }
        D a = 1, b = -v[0][0] - v[1][1],
          c = v[0][0] * v[1][1] - v[1][0] * v[0][1];
        D d = b * b - 4 * a * c;
        D x1 = (-b - sqrt(d)) / (2 * a);
        D x2 = abs(x1) <= 1e-11 ? -b : c / x1;
        // cerr<<x1<<" "<<x2<<endl;
        M p;
        bool fst = 0;
        {
            bool ok = 0;
            if (!ok) {
                p.v[0][0] = -v[0][1];
                p.v[0][1] = v[0][0] - x1;
                D tmp = p.v[0][0] * p.v[0][0] + p.v[0][1] * p.v[0][1];
                if (abs(tmp) >= 1e-11) {
                    p.v[0][0] /= sqrt(tmp);
                    p.v[0][1] /= sqrt(tmp);
                    ok = 1;
                }
            }
            if (!ok) {
                p.v[0][0] = v[1][1] - x1;
                p.v[0][1] = -v[1][0];
                D tmp = p.v[0][0] * p.v[0][0] + p.v[0][1] * p.v[0][1];
                if (abs(tmp) >= 1e-11) {
                    p.v[0][0] /= sqrt(tmp);
                    p.v[0][1] /= sqrt(tmp);
                    ok = 1;
                }
            }
            if (!ok) {
                p.v[0][0] = 1;
                p.v[0][1] = 0;
            }
        }
        p.v[1][0] = p.v[0][1];
        p.v[1][1] = -p.v[0][0];
        // rep(i,2){
        //     rep(j,2){
        //         cout<<p.v[i][j]<<" \n"[j];
        //     }
        // }
        M m;
        m.v = {array<D, 2>{std::cbrt(x1), 0l}, array<D, 2>{0l, std::cbrt(x2)}};
        // rep(i,2){
        //     rep(j,2){
        //         cout<<m.v[i][j]<<" \n"[j];
        //     }
        // }
        // M xx=transpose(p);
        // rep(i,2){
        //     rep(j,2){
        //         cout<<xx.v[i][j]<<" \n"[j];
        //     }
        // }
        return transpose(p) * m * (p);
    }
};

int main() {
    lint t;
    cin >> t;
    while (t--) {
        matrix x;
        rep(i, 2) rep(j, 2) cin >> x.v[i][j];
        rep(i, 2) rep(j, 2) x.v[i][j] *= 1000000;
        x = x.cbrt();
        rep(i, 2) {
            rep(j, 2) {
                cout << fixed << setprecision(5) << x.v[i][j] * 0.01
                     << " \n"[j];
            }
        }
    }
}
yukicoder

結果

テストケース

ソースコード
