Project Euler 363 Bézier Curves（几何+二分）

题目链接：

https://projecteuler.net/problem=363

题目：

A cubic Bézier curve is defined by four points: (P_0, P_1, P_2) and (P_3).

The curve is constructed as follows:

On the segments (P_0P_1), (P_1P_2) and (P_2P_3) the points (Q_0,Q_1) and (Q_2) are drawn such that

(P_0Q_0 / P_0P_1 = P_1Q_1 / P_1P_2 = P_2Q_2 / P_2P_3 = t (t in [0,1]).)

On the segments (Q_0Q_1) and (Q_1Q_2) the points (R_0) and (R_1) are drawn such that

(Q_0R_0 / Q_0Q_1 = Q_1R_1 / Q_1Q_2 = t) for the same value of (t).

On the segment (R_0R_1) the point B is drawn such that (R_0B / R_0R1 = t) for the same value of (t).

The Bézier curve defined by the points (P_0, P_1, P_2, P_3) is the locus of B as (Q_0) takes all possible positions on the segment (P_0P_1).

(Please note that for all points the value of t is the same.)

At this (external) web address you will find an applet that allows you to drag the points (P_0, P_1, P_2 and P_3) to see what the Bézier curve (green curve) defined by those points looks like. You can also drag the point (Q_0) along the segment (P_0P_1).

From the construction it is clear that the Bézier curve will be tangent to the segments (P_0P_1) in (P_0) and (P_2P_3) in (P_3).

A cubic Bézier curve with $P_0=(1,0), P_1=(1,v), P_2=(v,1) and P_3=(0,1) $is used to approximate a quarter circle.

The value v > 0 is chosen such that the area enclosed by the lines (OP_0, OP_3) and the curve is equal to π/4 (the area of the quarter circle).

By how many percent does the length of the curve differ from the length of the quarter circle?

That is, if L is the length of the curve, calculate 100 × (L − π/2)/(π/2)

Give your answer rounded to 10 digits behind the decimal point.

题解：

这道题超级有意思。

不查一下wiki根本不知道怎么下手...

https://en.wikipedia.org/wiki/Bézier_curve

因为Bézier curve上的4个点是：((1, 0), (1, v), (v, 1), (0, 1))

所以从wiki上得Specific cases可以知道，将4个点代入 (B(t)) 得到 (x(t))和 (y(t))：

+ - (B(t) = (1-t)^3 P_0 + 3(1-t)^2 t P_1 + 3(1-t) t^2 P_2 + t^3 P_3)

+ - (x(t) = (1 - t)^3 + 3(1-t)^2t + 3(1-t) t^2 v)

+ - (y(t) = 3(1-t)^2tv + 3(1-t) t^2 + t^3)

其中，(v) 未知，(0 <= t <= 1). (x(t), y(t)) 就是计算曲线上的坐标。

然后我们直接二分 (v) 就可以了。

得到 (v) 后可以直接计算曲线的长度(arc_length_of_a_curve)，也可以将曲线化为无数的小线段进行逼近。

代码：

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
const int maxn = 1e8;
const int mod = 1e9;
const double pi = acos(-1.0);

double x(double t, double v)
{
    return (1-t)*(1-t)*(1-t) + 3*(1-t)*(1-t)*t + 3*(1-t)*t*t*v;
}

double y(double t, double v)
{
    return 3*t*(1-t)*(1-t)*v + t*t*t + 3*t*t*(1-t);
}

double cal_distance(double a,double b,double c,double d)
{
  return (double)sqrt((a-c)*(a-c)+(b-d)*(b-d));
}

//https://en.wikipedia.org/wiki/B%C3%A9zier_curve
int main(int argc, char const *argv[]) {

  int dot =  500000;
  double left = 0.0, right = 1.0;
  double v = 0.0, area = 0.0;
  while (right - left > 1e-15) {
     v = (right + left) / 2;
     area = 0.0;
    for(int i = 0;i < dot;i++) {
      double now = (double) i / dot * 1.0;
      double next = (double) (i+1) / dot * 1.0;
      //cross product
      area += ( x(now,v) * y(next,v) - x(next,v) * y(now,v) ) / 2.0;
    }
    if(area > pi / 4.0) {
      right = v;
    }
    else {
      left = v;
    }
  }

  std::cout << "binary search finish !" << '
';
  std::cout << "left = " << left << '
';
  std::cout << "v = " << v << '
';
  std::cout << "area = " << area << '
';

  double L = 0.0;
  dot = 500000;
  for(int i = 0;i < dot;i++){
    double now = (double)i / dot * 1.0;
    double next = (double)(i + 1) / dot * 1.0;
    L += cal_distance(x(now, v), y(now, v), x(next, v), y(next, v));
  }
  printf("%.12f
",100.0*(L-pi/2.0)/pi*2.0);
  cerr << "Time elapsed: " << 1.0 * clock() / CLOCKS_PER_SEC << " s.
";
  return 0;
}