用jacob迭代解方程

来源：互联网发布：中国全面战争动员知乎编辑：程序博客网时间：2024/05/22 08:05

昨天晚上在学校OJ上面做题

http://acm.sjtu.edu.cn/OnlineJudge/problem/1139

题意简洁明了做法显然是高斯消元，因为排列可以看成个置换，把每个循环节的大小抽出来记成状态就可以了

这样状态总数是partitions(11) = 56个，可以接受

直接高斯消元WA了，据说有精度问题在群上咨询叉姐后得知要用jacob迭代

...照样WA了..

收获是学习了一种新的解方程的思路，而且还是挺优美的

首先戳进维基http://zh.wikipedia.org/wiki/%E9%9B%85%E5%8F%AF%E6%AF%94%E6%B3%95

可以清晰地看到这个迭代过程

当然.. 这个算法直接做到收敛是坑爹的复杂度奇高无比

观察那个迭代的式子，发现可以写成一个矩阵的形式（代码里面的jacob矩阵），然后就可以快速迭代若干次

因为它是收敛的，所以在次数很大的情况下误差一般就可以接受了

这个复杂度还是有点高的 .. 比普通的高斯消元多上一个快速幂的log

（求教怎么贴代码不吃tab..）

#include <stdio.h>
#include <stdlib.h>
#include <algorithm>
#include <iostream>
#include <vector>
#include <map>
#define eps 1e-9
#define ld long double

using namespace std;

struct Matrix
{
int n;
ld v[105][105];

Matrix () {}

Matrix (int _n)
{
n = _n;
for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++)
v[i][j] = 0;
}

friend Matrix operator * (Matrix &a, Matrix &b)
{
int i, j, k;
Matrix ret = Matrix(a.n);
for (i = 1; i <= a.n; i++)
for (j = 1; j <= a.n; j++)
for (k = 1; k <= a.n; k++)
ret.v[i][j] += a.v[i][k] * b.v[k][j];
return ret;
}
}jacob, j1, j2;

vector <int> Now, Vec[105];
vector <int>::iterator it;
map <vector <int>, int> Number;
int n, q, tot, i, j, k, l1, l2, l3, b[105], c[105], Hash[105];
ld a[105][105], ans[105], ans2[105], temp, gl;

void dfs(int i, int min, int s);
int Get(int a[]);
Matrix sumi(Matrix a, int b);

int main()
{
freopen("input.txt", "r", stdin); freopen("output.txt", "w", stdout);
scanf("%d %d", &n, &q);
dfs(1, 1, 0);
gl = 1 / (double)(n * (n - 1) * (n - 2));
for (i = 1; i <= tot; i++)
{
j = 0;
a[i][i] -= 1.0;
a[i][tot + 1] = -1.0;
for (it = Vec[i].begin(); it != Vec[i].end(); ++it)
{
for (k = 1; k < (*it); k++)
b[j + k] = j + k + 1;
b[j + k] = j + 1;
j += k;
}

for (l1 = 1; l1 < n - 1; l1++)
for (l2 = l1 + 1; l2 < n; l2++)
for (l3 = l2 + 1; l3 <= n; l3++)
{
for (j = 1; j <= n; j++)
c[j] = b[j];
a[i][Get(c)] += gl;

for (j = 1; j <= n; j++)
c[j] = b[j];
swap(c[l2], c[l3]);
a[i][Get(c)] += gl;

for (j = 1; j <= n; j++)
c[j] = b[j];
swap(c[l1], c[l2]);
a[i][Get(c)] += gl;

for (j = 1; j <= n; j++)
c[j] = b[j];
swap(c[l1], c[l2]), swap(c[l2], c[l3]);
a[i][Get(c)] += gl;

for (j = 1; j <= n; j++)
c[j] = b[j];
swap(c[l1], c[l2]), swap(c[l1], c[l3]);
a[i][Get(c)] += gl;

for (j = 1; j <= n; j++)
c[j] = b[j];
swap(c[l1], c[l3]);
a[i][Get(c)] += gl;
}
}

for (i = 1; i <= tot + 1; i++)
a[1][i] = 0.0;
a[1][1] = 1.0;

// jacob
j1 = j2 = Matrix(tot + 1);
for (i = 1; i <= tot; i++)
{
for (j = 1; j <= tot; j++)
if (i != j)
j1.v[j][i] = -a[i][j];
j1.v[j][i] = a[i][tot + 1];
}

for (i = 1; i <= tot; i++)
j2.v[i][i] = 1.0 / a[i][i];
j2.v[i][i] = j1.v[i][i] = 1.0;
jacob = j1 * j2;
jacob = sumi(jacob, 2000000000);
jacob = sumi(jacob, 2000000000);
for (i = 1; i <= tot; i++)
for (j = 1; j <= tot + 1; j++)
ans[i] += jacob.v[j][i];

for (i = 1; i <= q; i++)
{
for (j = 1; j <= n; j++)
scanf("%d", &c[j]);
int ret = Get(c);
cout << (long long)(ans[Get(c)] + eps) << endl;
}

return 0;
}

void dfs(int i, int min, int s)
{
int x;
if (s == n)
{
Now.clear();
for (x = 1; x < i; x++)
Now.push_back(b[x]);
Number[Now] = ++tot;
Vec[tot] = Now;
return;
}

for (x = min; s + x <= n; x++)
b[i] = x, dfs(i + 1, x, s + x);
}

int Get(int a[])
{
int i, j, s;
Now.clear();
for (i = 1; i <= n; i++)
Hash[i] = 0;
for (i = 1; i <= n; i++)
{
if (Hash[i] == 1)
continue;
s = 1, j = a[i], Hash[i] = 1;
while (j != i)
Hash[j] = 1, s++, j = a[j];
Now.push_back(s);
}

sort(Now.begin(), Now.end());
return Number[Now];
}

Matrix sumi(Matrix a, int b)
{
Matrix ret, c = a;
int tag = 0;
while (b > 0)
{
if (b & 1)
{
if (tag == 0)
tag = 1, ret = c;
else
ret = ret * c;
}

c = c * c;
b >>= 1;
}

return ret;
}