結果

問題 No.620 ぐるぐるぐるりん
コンテスト
ユーザー rickytheta
提出日時 2017-12-20 11:52:03
言語 C++11（廃止可能性あり）
(gcc 13.3.0)
結果
AC  
実行時間 18 ms / 1,000 ms
コード長 4,115 bytes
コンパイル時間 1,867 ms
コンパイル使用メモリ 162,600 KB
実行使用メモリ 6,944 KB
最終ジャッジ日時 2024-06-13 02:03:55
合計ジャッジ時間 4,044 ms
ジャッジサーバーID
（参考情報） 		judge1 / judge5
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
sample AC * 3
other AC * 28
権限があれば一括ダウンロードができます
コンパイルメッセージ
main.cpp: In function ‘void solve()’:
main.cpp:54:8: warning: ignoring return value of ‘int scanf(const char*, ...)’ declared with attribute ‘warn_unused_result’ [-Wunused-result]
   54 |   scanf("%d%Lf%Lf%Lf%Lf%Lf",&t,&p,&omega,&v,&gx,&gy);
      |   ~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
main.cpp: In function ‘int main()’:
main.cpp:123:8: warning: ignoring return value of ‘int scanf(const char*, ...)’ declared with attribute ‘warn_unused_result’ [-Wunused-result]
  123 |   scanf("%d",&n);
      |   ~~~~~^~~~~~~~~

ソースコード

diff # 

#include <bits/stdc++.h>

using namespace std;

typedef int _loop_int;
#define REP(i,n) for(_loop_int i=0;i<(_loop_int)(n);++i)
#define FOR(i,a,b) for(_loop_int i=(_loop_int)(a);i<(_loop_int)(b);++i)
#define FORR(i,a,b) for(_loop_int i=(_loop_int)(b)-1;i>=(_loop_int)(a);--i)

#define DEBUG(x) cout<<#x<<": "<<x<<endl
#define DEBUG_VEC(v) cout<<#v<<":";REP(i,v.size())cout<<" "<<v[i];cout<<endl

// ux, uy を変化させると最後の位置も線形に変化するので、2t個のベクトルを合わせて目的地へのベクトルにする問題になる
// タイミング i で、(ux,uy)=(1,0)とした時の変化ベクトルを(vx_{2i}, vy_{2i})とし、
//             (ux,uy)=(0,1)とした時の変化ベクトルを(vx_{2i+1}, vy_{2i+1})とする
// 何もしなくても動く分を (gx,gy) から引いたものを (ax,ay) とすれば、次の問題を解くことになる
// problem. sum_i (p_i * vx_i) == ax
//          sum_i (p_i * vy_i) == ay
//          sum_i p_i^2    <= p
// p は最適解ギリギリを出してくるだろうから、値を最小化する問題に書き換える
// minimize. sum p_i^2
// subject to. sum_i (p_i * vx_i) == ax
//             sum_i (p_i * vy_i) == ay
// これは二次計画問題で制約が線形な等式のみの場合
// 書き直すと、A を vx_i, vy_i を並べた 2×2t の行列、b を [ax; ay]、I を 2t×2t の単位行列とすれば、
// minimize (p * I * p^T)
// subject to. A * p == b
// となる
// http://www.me.titech.ac.jp/~mizu_lab/text/PDF-NLP/NLP1-QP-problem.pdf
// の4ページとかにあるようにこれは凸2次計画問題で、(5)の式を使って最適解を求められる
// 今回の場合、
// A^T y = x
// Ax = b
// を解くことになるが、代入すると AA^Ty = b なので、これを解いて y を求めて、x = A^T y で解に戻すということをすれば良い
// 必要なのは AA^T の逆行列計算ぐらいだが、これは 2×2 行列のため容易に求まり、他は単純な行列の積で求まる
// 誤差がやばいので long double さんに来てもらいました

#define double long double

int t;
double p;
double omega, v;
double gx, gy;

int m;
double A[2][35];
double b[35];
double AAT[2][2];
double AATinv[2][2];
double yy[2];
double xx[35];

void solve(){
  scanf("%d%Lf%Lf%Lf%Lf%Lf",&t,&p,&omega,&v,&gx,&gy);
  vector<double> vx, vy;
  REP(i,t){
    double x = 1.0, y = 0.0;
    FOR(j,i+1,t){
      double nx = x - omega*y + v*x;
      double ny = y + omega*x + v*y;
      x = nx;
      y = ny;
    }
    vx.push_back(x);
    vy.push_back(y);
    x = 0.0, y = 1.0;
    FOR(j,i+1,t){
      double nx = x - omega*y + v*x;
      double ny = y + omega*x + v*y;
      x = nx;
      y = ny;
    }
    vx.push_back(x);
    vy.push_back(y);
  }
  double sx = 1.0, sy = 0.0;
  REP(i,t){
    double nx = sx - omega*sy + v*sx;
    double ny = sy + omega*sx + v*sy;
    sx = nx;
    sy = ny;
  }
  double ax = gx-sx, ay = gy-sy;

  m = 2*t;
  // Ax = b
  REP(i,m)A[0][i] = vx[i];
  REP(i,m)A[1][i] = vy[i];
  b[0] = ax;
  b[1] = ay;
  // AAT = AA^T
  REP(i,2)REP(j,2)AAT[i][j] = 0.0;
  REP(i,2)REP(j,2)REP(k,m)AAT[i][j] += A[i][k] * A[j][k];
  // AAT y = b
  double detA = AAT[0][0] * AAT[1][1] - AAT[0][1] * AAT[1][0];
  AATinv[0][0] = AAT[1][1] / detA;
  AATinv[0][1] = -AAT[0][1] / detA;
  AATinv[1][0] = -AAT[1][0] / detA;
  AATinv[1][1] = AAT[0][0] / detA;
  // y = AATinv b
  yy[0] = AATinv[0][0] * b[0] + AATinv[0][1] * b[1];
  yy[1] = AATinv[1][0] * b[0] + AATinv[1][1] * b[1];
  // x = A^T y
  REP(i,m)xx[i] = A[0][i] * yy[0] + A[1][i] * yy[1];

  // simulate
  double simx = 1.0, simy = 0.0;
  double power = 0.0;
  REP(i,t){
    double ux = xx[2*i], uy = xx[2*i+1];
    printf("%.40Lf %.40Lf\n", ux, uy);
    power += ux*ux + uy*uy;
    double nx = simx - omega*simy + v*simx + ux;
    double ny = simy + omega*simx + v*simy + uy;
    simx = nx;
    simy = ny;
  }
  // printf("(%.5f, %.5f) is goal, (%.5f, %.5f) is answer, power is %.5f(<= %.5f)\n", gx, gy, simx, simy, power, p);
}

int main(){
  int n;
  scanf("%d",&n);
  while(n--)solve();
  return 0;
}
0