[ create a new paste ] login | about

Link: http://codepad.org/o7FsJVj3    [ raw code | fork ]

hecomi - C++, pasted on Jul 17:
#pragma once
#pragma warning(disable: 4996)

#include <iostream>
#include <iomanip>
#include <vector>
#include <functional>
#include <algorithm>
#include <cmath>
#include <boost/lambda/lambda.hpp>
#include <boost/lambda/bind.hpp>
#include <boost/lambda/if.hpp>
#include <boost/function.hpp>
#include <boost/timer.hpp>
#include <boost/multi_array.hpp>

namespace calc {
using namespace boost::lambda;
using namespace std;


const double ERROR_VALUE = DBL_MAX;

/*
@brief 等間隔ベクトルをセット
@param[out] y 代入結果を格納する変数
@param[in]  first 下限
@param[in]  end 上限
@param[in]  div セットする総データ数(下限~上限までの分割数)
*/
template <class T, template <class A, class Allocator = std::allocator<A> > class Container>
void makeVector(Container<T> &x, const T first, const T end, const int div)
{
	x.clear();
	if (div <= 0) {
		std::cout << "Error! div <= 0 (@makeVector)" << std::endl;
	}
	for (int i=0; i<=div; i++) {
		x.push_back(first + ((end - first)*static_cast<T>(i))/static_cast<T>(div));
	}
}

/*
@brief ファンクタを用いた代入
@param[out] y 代入結果を格納する変数
@param[in]  x 代入に用いる引数
@param[in]  func y=func(x)
*/
template <class T, template <class A, class Allocator = std::allocator<A> > class Container, typename Functor>
void setVector(Container<T> &y, const Container<T> &x, Functor func)
{
	y.clear();
	std::transform(x.begin(), x.end(), std::back_inserter(y), func);
}


/*
@brief ファンクタを用いた代入(2D)
@param[out] z 代入結果を格納する変数
@param[in]  x 代入に用いる引数
@param[in]  y 代入に用いる引数
@param[in]  func z=func(x,y)
*/
template <class T, template <class A, class Allocator = std::allocator<A> > class Container, typename Functor>
void setVector(boost::multi_array<T, 2> &z, const Container<T> &x, const Container<T> &y, Functor func)
{
	const int zSize = z.shape()[0]*z.shape()[1];
	if (zSize != x.size()*y.size()) {
		std::cout << "Error! @setVector" << std::endl;
		return;
	}

	for (int i=0; i<static_cast<int>(x.size()); i++) {
		for (int j=0; j<static_cast<int>(y.size()); j++) {
			z[i][j] = func(x[i], y[j]);
		}
	}
}

/*
@brief 線形補間クラス
*/
template <class T, template <class A, class Allocator = std::allocator<A> > class Container = std::vector>
class CLerp : public std::unary_function<T, T>
{
private:
	bool Valid;
	Container<T> X, Y;
	const int N;

public:
	/*
	@brief コンストラクタ.内部データを初期化.
	@param[in] x 単調増加x方向ベクトル
	@param[in] y 線形補間対称ベクトル
	*/
	CLerp(Container<T> x, Container<T> y) : Valid(false), N(x.size()-1) {
		// 要素コピー
		std::copy(x.begin(), x.end(), std::back_inserter(X));
		std::copy(y.begin(), y.end(), std::back_inserter(Y));

		// 要素数が2個以上か調べる
		if (X.size() < 2 || Y.size() < 2) {
			std::cout << "Error! The size of X or Y must be greater than 2. (@CLerp::Lerp)" << std::endl;
			return;
		}

		// 要素数が同じか調べる
		if (X.size() != Y.size()) {
			std::cout << "Error! The size of X and Y are different. (@CLerp::Lerp)" << std::endl;
			return;
		}

		// 単調増加か調べる
		for (int i=0; i<N-1; i++) {
			if (X[i] > X[i+1]) {
				std::cout << "Error! X must be monotonically increasing." << std::endl;
				return;
			}
		}

		Valid = true;
	}

	/*
	@brief 線形補間した結果を返す
	@param[in] x 位置
	*/
	T operator()(T x) const {
		// コンストラクタが正常でない場合終了
		if (!Valid) {
			return ERROR_VALUE;
		}

		// 最初の要素より小さかった場合,最初の2つの要素を線形補間
		if (x < X[0]) {
			return Y[0] + (Y[1] - Y[0])/(X[1] - X[0]) * (x - X[0]);
		}
		// 最初の要素より大きかった場合,最後の2つの要素を線形補間
		if (x > X[N]) {
			return Y[N] + (Y[N] - Y[N-1])/(X[N] - X[N-1]) * (x - X[N]);
		}
		// 範囲内の場合
		int cnt = 0, prev, next;
		bool flag = false;
		while (cnt < N) {
			if (x >= X[cnt] && x <= X[cnt+1]) {
				prev = cnt;
				next = cnt+1;
				flag = true;
				break;
			}
			cnt++;
		}
		return Y[prev] + (Y[next] - Y[prev])/(X[next] - X[prev]) * (x - X[prev]);
	}

};

/*
@brief 積分計算を行う
@param[in]  from 積分範囲(下)
@param[in]  to 積分範囲(上)
@param[in]  integrand 被積分関数ファンクタ
@param[in]  N 分割数
*/
template <class T, typename Functor>
T integral(const T from, const T to, Functor integrand, const int div = 100)
{
	// N
	const int N = div*2;
	// ⊿x
	const T dx = (to - from) * 2 / N;
	
	// x, y
	std::vector<T> x, y;
	makeVector(x, from, to, N);
	setVector(y, x, integrand);

	// 結果(シンプソンの積分)
	T res = 0;
	int i = 0;
	while (i < N) {
		res += dx/6 * (y[i] + 4*y[i+1] + y[i+2]);
		i+=2;
	}

	return res;
}

/*
@brief 与えられた2つの変数の範囲内にある根を二分法によって求める.
@param[in]  from 根の範囲(下)
@param[in]  to 根の範囲(上)
@param[in]  integrand 目的の方程式(0になる時の根)
@param[in]  accuracy 二分法の繰り返し回数
@param[in]  err 正確性
*/
template <class T, typename Functor>
T root(T from, T to, Functor equation, const int accuracy = 50, const double err = 1e-6)
{
	// accuracyをチェック
	if (accuracy <= 0) {
		std::cout << "Error! 'accuracy' must be greater than 1 (@root)." << std::endl;
		return ERROR_VALUE;
	}

	// errをチェック
	if (err <= 0) {
		std::cout << "Error! 'err' must be greater than 0 (@root)." << std::endl;
		return ERROR_VALUE;
	}

	// fromとtoが同じでないか
	if (from == to) {
		std::cout << "Error! 'from' is same as 'to'." << std::endl;
		return ERROR_VALUE;
	}

	// 計算の誤差
	const T maxErr = std::max(abs(equation(from)), abs(equation(to)))*err;

	// 間に解が無い場合は,範囲を広げて探索
	int cnt = 0;
	const int maxTry = 10;
	while (equation(from) * equation(to) >= 0) {
		// どちらかが解の場合,その値を返す
		if (abs(equation(from)) <= maxErr) {
			return from;
		}
		if (abs(equation(to)) <= maxErr) {
			return to;
		}
		if (cnt > maxTry) {
			std::cout << "Error! The solution to the equation cannot be found (@root)." << std::endl;
			return ERROR_VALUE;
		}
		T delta = to - from;
		from -= delta;
		to   += delta;
		cnt++;
	}

	T mid;
	for (int i=0; i<accuracy; i++) {
		// 中間値計算
		mid = (from + to)/2;

		// ピッタリの場合は値を返す
		if (abs(equation(mid)) <= maxErr) {
			break;
		}

		// 中間値を代入
		if (equation(from) * equation(mid) < 0) {
			to = mid;
		} else {
			from = mid;
		}
	}

	return mid;
}

/*
@brief boost::timerを利用して簡単に時間計測を行う
*/
class CTimer : public boost::timer
{
public:
	void operator()(const char* title)
	{
		cout << setw(20) << left << title;
		cout << ": " << elapsed() << " sec" << endl;
		restart();
	}
};


}


Create a new paste based on this one


Comments: