Динамические мат. функции в C++

Здравствуйте, Хабраюзеры.
Недавно я прочитал здесь статью об анонимных функциях в С++, и тут же у меня в голове возникла мысль: нужно срочно написать класс для работы с функциями, которые нам знакомы из математики. А именно, принимающими вещественный аргумент и возвращающими вещественное же значение. Нужно дать возможность обращаться с такими объектами максимально просто, не задумываясь о реализации.
И вот, как я это реализовал.

Проблема.

Все эти лямбда-выражения ведут себя порой довольно странно, по крайней мере для меня. Связано это, вероятно, с тем, что я не могу до конца понять, как устроен механизм создания этих выражений. Создавать большие функции я буду на основе уже имеющихся примитивными действиями, т.е. что-то вроде f(x) = a(x) + b(x). А это значит, что все лямбды, созданные даже лишь как промежуточные звенья в построении функций, должны сохраняться для обращения к ним. Сейчас попробую объяснить попонятнее.

Скажем, мы хотим, чтобы некая процедура принимала пару функций A и B и возвращала новое выражение, например, 5A + B. Наша процедура создаст лямбду 5A, затем создаст 5A + B, используя 5А и В. Полученное процедура вернет и завершится, в этот момент лямбда 5А пропадет из области видимости, и возвращенное выражение попросту не будет работать.

Решение.

Я решил создать одну глобальную коллекцию всех лямбд, все строящиеся выражения хранятся в ней, а в объектах будут присутствовать лишь указатели на элементы коллекции. Поясню кодом:

#include <xstddef>
#include <functional>
#include <list>

typedef std::tr1::function<double(double)> realfunc; // y = f(x) as in maths

class func_t
{
protected:
	static std::list<realfunc> all_functions; //коллекция всех функций

	realfunc *f; //указатель на элемент коллекции

public:
	func_t();
	func_t(const double);
	func_t(const realfunc&);
	func_t(const func_t&);

	~func_t() {};

	friend func_t operator+ (const func_t&, const func_t&);
	friend func_t operator- (const func_t&, const func_t&);
	friend func_t operator* (const func_t&, const func_t&);
	friend func_t operator/ (const func_t&, const func_t&);
	friend func_t operator^ (const func_t&, const func_t&);

	func_t operator() (const func_t&);

	double operator() (const double);

};

Для начала конструкторы.
Конструктор по умолчанию (без параметров) будет создавать функцию, возвращающую аргумент. Это будет точка отправления. f(x) = x.
Остальные понятны: второй создает функцию-константу — f(x) = c, третий превращает лямбду нужного типа в объект моего класса, последний — просто конструктор копирования.

Реализация конструкторов, по ней сразу будет видно, как устроены все методы класса:

func_t::func_t()
{
	f = &(*all_functions.begin());
}

func_t::func_t(const double c)
{
	func_t::all_functions.push_back(
		[=](double x)->double {return c;}
	);
	this->f = &all_functions.back();
}

func_t::func_t(const realfunc &realf)
{
	func_t::all_functions.push_back(realf);
	this->f = &all_functions.back();
}

func_t::func_t(const func_t &source)
{
	this->f = source.f;
}

Как видите, я создаю лямбду, толкаю ее в конец коллекции и возвращаю объект с указателем на нее.
Хочу сразу пояснить первый же конструктор. Как я уже говорил, создание «аргумента», т.е. функции f(x)=x является началом почти любой работы с моим классом, поэтому я решил особо выделить эту функцию и положил в первую же ячейку коллекции это выражение. И тогда при вызове конструктора по умолчанию объект всегда получает указатель на первый элемент коллекции.

Ах да, чуть не забыл, все конструкторы могут быть использованы для неявного преобразования, что и создает основное удобство использования.

Далее, операторы. С ними всё просто. Покажу три из них: первый реализует сложение, второй композицию функций, третий оценку.

func_t operator+ (const func_t &arg, const func_t &arg2)
{
	realfunc realf = [&](double x)->double {
		return (*arg.f)(x) + (*arg2.f)(x);
	};
	return func_t(realf);
}

func_t func_t::operator() (const func_t &arg)
{
	realfunc realf = [&](double x)->double {
		return (*f)((*arg.f)(x));
	};
	return func_t(realf);
}

double func_t::operator() (const double x)
{
	return (*f)(x);
}

Всё просто, правда? В первых двух опять создается нужная лямбда, а потом посылается в конструктор объекта. В третьем методе вообще халява =)
Поясню, что я везде использую передачу окружения по ссылкам (это [&] перед лямбдой) для доступа к аргументам метода.

В принципе, это всё. Теперь пара технических деталей. Я не очень силень в инициализации статических полей, поэтому пришлось вспоминать ужасно громозкие и некрасивые способы, которые я когда-то где-то подсмотрел.

protected:
	static class func_t_static_init_class
	{ public: func_t_static_init_class(); };
	static func_t_static_init_class func_t_static_init_obj;

...

//Static:
std::list<realfunc> func_t::all_functions = std::list<realfunc>();
func_t::func_t_static_init_class func_t::func_t_static_init_obj = func_t::func_t_static_init_class();

func_t::func_t_static_init_class::func_t_static_init_class()
{
	func_t::all_functions.push_back(
		[](double x)->double {return x;}
	);
}

Ну вы поняли, создаю список и пихаю в него первый элемент, о котором уже упомянал.
Прошу прощения за этот ужас, только учусь программировать.

Бонусы.

В принципе, вот и всё. Осталась пара вещей, которые я делал уже чисто для интереса (хотя, собственно, как и всё).

Во-первых, перегрузим парочку функций из cmath.

friend func_t sin(const func_t&);
friend func_t cos(const func_t&);
friend func_t tan(const func_t&);
friend func_t abs(const func_t&);

...

func_t sin(const func_t& arg)
{
	realfunc realf = [&](double x)->double {
		return sin((*arg.f)(x));
	};
	return func_t(realf);
}

Во-вторых, куда ж без производных и первообразных =)

static double delta_x;

func_t operator~ ();
func_t operator| (const double);

...

double func_t::delta_x = 0.01;

func_t func_t::operator~ ()
{
	realfunc realf = [&](double x)->double {
		return ((*f)(x + delta_x / 2) - (*f)(x - delta_x / 2)) / delta_x;
	};
	return func_t(realf);
}

func_t func_t::operator| (double first_lim)
{

	realfunc realf = [=](double x)->double {
		double l_first_lim = first_lim; //will move with this copy of first_lim
		double area = 0;
		bool reverse = x < first_lim; //first_lim > second_lim?
		if (reverse) {
			l_first_lim = x;
			x = first_lim;
		}
		double l_delta_x = delta_x; //step
		while (l_first_lim < x) { //move along the whole span
			if ((l_first_lim += l_delta_x) > x) //stepped too far?
				l_delta_x += x - l_first_lim; //the last l_delta_x may be shorter
			/* integral summ, the point is chosen between
			   the point for f(x) is chosen between l_first_lim and l_first_lim + l_delta_x */
			area += l_delta_x * (*f)(l_first_lim + l_delta_x / 2);
		}
		return area * (reverse?-1:1);
	};
	return func_t(realf);

}

Не буду вдаваться в объяснение данного кода, он скучный и является абсолютно наивной реализацией производной в точке и интеграла с переменным верхним пределом (в обоих случаях вместо предела используется фиксированная дельта).

Заключение.

Ну вот, собственно, и всё. Надеюсь, что было интересно. Не знаю, имеет ли что-то такое хоть какой-то смысл, но я попрактиковался в классах, всяких & и *, а для меня это главное =) Спасибо за внимание.

Упс! Еще кое-что.

Ну да, как это использовать.
Например вот так:

func_t f1 = cos(5 * func_t() + 8);

это создаст, как видно, функцию
f1(x) = cos(5x + 8)

или вот так:

funt_t x = func_t();
func_t f = x + f1(x / ~f1);

это f(x) = x + f1(x / f1`(x)) кэп мимо проходил

или в конце концов так:

realfunc g = [](double->double) {
    ...
    //что вашей душе угодно, хоть сокеты создавайте
    ...
}
func_t f3 = g;

Заключение номер 2.

Ну теперь точно всё. Еще раз спасибо за внимание!
Если что, полный и недоделанный код тут.