具体推导看书<<数值分析>>
code：
#include <iostream>
using namespace std;
const int MAXN = 100;
int n;
double x[MAXN], y[MAXN]; //下标从0..n
double alph[MAXN], beta[MAXN], a[MAXN], b[MAXN];
double h[MAXN];
double m[MAXN]; //各点的一阶导数;
inline double sqr(double pa) 
{
return pa * pa;
}
double sunc(double p, int i) {
return (1 + 2 * (p - x[i]) / (x[i + 1] - x[i])) * sqr((p - x[i + 1]) / (x[i + 1] - x[i])) * y[i]
+ (1 + 2 * (p - x[i + 1]) / (x[i] - x[i + 1])) * sqr((p - x[i]) / (x[i + 1] - x[i])) * y[i + 1]
+ (p - x[i]) * sqr((p - x[i + 1]) / (x[i] - x[i + 1])) * m[i]
+ (p - x[i + 1]) * sqr((p - x[i]) / (x[i + 1] - x[i])) * m[i + 1];
}
int main() {
int i, j;
double xx;
freopen("threeInsert.in", "r", stdin);
scanf("%d", &n);
for (i = 0; i <= n; i++) scanf("%lf%lf", &x[i], &y[i]);
// scanf("%lf%lf", &m[0], &m[n]);
for (i = 0; i <= n - 1; i++) h[i] = x[i + 1] - x[i];
//第一种边界条件
//alph[0] = 0; alph[n] = 1; beta[0] = 2 * m[0]; beta[n] = 2 * m[n];
//第二种边界条件
alph[0] = 1; alph[n] = 0; beta[0] = 3 * (y[1] - y[0]) / h[0]; beta[n] = 3 * (y[n] - y[n - 1] / h[n - 1]);
for (i = 1; i <= n - 1; i++) {
alph[i] = h[i - 1] / (h[i - 1] + h[i]);
beta[i] = 3 * ((1 - alph[i]) * (y[i] - y[i - 1]) / h[i - 1] + alph[i] * (y[i + 1] - y[i]) / h[i]);
}
a[0] = - alph[0] / 2; b[0] = beta[0] / 2;
for (i = 1; i <= n; i++) {
a[i] = - alph[i] / (2 + (1 - alph[i]) * a[i - 1]);
b[i] = (beta[i] - (1 - alph[i]) * b[i - 1]) / (2 + (1 - alph[i]) * a[i - 1]);
}
m[n + 1] = 0;
for (i = n; i >= 0; i--) {
m[i] = a[i] * m[i + 1] + b[i];
}
scanf("%lf", &xx);
for (i = 0; i < n; i++) {
if (xx >= x[i] && xx <= x[i + 1]) break;
}
printf("%lf\n", sunc(xx, i));
return 0;
} posted @
2007-10-20 13:07 豪 阅读(3425) |
评论 (4) |
编辑 收藏
Sailboat
Problem H: Sailboat
In the sailboat race, the contestant is requested to along with the prearrange path. Sailing ship's power comes from wind power and contestant's manpower. The wind power can completely used.
In a competition, the contestants are requested to along with a 1/4 circles with radius R, the sailboat will goto east from south. During this process, the wind direction is straight from west to the east with constant speed and power.
In order to maintain the travel direction, the athlete must adjust the sail to the vertical angle from movement direction in any time.
If the speed of sailboat is proportional to the power at movement direction, the proportional factor is k. Supposes the wind power is f, the athlete manpower is h, please given the time of sailboat from the beginning to the end.
Input
The first line of each case consists of 4 double number, that is radius of path: R, wind power: f,athlete manpower: h and proportional factor:k. In order to avoid the floating point error, you needn't output the answer directly. The next line is a integer n, the following n lines gives a double value which is candidate answer.
Output
For each candidate of each case, Only "Yes" or "No" should be printed. Output "Yes" if the relative error to your answer is less than 3%, otherwise "No". For example, if the model answer is 100, and the candidate is 98 or 102, you should output "Yes". Output one blank line between neighboring case
Sample Input
1.0 2.0 1.0 1.0
2
0.35
0.76
Sample Output
No
Yes
Problem Source: provided by skywind
#include <iostream>
#include <cmath>
using namespace std;
const int MAXN = 100;
const double PI = acos(-1.0);
double R, F, H, K, ans;
int n, cas;
double func(double x) {
return R / K / (H + F * cos(x));
}
double romberg(double a, double b, double EPS = 1e-6) {
double t[MAXN][MAXN] = {0}, tmp;
int i, j, k, k2, m, m4;
t[0][0] = (func(a) + func(b)) * (b - a) / 2;
k = 1; k2 = 1;
while (1) {
tmp = 0;
for (i = 1; i <= k2; i++) {
tmp += func(a + (2 * i - 1) * (b - a) / (2 * k2));
}
t[0][k] = (t[0][k - 1] + tmp * (b - a) / k2) / 2;
for (m = 1, m4=4; m <= k; m++, m4 *= 4) {
t[m][k - m] = (m4 * t[m - 1][k - m + 1] - t[m - 1][k - m]) / (m4 - 1);
}
if (fabs(t[k][0] - t[k - 1][0]) < EPS) break;
k++; k2 *= 2;
}
return t[k][0];
}
void solve() {
double tmp;
scanf("%lf", &tmp);
if (fabs(tmp - ans) / ans < 0.03) printf("Yes\n");
else printf("No\n");
}
int main() {
freopen("2457.in", "r", stdin);
while (scanf("%lf%lf%lf%lf%d", &R, &F, &H, &K, &n) != EOF) {
if (cas) printf("\n");
else cas++;
ans = romberg(0, PI/2);
while (n--) {
solve();
}
}
return 0;
}
posted @
2007-10-20 01:02 豪 阅读(670) |
评论 (0) |
编辑 收藏Easy Problem
Time limit:1000 ms Memory limit:65536 KB
Total Submit:1755 (462 users) Accepted Submit:366 (332 users)
Description
In this problem, you're to calculate the distance between a point
P(
xp,
yp,
zp) and a segment (
x1,
y1,
z1) − (
x2,
y2,
z2), in a 3D space, i.e. the minimal distance from P to any point
Q(
xq,
yq,
zq) on the segment (a segment is part of a line).
Input
The first line contains a single integer T (1 ≤
T ≤ 1000), the number of test cases. Each test case is a single line containing 9 integers
xp,
yp,
zp,
x1,
y1,
z1,
x2,
y2,
z2. These integers are all in [-1000,1000].
Output
For each test case, print the case number and the minimal distance, to two decimal places.
Sample Input
3
0 0 0 0 1 0 1 1 0
1 0 0 1 0 1 1 1 0
-1 -1 -1 0 1 0 -1 0 -1
Sample Output
Case 1: 1.00
Case 2: 0.71
Case 3: 1.00
Problem Source
The 32nd ACM-ICPC Beijing First Round Internet Contest
其实和二分差不多，划个函数曲线出来，分三段，比划一下就很容易理解了：）
#include <iostream>
#include <cmath>
using namespace std;
double dist(double l[], double r[]) {
return sqrt((l[0]-r[0])*(l[0]-r[0])+(l[1]-r[1])*(l[1]-r[1])+(l[2]-r[2])*(l[2]-r[2]));
}
int main() {
// freopen("1024.in", "r", stdin);
int n, cas=0;
double l[3], r[3], p[3], p1[3], p2[3], d1, d2;
scanf("%d", &n);
while (n--) {
scanf("%lf%lf%lf%lf%lf%lf%lf%lf%lf", &p[0], &p[1], &p[2], &l[0], &l[1], &l[2], &r[0], &r[1], &r[2]);
while (dist(l, r) > 1e-4) {
p1[0] = (l[0] + r[0]) / 2;
p1[1] = (l[1] + r[1]) / 2;
p1[2] = (l[2] + r[2]) / 2;
p2[0] = (r[0] + p1[0]) / 2;
p2[1] = (r[1] + p1[1]) / 2;
p2[2] = (r[2] + p1[2]) / 2;
d1 = dist(p1, p); d2 = dist(p2, p);
if (d2 >= d1) {
r[0] = p2[0]; r[1] = p2[1]; r[2] = p2[2];
} else {
l[0] = p1[0]; l[1] = p1[1]; l[2] = p1[2];
}
}
printf("Case %d: %.2lf\n", ++cas, dist(p,l));
}
}
posted @
2007-10-18 11:00 豪 阅读(1173) |
评论 (0) |
编辑 收藏
昨晚去图书馆看了《计算机图形学——OpenGL实现》关于Bresenham算法的另一种推导方式。
Bresenham最精妙之处在于通过方程变换，然后得到迭代方程，从而消除了浮点运算。
下面简单写写自己对中点法推导的理解：
记：W = bx - ax, H = by - ay
所以 （ax, ay）和（bx, by）的理想直线为：
-W*(y-ay) + H*(x-ax) = 0
记：函数 f(x, y) = -2*W*(y-ay) + 2*H*(x-ax);
f(x,y)有如下性质：
f(x, y) < 0, 那么(x, y)在直线上方
f(x, y) > 0, 那么(x, y)在直线下方
现考虑 点L(Px+1, Py), 点U(Px+1, Py+1), 则LU中点M(Px+1, Py+1/2) 有：
如果f(Mx, My) < 0, 则M在理想直线上方， 所以选择L
如果f(Mx, My) > 0, 则M在理想直线下方， 所以选择U
则：
f(Mx,My) = -2*w*(Py+1/2-ay) + 2*H*(Px+1-ax)
当 x从Px+1移动到Px+2时， 考虑f变化M'和M'':
M'：在前一步没有增加y, M' = (Px+2, Py+1/2)
M''：在前一步增加了y, M' = (Px+2, Py+3/2)
对于 M'：
f(M'x, M'y) = -2*w*(Py+1/2-ay) + 2*H*(Px+2-ax) = f(Mx, My) + 2 * H
对于 M''：
f(M''x, M''y) = -2*w*(Py+3/2-ay) + 2*H*(Px+2-ax) = f(Mx, My) - 2 * (W-H)
所以
对于下一个“测试量”都有一个常数增量：前一次没有增加y，增量为2*H，如果增加了y，则增量为-2*(W-H)
对于初始条件：x = ax, y = ay
M = (ax+1, ay+1/2);
f(Mx, My) = -2*W*(ay+1/2-ay) + 2*H(ax+1-ax) = 2*H-W
Code:
#include <stdlib.h>
#include <math.h>
#include <GL/glut.h>
void myInit() {
glClearColor(1.0, 1.0, 1.0, 0.0);
glColor3f(0.0, 0.0, 0.0);
//glPointSize(2.0);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
gluOrtho2D(0.0, 640.0, 0.0, 480.0);
}
void setPixel(int x, int y) {
glBegin(GL_POINTS);
glVertex2i(x, y);
glEnd();
}
void lineBres(int xs, int ys, int xe, int ye) {
int W = xe - xs, H = ye - ys, f = 2 * H - W, tH = 2 * H, tHW = 2 * (H - W);
int x, y;
if (xs > xe) {
x = xe;
y = ye;
xe = xs;
} else {
x = xs;
y = ys;
}
while (x <= xe) {
setPixel(x, y);
x++;
if (f<0) {
f += tH;
} else {
y++;
f += tHW;
}
}
}
void myDisplay() {
glClear(GL_COLOR_BUFFER_BIT);
lineBres(20, 10, 300, 180);
glFlush();
}
int main(int argc, char **argv) {
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_SINGLE|GLUT_RGB);
glutInitWindowSize(640, 480);
glutInitWindowPosition (100, 150);
glutCreateWindow("Bresenham画线");
glutDisplayFunc(myDisplay);
myInit();
glutMainLoop();
return 0;
}
posted @
2007-10-11 11:57 豪 阅读(1106) |
评论 (0) |
编辑 收藏
扩展了一点，有兴趣的可以去看看，MFC去，sigh~
http://www.cppblog.com/qywyh/articles/32740.html
posted @
2007-09-23 21:34 豪 阅读(1386) |
评论 (0) |
编辑 收藏