読者です 読者をやめる 読者になる 読者になる

AOJ 2220 The Number of the Real Roots of a Cubic Equation

問題文: The Number of the Real Roots of a Cubic Equation | Aizu Online Judge

誤差と戦うゲー

正直x*x*xとか書くのが面倒くさいからと言って、pow(x,3)とか書いてるとかなり痛い目にあう。

あとsqrtとかあまりよくないから自分のsqrtつくるべし。

3回目のtryでやっと通った。

解法: 三次方程式には解の公式があるので使ってみたり、二分探索で1つ解を求めて、それで式を割ったものを二次方程式の解の公式使って解いたりいろいろ。ただどれも誤差恐い。

今回は2つの方法で解きました。

カルダノの公式(三次方程式の解の公式)用いたやつ

(参考:C言語による最新アルゴリズム辞典?)

// 2220
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <iostream>
using namespace std;
const double PI = 2 * asin(1.0);


double cuberoot(double x){
	double s , prev;
	int positive;
	if(x == 0)  return 0;
	if(x>0)positive = 1; else {positive = 0 ; x = -x;}
	if(x>1) s = x ; else s = 1;
	do{
		prev = s ;  s = (x/(s*s)+2*s)/3;
	}while(s < prev);
	if(positive) return prev; else return -prev;
}
int M , P;

void check(double s){
	if(abs(s) < 1e-10) return ;
	if(s > 0) P++;
	else M++;
}
void cardano(double a, double b , double c , double d){
	double p , q , t , a3 , b3 , x1 , x2 , x3;
	b/=(3*a);
	c/=a;
	d/=a;
	p = b * b - c / 3;
	q = (b * (c-2*b*b)-d)/2;
	a = q * q - p * p * p;
	if(a == 0){
		q = cuberoot(q); x1 = 2 * q - b ; x2 = -q - b;
		//printf("x = %g, %g(重解)\n",x1,x2);
		check(x1);
		check(x2);
		check(x2);
	}else if( a > 0){
		if(q>0) a3 = cuberoot(q+sqrt(a));
		else a3 = cuberoot(q-sqrt(a));
		b3 = p / a3;
		x1 = a3 + b3 - b;
		x2 = - 0.5 * (a3 + b3) - b;
		x3 = fabs(a3-b3) * sqrt(3.0) / 2;
		//printf("x = %g; %g +- %g i\n",x1,x2,x3);
		check(x1);
	}else{
		a = sqrt(p) ; t = acos( q / (p * a) ); a *= 2;
		x1 = a * cos(t/3)-b;
		x2 = a * cos((t+2*PI)/3)-b;
		x3 = a * cos((t+4*PI)/3)-b;
		//printf("x = %g , %g , %g\n",x1,x2,x3);
		check(x1);
		check(x2);
		check(x3);
	}
}
int main(){ 
	int n;
	cin >> n;
	for(int i = 0 ; i < n ; i++){
		double a , b , c , d;
		cin >> a >> b >> c >> d;
		M = P = 0;
		cardano(a,b,c,d);
		cout << P << " " << M << endl;
	}
}

二分探索で1つ解見つけてから、その得た解を使って整式を割ったりした解法。

#include <iostream>
#include <cmath>
using namespace std;
const double EPS = 1e-10;
int m[3] = {0};

int is(double a){
	if(abs(a) < EPS)return 1;
	if(a < 0)return 0;
	if(a > 0)return 2;
}

double mysqrt(double x){
	double l = 0 , r = x;
	for(int i = 0 ; i < 64 ; i++){
		double m = (l+r) / 2.0;
		if(m * m < x) l = m;
		else r = m;
	}
	return l;
}

double solve(double a,double b,double c,double d){
	//cout << a << " " << b << " " << c << " " << d << endl;
	if(a == 0){
		if(c * c - 4 * b * d < -EPS)return 0.0;
		double one = ( -c + mysqrt(c*c-4*b*d)  )/(2*b);
		double two = ( -c - mysqrt(c*c-4*b*d)  )/(2*b);
		//cout << one << "<>" << two << endl;
		m[is(one)]++;
		m[is(two)]++;
		return 0.0;
	}
	
	double l = -1000.0 , r = 1000.0;
	
	if(a < 0){ a *= -1; b *= -1; c *= -1; d *= -1;}
	
	for(int i=0;i<64;i++){
		double x = (l+r)/2.0;
		double res = a*x*x*x+b*x*x+c*x+d;
		if(res>0){
			r = x;
		}else{
			l = x;
		}
	}
	
	return l;
}
int main(){
	int N; cin >> N;
	while(N--){
		int a,b,c,d;
		cin >> a >> b >> c >> d;
		m[0] = m[1] = m[2] = 0;
		
		double n = solve(a,b,c,d); // 三次方程式には必ず一つ以上実根が含まれているので、それを見つける。
		m[is(n)]++;
		
		solve(0,a,a*n+b,n*(a*n+b)+c);
		cout << m[2] << " " << m[0] << endl;
	}
}