具體推導看書<<數值分析>>
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 豪 閱讀(3603) |
評論 (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 豪 閱讀(704) |
評論 (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 豪 閱讀(1216) |
評論 (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 豪 閱讀(1166) |
評論 (0) |
編輯 收藏
擴展了一點，有興趣的可以去看看，MFC去，sigh~
http://www.shnenglu.com/qywyh/articles/32740.html
posted @
2007-09-23 21:34 豪 閱讀(1427) |
評論 (0) |
編輯 收藏