AOJ 2097 : Triangles - 土下座しながら探索中

問題概要：
1辺の長さが１の正三角形をn個、任意の点を中心として 360/n 度間隔で置いていく
この時n個の正三角形が覆う面積はどのくらいか？

解法：
正多角形から余計な部分を引いて面積を出した

上の図の黒い矢印の部分の三角形を正多角形から引く

n個の正三角形を置く亊でできる頂点の数をVとする
nが３の倍数の場合、これまでの正三角形とかぶるものがでてくるのでV == n となる
そうでない場合、どれもかぶらないので V == 3 * n となる

以下は　n = 2 の例

先ほどの黒い矢印がさす三角形の面積を求める
まず、θ２は正三角形の１つの角（６０度）の半分なので30度となる
θ３は、正V角形のとき(360/V)/2となる
点P1、P２は原点（０、０）から点（１、０）を６０度回転させて正三角形の座標を適当に作り、そこから三角形の外心を求め中心点をそれらの点から引くことで求める
Lはabs(P1-P2)/2となる
L2は正弦定理から、
L/sin(θ３) = L2/sin(θ２)
L2 = L * sin(θ２) / sin(θ3)
点P1とP2の距離をlenとする

求めたかった三角形の面積subは L * L2　（ = L * L2 / 2 * 2 ）となる
あとは正V角形の面積 ( = ( V * len * len ) / ( 4 * tan( PI / V ) ) ) から sub * V を引けば良い

ちなみにJAGの解説を読むとこんな長くやらなくても１０行くらいで解く亊ができる

実装時間：
４、５時間
最初 V をnが３の倍数かどうかで判定していた訳ではなく実際にsetにつめて数えていた
当然ながらまともに数える亊などできないのでnがおおきい時まともな答えがでていなかった
それを実装が下手で誤差がおおきくなりすぎているものと勘違いして数時間誤差を減らそうとしていたためこの様になった
（諦めてJAG解説を読んで本当の原因を知る）

コード :

double R;

Point getMP(Point a,Point b,Point c)
{
  double nm[3] = {norm(a),norm(b),norm(c)};
  double u = 0.5/(a.x*b.y-b.x*a.y+b.x*c.y-c.x*b.y+c.x*a.y-a.x*c.y);
  double x = u * (nm[0]*b.y-nm[1]*a.y+nm[1]*c.y-nm[2]*b.y+nm[2]*a.y-nm[0]*c.y);
  double y = u * (a.x*nm[1]-b.x*nm[0]+b.x*nm[2]-c.x*nm[1]+c.x*nm[0]-a.x*nm[2]);
  R = sqrt((x-a.x)*(x-a.x)+(y-a.y)*(y-a.y)); //半径
  return Point(x,y);
}

int n,V;

int main(){

  while(cin >> n,n){

    if( n == 1 )
      {
	cout << "0.433013" << endl;
	continue;
      }

    if(n % 3)V = 3 * n;
    else     V = n;

    double theta = 360.0 / V;
    double theta2 = ( ( 180.0 - theta ) * 0.5 ) - 30.0;
    assert( theta2 >= 0.0 );  

    Point sp = Point(1,0);
    Point np = rotate(sp,toRad(60.0));
    Point mp = getMP(Point(0,0),sp,np);

    /* solution - 1 
    double h = R  / ( sqrt(3.0) + ( 1.0 / tan( toRad(theta*0.5) ) ) );
    double area2 = ( h * R / 2 ) * V * 2;
    cout << setiosflags(ios::fixed) << setprecision(7) << area2  << endl;  
    continue;
    */

    // solution - 2 
    np = np - mp;
    Point np2 = rotate(np,toRad(theta));

    double len = abs(np-np2);
    double L = len * 0.5;
    double theta3 = 180.0 - ( 90.0 + theta2 );
    assert( theta3 >= 0 );
    double L2 = L * sin(toRad(theta2)) / sin(toRad(theta3));

    double sub = L * L2;

    double area = ( V * len * len ) / ( 4 * tan( M_PI / V ) );

    cout << setiosflags(ios::fixed) << setprecision(7) << area - sub * V << endl;


  }
  return 0;
}