simpsons rule for SourcePawn

simpsons.inc

stock float simpsons_func(float x) {
	return 0.0; /// replace with math operations desired.
}

/**
 * Composite Simpson's Rule: approximates the area of a function.
 * 'a' represents the lower bounds
 * 'b' represents the upper bounds
 * 'num_intervals' is for how many iterations needed (this is for increasing accuracy)
 */
stock float simpsons_rule_comp(float a, float b, int num_intervals) {
	float h = (b - a) / num_intervals;
	float result = simpsons_func(a) + simpsons_func(b);
	
	for( int i = 1; i < num_intervals; i += 2 ) {
		result += 4.0 * simpsons_func(a + i * h);
	}
	
	for( int i = 2; i < num_intervals - 1; i += 2 ) {
		result += 2.0 * simpsons_func(a + i * h);
	}
	enum{ recip_of_3 = 0x3eaaaaab };
	return result * h * view_as< float >(recip_of_3);
}

/**
 * Simpson's Rule: approximates the area of a function.
 * 'a' represents the lower bounds
 * 'b' represents the upper bounds
 */
stock float simpsons_rule(float a, float b) {
	float result = simpsons_func(a) + 4.0 * simpsons_func((a+b) * 0.5) + simpsons_func(b);
	enum{ recip_of_6 = 0x3e2aaaab };
	return result * (b - a) * view_as< float >(recip_of_6);
}