daily pastebin goal
57%
SHARE
TWEET

Code

a guest Sep 29th, 2013 9 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <iostream>
  2. #include <iomanip>
  3. #include <sstream>
  4. #include <fstream>
  5. #include <string>
  6. #include <cstdlib>
  7. #include <cstdio>
  8. #include <climits>
  9. #include <cstring>
  10. #include <cassert>
  11. #include <ctime>
  12. #include <cmath>
  13. #include <vector>
  14. #include <queue>
  15. #include <stack>
  16. #include <list>
  17. #include <set>
  18. #include <map>
  19. #include <bitset>
  20. #include <algorithm>
  21. #include <utility>
  22. #include <numeric>
  23. #include <functional>
  24.  
  25. #define forn(i, n) for (int i = 0; i < int(n); i++)
  26. #define forl(i, n) for (int i = 1; i <= int(n); i++)
  27. #define ford(i, n) for (int i = int(n) - 1; i >= 0; i--)
  28. #define fore(i, l, r) for (int i = int(l); i <= int(r); i++)
  29. #define correct(x, y, n, m) (0 <= (x) && (x) < (n) && 0 <= (y) && (y) < (m))
  30. #define all(a) (a).begin(), (a).end()
  31. #define sz(a) int((a).size())
  32. #define pb(a) push_back(a)
  33. #define mp(x, y) make_pair((x), (y))
  34. #define ft first
  35. #define sc second
  36. #define x first
  37. #define y second
  38. #define eprintf(...) fprintf(stderr, __VA_ARGS__)
  39.  
  40. using namespace std;
  41.  
  42. typedef long long li;
  43. typedef long double ld;
  44. typedef pair<int, int> pti;
  45.  
  46. template<typename X> inline X abs(const X& a) { return a < 0? -a: a; }
  47. template<typename X> inline X sqr(const X& a) { return a * a; }
  48.  
  49. const int INF = int(1e9);
  50. const li INF64 = li(1e18);
  51. const ld EPS = 1e-9, PI = 3.1415926535897932384626433832795;
  52. const ld inf = 1e100;
  53.  
  54. struct pt
  55. {
  56.         ld x, y;
  57.         pt() {}
  58.         pt(ld x, ld y): x(x), y(y) {}
  59. };
  60.  
  61. inline bool operator< (const pt& a, const pt& b)
  62. {
  63.         if (abs(a.x - b.x) > EPS) return a.x < b.x;
  64.         if (abs(a.y - b.y) > EPS) return a.y < b.y;
  65.         return false;
  66. }
  67.  
  68. inline bool operator== (const pt& a, const pt& b)
  69. {
  70.         if (abs(a.x - b.x) > EPS) return false;
  71.         if (abs(a.y - b.y) > EPS) return false;
  72.         return true;
  73. }
  74.  
  75. inline pt operator- (const pt& a, const pt& b) { return pt(a.x - b.x, a.y - b.y); }
  76. inline pt operator+ (const pt& a, const pt& b) { return pt(a.x + b.x, a.y + b.y); }
  77. inline pt operator* (const pt& a, const ld& b) { return pt(a.x * b, a.y * b); }
  78. inline pt operator/ (const pt& a, const ld& b) { return pt(a.x / b, a.y / b); }
  79.  
  80. inline ld len(const pt& v) { return sqrt(sqr(v.x) + sqr(v.y)); }
  81. inline ld dist(const pt& a, const pt& b) { return len(a - b); }
  82.  
  83. inline pt normal(const pt& v) { return len(v) < EPS ? v : v / len(v); }
  84. inline pt normal(const pt& v, const ld& l) { return normal(v) * l; }
  85.  
  86. inline ld dot(const pt& a, const pt&b) { return a.x * b.x + a.y * b.y; }
  87. inline ld cross(const pt& a, const pt&b) { return a.x * b.y - a.y * b.x; }
  88.  
  89. const int N = 100 + 3;
  90.  
  91. int n;
  92. pt a[N];
  93.  
  94. inline bool read()
  95. {
  96.         if (scanf("%d", &n) != 1)
  97.                 return false;
  98.                
  99.         forn(i, n)
  100.         {
  101.                 double x, y;
  102.                 assert(scanf("%lf%lf", &x, &y) == 2);
  103.                
  104.                 a[i] = pt(x, y);
  105.         }
  106.        
  107.         return true;
  108. }
  109.  
  110. int szr;
  111. ld r[2];
  112.  
  113. inline void solve(const ld& a, const ld& b, const ld& c)
  114. {
  115.         szr = 0;
  116.  
  117.         if (abs(a) > EPS)
  118.         {
  119.                 ld D = sqr(b) - 4 * a * c;
  120.                 if (D + EPS < 0) return;
  121.  
  122.                 D = sqrt(max(ld(0), D));               
  123.                 r[szr++] = (-b + D) / (2 * a);
  124.                 r[szr++] = (-b - D) / (2 * a);         
  125.                 return;
  126.         }
  127.        
  128.         if (abs(b) > EPS)
  129.         {
  130.                 r[szr++] = -c / b;
  131.                 return;
  132.         }
  133.        
  134.         if (abs(c) < EPS)
  135.         {
  136.                 r[szr++] = 0;
  137.                 return;
  138.         }
  139. }
  140.  
  141. inline ld solve(ld sx, ld sy, ld sxy, ld sx2, ld sy2, int n)
  142. {
  143.         if (n == 0) return 0;
  144.        
  145.         pt c(sx / n, sy / n);
  146.        
  147.         sxy = sxy + c.x * c.y * n - c.x * sy - c.y * sx;
  148.         sx2 = sx2 + sqr(c.x) * n - 2 * c.x * sx;
  149.         sy2 = sy2 + sqr(c.y) * n - 2 * c.y * sy;
  150.         sx = sy = 0;
  151.        
  152.         ld res = inf;
  153.        
  154.         forn(iter, 2)
  155.         {
  156.                 solve(sxy, (sy2 - sx2) * (iter == 0 ? 1 : -1), -sxy);
  157.                
  158.                 forn(i, szr)
  159.                 {
  160.                         const ld A = (iter == 0 ? r[i] : 1), B = (iter == 0 ? 1 : r[i]);                       
  161.                         res = min(res, (sqr(A) * sx2 + 2 * A * B * sxy + sqr(B) * sy2) / (sqr(A) + sqr(B)));
  162.                 }
  163.         }
  164.        
  165.         assert(res < inf / 2);
  166.         return res;
  167. }
  168.  
  169. pt pole, dir;
  170. ld A, B, C;
  171.  
  172. inline ld proj(const pt& a) { return dot(a - pole, dir); }
  173.  
  174. inline int sign(const pt& p)
  175. {
  176.         ld val = A * p.x + B * p.y + C;
  177.         return val + EPS < 0 ? -1 : val > EPS ? 1 : 0;
  178. }
  179.  
  180. inline bool cmp(const pt& a, const pt& b)
  181. {
  182.         ld v1 = proj(a), v2 = proj(b);
  183.         if (abs(v1 - v2) > EPS) return v1 < v2;
  184.         return false;
  185. }
  186.  
  187. pt p[N];
  188.  
  189. ld x[2][N];
  190. ld y[2][N];
  191. ld xy[2][N];
  192. ld x2[2][N];
  193. ld y2[2][N];
  194. int m[2][N];
  195.  
  196. inline ld sumX(int up, int l, int r)
  197. {
  198.         if (l > r) return 0;
  199.         return x[up][r] - (l == 0 ? 0 : x[up][l - 1]);
  200. }
  201.  
  202. inline ld sumY(int up, int l, int r)
  203. {
  204.         if (l > r) return 0;
  205.         return y[up][r] - (l == 0 ? 0 : y[up][l - 1]);
  206. }
  207.  
  208. inline ld sumXY(int up, int l, int r)
  209. {
  210.         if (l > r) return 0;
  211.         return xy[up][r] - (l == 0 ? 0 : xy[up][l - 1]);
  212. }
  213.  
  214. inline ld sumX2(int up, int l, int r)
  215. {
  216.         if (l > r) return 0;
  217.         return x2[up][r] - (l == 0 ? 0 : x2[up][l - 1]);
  218. }
  219.  
  220. inline ld sumY2(int up, int l, int r)
  221. {
  222.         if (l > r) return 0;
  223.         return y2[up][r] - (l == 0 ? 0 : y2[up][l - 1]);
  224. }
  225.  
  226. inline ld sumM(int up, int l, int r)
  227. {
  228.         if (l > r) return 0;
  229.         return m[up][r] - (l == 0 ? 0 : m[up][l - 1]);
  230. }
  231.  
  232. inline void clearSum()
  233. {
  234.         forn(i, 2)
  235.                 forn(j, n)
  236.                 {
  237.                         x[i][j] = 0;
  238.                         y[i][j] = 0;
  239.                         xy[i][j] = 0;
  240.                         x2[i][j] = 0;
  241.                         y2[i][j] = 0;
  242.                         m[i][j] = 0;
  243.                 }
  244. }
  245.  
  246. inline void makeSum()
  247. {
  248.         forn(i, 2)
  249.                 forl(j, n - 1)
  250.                 {
  251.                         x[i][j] += x[i][j - 1];
  252.                         y[i][j] += y[i][j - 1];
  253.                         xy[i][j] += xy[i][j - 1];
  254.                         x2[i][j] += x2[i][j - 1];
  255.                         y2[i][j] += y2[i][j - 1];
  256.                         m[i][j] += m[i][j - 1];        
  257.                 }
  258. }
  259.  
  260. int szv;
  261. pt v[N];
  262.  
  263. const int MASK = (1 << 4) + 3;
  264.  
  265. ld mx[MASK];
  266. ld my[MASK];
  267. ld mxy[MASK];
  268. ld mx2[MASK];
  269. ld my2[MASK];
  270. int mm[MASK];
  271.  
  272. int least[MASK];
  273.  
  274. inline void solve()
  275. {
  276.         ld SX = 0, SY = 0, SXY = 0, SX2 = 0, SY2 = 0;
  277.         int N = n;
  278.        
  279.         forn(i, n)
  280.         {
  281.                 SX += a[i].x;
  282.                 SY += a[i].y;
  283.                 SXY += a[i].x * a[i].y;
  284.                 SX2 += sqr(a[i].x);
  285.                 SY2 += sqr(a[i].y);
  286.         }
  287.        
  288.         ld ans = inf;
  289.        
  290.         forn(ii, n)
  291.                 forn(jj, ii)
  292.                 {
  293.                         pole = a[ii], dir = normal(a[jj] - a[ii]);
  294.                         A = -dir.y, B = dir.x, C = -(A * pole.x + B * pole.y);
  295.                        
  296.                         {
  297.                                 ld v = sqr(A) + sqr(B);
  298.                                 A /= v, B /= v, C /= v;
  299.                         }
  300.                        
  301.                         forn(i, n) p[i] = a[i];
  302.                         sort(p, p + n, cmp);
  303.                        
  304.                         clearSum();
  305.                        
  306.                         forn(i, n)
  307.                         {
  308.                                 int s = sign(p[i]);                            
  309.                                 if (s == 0) continue;
  310.                                
  311.                                 int up = s > 0;
  312.                                
  313.                                 x[up][i] = p[i].x;
  314.                                 y[up][i] = p[i].y;
  315.                                 xy[up][i] = p[i].x * p[i].y;
  316.                                 x2[up][i] = sqr(p[i].x);
  317.                                 y2[up][i] = sqr(p[i].y);
  318.                                 m[up][i] = 1;
  319.                         }
  320.                        
  321.                         makeSum();
  322.                        
  323.                         for (int i = 0, j = 0; i < n; i = j)
  324.                         {
  325.                                 szv = 0;
  326.                                 v[szv++] = a[ii], v[szv++] = a[jj];
  327.                                
  328.                                 while (abs(proj(p[i]) - proj(p[j])) < EPS)
  329.                                 {
  330.                                         v[szv++] = p[j];
  331.                                         j++;
  332.                                 }
  333.                                
  334.                                 sort(v, v + szv);
  335.                                 szv = int(unique(v, v + szv) - v);
  336.                                
  337.                                 assert(szv >= 2 && szv <= 4);
  338.                                
  339.                                 mx[0] = my[0] = mxy[0] = mx2[0] = my2[0] = mm[0] = 0;
  340.                                
  341.                                 forl(mask, (1 << szv) - 1)
  342.                                 {
  343.                                         int idx = least[mask];
  344.                                         int pmask = mask ^ (1 << idx);
  345.                                        
  346.                                         mx[mask] = mx[pmask] + v[idx].x;
  347.                                         my[mask] = my[pmask] + v[idx].y;                                       
  348.                                         mxy[mask] = mxy[pmask] + v[idx].x * v[idx].y;
  349.                                         mx2[mask] = mx2[pmask] + sqr(v[idx].x);
  350.                                         my2[mask] = my2[pmask] + sqr(v[idx].y);
  351.                                         mm[mask] = mm[pmask] + 1;
  352.                                 }
  353.                                
  354.                                 ld curx = sumX(0, 0, i - 1) + sumX(1, j, n - 1);
  355.                                 ld cury = sumY(0, 0, i - 1) + sumY(1, j, n - 1);
  356.                                 ld curxy = sumXY(0, 0, i - 1) + sumXY(1, j, n - 1);
  357.                                 ld curx2 = sumX2(0, 0, i - 1) + sumX2(1, j, n - 1);
  358.                                 ld cury2 = sumY2(0, 0, i - 1) + sumY2(1, j, n - 1);
  359.                                 int curn = sumM(0, 0, i - 1) + sumM(1, j, n - 1);
  360.                                
  361.                                 forn(mask, 1 << szv)
  362.                                 {
  363.                                         curx += mx[mask];
  364.                                         cury += my[mask];
  365.                                         curxy += mxy[mask];
  366.                                         curx2 += mx2[mask];
  367.                                         cury2 += my2[mask];
  368.                                         curn += mm[mask];
  369.                                        
  370.                                         ans = min(ans, solve(curx, cury, curxy, curx2, cury2, curn) + solve(SX - curx, SY - cury, SXY - curxy, SX2 - curx2, SY2 - cury2, N - curn));
  371.                                        
  372.                                         curx -= mx[mask];
  373.                                         cury -= my[mask];
  374.                                         curxy -= mxy[mask];
  375.                                         curx2 -= mx2[mask];
  376.                                         cury2 -= my2[mask];
  377.                                         curn -= mm[mask];
  378.                                 }
  379.                         }
  380.                 }
  381.                
  382.         assert(ans < inf / 2);
  383.         cout << ans / n << endl;
  384. }
  385.  
  386. int main()
  387. {
  388. #ifdef SU2_PROJ
  389.     assert(freopen("input.txt", "rt", stdin));
  390.     //assert(freopen("output.txt", "wt", stdout));
  391. #endif
  392.  
  393.         forl(i, MASK - 1)
  394.         {
  395.                 forn(j, 4)
  396.                 {
  397.                         if (i & (1 << j))
  398.                         {
  399.                                 least[i] = j;
  400.                                 break;
  401.                         }
  402.                 }
  403.         }
  404.    
  405.     cout << setprecision(10) << fixed;
  406.     cerr << setprecision(10) << fixed;
  407.    
  408.     assert(read());
  409.     solve();
  410.        
  411.     return 0;
  412. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top