Advertisement
Guest User

Untitled

a guest
Dec 2nd, 2015
625
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 1.62 KB | None | 0 0
  1. double CalcHalf(LL r1, LL r2, LL dx, LL dy)
  2. {
  3.     LL l2 = dx * dx + dy*dy;
  4.     double l = sqrt( static_cast<double>(l2) );
  5.     double x = static_cast<double>(r1*r1 - r2*r2 + l2) / (2.0 * l);
  6.     double angle = 0.0;
  7.    
  8.     if (sqrt(2.0) * x < static_cast<double>(r1))
  9.     {
  10.         angle = acos( x / r1 );
  11.     }
  12.     else
  13.     {
  14.         double S = sqrt( static_cast<double>( (r1+r2)*(r1+r2) - l2 ) * static_cast<double>(l2 - (r1-r2)*(r1-r2)) ) / 4.0;
  15.         double y = 2.0 * S / l;
  16.         if (y / r1 < 0.1)
  17.         {
  18.             // SPECIAL CASE FOR SMALL ANGLES
  19.             // asin(x)= x + 1/6 x^3 + 3/40 x^5 + 5/112 x^7 + 35/1152 x^9 + 63/2816 x^11 + 231/13312 x^13 + 143/10240 x^15 +...
  20.             double c[12] = { 0, 0, 0, 1.0 / 6.0, 0.0, 3.0 / 40.0, 0.0, 5.0 / 112.0, 0.0, 35.0 / 1152.0, 0.0, 63.0 / 2816.0 };
  21.             double v[12] = { 0 };
  22.             double vSin = 1.0;
  23.             double sin2A = 2.0 * (y / r1) * (x/r1);
  24.             for (size_t i = 0; i < 12; ++i)
  25.             {
  26.                 v[i] = vSin * c[i];
  27.                 vSin *= sin2A;
  28.             }
  29.  
  30.             double s = 0.0;
  31.             for (int i = 11; i >= 0; --i)
  32.                 s += v[i];
  33.  
  34.             return static_cast<double>(r1*r1) * s / 2.0;
  35.         }
  36.  
  37.         angle = asin(y / r1);
  38.     }
  39.  
  40.     return static_cast<double>(r1*r1) * ( 2 * angle - sin(2*angle) ) / 2.0 ;
  41. }
  42.  
  43. void SolveD()
  44. {
  45.     LL x1, y1, r1;
  46.     LL x2, y2, r2;
  47.  
  48.     cin >> x1 >> y1 >> r1;
  49.     cin >> x2 >> y2 >> r2;
  50.  
  51.     x2 -= x1;
  52.     y2 -= y1;
  53.     x1 = y1 = 0;
  54.  
  55.     if ((r1 + r2)*(r1 + r2) <= x2*x2 + y2*y2)
  56.     {
  57.         cout << 0;
  58.     }
  59.     else if (x2*x2 + y2*y2 <= (r2 - r1) * (r2 - r1))
  60.     {
  61.         double thePi = 3.14159265358979;
  62.         cout << std::setprecision(17) << min(r1, r2) * thePi * min(r1, r2);
  63.     }
  64.     else
  65.     {
  66.         cout << std::setprecision(17) << CalcHalf(r1, r2, x2, y2) + CalcHalf(r2, r1, x2, y2);
  67.     }
  68. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement