結果
| 問題 | No.2447 行列累乗根 |
| ユーザー |
 hotman78
|
| 提出日時 | 2023-04-09 21:53:20 |
| 言語 | C++17 (gcc 12.3.0 + boost 1.83.0) |
| 結果 |
AC
|
| 実行時間 | 321 ms / 2,000 ms |
| コード長 | 19,193 bytes |
| コンパイル時間 | 4,561 ms |
| コンパイル使用メモリ | 251,840 KB |
| 実行使用メモリ | 6,824 KB |
| 最終ジャッジ日時 | 2024-11-15 03:51:43 |
| 合計ジャッジ時間 | 12,065 ms |
|
ジャッジサーバーID （参考情報） |
judge2 / judge3 |
（要ログイン）
テストケース
テストケース表示| 入力 | 結果 | 実行時間 実行使用メモリ |
|---|---|---|
| testcase_00 | AC | 2 ms
6,816 KB |
| testcase_01 | AC | 2 ms
6,816 KB |
| testcase_02 | AC | 2 ms
6,820 KB |
| testcase_03 | AC | 2 ms
6,824 KB |
| testcase_04 | AC | 3 ms
6,820 KB |
| testcase_05 | AC | 2 ms
6,820 KB |
| testcase_06 | AC | 2 ms
6,816 KB |
| testcase_07 | AC | 2 ms
6,816 KB |
| testcase_08 | AC | 2 ms
6,816 KB |
| testcase_09 | AC | 3 ms
6,816 KB |
| testcase_10 | AC | 2 ms
6,816 KB |
| testcase_11 | AC | 2 ms
6,820 KB |
| testcase_12 | AC | 2 ms
6,816 KB |
| testcase_13 | AC | 2 ms
6,820 KB |
| testcase_14 | AC | 2 ms
6,816 KB |
| testcase_15 | AC | 3 ms
6,816 KB |
| testcase_16 | AC | 2 ms
6,816 KB |
| testcase_17 | AC | 2 ms
6,816 KB |
| testcase_18 | AC | 5 ms
6,820 KB |
| testcase_19 | AC | 35 ms
6,816 KB |
| testcase_20 | AC | 313 ms
6,816 KB |
| testcase_21 | AC | 307 ms
6,820 KB |
| testcase_22 | AC | 309 ms
6,816 KB |
| testcase_23 | AC | 312 ms
6,816 KB |
| testcase_24 | AC | 297 ms
6,820 KB |
| testcase_25 | AC | 302 ms
6,820 KB |
| testcase_26 | AC | 298 ms
6,820 KB |
| testcase_27 | AC | 321 ms
6,820 KB |
ソースコード
// author: hotman78
// date: 2023/04/09-21:52:50
// --- 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])<=100){
// 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]) <= 100) {
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];
}
}
}
}