JMK no matter what

오늘부터 spoj 에서 하루에 하나씩 문제를 풀기로 맘먹었다. 문제 링크

원 출처는 CERC 2000. f() 가 polynomial 이고, f(1), f(2), .. f(n) 이 주어질 때 f(n+1), f(n+2), ... f(n+c) 까지를 계산하는 문제. f() 는 물론 최소 차수를 갖도록 한다.

이것.. A_ij = (j+1)^i 인 행렬 (i, j 는 0부터) 을 만들고 f(x) 를 포함하는 열벡터를 P 라고 하면 우리가 원하는 차수는 A^-1 * P 가 된다. 이걸로 풀어야겠구나 하고 30분 꼬박 걸려서 행렬 라이브러리를 만들고 디버깅했으나 TLE. 생각해보니 n차 다항식 수열의 계차수열은 n-1차... 5분만에 재귀호출로 풀었다. 주여. 소스코드

lang:cpp
#include<iostream>
#include<cstring>
#include<algorithm>
#include<sstream>
#include<string>
#include<vector>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<fstream>
#include<cassert>
#include<numeric>
#include<set>
#include<map>
#include<queue>
#include<list>
#include<deque>
using namespace std;

#define FOR(i,a,b) for(int i = (a); i < (b); ++i)
#define REP(i,n) FOR(i,0,n)
#define FORE(it,x) for(typeof(x.begin()) it=x.begin();it!=x.end();++it)
#define pb push_back
#define all(x) (x).begin(),(x).end()
#define CLEAR(x,with) memset(x,with,sizeof(x))
#define sz size()
typedef vector<int> VI;
typedef vector<VI> VVI;
typedef long long ll;

VI complete(const VI& seq, int more)
{
    if(seq.sz == 1) return VI(more, seq[0]);
    if(count(all(seq), seq[0]) == seq.sz) return VI(more, seq[0]);
    VI diff(seq.sz-1);
    REP(i,diff.sz) diff[i] = seq[i+1] - seq[i];
    VI diff2 = complete(diff, more);
    VI ret(1, seq.back() + diff2[0]);
    FOR(i,1,more)
        ret.pb(ret.back() + diff2[i]);
    return ret;
}

int main()
{
    int cases;
    scanf("%d", &cases);
    for(int cc = 0; cc < cases; ++cc)
    {
        int n, c;
        scanf("%d %d", &n, &c);
        VI given(n);
        REP(i,n) scanf("%d", &given[i]);
        VI res = complete(given, c);
        REP(i,c)
        {
            if(i) printf(" ");
            printf("%d", res[i]);
        }
        printf("\n");

    }
}

불쌍하니 행렬 쓰는 소스 코드도... 사실 행렬에 너무 큰 숫자들이 들어가서 numerical stability 때문에 안될걸 알면서도 내 보았다..;;

lang:cpp
#include<iomanip>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<sstream>
#include<string>
#include<vector>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<fstream>
#include<cassert>
#include<numeric>
#include<set>
#include<map>
#include<queue>
#include<list>
#include<deque>
using namespace std;

#define FOR(i,a,b) for(int i = (a); i < (b); ++i)
#define REP(i,n) FOR(i,0,n)
#define FORE(it,x) for(typeof(x.begin()) it=x.begin();it!=x.end();++it)
#define pb push_back
#define all(x) (x).begin(),(x).end()
#define CLEAR(x,with) memset(x,with,sizeof(x))
#define sz size()
typedef vector<int> VI;
typedef vector<VI> VVI;
typedef long long ll;

typedef vector<double> row;
typedef vector<row> matrix;

row operator * (const row& a, double b)
{
    row c(a.size()); for(int i = 0; i < a.size(); ++i) c[i] = a[i] * b;
    return c;
}

row operator + (const row& a, const row& b)
{
    row c(a.size()); for(int i = 0; i < a.size(); ++i) c[i] = a[i] + b[i];
    return c;
}

matrix operator * (const matrix& a, const matrix& b)
{
    matrix c(a.size(), row(b[0].size(), 0));
    for(int i = 0; i < a.size(); ++i)
        for(int j = 0; j < b[0].size(); ++j)
            for(int k = 0; k < a[0].size(); ++k)
                c[i][j] += a[i][k] * b[k][j];
    return c;
}

matrix identity(int s)
{
    matrix r(s, row(s, 0));
    for(int i = 0; i < s; ++i) r[i][i] = 1;
    return r;
}

ostream& operator << (ostream& os, const matrix& m)
{
    os << setprecision(2);
    for(int i = 0; i < m.size(); ++i)
    {
        for(int j = 0; j < m[i].size(); ++j)
        {
            os.width(8);
            os<< m[i][j];
        }
        os << endl;
    }
    return os;
}

matrix inv(matrix a)
{
    int n = a.size();
    matrix b = identity(n);

    for(int j = 0; j < n; ++j)
    {
        int maxRow = j;
        for(int i = j + 1; i < n; ++i)
            if(fabs(a[i][j]) > fabs(a[maxRow][j]))
                maxRow = j;
        if(fabs(a[maxRow][j]) < 1e-5) throw 1;
        if(maxRow != j)
        {
            swap(a[j], a[maxRow]);
            swap(b[j], b[maxRow]);
        }
        double div = a[j][j];
        a[j] = a[j] * (1.0 / div); b[j] = b[j] * (1.0 / div);
        for(int i = 0; i < n; ++i)
            if(i != j)
            {
                b[i] = b[i] + (b[j] * -a[i][j]);
                a[i] = a[i] + (a[j] * -a[i][j]);
            }
    }
    return b;
}

double eval(const matrix& coef, int ni)
{
    double v = 0;
    for(int j = coef.size()-1; j >= 0; --j)
        v = v * ni + coef[j][0];
    return v;
}

bool poly(const matrix& P, int deg, matrix& ret)
{
    // deg 0: constant function, 1 numbers
    // deg 1: linear function, 2 numbers
    // use first n numbers.
    int n = deg + 1;
        matrix Q(P.begin(), P.begin() + n);
    matrix A(n, row(n));
    // A[i][j] = (i+1)^j
    for(int i = 0; i < n; ++i)
        A[0][i] = A[i][0] = 1;
    for(int i = 1; i < n; ++i)
        for(int j = 1; j < n; ++j)
            A[i][j] = A[i][j-1] * (i+1);
    //cout << "Trying for deg = " << deg << endl;
    //cout << A << endl;
    try
    {
        // A*D=Q => D=A^-1 * Q
        matrix invA = inv(A);
        ret = invA * Q;
    }
    catch(int)
    {
        return false;
    }

    for(int i = n+1; i <= P.size(); ++i)
    {
        if(fabs(eval(ret, i) - P[i-1][0]) > 1e-7)
            return false;
    }
    return true;
}

int main()
{
    int cases;
    scanf("%d", &cases);
    for(int cc = 0; cc < cases; ++cc)
    {
        int n, c;
        scanf("%d %d", &n, &c);
        matrix P(n, row(1));
        for(int i = 0; i < n; ++i)
            scanf("%lf", &P[i][0]);

        // check for degree -1

        int deg = 0;
        matrix coef;
        while(!poly(P, deg, coef))
            ++deg;
        for(int i = 0; i < c; ++i)
        {
            int ni = n + i + 1;
            if(i) printf(" ");
            printf("%.0lf", eval(coef, ni));
        }
        printf("\n");
    }
}