HNUOJ

BridgeTime Limit: 1000ms, Special Time Limit:2500ms,Memory Limit:32768KBTotal submit users: 153, Accepted users:110Problem 10046 : No special judgementProblem description

A suspension bridge suspends the roadway from huge main cables, which extend from one end of the bridge to the other. These cables rest on top of high towers and are secured at each end by anchorages. The towers enable the main cables to be draped over long distances.

Suppose that the maximum distance between two neighboring towers is D, and that the distance from the top of a tower to the roadway is H. Also suppose that the shape of a cable between any two neighboring towers is the same symmetric parabola (as shown in the figure). Now given B, the length of the bridge and L, the total length of the cables, you are asked to calculate the distance between the roadway and the lowest point of the cable, with minimum number of towers built (Assume that there are always two towers built at the two ends of a bridge).

Input

Standard input will contain multiple test cases. The first line of the input is a single integer T (1 <= T <= 10) which is the number of test cases. T test cases follow, each preceded by a single blank line.

For each test case, 4 positive integers are given on a single line. D - the maximum distance between two neighboring towers; H - the distance from the top of a tower to the roadway; B - the length of the bridge; and L - the total length of the cables.

It is guaranteed that B <= L. The cable will always be above the roadway.

Output

Results should be directed to standard output. Start each case with "Case #:" on a single line, where # is the case number starting from 1. Two consecutive cases should be separated by a single blank line.

For each test case, print the distance between the roadway and the lowest point of the cable, as is described in the problem. The value must be accurate up to two decimal places.

Sample Input 2 20 101 400 4042 1 2 3 4 Sample Output Case 1: 1.00 Case 2: 1.60 Judge Tips

方法一：

桥的间隔数为n = ceil(B/D)，每段绳子的长度为L / n，相邻两塔之间的距离为 B / n

主要问题还是在于已知抛物线的开口宽度w 和 抛物线的高度h 求抛物线的长度

弧长积分公式为：

设抛物线方程为f(x) = ax²，则这段抛物线弧长为

查积分表或者自己分部积分算一下： 

二分抛物线高度x，使得每段抛物线长度为L / n，所求答案为H - x

注：

C语言的math.h头文件中有ceil和floor两个函数；

double ceil(double x); float ceilf(float x); long double ceill(long double x); double floor(double x); float floorf(float x); long double floorl(long double x); 上面一个是把一个浮点数向上取整，下面一个是向下取整。

代码：

#include <cstdio> #include <cmath> inline double F(double a, double x) {//sqrt(a^2+x^2)的原函数 double a2 = a*a, x2 = x*x; double s = sqrt(a2+x2); return (x*s + a2*log(x+s))/2;//∫√(1+x^2 ) dx=x/2 √(1+x^2 )+1/2 log_e(x+√(1+x^2 ))+C } double length(double w, double h) {//宽为w，高为h的抛物线的长度 double a = 4*h/w/w; double b = 0.5/a; return 4*a*(F(b, w/2) - F(b, 0)); } int main() { //freopen("in.txt", "r", stdin); int T; scanf("%d", &T); for(int kase = 1; kase <= T; kase++) { int D, H, B, L; scanf("%d%d%d%d", &D, &H, &B, &L); int n = (B-1)/D + 1; //间隔数 double d = (double)B / n; //间隔 double l = (double)L / n; //每段绳长 double Left = 0, Right = H; while(Right - Left > 1e-5) {//二分求抛物线高度 double mid = (Right + Left) / 2; if(length(d, mid) > l) Right = mid; else Left = mid; } if(kase > 1) puts(""); printf("Case %d:\n%.2f\n", kase, H-Left); } return 0; } 方法二：

自适应Simpson算法

Simpson积分公式

原理是用二次曲线逼近原函数

引理：在平面直角坐标系中，由任意三点(x1, y1), (x2, y2), (x3, y3)（x1<x2<x3，x2 = (x1 + x3)/2）确定的抛物线y= f(x)在[x1, x3]的

定积分为

 

设f(x) = ax2 + bx + c，于是

左边 = (ax33/3 + bx32/2 + cx3) – (ax13/3 + bx12/2 + cx3)

    =1/6 × [ 2a(x33 – x13) + 3b(x32 – x12) + 6c(x3 – x1) ]

    =1/6 × (x3 – x1) ×[ 2a(x32 + x3x1 + x12) + 2b(x3 + x1) + 6c ]

    =1/6 (x3 – x1)

     ×[ (ax32 + bx3 + c) + (ax12 + bx1 + c) + a(x3+x1)2 + 2b(x3+x1) + 4c ]

    =1/6 (x3 – x1) [ y1 + y3 + a(x3+x1)2 + 2b(x3+x1) + 4c ]

注意到x2 = (x1 + x3)/2，所以

左边 = 1/6 (x3 – x1) [ y1 + y3 + 4ax22 + 4bx2+ 4c ]

    = 1/6 (x3 – x1) [ y1 + y3 + 4y2 ] = 右边

那么，如果我们把要积分的区域[L, R]划分成n份（n为偶数），份与份之间的分割线的横坐标为x0~xn，其中x0 = L, xn = R，对应的函数值分别为y0~yn，那么，对第0,1,2条分割线；第2,3,4条分割线；第4,5,6条分割线……之间的区域各自进行二次曲线拟合，即应用上面的公式，把结果累加起来，可以得到下面的Simpson积分公式：

  

有时，如果对划分的份数不太确定，可以使用一种叫做“自适应Simpson”的方法。对要积分的区域[L, R]分成两半[L, mid]和[mid, R]，如果用二次曲线对(L, mid, R)求出的积分值，和对(L, (L+mid)/2, mid)和(mid, (mid+R)/2, R)分别积分在加起来的值，相差不超过某个精确度，那么就可以停止划分。

（要注意误差的叠加效应，当划分份数比较多时，每一份的误差虽然很小，但叠加在一起可能很大，所以精确度尽量设小）

如果我们能确定被积函数（就是覆盖长度或面积关于扫描位置的函数）是个一次多项式或者二次多项式，但是我们不想把函数式求出来，那也没关系，用上面的方法积分，因为用了二次曲线拟合，答案还是准确的。

要注意的一点是，被积函数在我们想积分的一段上应该是有连续的取值的，就是说，我们积分的扫描线范围[L, R]，在[L, R]内不能出现的某条扫描线，它没有穿过任何图形。（这样很有可能L, R, mid都没有穿过图形，求出的答案为0）所以一开始要把有图形的“竖线区”预处理出来。预处理的方法是，每个图形都有一个“被穿过”段，把这些段进行合并就可以了。

自适应积分法的基本步骤如下：

(1) 将积分区间[a,b]分成两个相等的1级子区间[a,a+0.5h]和[a+0.5h,a+h],且h=b-a;

(2) 在上述两个1级子区间上用辛普森积分得到积分I_{a,a+0.5h}^{(1)}和I_{a+0.5h,a+h}^{(1)};

(3) 将子区间[a,a+0.5h]分成两个相等的2级子区间[a,a+0.5^2*h]和[a+0.5^2*h,a+0.5h];

(4)  采用辛普森积分计算得到:I_{a,a+0.5h}^{(2)}=I_{a,a+0.5^2*h}^{(1)}+I_{a+0.5^2*h,a+0.5h}^{(1)};

(5) 比较I_{a,a+0.5h}^{(2)}和I_{a,a+0.5h}^{(1)}, 如果|I_{a,a+0.5h}^{(1)} - I_{a,a+0.5h}^{(2)}|<10*0.5*epsi，

其中epsi为整体积分所需要的精度，则认为子区间[a,a+0.5h]上的积分I_{a,a+0.5h}^{(1)}

已达到所需精度, 不需要再细分;

否则就需要再细分, 对每个2级子区间做同样的判断.1级子区间[a+0.5h,a+h]的操作过程完全与上面相同

Simpson算法代码：

double F(double x) { //Simpson公式用到的函数 } double simpson(double a, double b)//三点Simpson法，这里要求F是一个全局函数 { double c = a + (b - a) / 2; return (F(a) + 4 * F(c) + F(b))*(b - a) / 6; } double asr(double a, double b, double eps, double A)//自适应Simpson公式（递归过程）。已知整个区间[a,b]上的三点Simpson值A { double c = a + (b - a) / 2; double L = simpson(a, c), R = simpson(c, b); if (fabs(L + R - A) <= 15 * eps)return L + R + (L + R - A) / 15.0; return asr(a, c, eps / 2, L) + asr(c, b, eps / 2, R); } double asr(double a, double b, double eps)//自适应Simpson公式（主过程） { return asr(a, b, eps, simpson(a, b)); }

最终解决问题代码：

#include<stdio.h> #include<stdlib.h> #include<math.h> double w,h,_a;//_a:设 F(x)=_a*x^2 double F(double x) { return sqrt(1+4*_a*_a*x*x); } double simpson(double a,double b)//三点simpson法 { double c=(a+b)/2; return (b-a)/6*(F(a)+4*F(c)+F(b)); } double asr(double a,double b,double eps,double A)//递归过程 { double c=(a+b)/2,L=simpson(a,c),R=simpson(c,b); if(fabs(L+R-A)<=15*eps) return L+R+(L+R-A)/15; return asr(a,c,eps/2,L)+asr(c,b,eps/2,R); } int main() { double H,L,left,right; int T,D,B,n,i; scanf("%d",&T); for(i=1;i<=T;i++) { scanf("%d%lf%d%lf",&D,&H,&B,&L); n=(B+D-1)/D;//平均分成的段数 L/=(double)n;//每段绳索长 w=(double)B/(double)n;//每段宽度 left=0; right=H; while(right-left>1e-4)//二分桥的"深度"h==H-y { h=(left+right)/2; _a=4*h/(w*w); if(2*asr(0,w/2,1e-4,simpson(0,w/2))<=L) left=h; else right=h; } if(i>1) printf("\n"); printf("Case %d:\n%.2lf\n",i,H-h); } return 0; }