「ライブラリ」の編集履歴（バックアップ）一覧はこちら

「ライブラリ」(2006/06/04 (日) 03:40:03) の最新版変更点

追加された行は緑色になります。

削除された行は赤色になります。

*平面幾何計算
 class Double2D{
     double x;
     double y;
     static final double EPS = 1.0E-9;
     Double2D(double x,double y){
         this.x = x;
         this.y = y;
     }
     @Override
     public boolean equals(Object o){
         Double2D p = (Double2D)o;
         return x==p.x&&y==p.y;
     }
     public double squareDist(Double2D p){
         return (x-p.x)*(x-p.x) + (y-p.y)*(y-p.y);
     }
     public double dist(Double2D p){
         return Math.sqrt(this.squareDist(p));
     }
     // thisとpを通る半径rの円の中心
     public Double2D[] centerOfPassingCircle1(Double2D p,double r){
         if(this.dist(p)>r) return null; // 境界条件はテキトー
         double a = p.x - x;
         double b = p.y - y;
         double d = this.dist(p);
         double qx0 = x / 2 + p.x / 2 + Math.sqrt(r - d*d/4)*(-b)/d;
         double qy0 = y / 2 + p.y / 2 + Math.sqrt(r - d*d/4)*a/d;
         double qx1 = x / 2 + p.x / 2 - Math.sqrt(r - d*d/4)*(-b)/d;
         double qy1 = y / 2 + p.y / 2 - Math.sqrt(r - d*d/4)*a/d;
         Double2D[] ret = new Double2D[2];
         ret[0] = new Double2D(qx0,qy0);
         ret[1] = new Double2D(qx1,qy1);
         return ret;
     }
     // 三点の作る三角形の面積
     public static double squareOfTriangle(Double2D p1,Double2D p2,Double2D p3){
         double x1 = p1.x-p2.x;
         double x2 = p3.x-p2.x;
         double y1 = p1.y-p2.y;
         double y2 = p3.y-p2.y;
         double s = (x1*y2-x2*y1)/2;
         return Math.abs(s);
     }
     // 多角形の内部にあるかどうか 自己ループはないものとする
     public static boolean inPolygon(ArrayList<Double2D> poly){
         double acc = 0.0;
         int n = poly.size();
         for(int i=0;i<n-1;i++){
             acc += arg(poly.get(i),poly.get(i+1));
         }
         acc += arg(poly.get(n-1),poly.get(0));
         
         if(Math.abs(acc)<EPS) return false;
         else return true;
     }
     // -Pi < theta < Pi の範囲で偏角を返す　複素数でいう Arg(p2/p1)
     public static double arg(Double2D p1,Double2D p2){
         double ax = p1.x;
         double ay = p1.y;
         double bx = p2.x;
         double by = p2.y;
         double theta = Math.atan2((ax*by-ay*bx)/(ax*ax+ay*ay),(ax*bx+ay*by)/(ax*ax+ay*ay));
         return theta;
     }
     // 外接円の中心 Online Judge 検証なし 簡易検証あり
     public static Double2D centerOfCircumscribingCircle(Double2D p1, Double2D p2, Double2D p3){
         double x1 = p1.x, y1 = p1.y, x2 = p2.x, y2 = p2.y, x3 = p3.x, y3 = p3.y;
         double det = (x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2));
         if(Math.abs(det)<EPS) return null;
         double sx = 
             -y1*(x2*x2+y2*y2-x3*x3-y3*y3)
             -y2*(x3*x3+y3*y3-x1*x1-y1*y1)
             -y3*(x1*x1+y1*y1-x2*x2-y2*y2);
         double sy =
              x1*(x2*x2+y2*y2-x3*x3-y3*y3)
             +x2*(x3*x3+y3*y3-x1*x1-y1*y1)
             +x3*(x1*x1+y1*y1-x2*x2-y2*y2);
         return new Double2D(0.5*sx/det, 0.5*sy/det);
     }
 }
 
 class Line{
     double a;
     double b;
     double c;
     static final double EPS = 1.0E-9;
     Line(double a, double b, double c) {
         this.a = a;
         this.b = b;
         this.c = c;
     }
     Line(Double2D p1,Double2D p2){
         this.a = p1.y-p2.y;
         this.b = -p1.x+p2.x;
         this.c = -p1.x*(p1.y-p2.y)+p1.y*(p1.x-p2.x);
     }
     // 線分と線分の交点 ない場合は null
     // 2線分が重なりをもつときは考えていない
     public Double2D crossPoint(Line l){
         // 平行
         if(Math.abs(this.a*l.b-this.b*l.a)<EPS){
             return null;
         }
         else{
             double x = (-l.b*this.c+this.b*l.c)/(this.a*l.b-this.b*l.a);
             double y = (l.a*this.c-this.a*l.c)/(this.a*l.b-this.b*l.a);
             return new Double2D(x,y);            
         }
     }
 }
 
 class Circle{
     double x;
     double y;
     double r;
     static final double EPS = 1.0E-9;
     Circle(double x, double y, double r) {
         this.x = x;
         this.y = y;
         this.r = r;
     }
     public double centerDist(Circle c){
         return Math.sqrt((x-c.x)*(x-c.x)+(y-c.y)*(y-c.y));
     }
     // 2円の交点
     public Double2D[] crossPoint1(Circle c){
         double d = centerDist(c);
         if(r+c.r > d) return null; // 境界条件はテキトー
         double d1 = (d*d + r*r - c.r*c.r)/(2*d);
         double ex = x + (c.x-x)*d1/d;
         double ey = y + (c.y-y)*d1/d;
         double l = Math.sqrt(r*r-d1*d1);
         double lx = -(c.y-y)*l/d;
         double ly = (c.x-x)*l/d;
         Double2D[] ret = new Double2D[2];
         ret[0] = new Double2D(ex+lx, ey+ly);
         ret[1] = new Double2D(ex-lx, ey-ly);
         return ret;
     }
     // 1点を通る接線
     public Line[] tangentLine(Double2D p){
         Line[] ret = new Line[2];
         if(Math.abs((x-p.x)*(x-p.x)-r*r)<EPS){
             ret[0] = new Line(1, 0, -p.x);
             double m = ((y-p.y)*(y-p.y)-r*r)/(2*(x-p.x)*(y-p.y));
             ret[1] = new Line(m, -1, -m*p.x+p.y);
         }
         else{
             double m1 = ((x-p.x)*(y-p.y)+r*Math.sqrt((x-p.x)*(x-p.x)+(y-p.y)*(y-p.y)-r*r))/((x-p.x)*(x-p.x)-r*r);
             double m2 = ((x-p.x)*(y-p.y)-r*Math.sqrt((x-p.x)*(x-p.x)+(y-p.y)*(y-p.y)-r*r))/((x-p.x)*(x-p.x)-r*r);
             ret[0] = new Line(m1, -1, -m1*p.x+p.y);
             ret[1] = new Line(m2, -1, -m2*p.x+p.y);
         }
         return ret;
     }
 }
 
 class LineSegment{
     Double2D p1;
     Double2D p2;
     static final double EPS = 1.0E-9;
     LineSegment(double x1, double y1, double x2, double y2){
         p1 = new Double2D(x1, y1);
         p2 = new Double2D(x2, y2);
     }
     // 2線分が重なりをもつときは考えていない
     public Double2D crossPoint(LineSegment ls){
         Double2D p3 = ls.p1;
         Double2D p4 = ls.p2;
         if(fcompare(Math.max(p1.x, p2.x), Math.min(p3.x, p4.x))<0 ||
            fcompare(Math.max(p3.x, p4.x), Math.min(p1.x, p2.x))<0 ||
            fcompare(Math.max(p1.y, p2.y), Math.min(p3.y, p4.y))<0 ||
            fcompare(Math.max(p3.y, p4.y), Math.min(p1.y, p2.y))<0) return null;
         
         double dx00 = p2.x - p1.x;
         double dy00 = p2.y - p1.y;
         double dx01 = p3.x - p1.x;
         double dy01 = p3.y - p1.y;
         double dx02 = p4.x - p1.x;
         double dy02 = p4.y - p1.y;
         
         double dx10 = p4.x - p3.x;
         double dy10 = p4.y - p3.y;
         double dx11 = p1.x - p3.x;
         double dy11 = p1.y - p3.y;
         double dx12 = p2.x - p3.x;
         double dy12 = p2.y - p3.y;
 
         // 平行
         if(Math.abs(dx00*dy10-dy00*dx10)<EPS){
             if(p1.x==p3.x&&p1.y==p3.y) return p1;
             if(p1.x==p4.x&&p1.y==p4.y) return p1;
             if(p2.x==p3.x&&p2.y==p3.y) return p2;
             if(p2.x==p4.x&&p2.y==p4.y) return p2;
         }
         
         if(fcompare((dx00*dy01-dx01*dy00)*(dx00*dy02-dx02*dy00), 0)<=0 &&
            fcompare((dx10*dy11-dx11*dy10)*(dx10*dy12-dx12*dy10), 0)<=0){
             Line l1 = new Line(p1, p2);
             Line l2 = new Line(p3, p4);
             return l1.crossPoint(l2);
         }
         else return null;
     }
     private int fcompare(double a, double b){
         if(Math.abs(a-b)<EPS) return 0;
         if(a>b) return 1;
         else return -1;
     }
 }
*空間幾何計算
 class Double3D{
     double x;
     double y;
     double z;
     Double3D(double x, double y, double z) {
         this.x = x;
         this.y = y;
         this.z = z;
     }
     public double dist(Double3D p){
         double dx = p.x-this.x;
         double dy = p.y-this.y;
         double dz = p.z-this.z;
         return Math.sqrt(dx*dx+dy*dy+dz*dz);
     }
     //determinant of matrix |a1 a2 a3| where a1,a2,a3 are column vectors
     public static double det(Double3D a1,Double3D a2,Double3D a3){
         double d1 = a1.x*(a2.y*a3.z-a2.z*a3.y);
         double d2 = a1.y*(a2.x*a3.z-a2.z*a3.x);
         double d3 = a1.z*(a2.x*a3.y-a2.y*a3.x);
         return d1-d2+d3;
     }
     /* solution of linear equation
      * where coefficient matrix is (a1 a2 a3) and c is constant vector
      * when the solution is not unique, return null
      */
     public static Double3D solveLinearEquation(Double3D a1,Double3D a2,Double3D a3,Double3D c){
         double d = det(a1,a2,a3);
         if(Math.abs(d)<1E-10) return null;
         double o1 = det(c,a2,a3);
         double o2 = det(a1,c,a3);
         double o3 = det(a1,a2,c);
         return new Double3D(o1/d,o2/d,o3/d);
     }
     /* center of circumcircle of Triangle pqr
      * when p,q,r are on a line, return null
      */
     public static Double3D centerOfCircumCircle(Double3D p,Double3D q,Double3D r){
         double a11 = 2*(p.x-q.x);
         double a12 = 2*(p.y-q.y);
         double a13 = 2*(p.z-q.z);
         double c1 = p.x*p.x+p.y*p.y+p.z*p.z-q.x*q.x-q.y*q.y-q.z*q.z;
         double a21 = 2*(p.x-r.x);
         double a22 = 2*(p.y-r.y);
         double a23 = 2*(p.z-r.z);
         double c2 = p.x*p.x+p.y*p.y+p.z*p.z-r.x*r.x-r.y*r.y-r.z*r.z;
         double a31 = p.y*q.z+q.y*r.z+r.y*p.z-p.z*q.y-q.z*r.y-r.z*p.y;
         double a32 = p.z*q.x+q.z*r.x+r.z*p.x-p.x*q.z-q.x*r.z-r.x*p.z;
         double a33 = p.x*q.y+q.x*r.y+r.x*p.y-p.y*q.x-q.y*r.x-r.y*p.x;
         double c3 = det(p,q,r);
         Double3D a1 = new Double3D(a11,a21,a31);
         Double3D a2 = new Double3D(a12,a22,a32);
         Double3D a3 = new Double3D(a13,a23,a33);
         Double3D c = new Double3D(c1,c2,c3);
         return solveLinearEquation(a1,a2,a3,c);
     }
     /* center of circumsphere of tetrahedron pqrs
      * when p,q,r,s are on a plane, return null
      */
     public static Double3D centerOfCircumSphere(Double3D p,Double3D q,Double3D r,Double3D s){
         double a11 = 2*(p.x-q.x);
         double a12 = 2*(p.y-q.y);
         double a13 = 2*(p.z-q.z);
         double c1 = p.x*p.x+p.y*p.y+p.z*p.z-q.x*q.x-q.y*q.y-q.z*q.z;
         double a21 = 2*(p.x-r.x);
         double a22 = 2*(p.y-r.y);
         double a23 = 2*(p.z-r.z);
         double c2 = p.x*p.x+p.y*p.y+p.z*p.z-r.x*r.x-r.y*r.y-r.z*r.z;
         double a31 = 2*(p.x-s.x);
         double a32 = 2*(p.y-s.y);
         double a33 = 2*(p.z-s.z);
         double c3 = p.x*p.x+p.y*p.y+p.z*p.z-s.x*s.x-s.y*s.y-s.z*s.z;
         Double3D a1 = new Double3D(a11,a21,a31);
         Double3D a2 = new Double3D(a12,a22,a32);
         Double3D a3 = new Double3D(a13,a23,a33);
         Double3D c = new Double3D(c1,c2,c3);
         return solveLinearEquation(a1,a2,a3,c);
     }
     public static Double3D midpoint(Double3D p,Double3D q){
         return new Double3D((p.x+q.x)/2,(p.y+q.y)/2,(p.z+q.z)/2);
     }
 }

*高速Kruskal with 高速MFSet
O(NlogN)のMFsetを用いたKruskalのアルゴリズム．
 class KruskalSolver{
     int n;
     double[][] dist;
     double len;
     private PriorityQueue<Edge> q;
     KruskalSolver(int n,double dist[][]){
         this.n = n;
         this.dist = dist;
         this.q = new PriorityQueue<Edge>();
         this.len = 0.0;
         
         for(int i=0;i<n;i++){
             for(int j=i+1;j<n;j++){
                 if(dist[i][j]<Double.MAX_VALUE){
                     q.offer(new Edge(i, j, dist[i][j]));
                 }
             }
         }
         
         //Kruskal's algorithm
         MFSet mfs = new MFSet(n);
         while(!q.isEmpty()){
             Edge e = q.poll();
             int i = e.p1;
             int j = e.p2;
             if(mfs.find(i)!=mfs.find(j)){
                 mfs.merge(mfs.find(i),mfs.find(j));
                 len += e.length;
             }
         }
     }
 }
 
 class Edge implements Comparable<Edge>{
     int p1;
     int p2;
     double length;
     Edge(int i,int j,double l){
         p1 = i;
         p2 = j;
         length = l;
     }
     public int compareTo(Edge e){
         if(length>e.length) return 1;
         else if(length==e.length) return 0;
         else return -1;
     }
 }
 
 class MFSet{
     int n;
     private int[] setsize;
     private int[] firstelem;
     private int[] set;
     private int[] next;
     MFSet(int n){
         this.n = n;
         setsize = new int[n];
         Arrays.fill(setsize, 1);
         firstelem = new int[n];
         for(int i=0;i<n;i++) firstelem[i] = i;
         set = new int[n];
         for(int i=0;i<n;i++) set[i] = i;
         next = new int[n];
         for(int i=0;i<n;i++) next[i] = -1;
     }
     int find(int x){
         return set[x];
     }
     void merge(int a, int b){
         if(setsize[a]<setsize[b]){
             merge(b, a);
             return;
         }
         int i = firstelem[b];
         while(next[i]>=0){
             set[i] = a;
             i = next[i];
         }
         set[i] = a;
         next[i] = firstelem[a];
         firstelem[a] = firstelem[b];
         setsize[a] = setsize[a] + setsize[b];
         setsize[b] = 0;
         firstelem[b] = -1;
     }
 }
*Kruskal
O(N^2)のMFSetを用いたKruskalのアルゴリズム．
 class KruskalSolver{
     int n;
     double[][] dist;
     double len;
     private PriorityQueue<Edge> q;
     int[] set;
     KruskalSolver(int n,double dist[][]){
         this.n = n;
         this.dist = dist;
         this.len = 0.0;
         this.q = new PriorityQueue<Edge>();
         this.set = new int[n];
         for(int i=0;i<n;i++) set[i] = i;
         
         for(int i=0;i<n;i++){
             for(int j=i+1;j<n;j++){
                 if(dist[i][j]<Double.MAX_VALUE){
                     q.offer(new Edge(i, j, dist[i][j]));
                 }
             }
         }
         
         //Kruskal's algorithm
         while(!q.isEmpty()){
             Edge e = q.poll();
             int i = e.p1;
             int j = e.p2;
             int si = set[i];
             int sj = set[j];
             if(si!=sj){
                 len += e.length;
                 for(int k=0;k<n;k++){
                     if(set[k]==sj) set[k] = si;
                 }
             }
         }
     }
 }
 
 class Edge implements Comparable<Edge>{
     int p1;
     int p2;
     double length;
     Edge(int i,int j,double l){
         p1 = i;
         p2 = j;
         length = l;
     }
     public int compareTo(Edge e){
         if(length>e.length) return 1;
         else if(length==e.length) return 0;
         else return -1;
     }
 }
*Floyd
全点対間最短路．
 class FloydSolver{
     int n;
     double[][] cost;
     int[][] path;
     FloydSolver(int n, double[][] dist) {
         this.n = n;
         cost = new double[n][n];
         path = new int[n][n];
         for(int i=0;i<n;i++){
             for(int j=0;j<n;j++){
                 cost[i][j] = dist[i][j];
                 path[i][j] = -1;
             }
         }
         for(int k=0;k<n;k++){
             for(int i=0;i<n;i++){
                 for(int j=0;j<n;j++){
                     if(cost[i][k]==Double.MAX_VALUE ||
                        cost[k][j]==Double.MAX_VALUE) continue;
                     if(cost[i][k]+cost[k][j]<cost[i][j]){
                         cost[i][j] = cost[i][k] + cost[k][j];
                         path[i][j] = k;
                     }
                 }
             }
         }
     }
     public double pathLength(int v1, int v2){
         return cost[v1][v2];
     }
     public ArrayList<Integer> path(int v1, int v2){
         if(cost[v1][v2]==Double.MAX_VALUE) return null;
         ArrayList<Integer> ret = new ArrayList<Integer>();
         ret.add(v1);
         path(v1, v2, ret);
         ret.add(v2);
         return ret;
     }
     private void path(int v1, int v2, ArrayList<Integer> ret){
         int k = path[v1][v2];
         if(k==-1) return;
         else{
             path(v1, k, ret);
             ret.add(k);
             path(k, v2, ret);
         }
     }
 }

*平面幾何計算
 class Double2D{
     double x;
     double y;
     static final double EPS = 1.0E-9;
     Double2D(double x,double y){
         this.x = x;
         this.y = y;
     }
     @Override
     public boolean equals(Object o){
         Double2D p = (Double2D)o;
         return x==p.x&&y==p.y;
     }
     public double squareDist(Double2D p){
         return (x-p.x)*(x-p.x) + (y-p.y)*(y-p.y);
     }
     public double dist(Double2D p){
         return Math.sqrt(this.squareDist(p));
     }
     // thisとpを通る半径rの円の中心
     public Double2D[] centerOfPassingCircle1(Double2D p,double r){
         if(this.dist(p)>r) return null; // 境界条件はテキトー
         double a = p.x - x;
         double b = p.y - y;
         double d = this.dist(p);
         double qx0 = x / 2 + p.x / 2 + Math.sqrt(r - d*d/4)*(-b)/d;
         double qy0 = y / 2 + p.y / 2 + Math.sqrt(r - d*d/4)*a/d;
         double qx1 = x / 2 + p.x / 2 - Math.sqrt(r - d*d/4)*(-b)/d;
         double qy1 = y / 2 + p.y / 2 - Math.sqrt(r - d*d/4)*a/d;
         Double2D[] ret = new Double2D[2];
         ret[0] = new Double2D(qx0,qy0);
         ret[1] = new Double2D(qx1,qy1);
         return ret;
     }
     // 三点の作る三角形の面積
     public static double squareOfTriangle(Double2D p1,Double2D p2,Double2D p3){
         double x1 = p1.x-p2.x;
         double x2 = p3.x-p2.x;
         double y1 = p1.y-p2.y;
         double y2 = p3.y-p2.y;
         double s = (x1*y2-x2*y1)/2;
         return Math.abs(s);
     }
     // 多角形の内部にあるかどうか 自己ループはないものとする
     public static boolean inPolygon(ArrayList<Double2D> poly){
         double acc = 0.0;
         int n = poly.size();
         for(int i=0;i<n-1;i++){
             acc += arg(poly.get(i),poly.get(i+1));
         }
         acc += arg(poly.get(n-1),poly.get(0));
         
         if(Math.abs(acc)<EPS) return false;
         else return true;
     }
     // -Pi < theta < Pi の範囲で偏角を返す　複素数でいう Arg(p2/p1)
     public static double arg(Double2D p1,Double2D p2){
         double ax = p1.x;
         double ay = p1.y;
         double bx = p2.x;
         double by = p2.y;
         double theta = Math.atan2((ax*by-ay*bx)/(ax*ax+ay*ay),(ax*bx+ay*by)/(ax*ax+ay*ay));
         return theta;
     }
     // 外接円の中心 Online Judge 検証なし 簡易検証あり
     public static Double2D centerOfCircumscribingCircle(Double2D p1, Double2D p2, Double2D p3){
         double x1 = p1.x, y1 = p1.y, x2 = p2.x, y2 = p2.y, x3 = p3.x, y3 = p3.y;
         double det = (x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2));
         if(Math.abs(det)<EPS) return null;
         double sx = 
             -y1*(x2*x2+y2*y2-x3*x3-y3*y3)
             -y2*(x3*x3+y3*y3-x1*x1-y1*y1)
             -y3*(x1*x1+y1*y1-x2*x2-y2*y2);
         double sy =
              x1*(x2*x2+y2*y2-x3*x3-y3*y3)
             +x2*(x3*x3+y3*y3-x1*x1-y1*y1)
             +x3*(x1*x1+y1*y1-x2*x2-y2*y2);
         return new Double2D(0.5*sx/det, 0.5*sy/det);
     }
 }
 
 class Line{
     double a;
     double b;
     double c;
     static final double EPS = 1.0E-9;
     Line(double a, double b, double c) {
         this.a = a;
         this.b = b;
         this.c = c;
     }
     Line(Double2D p1,Double2D p2){
         this.a = p1.y-p2.y;
         this.b = -p1.x+p2.x;
         this.c = -p1.x*(p1.y-p2.y)+p1.y*(p1.x-p2.x);
     }
     // 線分と線分の交点 ない場合は null
     // 2線分が重なりをもつときは考えていない
     public Double2D crossPoint(Line l){
         // 平行
         if(Math.abs(this.a*l.b-this.b*l.a)<EPS){
             return null;
         }
         else{
             double x = (-l.b*this.c+this.b*l.c)/(this.a*l.b-this.b*l.a);
             double y = (l.a*this.c-this.a*l.c)/(this.a*l.b-this.b*l.a);
             return new Double2D(x,y);            
         }
     }
 }
 
 class Circle{
     double x;
     double y;
     double r;
     static final double EPS = 1.0E-9;
     Circle(double x, double y, double r) {
         this.x = x;
         this.y = y;
         this.r = r;
     }
     public double centerDist(Circle c){
         return Math.sqrt((x-c.x)*(x-c.x)+(y-c.y)*(y-c.y));
     }
     // 2円の交点
     public Double2D[] crossPoint1(Circle c){
         double d = centerDist(c);
         if(r+c.r > d) return null; // 境界条件はテキトー
         double d1 = (d*d + r*r - c.r*c.r)/(2*d);
         double ex = x + (c.x-x)*d1/d;
         double ey = y + (c.y-y)*d1/d;
         double l = Math.sqrt(r*r-d1*d1);
         double lx = -(c.y-y)*l/d;
         double ly = (c.x-x)*l/d;
         Double2D[] ret = new Double2D[2];
         ret[0] = new Double2D(ex+lx, ey+ly);
         ret[1] = new Double2D(ex-lx, ey-ly);
         return ret;
     }
     // 1点を通る接線
     public Line[] tangentLine(Double2D p){
         Line[] ret = new Line[2];
         if(Math.abs((x-p.x)*(x-p.x)-r*r)<EPS){
             ret[0] = new Line(1, 0, -p.x);
             double m = ((y-p.y)*(y-p.y)-r*r)/(2*(x-p.x)*(y-p.y));
             ret[1] = new Line(m, -1, -m*p.x+p.y);
         }
         else{
             double m1 = ((x-p.x)*(y-p.y)+r*Math.sqrt((x-p.x)*(x-p.x)+(y-p.y)*(y-p.y)-r*r))/((x-p.x)*(x-p.x)-r*r);
             double m2 = ((x-p.x)*(y-p.y)-r*Math.sqrt((x-p.x)*(x-p.x)+(y-p.y)*(y-p.y)-r*r))/((x-p.x)*(x-p.x)-r*r);
             ret[0] = new Line(m1, -1, -m1*p.x+p.y);
             ret[1] = new Line(m2, -1, -m2*p.x+p.y);
         }
         return ret;
     }
 }
 
 class LineSegment{
     Double2D p1;
     Double2D p2;
     static final double EPS = 1.0E-9;
     LineSegment(double x1, double y1, double x2, double y2){
         p1 = new Double2D(x1, y1);
         p2 = new Double2D(x2, y2);
     }
     // 2線分が重なりをもつときは考えていない
     public Double2D crossPoint(LineSegment ls){
         Double2D p3 = ls.p1;
         Double2D p4 = ls.p2;
         if(fcompare(Math.max(p1.x, p2.x), Math.min(p3.x, p4.x))<0 ||
            fcompare(Math.max(p3.x, p4.x), Math.min(p1.x, p2.x))<0 ||
            fcompare(Math.max(p1.y, p2.y), Math.min(p3.y, p4.y))<0 ||
            fcompare(Math.max(p3.y, p4.y), Math.min(p1.y, p2.y))<0) return null;
         
         double dx00 = p2.x - p1.x;
         double dy00 = p2.y - p1.y;
         double dx01 = p3.x - p1.x;
         double dy01 = p3.y - p1.y;
         double dx02 = p4.x - p1.x;
         double dy02 = p4.y - p1.y;
         
         double dx10 = p4.x - p3.x;
         double dy10 = p4.y - p3.y;
         double dx11 = p1.x - p3.x;
         double dy11 = p1.y - p3.y;
         double dx12 = p2.x - p3.x;
         double dy12 = p2.y - p3.y;
 
         // 平行
         if(Math.abs(dx00*dy10-dy00*dx10)<EPS){
             if(p1.x==p3.x&&p1.y==p3.y) return p1;
             if(p1.x==p4.x&&p1.y==p4.y) return p1;
             if(p2.x==p3.x&&p2.y==p3.y) return p2;
             if(p2.x==p4.x&&p2.y==p4.y) return p2;
         }
         
         if(fcompare((dx00*dy01-dx01*dy00)*(dx00*dy02-dx02*dy00), 0)<=0 &&
            fcompare((dx10*dy11-dx11*dy10)*(dx10*dy12-dx12*dy10), 0)<=0){
             Line l1 = new Line(p1, p2);
             Line l2 = new Line(p3, p4);
             return l1.crossPoint(l2);
         }
         else return null;
     }
     private int fcompare(double a, double b){
         if(Math.abs(a-b)<EPS) return 0;
         if(a>b) return 1;
         else return -1;
     }
 }
*空間幾何計算
 class Double3D{
     double x;
     double y;
     double z;
     Double3D(double x, double y, double z) {
         this.x = x;
         this.y = y;
         this.z = z;
     }
     public double dist(Double3D p){
         double dx = p.x-this.x;
         double dy = p.y-this.y;
         double dz = p.z-this.z;
         return Math.sqrt(dx*dx+dy*dy+dz*dz);
     }
     //determinant of matrix |a1 a2 a3| where a1,a2,a3 are column vectors
     public static double det(Double3D a1,Double3D a2,Double3D a3){
         double d1 = a1.x*(a2.y*a3.z-a2.z*a3.y);
         double d2 = a1.y*(a2.x*a3.z-a2.z*a3.x);
         double d3 = a1.z*(a2.x*a3.y-a2.y*a3.x);
         return d1-d2+d3;
     }
     /* solution of linear equation
      * where coefficient matrix is (a1 a2 a3) and c is constant vector
      * when the solution is not unique, return null
      */
     public static Double3D solveLinearEquation(Double3D a1,Double3D a2,Double3D a3,Double3D c){
         double d = det(a1,a2,a3);
         if(Math.abs(d)<1E-10) return null;
         double o1 = det(c,a2,a3);
         double o2 = det(a1,c,a3);
         double o3 = det(a1,a2,c);
         return new Double3D(o1/d,o2/d,o3/d);
     }
     /* center of circumcircle of Triangle pqr
      * when p,q,r are on a line, return null
      */
     public static Double3D centerOfCircumCircle(Double3D p,Double3D q,Double3D r){
         double a11 = 2*(p.x-q.x);
         double a12 = 2*(p.y-q.y);
         double a13 = 2*(p.z-q.z);
         double c1 = p.x*p.x+p.y*p.y+p.z*p.z-q.x*q.x-q.y*q.y-q.z*q.z;
         double a21 = 2*(p.x-r.x);
         double a22 = 2*(p.y-r.y);
         double a23 = 2*(p.z-r.z);
         double c2 = p.x*p.x+p.y*p.y+p.z*p.z-r.x*r.x-r.y*r.y-r.z*r.z;
         double a31 = p.y*q.z+q.y*r.z+r.y*p.z-p.z*q.y-q.z*r.y-r.z*p.y;
         double a32 = p.z*q.x+q.z*r.x+r.z*p.x-p.x*q.z-q.x*r.z-r.x*p.z;
         double a33 = p.x*q.y+q.x*r.y+r.x*p.y-p.y*q.x-q.y*r.x-r.y*p.x;
         double c3 = det(p,q,r);
         Double3D a1 = new Double3D(a11,a21,a31);
         Double3D a2 = new Double3D(a12,a22,a32);
         Double3D a3 = new Double3D(a13,a23,a33);
         Double3D c = new Double3D(c1,c2,c3);
         return solveLinearEquation(a1,a2,a3,c);
     }
     /* center of circumsphere of tetrahedron pqrs
      * when p,q,r,s are on a plane, return null
      */
     public static Double3D centerOfCircumSphere(Double3D p,Double3D q,Double3D r,Double3D s){
         double a11 = 2*(p.x-q.x);
         double a12 = 2*(p.y-q.y);
         double a13 = 2*(p.z-q.z);
         double c1 = p.x*p.x+p.y*p.y+p.z*p.z-q.x*q.x-q.y*q.y-q.z*q.z;
         double a21 = 2*(p.x-r.x);
         double a22 = 2*(p.y-r.y);
         double a23 = 2*(p.z-r.z);
         double c2 = p.x*p.x+p.y*p.y+p.z*p.z-r.x*r.x-r.y*r.y-r.z*r.z;
         double a31 = 2*(p.x-s.x);
         double a32 = 2*(p.y-s.y);
         double a33 = 2*(p.z-s.z);
         double c3 = p.x*p.x+p.y*p.y+p.z*p.z-s.x*s.x-s.y*s.y-s.z*s.z;
         Double3D a1 = new Double3D(a11,a21,a31);
         Double3D a2 = new Double3D(a12,a22,a32);
         Double3D a3 = new Double3D(a13,a23,a33);
         Double3D c = new Double3D(c1,c2,c3);
         return solveLinearEquation(a1,a2,a3,c);
     }
     public static Double3D midpoint(Double3D p,Double3D q){
         return new Double3D((p.x+q.x)/2,(p.y+q.y)/2,(p.z+q.z)/2);
     }
 }

*高速Kruskal with 高速MFSet
O(NlogN)のMFsetを用いたKruskalのアルゴリズム．
 class KruskalSolver{
     int n;
     double[][] dist;
     double len;
     private PriorityQueue<Edge> q;
     KruskalSolver(int n,double dist[][]){
         this.n = n;
         this.dist = dist;
         this.q = new PriorityQueue<Edge>();
         this.len = 0.0;
         
         for(int i=0;i<n;i++){
             for(int j=i+1;j<n;j++){
                 if(dist[i][j]<Double.MAX_VALUE){
                     q.offer(new Edge(i, j, dist[i][j]));
                 }
             }
         }
         
         //Kruskal's algorithm
         MFSet mfs = new MFSet(n);
         while(!q.isEmpty()){
             Edge e = q.poll();
             int i = e.p1;
             int j = e.p2;
             if(mfs.find(i)!=mfs.find(j)){
                 mfs.merge(mfs.find(i),mfs.find(j));
                 len += e.length;
             }
         }
     }
 }
 
 class Edge implements Comparable<Edge>{
     int p1;
     int p2;
     double length;
     Edge(int i,int j,double l){
         p1 = i;
         p2 = j;
         length = l;
     }
     public int compareTo(Edge e){
         if(length>e.length) return 1;
         else if(length==e.length) return 0;
         else return -1;
     }
 }
 
 class MFSet{
     int n;
     private int[] setsize;
     private int[] firstelem;
     private int[] set;
     private int[] next;
     MFSet(int n){
         this.n = n;
         setsize = new int[n];
         Arrays.fill(setsize, 1);
         firstelem = new int[n];
         for(int i=0;i<n;i++) firstelem[i] = i;
         set = new int[n];
         for(int i=0;i<n;i++) set[i] = i;
         next = new int[n];
         for(int i=0;i<n;i++) next[i] = -1;
     }
     int find(int x){
         return set[x];
     }
     void merge(int a, int b){
         if(setsize[a]<setsize[b]){
             merge(b, a);
             return;
         }
         int i = firstelem[b];
         while(next[i]>=0){
             set[i] = a;
             i = next[i];
         }
         set[i] = a;
         next[i] = firstelem[a];
         firstelem[a] = firstelem[b];
         setsize[a] = setsize[a] + setsize[b];
         setsize[b] = 0;
         firstelem[b] = -1;
     }
 }
*Kruskal
O(N^2)のMFSetを用いたKruskalのアルゴリズム．
 class KruskalSolver{
     int n;
     double[][] dist;
     double len;
     private PriorityQueue<Edge> q;
     int[] set;
     KruskalSolver(int n,double dist[][]){
         this.n = n;
         this.dist = dist;
         this.len = 0.0;
         this.q = new PriorityQueue<Edge>();
         this.set = new int[n];
         for(int i=0;i<n;i++) set[i] = i;
         
         for(int i=0;i<n;i++){
             for(int j=i+1;j<n;j++){
                 if(dist[i][j]<Double.MAX_VALUE){
                     q.offer(new Edge(i, j, dist[i][j]));
                 }
             }
         }
         
         //Kruskal's algorithm
         while(!q.isEmpty()){
             Edge e = q.poll();
             int i = e.p1;
             int j = e.p2;
             int si = set[i];
             int sj = set[j];
             if(si!=sj){
                 len += e.length;
                 for(int k=0;k<n;k++){
                     if(set[k]==sj) set[k] = si;
                 }
             }
         }
     }
 }
 
 class Edge implements Comparable<Edge>{
     int p1;
     int p2;
     double length;
     Edge(int i,int j,double l){
         p1 = i;
         p2 = j;
         length = l;
     }
     public int compareTo(Edge e){
         if(length>e.length) return 1;
         else if(length==e.length) return 0;
         else return -1;
     }
 }
*Floyd
全点対間最短路．
 class FloydSolver{
     int n;
     double[][] cost;
     int[][] path;
     FloydSolver(int n, double[][] dist) {
         this.n = n;
         cost = new double[n][n];
         path = new int[n][n];
         for(int i=0;i<n;i++){
             for(int j=0;j<n;j++){
                 cost[i][j] = dist[i][j];
                 path[i][j] = -1;
             }
         }
         for(int k=0;k<n;k++){
             for(int i=0;i<n;i++){
                 for(int j=0;j<n;j++){
                     if(cost[i][k]==Double.MAX_VALUE ||
                        cost[k][j]==Double.MAX_VALUE) continue;
                     if(cost[i][k]+cost[k][j]<cost[i][j]){
                         cost[i][j] = cost[i][k] + cost[k][j];
                         path[i][j] = k;
                     }
                 }
             }
         }
     }
     public double pathLength(int v1, int v2){
         return cost[v1][v2];
     }
     public ArrayList<Integer> path(int v1, int v2){
         if(cost[v1][v2]==Double.MAX_VALUE) return null;
         ArrayList<Integer> ret = new ArrayList<Integer>();
         ret.add(v1);
         path(v1, v2, ret);
         ret.add(v2);
         return ret;
     }
     private void path(int v1, int v2, ArrayList<Integer> ret){
         int k = path[v1][v2];
         if(k==-1) return;
         else{
             path(v1, k, ret);
             ret.add(k);
             path(k, v2, ret);
         }
     }
 }
*最大流問題
 class MaxFlowSolver {
     int size;//number of nodes
     int source;//nodeID of source
     int sink;//nodeID of sink
     int[][] network;
     int[] prev;
 
     MaxFlowSolver(int n, int s, int t, int[][] net) {
         size = n;
         source = s;
         sink = t;
         network = net;
         prev = new int[size];
     }
 
     public int solve() {
         int flow = 0;
         while (searchPath()) {
             int flowinc = Integer.MAX_VALUE;
 
             // calculate how much flow increases
             int v = sink;
             while (prev[v]<Integer.MAX_VALUE) {
                 flowinc = Math.min(flowinc, network[prev[v]][v]);
                 v = prev[v];
             }
 
             //flow increases
             flow += flowinc;
 
             //update network array
             v = sink;
             while(prev[v]<Integer.MAX_VALUE){
                 network[prev[v]][v] -= flowinc;
                 network[v][prev[v]] += flowinc;
                 v = prev[v];
             }
         }
 
         return flow;
     }
 
     //search shortest path from source to sink by BFS
     private boolean searchPath() {
         Arrays.fill(prev, -1);
 
         Queue<Integer> q = new LinkedList<Integer>();
         q.offer(source);
         prev[source] = Integer.MAX_VALUE;
 
         L: while (!q.isEmpty()) {
             int u = q.poll();
             for (int i = 0; i < size; i++) {
                 if (prev[i] < 0 && network[u][i] > 0) {
                     prev[i] = u;
 
                     //discoverd path to sink
                     if (i == sink) break L;
 
                     q.offer(i);
                 }
             }
         }
 
         return prev[sink]>=0;
     }
 }