NEERC/asteroids/asteroids_petr.java

   1 import java.io.BufferedReader;
   2 import java.io.FileReader;
   3 import java.io.PrintWriter;
   4 import java.io.IOException;
   5 import java.util.*;
   6
   7 public class asteroids_petr implements Runnable {
   8     static final int MAX_POWER = 3;
   9
  10     static class Point {
  11         double x;
  12         double y;
  13         double z;
  14
  15         Point(double x, double y, double z) {
  16             this.x = x;
  17             this.y = y;
  18             this.z = z;
  19         }
  20
  21         void swapXY() {
  22             double t = x;
  23             x = y;
  24             y = t;
  25         }
  26
  27         void swapXZ() {
  28             double t = x;
  29             x = z;
  30             z = t;
  31         }
  32
  33         public void rotateXY(double v) {
  34             double nx = x * Math.cos(v) - y * Math.sin(v);
  35             double ny = x * Math.sin(v) + y * Math.cos(v);
  36             x = nx;
  37             y = ny;
  38         }
  39     }
  40
  41     static class Point2D implements Comparable<Point2D> {
  42         double x;
  43         double y;
  44         double angle;
  45
  46         Point2D(double x, double y) {
  47             this.x = x;
  48             this.y = y;
  49         }
  50
  51         public int compareTo(Point2D point2D) {
  52             return Double.compare(angle, point2D.angle);
  53         }
  54     }
  55
  56     static int count = 0;
  57
  58     static class Plane {
  59         double a;
  60         double b;
  61         double c;
  62         double d;
  63
  64         public Plane(Point p1, Point p2, Point p3) {
  65             a = (p2.z - p1.z) * (p3.y - p1.y) - (p2.y - p1.y) * (p3.z - p1.z);
  66             b = (p2.x - p1.x) * (p3.z - p1.z) - (p2.z - p1.z) * (p3.x - p1.x);
  67             c = (p2.y - p1.y) * (p3.x - p1.x) - (p2.x - p1.x) * (p3.y - p1.y);
  68             double z = Math.sqrt(a * a + b * b + c * c);
  69             a /= z;
  70             b /= z;
  71             c /= z;
  72             d = -(a * p1.x + b * p1.y + c * p1.z);
  73         }
  74
  75         public int side(Point p) {
  76             double z = sideDouble(p);
  77             if (Math.abs(z) < 1e-11)
  78                 return 0;
  79             else if (z < 0)
  80                 return -1;
  81             else
  82                 return 1;
  83         }
  84
  85         public double sideDouble(Point p) {
  86             return a * p.x + b * p.y + c * p.z + d;
  87         }
  88     }
  89
  90     static class Polyhedron {
  91         Point[] points;
  92
  93         void swapXY() {
  94             for (Point p : points)
  95                 p.swapXY();
  96         }
  97
  98         void swapXZ() {
  99             for (Point p : points)
 100                 p.swapXZ();
 101         }
 102
 103         double integrate(double[] what) {
 104             Set<Double> allXs = new TreeSet<Double>();
 105             for (Point p : points)
 106                 allXs.add(p.x);
 107             Double prev = null;
 108             double sum = 0;
 109             for (Double cur : allXs) {
 110                 if (prev != null) {
 111                     double[] areas = new double[MAX_POWER + 1];
 112                     double[] xs = new double[MAX_POWER + 1];
 113                     for (int i = 0 ; i <= MAX_POWER; ++i) {
 114                         xs[i] = (prev * (MAX_POWER - i) + cur * i) / MAX_POWER;
 115                         areas[i] = getArea(xs[i]) * evalPolynomial(what, xs[i]);
 116                     }
 117                     sum += (cur - prev) / 8.0 * (areas[0] + 3 * areas[1] + 3 * areas[2] + areas[3]);
 118                 }
 119                 prev = cur;
 120             }
 121             return sum;
 122         }
 123
 124         private double evalPolynomial(double[] koef, double x) {
 125             double res = 0;
 126             for (int i = koef.length - 1; i >= 0; --i)
 127                 res = res * x + koef[i];
 128             return res;
 129         }
 130
 131         static Random random = new Random();
 132
 133         private double getArea(double x) {
 134             List<Point2D> p2d = new ArrayList<Point2D>();
 135             for (Point p1 : points) {
 136                 if (p1.x == x) {
 137                     p2d.add(new Point2D(p1.y, p1.z));
 138                     continue;
 139                 }
 140                 if (p1.x < x)
 141                     for (Point p2 : points) {
 142                         if (p2.x == x) {
 143                             continue;
 144                         }
 145                         if (p2.x > x) {
 146                             p2d.add(new Point2D(p1.y + (p2.y - p1.y) * (x - p1.x) / (p2.x - p1.x), p1.z + (p2.z - p1.z) * (x - p1.x) / (p2.x - p1.x)));
 147                         }
 148                     }
 149             }
 150             if (p2d.isEmpty())
 151                 return 0.0;
 152             Point2D corner = new Point2D(0, 0);
 153             double totalWeight = 0.0;
 154             for (Point2D p : p2d) {
 155                 double weight = random.nextDouble();
 156                 corner.x += weight * p.x;
 157                 corner.y += weight * p.y;
 158                 totalWeight += weight;
 159             }
 160             double dx = -corner.x / totalWeight;
 161             double dy = -corner.y / totalWeight;
 162             for (Point2D p : p2d) {
 163                 p.x += dx;
 164                 p.y += dy;
 165                 p.angle = Math.atan2(p.y, p.x);
 166             }
 167             Collections.sort(p2d);
 168             corner = p2d.get(0);
 169             for (Point2D p : p2d) {
 170                 if (p.x < corner.x || p.x == corner.x && p.y < corner.y)
 171                     corner = p;
 172             }
 173             int start =0;
 174             for (int i = 0; i < p2d.size(); ++i)
 175                 if (p2d.get(i) == corner)
 176                     start = i;
 177             List<Point2D> hull = new ArrayList<Point2D>();
 178             for (int i = 0; i <= p2d.size(); ++i) {
 179                 hull.add(p2d.get((start + i) % p2d.size()));
 180                 while (hull.size() >= 3) {
 181                     Point2D a = hull.get(hull.size() - 3);
 182                     Point2D b = hull.get(hull.size() - 2);
 183                     Point2D c = hull.get(hull.size() - 1);
 184                     /*if ((b.x - a.x) * (c.y - b.y) - (b.y - a.y) * (c.x - b.x) == 0 && (b.x - a.x) * (c.x - b.x) + (b.y - a.y) * (c.y - b.y) < 0)
 185                         throw new RuntimeException();*/
 186                     if ((b.x - a.x) * (c.y - b.y) - (b.y - a.y) * (c.x - b.x) <= 0)
 187                         hull.remove(hull.size() - 2);
 188                     else
 189                         break;
 190                 }
 191             }
 192             double area = 0;
 193             for (int i = 0; i < hull.size() - 1; ++i) {
 194                 Point2D a = hull.get(i);
 195                 Point2D b = hull.get(i + 1);
 196                 area += (a.y + b.y) * (a.x - b.x);
 197                 if (i < hull.size() - 2) {
 198                     Point2D c = hull.get(i + 2);
 199                     if ((b.x - a.x) * (c.y - b.y) - (b.y - a.y) * (c.x - b.x) <= 0)
 200                         throw new RuntimeException();
 201                 }
 202             }
 203             if (area < -1e-12)
 204                 throw new RuntimeException();
 205             return area / 2.0;
 206         }
 207
 208         public double getDistance(Point center) {
 209             double res = Double.MAX_VALUE;
 210             for (int i1 = 0; i1 < points.length; ++i1)
 211                 for (int i2 = i1 + 1; i2 < points.length; ++i2)
 212                     for (int i3 = i2 + 1; i3 < points.length; ++i3) {
 213                         Plane plane = new Plane(points[i1], points[i2], points[i3]);
 214                         boolean[] was = new boolean[3];
 215                         for (Point p : points)
 216                             was[plane.side(p) + 1] = true;
 217                         if (was[0] && was[2])
 218                             continue;
 219                         res = Math.min(res, Math.abs(plane.sideDouble(center)));
 220                     }
 221             return res;
 222         }
 223
 224         public void rotateXY(double v) {
 225             for (Point p : points)
 226                 p.rotateXY(v);
 227         }
 228     }
 229
 230     Polyhedron parsePolyhedron() throws IOException {
 231         int n = nextInt();
 232         Polyhedron p = new Polyhedron();
 233         p.points = new Point[n];
 234         for (int i = 0; i < n; ++i) {
 235             int x = nextInt();
 236             int y = nextInt();
 237             int z = nextInt();
 238             p.points[i] = new Point(x, y, z);
 239         }
 240         return p;
 241     }
 242
 243     private void solve() throws IOException {
 244         Polyhedron[] polyhedrons = new Polyhedron[2];
 245         for (int i = 0; i < 2; ++i) {
 246             polyhedrons[i] = parsePolyhedron();
 247             /*polyhedrons[i].rotateXY(0.54326);
 248             polyhedrons[i].swapXZ();
 249             polyhedrons[i].rotateXY(0.1324543);
 250             polyhedrons[i].swapXZ();
 251             polyhedrons[i].rotateXY(0.5416512);*/
 252         }
 253         double res = 0.0;
 254         for (Polyhedron polyhedron : polyhedrons) {
 255             double volume = polyhedron.integrate(new double[]{1});
 256             double xc = polyhedron.integrate(new double[]{0, 1}) / volume;
 257             polyhedron.swapXY();
 258             double yc = polyhedron.integrate(new double[]{0, 1}) / volume;
 259             polyhedron.swapXY();
 260             polyhedron.swapXZ();
 261             double zc = polyhedron.integrate(new double[]{0, 1}) / volume;
 262             polyhedron.swapXZ();
 263             System.out.println("Center of mass: (" + xc + ", " + yc + ", " + zc + ")");
 264             res += polyhedron.getDistance(new Point(xc, yc, zc));
 265         }
 266         writer.println(res);
 267     }
 268
 269
 270     public static void main(String[] args) throws InterruptedException {
 271         Thread thread = new Thread(new asteroids_petr());
 272         thread.start();
 273         thread.join();
 274         System.err.println(count);
 275     }
 276
 277     static final String TASK_ID = "asteroids";
 278     static final String IN_FILE = TASK_ID + ".in";
 279     static final String OUT_FILE = TASK_ID + ".out";
 280     BufferedReader reader;
 281     StringTokenizer tokenizer;
 282     PrintWriter writer;
 283
 284     public void run() {
 285         try {
 286             reader = new BufferedReader(new FileReader(IN_FILE));
 287             tokenizer = null;
 288             writer = new PrintWriter(OUT_FILE);
 289             solve();
 290             reader.close();
 291             writer.close();
 292         } catch (Exception e) {
 293             e.printStackTrace();
 294             System.exit(1);
 295         }
 296     }
 297
 298     int nextInt() throws IOException {
 299         return Integer.parseInt(nextToken());
 300     }
 301
 302     long nextLong() throws IOException {
 303         return Long.parseLong(nextToken());
 304     }
 305
 306     double nextDouble() throws IOException {
 307         return Double.parseDouble(nextToken());
 308     }
 309
 310     String nextToken() throws IOException {
 311         while (tokenizer == null || !tokenizer.hasMoreTokens()) {
 312             tokenizer = new StringTokenizer(reader.readLine());
 313         }
 314         return tokenizer.nextToken();
 315     }
 316 }