写了几个高斯消元的题,试一下模板
hdu3976
题意:给定50个节点,2000个电阻,求1到n的等效电阻
这题要用 节点电位法。
在1和n之间加一个1A的电流源,将n接地,这样求出的1号点电位数值上等于等效电阻。
#include <cstdio>
#include <iostream>
typedef long long ll;
using namespace std;
void in(int &x){
char ch; int minus = 0;
while (ch=getchar(), (ch<'0'||ch>'9') && ch!='-');
if (ch == '-') minus = 1, x = 0;
else x = ch-'0';
while (ch=getchar(), ch>='0'&&ch<='9') x = x*10+ch-'0';
if (minus) x = -x;
}
void in(ll &x){
char ch; int minus = 0;
while (ch=getchar(), (ch<'0'||ch>'9') && ch!='-');
if (ch == '-') minus = 1, x = 0;
else x = ch-'0';
while (ch=getchar(), ch>='0'&&ch<='9') x = x*10+ch-'0';
if (minus) x = -x;
}
void out(int x){
char hc[30];int len, minus=0;
if (x<0) minus = 1, x = -x;
len = 0; hc[len++] = x%10+'0';
while (x/=10) hc[len++] = x%10+'0';
if (minus) putchar('-');
for (int i=len-1; i>=0; i--) putchar(hc[i]);
}
void out(ll x){
char hc[30];int len, minus=0;
if (x<0) minus = 1, x = -x;
len = 0; hc[len++] = x%10+'0';
while (x/=10) hc[len++] = x%10+'0';
if (minus) putchar('-');
for (int i=len-1; i>=0; i--) putchar(hc[i]);
}
#include <cstring>
#include <cmath>
const int N = 50+5;
const double eps = 1e-6;
double g[N][N];
int n, m;
void solve(int n, double g[][N]){
int now = 0;
for (int c=0; c<n; c++){
int k = -1;
for (int i=now; i<n; i++) if (fabs(g[i][c])>eps)
if (k==-1 || fabs(g[i][c])>fabs(g[k][c])) k = i;
if (k == -1) continue;
if (k != now)
for (int j=c; j<=n; j++)
swap(g[now][j], g[k][j]);
for (int i=0; i<n; i++) if (i!=now && fabs(g[i][c])>eps){
double tmp = g[i][c]/g[now][c];
for (int j=c; j<=n; j++)
g[i][j] -= tmp*g[now][j];
}
now ++;
}
for (int i=0; i<n; i++) if (fabs(g[i][i])>eps)
g[i][n] /= g[i][i];
}
int main(){
#ifndef ONLINE_JUDGE
freopen("input.txt", "r", stdin);
#endif
int test;
in(test);
for (int T=1; T<=test; T++){
in(n); in(m);
memset(g, 0, sizeof(g));
int a, b, c;
for (int i=1; i<=m; i++){
in(a); in(b); in(c); if (!c) continue;
a--; b--; double tmp = 1.0/c;
g[a][a] += tmp;
g[b][b] += tmp;
g[a][b] -= tmp;
g[b][a] -= tmp;
}
for (int i=0; i<n; i++) g[i][n] = 0.0;
g[0][n] = 1.0; g[n-1][n] = -1.0;
solve(n, g);
printf("Case #%d: %.2lf\n", T, g[0][n]-g[n-1][n]);
}
return 0;
}
hdu5006
题意不变,求s到t之间的等效电阻,但是最多有10000个点,电阻的值只能随机取0或1
两节点之间电阻如果为0说明可以看做一个点,所以实际上点的数量非常少。
另外如果s和t在电路中不联通,则方程无解。
#include <cstdio>
#include <iostream>
typedef long long ll;
using namespace std;
void in(int &x){
char ch; int minus = 0;
while (ch=getchar(), (ch<'0'||ch>'9') && ch!='-');
if (ch == '-') minus = 1, x = 0;
else x = ch-'0';
while (ch=getchar(), ch>='0'&&ch<='9') x = x*10+ch-'0';
if (minus) x = -x;
}
void in(ll &x){
char ch; int minus = 0;
while (ch=getchar(), (ch<'0'||ch>'9') && ch!='-');
if (ch == '-') minus = 1, x = 0;
else x = ch-'0';
while (ch=getchar(), ch>='0'&&ch<='9') x = x*10+ch-'0';
if (minus) x = -x;
}
void out(int x){
char hc[30];int len, minus=0;
if (x<0) minus = 1, x = -x;
len = 0; hc[len++] = x%10+'0';
while (x/=10) hc[len++] = x%10+'0';
if (minus) putchar('-');
for (int i=len-1; i>=0; i--) putchar(hc[i]);
}
void out(ll x){
char hc[30];int len, minus=0;
if (x<0) minus = 1, x = -x;
len = 0; hc[len++] = x%10+'0';
while (x/=10) hc[len++] = x%10+'0';
if (minus) putchar('-');
for (int i=len-1; i>=0; i--) putchar(hc[i]);
}
#include <cstring>
#include <cmath>
#include <cassert>
const int N = 2500;
const int M = 40000+10;
const double eps = 1e-7;
int n, m, st, ed;
double g[N][N], ans[N];
int fa[N*4+10], u[M], v[M], c[M], ms[M/4];
int find(int x){
if (fa[x] == x) return x;
return fa[x] = find(fa[x]);
}
bool solve(int n, double g[][N], double ans[]){
int now = 0;
for (int c=0; c<n; c++){
int r = now;
for (int i=now; i<n; i++) if (fabs(g[i][c]) > fabs(g[r][c]))
r = i;
if (r != now)
for (int j=c; j<=n; j++) swap(g[r][j], g[now][j]);
if (fabs(g[now][c]) < eps) continue;
for (int i=0; i<n; i++) if (i != now && fabs(g[i][c])>eps){
double tmp = g[i][c]/g[now][c];
for (int j=c; j<=n; j++)
g[i][j] -= tmp*g[now][j];
}
now ++;
}
for (int j=0; j<n; j++)
for (int i=0; i<n; i++) if (fabs(g[i][j]) > eps){
ans[j] = g[i][n]/g[i][j];
break;
}
for (int i=0; i<n; i++) if (fabs(g[i][n])>eps){
int flag = 0;
for (int j=i; j<n; j++) if (fabs(g[i][j])>eps){
flag = 1; break;
}
if (!flag) return false;
}
return true;
}
int main(){
#ifndef ONLINE_JUDGE
freopen("input.txt", "r", stdin);
#endif
int test, start, finish;
in(test);
while (test--){
in(n); in(m); in(start); in(finish);
memset(ms, 0xff, sizeof(ms));
memset(g, 0, sizeof(g));
for (int i=1; i<=n; i++) fa[i] = i;
for (int i=1; i<=m; i++){
in(u[i]); in(v[i]); in(c[i]);
if (!c[i]){
int l = find(u[i]), r = find(v[i]);
if (l!=r) fa[l] = r;
}
}
int cnt = 0;
for (int i=1; i<=n; i++)
if (find(i) == i) ms[i] = cnt++;
else ms[i] = -1;
for (int i=1; i<=m; i++) if (c[i]){
int l = find(u[i]), r = find(v[i]);
if (l != r){
l = ms[l]; r = ms[r];
g[l][l] += 1; g[r][r] += 1;
g[l][r] -= 1; g[r][l] -= 1;
}
}
st = ms[find(start)]; ed = ms[find(finish)];
for (int i=0; i<cnt; i++) g[i][ed] = 0;
g[st][cnt] = 1; g[ed][cnt] = -1;
if (st == ed) printf("%.6lf\n", 0.0);//important
else
if (solve(cnt, g, ans)) printf("%.6lf\n", ans[st]);
else puts("inf");
}
return 0;
}
hdu4870
概率dp,50分看成1,1000分就成了20
dp[i][j] 表示一个账号为i 另一个为j 时的期望数
转移方程 dp[i][j] = p*dp[i+1][j] + (1-p)*dp[i-2][j] + 1
方程既跟后面的状态有关,又跟前面的状态有关,所以不能直接递推
这种题就要用高斯消元。
ps:写这题的时候有两个地方卡了很久
1、我用 q 代替 1-p ,导致了精度损失。
2、原来写高斯消元我会取绝对值最大的主元进行消元,据说可以减少精度损失。
但由于未知原因,这一题加了这个优化反而会损失精度。
#include <cstdio>
#include <cstring>
#include <cmath>
#include <cassert>
#include <algorithm>
using namespace std;
const double eps = 1e-9;
const int LIM = 210;
const int L = 20;
int n;
double p, q;
double g[LIM+5][LIM+5];//, gg[LIM+5][LIM+5];
int id[L+5][L+5], pre[LIM+5], next[LIM+5];
double solve(int n, double g[][LIM+5], double ans[]){
int now = 0, r;
for (int c=0; c<n; c++){
int r = now;
for (int i=now; i<n; i++) if (fabs(g[i][c]) > eps && fabs(g[i][c]) < fabs(g[r][c]))
r = i;
for (r=now; r<n; r++) if (fabs(g[r][c])>eps) break;
if (r != now)
for (int j=c; j<=n; j++) swap(g[r][j], g[now][j]);
if (fabs(g[now][c]) < eps) continue;
for (int i=0; i<n; i++) if (i != now && fabs(g[i][c])>eps){
double tmp = g[i][c]/g[now][c];
for (int j=c; j<=n; j++)
g[i][j] -= tmp*g[now][j];
}
now ++;
}
return g[0][n]/g[0][0];
}
void calc(int r, int c){
int tr = max(r-2, 0), tc = c;
pre[id[r][c]] = id[tr][tc];
tr = r+1; if (tr>tc) swap(tr, tc);
if (tc == L) return;
next[id[r][c]] = id[tr][tc];
}
int main(){
#ifndef ONLINE_JUDGE
freopen("input.txt", "r", stdin);
#endif
memset(pre, 0xff, sizeof(pre));
memset(next, 0xff, sizeof(next));
for (int i=0, cnt = 0; i<L; i++)
for (int j=i; j<L; j++)
id[i][j] = cnt++;
for (int i=0; i<L; i++)
for (int j=i; j<L; j++){
calc(i, j);
}
while (scanf("%lf", &p) != EOF){
memset(g, 0, sizeof(g));
q = 1.0-p;
for (int i=0; i<L; i++)
for (int j=i; j<L; j++){
int k = id[i][j];
g[k][k] += 1;
if (pre[k] != -1) g[k][pre[k]] -= (1-p);
if (next[k] != -1) g[k][next[k]] -= p;
g[k][LIM] += 1;
}
double ans[LIM+10];
printf("%.6lf\n", solve(LIM, g, ans));
}
return 0;
}
这一题还有另一种思路
原来状态方程是dp[i][j] = p*dp[i+1][j] + (1-p)*dp[i-2][j] + 1 ,跟前面、后面的状态都有关。
其实两个账号可以分开来看,一个账号你要达到20, 另一个则要达到19
所以可以用dp[i]表示从i到i+1所需期望步数,
转移方程为 dp[i] = p*1+(1-p)*(dp[i-2]+dp[i-1]+dp[i]),
这样就只涉及前面的状态了,可以直接o(n)的dp
#include <cstdio>
using namespace std;
const int N = 21;
double p;
double dp[N], s[N];
int main(){
double p;
while (scanf("%lf", &p) != EOF){
dp[0] = 1/p;
dp[1] = 1/p/p;
for (int i=2; i<=20; i++){
dp[i] = (1/p-1)*(dp[i-2]+dp[i-1]) + 1/p;
}
s[0] = dp[0];
for (int i=1; i<=20; i++)
s[i] = s[i-1]+dp[i];
printf("%.6lf\n", s[19]+s[18]);
}
return 0;
}