1
25
2010
0

数值积分算法Simpson

http://en.wikipedia.org/wiki/Adaptive_Simpson%27s_method

 

double simpson_rec( double a, double b, double S, double fa, double fb, double fc,
		int bottom, double epsilon){
	double c = (a+b)/2, h = b-a;
	double d = (a+c)/2, e = (c+b)/2;
	double fd = f(d), fe = f(e);
	double SL = (h/12)*(fa + 4*fd+fc);
	double SR = (h/12)*(fc + 4*fe+fb);
	double S2 = SL + SR;
	if( bottom <= 0 || fabs(S2 - S ) <= 15*epsilon){
		return S2 + (S2-S)/15;
	}
	return 	simpson_rec(a,c,SL,fa,fc,fd, bottom-1, epsilon/2) + 
		simpson_rec(c,b,SR,fc,fb,fe, bottom-1, epsilon/2);
}

double simpson( double a, double b, double epsilon, int rec){
	double c = (a+b)/2, h = b-a;
	double fa = f(a);
	double fb = f(b);
	double fc = f(c);
	double S = (h/6) * (fa+4*fc+fb);
	return simpson_rec(a,b,S,fa,fb,fc,rec,epsilon);
}

 

 

Category: 未分类 | Tags:
1
23
2010
0

大整数模板

一个简单的大整数模板,使用了vector来实现,效率应该会灰常低,不过估计打起来会比较快吧,至少在Topcoder上是有用的


#define ll long long
#define SZ size()

ll BASE = 1000000000;
struct BigInteger{
vector<ll> num;
void operator+=(const BigInteger& other){
vector<ll>& a = num;
const vector<ll>& b = other.num; // const CE!
a.resize( max(a.size(), b.size()) + 1, 0);
Rep(i,b.size()) if( (a[i] += b[i]) >= BASE ) a[i] -= BASE, a[i+1]++;
for( int i = b.size(); a[i] >= BASE ; ++i) a[i] -= BASE, a[i+1]++;
if( a.SZ >= 2 && a.back() == 0 ) a.resize( a.size() - 1);
}

void operator-=(const BigInteger& other){
                // 默认this >=
vector<ll>& a = num;
const vector<ll>& b = other.num; // const CE!
a.resize( max(a.size(), b.size()) + 1, 0);
Rep(i,b.size()) if( (a[i] -= b[i]) < 0 ) a[i] += BASE, a[i+1]--;
for( int i = b.size(); a[i] < 0; ++i) a[i] += BASE, a[i+1]--;
if( a.SZ >= 2 && !a.back()) a.resize( a.size() - 1);
}

void print(){
printf("%lld", num.back());
                                                // remember to change length
for(int i = num.size() - 2; i>= 0; --i)printf("%09lld", num[i]);        
        }         void scan(){
char str[1000];
scanf("%s", str);
int n = strlen(str);
num.resize((n+3)/4, 0); // ceil
for(int i = 0; i < n; ++i) num[(n-1-i)/4] = num[(n-1-i)/4]*10 + str[i] - '0';
}
void init(int a){
num.clear();
if( a < 0 ) a = 0;
for(; a; a/=BASE) num.push_back( a % BASE );
if( num.size() == 0 ) num.push_back(0);
}                 // return the remainder, b is int int operator/=(int b){
int c = 0;
for(int i = num.size() - 1; i>=0; --i){
c = c * BASE + num[i];
num[i] = c / b;
c %= b;
}
while( num.size() >= 2 && num.back() == 0 ) num.resize( num.size() - 1 );
return c;
}
BigInteger(int a = 0){
init(a);
}
};
BigInteger operator*(const BigInteger& _a, const BigInteger& _b){
BigInteger _c;
const vector<ll>& a = _a.num; // const CE!
const vector<ll>& b = _b.num; // const CE!
vector<ll>& c = _c.num;
c = vector<ll>(a.size() + b.size(), 0);
Rep(i,a.size())Rep(j,b.size()){
if( (c[i+j] += a[i] * b[j]) >= BASE )c[i+j+1] += c[i+j] / BASE, c[i+j] %= BASE;
}
while( c.SZ >= 2 && !c.back() )c.resize( c.size() - 1); Rep(i,c.SZ)assert(0<=c[i] && c[i] < BASE ); // 对调试很有帮助的断言
return _c;
}
bool operator==(const BigInteger& a, const BigInteger& b){return a.num == b.num;}
bool operator < (const BigInteger& a, const BigInteger& b ){
if( a.num.size() == b.num.size() )                return lexicographical_compare(a.num.rbegin(), a.num.rend(), b.num.rbegin(), b.num.rend()) ;
return a.num.size() < b.num.size();
}

 

 

Category: 未分类 | Tags:
1
21
2010
0

个人代码集合地

 

循环数组的处理

 

int mk[1024];
ll arr[2048];

	int period, pre;
	arr[0] = A;
	for(int i = 1; ; ++i){
		arr[i] = a * arr[i-1] * arr[i-1] % M + b * arr[i-1] + c;
		arr[i] %= M;
		if( mk[ arr[i] ] ){
			pre = mk[arr[i]];
			period = i - pre;
			break;
		}
		mk[ arr[i] ] = i;
	}

 

归并排序,顺便求逆序数

 

void merge_sort( int* v, int n){
	if( n < 20 ) Rep(__, n) Rep(i,n-1)if( v[i] > v[i+1] ) swap( v[i], v[i+1]), cnt++;
	if( n < 20 ) return;
	merge_sort(v, n/2);
	merge_sort(v+n/2, (n+1)/2);
	int l = 0, r = n/2;
	for(int i = 0; i<n ; ++i){
		if(r == n || l != n/2 && v[l] <= v[r] ) tmp[i] = v[l++] ;
		else tmp[i] = v[r], cnt += r - i, r++;
	}
//	Rep(i,n) cout << tmp[i] << " "; cout << endl;
	copy( tmp, tmp + n, v);
}

 

两圆相交求并

 

double inter_circle( double x1, double y1, double r1, double x2, double y2, double r2){
	double d = sqrt(sq(x1-x2)+sq(x2-y2));
	if( d >= r1 + r2 ) return pi*(r1*r1+r2*r2);
	if( d <= abs(r1-r2) || abs(d) < 1e-8) return pi*max( r1*r1, r2*r2);
	double t1 = acos(( d*d + r1*r1 - r2*r2) / 2 / d / r1);
	double t2 = acos(( d*d + r2*r2 - r1*r1) / 2 / d / r2);
	double area = r1*r1*pi + r2*r2*pi;
	area -= t1 * r1 * r1;
	area -= t2 * r2 * r2;
	double s = (r1 + r2 + d)/2;
	area += 2 * sqrt( s * (s-r1) * (s-r2) * (s-d));
	return area;
}

 

两球相交求并

 

double inter_ball( double r1,double r2, double d){
	if( d >= r1 + r2 ) return 0;
	if( d <= abs(r1-r2) || abs(d) < 1e-8) return pi*min( r1*r1, r2*r2);
	double p = (r1 + r2 + d)/2;
	double tri = sqrt( p * (p-r1) * (p-r2) * (p-d));
	double r = tri / d * 2; 
	double h1 = sqrt( sq(r1-r) );
	double h2 = sqrt( sq(r2-r) );
	double s = pi*r*r;
	double area = h1 * s / 3 + h2 * s / 3;
	double x1 = (d*d+r1*r1-r2*r2)/2/d - d;
	double x2 = (d*d+r2*r2-r1*r1)/2/d - d;
	area = pi * ( r1*r1*(r1-x1) - ( r1*r1*r1 - x1*x1*x1)/3);
	area += pi * ( r2*r2*(r2-x2) - ( r2*r2*r2 - x2*x2*x2)/3);
	return area;

 

 

 

Category: 未分类 | Tags:

Host by is-Programmer.com | Power by Chito 1.3.3 beta | Theme: Aeros 2.0 by TheBuckmaker.com