AOJ-1146, POJ-3011: Secrets in Shadows (影の秘密)

ICPC AOJ PKU

keyword

幾何　平面走査法　誤差　数値計算　Java

問題概要

中心が整数座標 in [1,30]*[1,30]、半径1の円柱がある。影の幅が最も大きくなるときの値と最も小さくなるときの角度を誤差EPS=1e-10以下で求める問題。

解法

まず、角度を決めたときの影の幅は比較的簡単に求めることができる。

円の中心から軸Xへ下ろした垂線の足と原点間の距離は内積を使って簡単に計算できる。後はそこを中心とする長さ2の線分が被覆する長さを求めてやればよい。この部分の計算量はソートが必要なのでO(N*log N)。
で、問題は最大の幅を持つ角度を求めてやることである。単純にやろうとすると次のような方法が考えられる。区間[0,π]を適当に分割して最大値を求め、最大値をとる値の周辺をさらに分割して調べ…という方法である。幅の関数はどう考えても連続かつ区分的に滑らかだし、入力が非常に小さな整数座標しかないのでそこまで病的なテストケースは作れないだろうから上手くいきそうな気がする。で、その方針で書いたコードが次の様なの。

import java.util.Arrays;
import java.util.Scanner;
import static java.lang.Math.*;

public class Main{

	int N;
	final int MESH = 1200;
	final double diameter = 2.0;
	final double EPS = 1e-12;
	Point[] ps;

	void run(){
		Scanner in = new Scanner(System.in);
		for(;;){
			N = in.nextInt();
			if(N==0) return ;
			ps = new Point[N];
			for(int i=0; i<N; i++){
				double x = in.nextDouble(), y = in.nextDouble();
				ps[i] = new Point(x,y);
			}
			System.out.printf("%.15f\n", solve(true));
			System.out.printf("%.15f\n", solve(false));
		}
	}

	double calc(double theta){
		double[] xs = new double[N];
		for(int i=0; i<N; i++){
			xs[i] = ps[i].x*sin(theta) + ps[i].y*cos(theta);
		}
		return rangeMerge(xs);
	}

	double rangeMerge(double[] xs){
		Arrays.sort(xs);
		double ret = diameter;
		for(int i=0; i<N-1; i++){
			ret += diameter - max(0.0, (xs[i] + diameter) - xs[i+1]);
		}
		return ret;
	}

	double solve(boolean findMin){
		double lb = 0.0, ub = PI, bestValue = findMin?Double.MAX_VALUE:-Double.MAX_VALUE, center = 0.0;
		while(ub-lb>EPS){
			double width = (ub - lb)/MESH;
			for(int i=0; i<MESH; i++){
				double result = calc(lb + width*i);
				if(findMin ^ (result > bestValue)){
					center = lb + width*i;
					bestValue = result;
				}
			}
			lb = center - 2.0*width;
			ub = center + 2.0*width;
		}
		return center;
	}

	class Point{
		double x, y;
		Point(double _x, double _y){
			this.x = _x; this.y = _y;
		}
	}

	public static void main(String args[]){
		new Main().run();
	}
}

これで答えっぽいものが出るけど誤差が微妙に取れない。原因は丸め誤差。浮動小数点を用いた計算なんてx*sin(theta)+y*cos(theta)の部分と足し算する部分くらいなのにそれでも丸め誤差から逃れられない。JavaでBigDecimalを使えばいいような気もするけど三角関数が無い(GMPにすらない)。上のコードで丸め誤差を回避する様々なテクニック(値は小さい順に足すとか、区間を併合してから幅をとるとか)を試したけどいずれも駄目だった。
本質的な話として、最大値の周りでは微分係数が0に近い可能性があるので、f(x+ε)≒f(x)+f'(x)*εの右辺2項目が桁落ちしてしまう(なので、long doubleにする位じゃ全然精度が足りない)。なので、最大値は十分な精度で求められるけど、最大値を与える値の精度は良くならない(三分探索ってどうやってこれを解決してるんだっけ)。
探索っぽく解くのは丸め誤差から抜けられないのでお約束の平面走査法で考える。区間[0,π]における特別な点を全部列挙すれば、真の解の近傍でうろうろすることもない。角度を0からπまで動かしたとき、どのタイミングで最大値をとりうるかというと、

  • X軸上の2点がすれ違う。
  • X軸上の2点間の距離が直径と等しくなる。
  • 連続なところで、関数の微分係数が0になる。

の3種類しかない。
上2つは愚直に全探索して調べる。三角関数の方程式を解くだけ。
3つ目の「連続な所」というのは上2つで出てきた値に挟まれた部分。関数の微分係数とかいうと面倒くさいようだけど、点がX軸上に現れる順序と隣接する2点間の距離が分かればa*sin(Θ)+b*cos(Θ)=0の形に落としてやるので結局三角関数の方程式を解くだけ。後は高校数学の範囲。

感想

幾何に見せかけた数値計算に見せかけた高校数学の問題だった。任意精度小数計算で三角関数が見つからなくて絶望していたけどWolfram alpha先生が何とかしてくれた。

import java.util.Arrays;
import java.util.Comparator;
import java.util.HashSet;
import java.util.Scanner;

import static java.lang.Math.*;

public class Main{

	int N;
	final int MESH = 1200;
	final double diameter = 2.0;
	final double EPS = 1e-13;
	Point[] ps;

	void run(){
		Scanner in = new Scanner(System.in);
		for(;;){
			N = in.nextInt();
			if(N==0) return ;
			ps = new Point[N];
			for(int i=0; i<N; i++){
				double x = in.nextDouble(), y = in.nextDouble();
				ps[i] = new Point(x,y);
			}
			solve();
		}
	}

	boolean inRange(double theta){
		return 0<=theta && theta<=PI;
	}

	double toRange(double theta){
		if(inRange(theta)){
			return theta;
		}
		else if(inRange(theta + PI)){
			return theta + PI;
		}
		else if(inRange(theta - PI)){
			return theta - PI;
		}
		else if(inRange(theta - 2*PI)){
			return theta - 2*PI;
		}
		else if(inRange(theta + 2*PI)){
			return theta + 2*PI;
		}
		return 0.0;
	}

	double[] findThreshold(){
		HashSet<Double> threshold = new HashSet<Double>();
		threshold.add(0.0);
		threshold.add(PI/2);
		threshold.add(PI);
		for(int i=0; i<N; i++){
			for(int j=0; j<N; j++)if(i!=j){
				double dy = ps[i].y - ps[j].y, dx = ps[i].x - ps[j].x;
				threshold.add(toRange(atan(dx/dy)));
				threshold.add(toRange(atan(-dx/dy)));
				threshold.add(toRange(atan(dy/dx)));
				threshold.add(toRange(atan(-dy/dx)));
			}
		}
		for(int i=0; i<N; i++){
			for(int j=0; j<N; j++)if(i!=j){
				double dy = ps[i].y - ps[j].y, dx = ps[i].x - ps[j].x;
				double a = hypot(dy, dx);
				if(a < diameter) continue;
				else{
					threshold.add(toRange(acos(diameter/a) + atan2(dy,dx)));
					threshold.add(toRange(acos(diameter/a) - atan2(dy,dx)));
					threshold.add(toRange(acos(diameter/a) + atan2(dx,dy)));
					threshold.add(toRange(acos(diameter/a) - atan2(dx,dy)));
					threshold.add(toRange(acos(-diameter/a) + atan2(dy,dx)));
					threshold.add(toRange(acos(-diameter/a) - atan2(dy,dx)));
					threshold.add(toRange(acos(-diameter/a) + atan2(dx,dy)));
					threshold.add(toRange(acos(-diameter/a) - atan2(dx,dy)));
				}
			}
		}
		Double[] tmp = threshold.toArray(new Double[0]);
		double[] tmpThreshold = new double[tmp.length];
		for(int i=0; i<tmp.length; i++){
			tmpThreshold[i] = tmp[i];
		}
		Arrays.sort(tmpThreshold);
		for(int i=0; i<tmpThreshold.length-1; i++){
			double theta = 0.5*(tmpThreshold[i+1] + tmpThreshold[i]);
			Pair[] vs = new Pair[N];
			for(int j=0; j<N; j++){
				vs[j] = new Pair(ps[j].x*sin(theta) + ps[j].y*cos(theta), j);
			}
			Arrays.sort(vs, new Cmp());
			double xx = 0, yy = 0;
			for(int j=1; j<N; j++){
				if(vs[j-1].value + diameter > vs[j].value){
					xx += ps[vs[j].index].x - ps[vs[j-1].index].x;
					yy += ps[vs[j].index].y - ps[vs[j-1].index].y;
				}
			}
			threshold.add(toRange(atan2(yy,xx)));
			threshold.add(toRange(atan2(xx,yy)));
			threshold.add(toRange(atan2(-yy,xx)));
			threshold.add(toRange(atan2(-xx,yy)));
		}
		if(false)for(int i=1; i<=30; i++){
			for(int j=1; j<=30; j++){
				threshold.add(atan2(i,j));
				threshold.add(toRange(atan2(-i,j)));
			}
		}
		tmp = threshold.toArray(new Double[0]);
		tmpThreshold = new double[tmp.length];
		for(int i=0; i<tmp.length; i++){
			tmpThreshold[i] = tmp[i];
		}
		Arrays.sort(tmpThreshold);
		return tmpThreshold;
	}

	double calc(double theta){
		double[] xs = new double[N];
		for(int i=0; i<N; i++){
			xs[i] = ps[i].x*sin(theta) + ps[i].y*cos(theta);
		}
		return rangeMerge(xs);
	}

	double rangeMerge(double[] xs){
		Arrays.sort(xs);
		double ret = diameter;
		for(int i=0; i<N-1; i++){
			ret += diameter - max(0.0, (xs[i] + diameter) - xs[i+1]);
		}
		return ret;
	}

	void solve(){
		double[] threshold = findThreshold();
		double minValue = Double.MAX_VALUE, maxValue = -Double.MAX_VALUE;
		double minAns = -1, maxAns = -1;
		for(double th:threshold){
			double result = calc(th);
			if(result > maxValue){
				maxValue = result;
				maxAns = th;
			}
			if(result < minValue){
				minValue = result;
				minAns = th;
			}
		}
		System.out.printf("%.13f\n", minAns + EPS);
		System.out.printf("%.13f\n", maxAns + EPS);
	}

	class Point{
		double x, y;
		Point(double _x, double _y){
			this.x = _x; this.y = _y;
		}
	}

	class Pair{
		double value;
		int index;
		Pair(double v, int i){
			this.value = v;
			this.index = i;
		}
	}

	class Cmp implements Comparator<Pair>{
		public int compare(Pair a, Pair b){
			double dif = a.value - b.value;
			return dif>0?1:dif<0?-1:0;
		}
	}

	public static void main(String args[]){
		new Main().run();
	}
}