Статья: "Самый быстрый калькулятор.Транслятор математических выражений в вашей программе."

Статья

Самый быстрый калькулятор.Транслятор математических выражений в вашей программе.

О чем эта статья?

В данной статье описывается механизм работы авторского класса C++ (далее просто класса) для трансляции математических выражений, составленных из математических функций, представленных в строковом виде (например, sin(X)/X+4.78 ). Здесь описывается динамическое создание макрофункции из машинных команд (далее дампа) по строковому выражению функции (далее просто функции).
Описываемый здесь класс способен:
a) Создавать дамп по заданной функции. (Например, по функции sin(X) будет сгенерирован, скажем, такой дамп: D905BC005100D9FEC3).
b) Исполнять дамп по заданной переменной. (Если передать управление дампу D905BC005100D9FEC3, то сопроцессор и процессор выполнят следущие операции:
D905BC005100 загрузит вещественное число, находящееся по адресу BC005100 (команда FLD);
D9FE вычислит синус загруженного числа (команда FSIN);
C3 возвратит управление туда, откуда был вызван дамп (команда RET).
Записав нужную переменную по адресу BC005100 и передав управление дампу, можно заставить сопроцессор выполнять сгенерированную функцию - код с разными с разными вещественными аргументами. При возврате в точку, откуда вызывался дамп, используя команду FSTP, мы можем посмотреть, что же подсчитал сопроцессор.
Сгенерированный код дамп содержит МИНИМАЛЬНЫЙ набор машинных команд, необходимых для вычислений. Дамп не содержит операций манипуляций с тестовой строкой функцией.
Это все и обьясняет быстродействие вычисления функций, представленных в строковом виде.

c) Вычислять значение функции в режиме интерпретации для заданной переменной.

Кратко описывается сопроцессор и работа сопроцессора при вычислении некоторых функций как моих, так и функции библиотеки MSVCRT (MS Visual C++ Run time library).
Изложение статьи, скорее, ориентировано на не очень искушенных программистов, хотя идеи класса будут, я думаю, любопытны остальным.

К статье прилагается описанный здесь класс; Демо проект MS VC++ 6.0, реализующий класс; Программа рекурсивного поиска файлов (проект MS VC++ 6.0 ); Программа построения функции любого вида, с использованием обьекта листа EXCEL, написанная на MS VB 6.0.

Зачем это нужно?

Описание класса поможет вам разобраться, что происходит внутри программы и сопроцессора при использовании математических функций.
Использование класса позволит снизить время вычисления функций, представленных в строковом виде. А при многократном вычислении заданной функции по различным аргументам позволит ЗНАЧИТЕЛЬНО снизить время вычисления. Сравниние скоростей вычисления описано в разделе А Скорость?.
Области возможного применения класса:
a) Моделирование процессов. При моделировании обычно требуется производить большие объемы вычислений, выводить графическое отображение зависимостей и часто модифицировать алгоритмы вычислений. Пример: написать программу построения функции любого вида, используя описываемый здесь класс, очень просто. А скорость подсчета некоторых функций, задаваемых пользователем программы в строковом выражении, будет даже выше, чем при построении тех же функций по зашитому, статичному алгоритму, сгенерированному по инструкциям высокоуровневого языка.
b) Конфигурирование и интеллектуальный учет.
Тонкое конфигурирование и настройка программ, критичных к скорости выполенния (например ядра игры) и бизнес-приложений (где важно удобство настроек) пользователем.
c) Внутренние - скриптовые языки программирования, встраиваемые внутри программ.
Хороший пример - встроенный язык программирования, встроенный WINAPM - для задания пользователем функций построения визуализированного звука. Скорость вычисления задаваемых пользователем функций можно значильно повысить.

   Технические возможности : класс позволяет создавать дампы по строковым выражениям функций, которые содержат следующие операторы и функции:

Функции:

sin(arg) - синус аргумента
   cos(arg) - косинус аргумента
   tg(arg) - тангенс аргумента
   ctg(arg) - котангенс аргумента
   asin(val) - обратный синус аргумента (в радианах)
   acos(val) - обратный косинус аргумента (в радианах)
   atg(val) - обратный тангенс аргумента (в радианах)
   изменение знака аргумента
   exp(val) - число е в степени val
   ln(val) - натуратьный логарифм аргумента
   round(val) округление вещественного числа к целому

   Операторы:

   * умножение
   / деление
   + сложение
   - вычитание
   ^ возведение в степень

ОГРАНИЧЕНИЯ КЛАССА:

ВНИМАНИЕ! Приоритетность всех операторов одинаковая. Приоритетность операций задается пользователем посредстром скобок. Например, функцию sin(X)+X*X нужно записывать так: sin(X)+(X*X), иначе класс примет ее за функцию (sin(X)+X)*X. (Это ограничение я ввожу, чтобы избежать возможность генерирования кода, который может повлечь за собой переполнение стека сопроцессора. Для рекурсивного метода разбора строкового выражения функции мне не удалось решить эту нехорошую возможность менее жестским ограничением).
Класс по умолчанию резервирет немного памяти под дамп. Если придется сильно нагружать класс, то придется изменять размер дампа.
Класс и реализованные алгоритмы (экспонента, степенные функции, обратные тригонометрические функции) еще не проходили серьезного тестирования.

Кратко о сопроцессоре.

Рабочая среда сопроцессора (далее FPU) состоит из:
   Аппаратуры
   Окружения
   Системы команд
   Интерфейса ввода-вывода.

Нас будет интересовать окружение и система команд FPU, из окружения рассмотрим только регистр состояния и стек FPU. Регистер состояния важен для нас, т.к в нем хранится описание ошибки.

FPU содержит восемь 80-битовых регистров ST7..ST0, предназначенных для хранения промежуточных результатов вычислений.

Регистр состояния FPU имеет формат:

IE	Недействительная операция. Возникает при попытке выполнения запрещенных математических операций или обращение к свободному регистру.
DE	Денормализованный результат. Возникает в случае, когда результат не удается представиль в нормализованной форме..
ZE	Деление на нуль..
OE	Если результат выполнения какой-то операции слишком велик и не может быть представлен в формате приемника результата..
UE	Антипереполнение. Возникает тогда, когда результат слишком мал для его представления в формате приемника результата операции, но отличен от нуля..
PE	Неточный результат. Возникает, скорее всего, в случае, когда в результате получается вещественное число, период которого больше разрядности мантисы - число округляется..
XX	Не используется.
ES	Флаг суммарной ошибки, который устанавливается в 1 при возникновении незамаскированного особого случая..
C0	Код условия. Он определяется по результату выполнения команд сравнения FCOM и команды нахождения остатка FPREM..
C1	Код условия.
C2	Код условия.
ST	Содержит номер X регистра STX, являющегося вершиной стека..
C3	Код условия.
B	Бит занятости. Он устанавлен, когда FPU выполняет команду или когда происходит прерывание от FPU..

Все команды сопроцессора можно разделить на несколько групп:
   команды пересылки данных
   арифметические команды
   команды сравнений чисел
   трансцендентные команды
   управляющие команды

Приведу краткое описание тех лишь команд процессора, которые используются в классе:

F2XM1      - вычисление степени x по основанию 2 минус 1
FABS      - модуль x
FADD      - сложение
FCHS      - изменение знака x
FCOS      - косинус x
FDIV      - деление
FMUL      - умножение
FPATAN   - арктангенс частного
FPREM      - вычисление остатка от деления
FPTAN      - вычисление тангенса
FSIN      - вычисление синуса
FSINCOS   - вычисление синуса и косинуса
FSQRT      - вычисление квадратного корня из
FSUB      - вычитание
FYL2X      - вычисление y*log2(x)
FYL2XP1   - вычисление y*log2(x+1)
FFLD      - загрузка переменной
FLDPI      - загрузка числа пи
FLD1      - загрузка единицы
FLDl2T   - загрузка числа ln(10)/ln(2)
FLDL2E   - загрузка числа 1/ln(2)
FLDLG2   - загрузка числа log(2)/ln(10)
FLDLN2   - загрузка числа ln(2)
FLDZ      - загрузка нуля
FXCH      - обмен
FSCALE   - вычисление x*2^y
FSTP      - выгрузка результата и освобождение
FRNDINT   - округление до целого числа
FFREE      - освобождение регистра стека сопроцессора
FINCSTP   - увеличение указателя стека на единицу

Работа с FPU на примере.

У сопроцессора нет команды вычисления обратного синуса. Чтобы его считать, приходится реализовывать процедуры, делающие это. В качестве примера приведу процедуру подсчета обратного синуса по формуле acrsin(x) = arctg(x/sqrt(1-X*X)), реализованную мной.
Запустите демо-программу. Выберите Файл + Тесты+1 Вычисление обратного синуса.
В статус-окне появится Тест 1: arcsin(6.543210E-001) = 7.132844E-001.
Вот процедура, которую вы запускаете:

Код:

void	ArcSinTst() {
 float x,Result;

x = 0.654321;
__asm {
fld	x	// загружаем аргумент
fld1		// загружаем 1
fld 	st(1)	// дублируем аргумент
;состояние стека FPU на этот момент  кол. 1 рисунка

fmul	st(0),st(0)	// умножаем
;состояние стека FPU на этот момент  кол. 2 рисунка

fsub		// вычитаем
;состояние стека FPU на этот момент  кол. 3 рисунка

fsqrt		// извлекаем корень
;состояние стека на этот момент  кол. 4 рисунка

fpatan		// вычисляем обратный тангенс
;состояние стека на этот момент  кол. 5 рисунка

fstp	Result	// выгружаем значение из верхушки стека
}


sprintf(szMathStr,"Тест 1: arcsin(%E) = %E",ttt,ttr);
SetDlgItemText(hWnd,IDC_RESULT,szMathStr);
}

Как это работает?

Класс может работать в двух режимах: ACTION_CALCULATE и ACTION_BUILD. Описание этих режимов приведенно ниже.

Список методов и переменных, который надо знать, чтобы формально использовать класс:

Аргумент функции, для которого требуется вычислить функцию. Эта переменная используется только при работе класса в режиме ACTION_CALCULATE.
float         CalcX;

Функция для сброса управляюцих указателей, для подготовки работы калькулятора в режиме ACTION_BUILD.
void   Ce()

Добавление вещественного числа в стек переменных дампа
void   Dump_SP_PushVar(
float flVar      // вещественное число, которое бедет добавлено
);

Запись слова в стек переменных дампа.
void   Dump_P_PushWord(WORD wCode);

Запуск дампа - вычисление дампа функции по переменной, находящейся в верхушке стека переменных дампа.
float         Dump_Execute();

Функция для запуска начала действия с функцией. Действие характеризуется переменной ActType.
float   Start (
char *ptFunction,   // указатель на функцию
int ActType      // характер действия
);

char *ptFunction
указатель на функцию. Имейте в виду, что функция не копируется во внутренний буфер. Непродумманое изменение строкового выражения функции в момент работы класса в режиме ACTION_BUILD или ACTION_CALCULATE может повлечь за собой неправильный результат работы или вызвать исключения класса.

int ActType;
   Может принимать определенные значения:

ACTION_CALCULATE
формирование дампа по указателю на строковое выражение функции ptFunction. Следите за указателем lpDump! По умолчанию под дамп резервируется 192 байт, что вполне достаточно для реализации довольно сложных функций (на каждую функцию затрачивается 2 байта (есть исключения) в дампе и 2 байта на каждую постоянную в области дампа со смещением DEF_STACK_OFFSET, равным по-умолчанию 176 байтам). В случае переполнения дампа в область переменных может попасть код. В результате получим, что функция будет выполняться не с теми постоянными, какими мы их закладывали. Это в лучшем случае. При других раскладах класс может спровоцировать GPF.

ACTION_BUILD
Вычисление значения функции без создания дампа методом рекурсивного разбора функции.

Исключения, возбуждаемые классом:

Для режима ACTION_CALCULATE
calc::SyntaxError      - синтаксичеккая ошибка

Для режима ACTION_BUILD
сalc::divzero      - деление на ноль
calc::SyntaxError      - синтаксичеккая ошибка
calc::MathError      - математическия ошибка

Возвращаемое значение:
В случае ACTION_CALCULATE возвращаемое значение представляет собой значение функции ptFunction для аргументта CalcX.

В случае ACTION_BUILD возвращаемое значение - NULL .

Пример формального использования класса

В режиме ACTION_CALCULATE

Код:

calc		calculator;
float		Result;
char		szMathStr[100];

	try {
calculator.CalcX = 1.78e+1;	// Устанавливаем аргумент X = 17,8
Result = calculator.Start(sin(X)/X,ACTION_CALCULATE);//  Вычисляем sin(17,8)/17,8
}	
// В случае возникновения исключений - обрабатываем их
catch (calc::divzero) {
MessageBox(NULL,"Деление на ноль.",NULL,MB_ICONERROR);
}
catch(calc::SyntaxError e) {
wsprintf(szMathStr,"Ошибка синтаксиса:%s.",e.p);
MessageBox(NULL,szMathStr,NULL,MB_ICONERROR);
}
catch(calc::MathError e) {
wsprintf(szMathStr,"Математическая ошибка:%s.",e.p);
MessageBox(NULL,szMathStr,NULL,MB_ICONERROR);
};

В режиме ACTION_BUILD

Код:

calc		calculator;
float		Result;
char		szMathStr[100];

// Создаем дамп функции acos(sin(X) * X/2) :
try {
		calculator.Start(acos(sin(X) * X/2),ACTION_BUILD);
	}
	catch(calc::SyntaxError e) {
		wsprintf(szMathStr,"Ошибка синтаксиса:%s.",e.p);
		MessageBox(NULL,szMathStr,NULL,MB_ICONERROR);
	}

	// вычисляем полученный дамп для значения аргумента 1e+0
calculator.Ce();				// Сбрасываем калькулятор
	calculator.Dump_SP_PushVar(1e+0);	// Сохраняем переменную
	Result = calculator.Dump_Execute();	// Получаем результат
	
	// вычисляем полученный дамп для значения аргумента 1.1e+0
calculator.Ce();				// Сбрасываем калькулятор
	calculator.Dump_SP_PushVar(1e+0);	// Сохраняем переменную
	Result = calculator.Dump_Execute();	// Получаем результат

Как я уже сказал, создание дампа функции или вычисление ее без создания дампа реализованно по народному алгоритму - рекурсии.
Если вы не посвящены в механику работы этого метода, я скажу, что основной идеей является многократрое использование переменных процедуры посредством вложенных вызовов рекурсивных процедур. Этот метод позволяет разумно выделять память, необходимую для обработки. (Как известно, память под локальные переменные процедуры выделяется в стеке, стандартный размер которого не все программы способны использовать в полной мере). Все алгоритмы, реализованные на основе рекурсии, можно реализовать и без нее. Для этого, обычно, требуется создание кучи, размер которой еще и мучительно выбирать (хочется, чтобы он был как можно меньше, и программа занимала меньше памяти, но и хочется, чтобы этот буфер не переполнялся), плюс процессору придется тратить время на администрирование кучи (против лишних инструкций CALL и инструкций освобождения стека в методе рекурсии). Что лучше - неизвестно.
Я знаю три метода, по которым можно реализовать синтаксический разбор строки. Но рекурсивный метод, на мой взгляд, наглядней.

Разобраться с работой метода рекурсии поможет простая программа, находящаяся в катологе tree_recurs, рекурсивного поика файлов в заданной дирректории и добавление их имен в элемент управления TreeView. В ней рекурсивной процедурой является небольшая процедура DoRecursion. Если вы хотите разобраться, как работает метед рекурсии для разбора строкового выражения функции в классе, я рекормендую скачала изучить работу этой программы.

Код:

// Процедура рекурсивного поиска файлов и добавления их в
// элемент управления TreeView
void mtree::DoRecursion(HTREEITEM stWorkItem) {
 HANDLE hFind;					
 HTREEITEM WorkItem;			
	if (FirstRecursionCall) {
		WorkItem = InsTreeItem((HTREEITEM)TVI_ROOT, WorkPath, 
			(HTREEITEM)TVI_FIRST, 0,0);
		FirstRecursionCall = false;
	} else {
		WorkItem = stWorkItem;
		strcat(WorkPath,fd.cFileName);
		strcat(WorkPath,"\");
	}
	strcat(WorkPath,Mask);
	hFind = FindFirstFile(WorkPath, &fd);
	pt = WorkPath + strlen(WorkPath) - strlen(Mask);
	*pt=0;
	if (hFind == INVALID_HANDLE_VALUE) goto nops;
	goto insert;
	while (FindNextFile(hFind, &fd)) {
insert: if ((fd.dwFileAttributes == FILE_ATTRIBUTE_DIRECTORY) &&
			((strcmp(fd.cFileName ,".") == NULL) ||
			(strcmp(fd.cFileName,"..") == NULL)
			)) continue;

		WorkItemTem = InsTreeItem(WorkItem,fd.cFileName,(HTREEITEM)TVI_SORT, 0,0);
	
		if (fd.dwFileAttributes == FILE_ATTRIBUTE_DIRECTORY)
			DoRecursion(WorkItemTem);
	}
nops: if (!FirstRecursionCall) {
		pt = WorkPath + strlen(WorkPath) - 2;
		while ( ((*pt) != '\') && (pt != WorkPath)) pt--;
		pt++;
		*pt = NULL;
	}
}

Добавьте свое!

Для добавления новой макрофункции (например, гиперболического синуса) вам понадобится добавить if(strcmp(TokenValue,"sh")==0) { } в блок case FUNC: { } метода calc::Dig().
При составлении макрофункции учитывайте следуещее:

Аргумент для посчета вашего макроса уже загружен в ST(0).
По ходу выполнения вашего макроса вам предоставлеется все восемь регистров ST(7)..ST(0).
После завершения работы вашего макроса результат вычислений должен быть представлен в ST(0), а все другие стековые регистры FPU должны быть свободны.
Понадобится представление вашего макроса в машинном коде.

Есть возможность использовать сложные макрофункции, написанные даже на языке высокого уровня, например, через таблицу экспорта.

Для этого следует описать макрос:

Код:

MACRO_EXPORT void	MacroTest_ch() {
 Float myval,result;
	MessageBox(NULL,"called from user macro","MessageBox()",0);

	__asm{FSTP myval}				// выгружаем аргумент из стека
result = (exp(myval)-exp(-myval))/2;		// вычисляем
	__asm{FLD result}				// загружаем результат в стек
	return;
}

В блоке if(strcmp(TokenValue,"MacroTest_ch")==0) { } необходимо подставлять код процедуры запуска макрофункции в дамп - CALL fpTestMacro,

А при инициализации класса - определять адрес макрофункции так:
FARPROC fpTestMacro = GetProcAddress(GetModuleHandle(NULL),"MacroTest");
макрофункции должны быть помечены для эспорта дирректирой MACRO_EXPORT.

Что внутри?

Ниже приведен код для вычисления функции возведения в степень - exp(). Здесь приведена только часть процедуры по выполнению, код который реально выполняется при использовании функции exp в C++ программе, скомпилированной, например, в среде MS Visual C++ 6.0, в штатном режиме (т.е. чаще всего).

Для простоты восприятия я оформил этот код в виде процедур, а не сплошным текстом, так, как это примерно выглядит на том языке, на котором это написанно. Назначение каждой процедуры написанно у нее после заголовка. Необходимость вводить такую функцию возникает из-за скудного набора команд FPU работы со степенью. Их всего две -

F2XM1 (возведение 2 в вещественную степень x, где x принадлежит диапазону 1<x<+1) и
FSCALE (возведение 2 в степень x, где x представленно в виде целого числа).

Чтобы понять, как может вычисляться экспонента, приведу математическое преобразование экспоненты:

Round(x/ln(2)) - целое число.

(x/ln(2) - Round(x/ln(2))) - вещественное число, принадлежащее (-1;+1).

Теперь, я думаю, стало понятно, как подсчитать экспоненту командами F2XM1 и FSCALE. Реализованная мной процедура подсчета экспоненты приведена ниже.

А Скорость?

Чтобы не быть голословным, приведу скорость работы класса по сравнению с достойным алгоритмом. В качестве этого достойного алгоритма возьму алгоритм вычислений, основанный на рекурсивной интерпретации функции, который характеризуется неплохим быстродействием, но реализуется проще.

В качестве недостойного алгоритма приведу алгоритм использования обьекта листа EXCEL для подсчета функции любого вида. Этот алгоритм реализован в программе GraphIX.
Где вы сами (на глаз) можете оценить скорость выполнения (программа GraphIX->файл->функции->стандартный вид).

На рисунке представленна зависимость времени выполения от возрастающей вычислительной нагрузки. Такую нагрузку дает, например, такая группа функций:

sin(X)+sin(X)
sin(X)+sin(X)+sin(X)

sin(X)+sin(X)+sin(X)+...+sin(X)

Если посмотреть на дампы, сгенерированные классом, (приведены в таблице), можно обьяснить линейную зависимость времени выполнения от возрастающей вычислительной нагрузки для нашего случая:

Дамп функции sin(X)+sin(X)+sin(X) отличается от дампа функции sin(X)+sin(X) тем же, чем дамп функции sin(X)+sin(X)+sin(X)+sin(X) отличается от дампа функции sin(X)+sin(X)+sin(X) - кодом D905BC005100D9FEDEC1.

Дамп для функции sin(X)+sin(X)
---------------------------------
D905BC005100 ;
загрузка XD9FE ;
вычисление sinD905BC005100 ;
загрузка XD9FE ;
вычисление sinDEC1 ;
сложениеC3 ;
возврат

Дамп для функций sin(X)+sin(X)+sin(X)
-----------------------------------------
D905BC005100 ;
загрузка XD9FE ;
вычисление sinD905BC005100 ;
загрузка XD9FE ;
вычисление sinDEC1 ;
сложениеD905BC005100 ;
загрузка XD9FE ;
вычисление sinDEC1 ;
сложениеC3 ;
возврат

Дамп для функций sin(X)+sin(X)+sin(X)+sin(X)
-----------------------------------------------
D905BC005100 ;
загрузка XD9FE ;
вычисление sinD905BC005100 ;
загрузка XD9FE ;
вычисление sinDEC1 ;
сложениеD905BC005100 ;
загрузка XD9FE ;
вычисление sinDEC1 ;
сложениеD905BC005100 ;
загрузка XD9FE ;
вычисление sinDEC1 ;
сложениеC3 ;
возврат

	Тестовая функция	Время выполнения дампа Этой функции (1000000 раз),сек.	Время вычисления функции методом интерпретации текстовой строки функции 100000 раз, сек.	Убыстрение, раз.
1	X	0.220	0.110	5
2	(exp(X)-exp(-X))/2	1.400	1.540	11
3	((exp(X) +0.4 )^ln(X))/(sin(X)+1)	1.860	2.960	16

Сравнение производилось на машине
Intel Сeleron 366MHZ,
В операционной системе MS Windows 98SE

PS: Код к статье качать отсюда

Версия для печати

Обсудить на форуме