Математические выражения в .NET (разбор, дифференцирование, упрощение, дроби, компиляция)

    Мне со школьных времен был интересен алгоритм вывода аналитических производных и упрощения выражений. Данная задача была актуальна впоследствии и в вузе. Тогда-то я реализовал ее, только получилось все не так, как хотелось: вместо кода IL у меня просто генерировался C# код в текстовом виде, сборки не выгружались, ну и к тому же не было возможности вывода производных в аналитическом виде. Однако потом я решил все-таки реализовать такую библиотеку, так как интерес остался. Стоит отметить, что таких библиотек в интернете большое количество, но нигде я не нашел именно этапа компиляции выражений в IL код, т.е. по сути везде выполняется интерпретация, которая не столь эффективна, в отличие от компиляции. Ну и к тому же я это разрабатывал чисто для себя, для изучения новых технологий, особо не надеясь, что результат моих трудов может где-нибудь потребоваться. Для нетерпеливых: исходники, программа.

    Используемые программы и библиотеки


    1. GOLD Parsing System — IDE для написания грамматик и генерации кода лексеров и парсеров под различные языки (C, C#, Java, JavaScript, Objective-C, Perl, Python, Ruby и др.). Основана на LALR парсинге.
    2. Visual Studio 2010
    3. GOLD.Engine — сборка под .NET, подключаемая для взаимодействия со сгенерированными таблицами.
    4. NUnit — Открытая среда юнит-тестирования приложений для .NET.
    5. ILSpy — OpenSource дизассемблер под .NET.

    Этапы, на которые я разбил весь процесс:
    1. Построение дерева выражения
    2. Вычисление аналитической производной
    3. Упрощение (симплификация) выражения
    4. Обработка рациональных дробей
    5. Компиляция выражения


    Построение дерева выражения


    Парсер-генератор GOLD я выбрал, поскольку уже имею опыт с ANTLR и захотелось чего-то нового. Ну а его преимущества, по сравнению с остальными, вы можете увидеть в данной таблице. Как видим, если сравнивать с ANTLR, то GOLD основан на алгоритме LALR, а не LL. А это значит, что в теории, сгенерированный парсер быстрее и мощнее, но с другой стороны, его нельзя отлаживать (в GOLD вообще подгружается файл в бинарном формате), а если бы и можно было бы, то это было бы совсем не очевидно и неудобно. Другим большим отличием является форма задания грамматики: BNF, а не EBNF. А это означает, что грамматика, записанная в данной форме, имеет немного больший размер из-за угловых скобок, знака определения и отсутствия условных вхождений и повторений (подробней в wiki. И, например, такое правило в EBNF, как

    ExpressionList = Expression {  ',' Expression }
    

    будет переписано в такой форме:

    <ExpressionList> ::= <ExpressionList> ',' <Expression> | <Expression>
    

    Итак, окончательная грамматика математических выражений имеет следующий вид:

    Грамматика математических выражений
    "Name"    = 'Mathematics expressions'
    "Author"  = 'Ivan Kochurkin'
    "Version" = '1.0'
    "About"   = ''
    
    "Case Sensitive" = False 
    "Start Symbol"  = <Statements>
    
    Id = {Letter}{AlphaNumeric}*
    
    Number1 = {Digit}+('.'{Digit}*('('{Digit}*')')?)?
    Number2 = '.'{Digit}*('('{Digit}*')')?
    
    AddLiteral = '+' | '-'
    MultLiteral = '*' | '/'
                                        
    <Statements> ::= <Statements> <Devider> <Statement>
                  | <Statements> <Devider>
                  | <Statement>              
    
    <Devider> ::= ';' | '.'
               
    <Statement> ::= <Expression> '=' <Expression>                           
                 | <Expression>
                                                             
    <Expression> ::= <FuncDef>
                  | <Addition>
     
    <FuncDef> ::= Id '(' <ExpressionList> ')'
               | Id '' '(' <ExpressionList> ')'
               | Id '(' <ExpressionList> ')' ''
                  
    <ExpressionList> ::= <ExpressionList> ',' <Expression>
                |  <Expression>
            
    <Addition> ::= <Addition> AddLiteral <Multiplication>
                | <Addition> AddLiteral <FuncDef>
                | <FuncDef> AddLiteral <Multiplication>
                | <FuncDef> AddLiteral <FuncDef> 
                | <Multiplication>           
    
    <Multiplication> ::= <Multiplication> MultLiteral <Exponentiation>
                    | <Multiplication> MultLiteral <FuncDef>  
                    | <FuncDef> MultLiteral <Exponentiation>
                    | <FuncDef> MultLiteral <FuncDef>
                    | <Exponentiation>
                    
    <Exponentiation> ::= <Exponentiation> '^' <Negation>
                    | <Exponentiation> '^' <FuncDef>
                    | <FuncDef> '^' <Negation>
                    | <FuncDef> '^' <FuncDef> 
                    | <Negation>
    
    <Negation> ::= AddLiteral <Value>
                    | AddLiteral <FuncDef>
                    | <Value>
    
    <Value> ::= Id
              | Number1
              | Number2
              | '(' <Expression> ')'
              | '|' <Expression> '|'
              | '(' <Expression> ')' ''
              | '|' <Expression> '|' ''
              | Id ''
    


    В грамматике вроде бы все очевидно, за исключением лексемы: Number1 = {Digit}+('.'{Digit}*('('{Digit}*')')?)? Данная конструкция позволяет парсить строки следующего вида: 0.1234(56), что означает рациональную дробь: 61111/495000. Про преобразование такой строки в дробь я расскажу позже.

    Итак, после того как скомпилированные грамматические таблицы (Compiled Grammar Table File) сгенерированы, и создан каркас класса парсера, необходимо модифицировать последний для того, чтобы строилось подходящее AST дерево. Каркас класса парсера в данном случае это .cs файл, в котором написан цикл обхода по всем всем правилам грамматики, обработка исключений в случае неверной последовательности лексем и других ошибок (кстати говоря в GOLD есть тоже разные генераторы таких каркасов, и я выбрал Cook .NET). Итак, в нашем случае «подходящее» обозначает дерево, состоящее из узлов, которые могут иметь типы:

    • Calculated — узел, представляющий рассчитанную константу в приемлемом для компилятора формате double, например «0.714285714285714», «0.122222222222222». Для каждой такой константы создается новый объект, даже если они одинаковые.
    • Value — узел, представляющий рассчитанную константу в формате рациональной дроби, т.е. 1 представляется как «1/1», «0.1(2)» — как «11/90», «0.1234» — как «617/5000». Для каждой такой константы создается новый объект, даже если они одинаковые.
    • Constant — узел, представляющий неопределенную константу, например «a», «b» и др. Если в выражении встречаются две таких константы, то они указывают на один узел.
    • Variable — узел, представляющий переменную, например «x», «y» и др. Если в выражении определенная переменная встречается несколько раз, то для нее создается только один объект, так же как и для неопределенной константы. Важно понимать, что разграничение константы и переменной важно только на этапе вывода аналитическое производной, а в других случаях это знание не обязательно.
    • Function — узел, представляющий функцию, например "+", «sin», "^" и другие. Единственный узел, который может иметь потомков.

    Все типы узлов наглядно представлены на рисунке:


    Из грамматики видно, что узлы полученного дерева либо не имеют детей (например, значения или переменные), либо имеют от одного до двух детей в случае стандартных функций (например, cos, сложение, возведение в степень). Все остальные функции могут иметь бОльшее количество аргументов, но они не рассматривались.

    Таким образом, все узлы (а точнее узлы с функциями) имеют от 0 до 2 детей включительно. Это верно с теоретической точки зрения. Но на самом деле пришлось сделать так, чтобы функции "+" и "*" имели неограниченное количество детей для того, чтобы облегчить задачу симплификации, которая будет рассмотрена позже. Кстати, от таких бинарных операций, как "-" и "/" тоже пришлось отказаться по таким же причинам (они были заменены на сложение с отрицанием правой части и умножением с обращением правой части).

    Итак, на этапе разбора все узлы для каждого правила попадают или извлекаются из буфера, который называется Nodes. Таким образом, в конце всего процесса в этом буфере остается одна или несколько функций с правой и левой частью. Также в процессе парсинга для того, чтобы функции сложения и умножения сразу создавались мультинарными, использовались дополнительные буферы ArgsCount и ArgsFuncTypes, хранящие текущее количество аргументов в функции и тип текущей функции соответственно. Таким образом, например для функции умножения, используется данный код:

    Обработка первого множителя:

    // <Multiplication> ::= <Exponentiation>
    if (MultiplicationMultiChilds)
    	PushFunc(KnownMathFunctionType.Mult);
    

    Обработка i-того множителя:

    // <Multiplication> ::= <Multiplication> MultLiteral <Exponentiation>
    // <Multiplication> ::= <Multiplication> MultLiteral <FuncDef>
    if (KnownMathFunction.BinaryNamesFuncs[r[1].Data.ToString()] == KnownMathFunctionType.Div)
    		Nodes.Push(new FuncNode(KnownMathFunctionType.Exp, new MathFuncNode[] { Nodes.Pop(), new ValueNode(-1) }));
    	ArgsCount[ArgsCount.Count - 1]++;
    

    Обработка последнего множителя:

    // <Addition> ::= <Multiplication>
    if (MultiplicationMultiChilds)
    	PushOrRemoveFunc(KnownMathFunctionType.Mult);
    if (AdditionMultiChilds)
    	PushFunc(KnownMathFunctionType.Add);
    


    Из этого кода видно, что, например, если никакого умножения и нет, то последнюю функции умножения с нулевым количеством элементов нужно удалить с помощью PushOrRemoveFunc. Аналогичные действия нужно провести и для сложения.

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

    Вычисление аналитической производной


    На данном этапе нужно преобразовать функцию в ее производную.

    Как известно из первого этапа, в созданном дереве выражения существуют только четыре типа узлов (пятый появляется позже). Для них и определим производную:

    • Value' = 0
    • Constant' = 0
    • Variable' = 1
    • Function' = Derivatives[Function]

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

    Список производных
    (f(x) ^ g(x))' = f(x) ^ g(x) * (f(x)' * g(x) / f(x) + g(x)' * ln(f(x)));");
    
    neg(f(x))' = neg(f(x)');");
    
    sin(f(x))' = cos(f(x)) * f(x)';
    cos(f(x))' = -sin(f(x)) * f(x)';
    tan(f(x))' = f(x)' / cos(f(x)) ^ 2;
    cot(f(x))' = -f(x)' / sin(f(x)) ^ 2;
    
    arcsin(f(x))' = f(x)' / sqrt(1 - f(x) ^ 2);
    arccos(f(x))' = -f(x)' / sqrt(1 - f(x) ^ 2);
    arctan(f(x))' = f(x)' / (1 + f(x) ^ 2);
    arccot(f(x))' = -f(x)' / (1 + f(x) ^ 2);
    
    sinh(f(x))' = f(x)' * cosh(x);
    cosh(f(x))' = f(x)' * sinh(x);
    
    arcsinh(f(x))' = f(x)' / sqrt(f(x) ^ 2 + 1);
    arcosh(f(x))' = f(x)' / sqrt(f(x) ^ 2 - 1);
    
    ln(f(x))' = f(x)' / f(x);
    log(f(x), g(x))' = g'(x)/(g(x)*ln(f(x))) - (f'(x)*ln(g(x)))/(f(x)*ln(f(x))^2);
    

    Стоит отметить, что все производные не записано жестко в коде, а тоже могут вводиться и парситься динамически.

    В этом списке нет сложения и умножения потому что, как упоминалось выше, это функции с несколькими аргументами. А для того, чтобы парсить такие правила, необходимо было бы еще написать много кода. По этим же причинам, здесь нет композиции функций, т.е. (f(g(x)))' = f(g(x))' * g(x)', а вместо этого все функции представлены как композиции функций.

    Также, если для функции не найдена подстановка в Derivatives (т.е. неопределенная функция), то она просто заменяется на функцию со штрихом: т.е. f(x) превратится в f(x)'.

    После того, как аналитическая производная успешно была получена, возникает проблема, которая заключается в том, что в получившееся выражении возникает очень много «мусора», такого как a + 0, a * 1, a * a^-1 и пр., а также получившееся выражение можно вычислить более оптимальным образом. Например, даже для простого выражения получится некрасивое выражение:

    (x^2 + 2)' = 0 + 2 * 1 * x ^ 1
    

    Для устранения таких таких недостатков используется упрощение.

    Упрощение (симплификация) выражения


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

    В начале топика я затронул тему про то, почему сложение и умножение представлено в виде n-арных функций. Представим ситуацию, изображенную на картинке ниже. Здесь видим, что a и -a сокращаются. Но как это сделать в случае бинарных функций? Для этого нужно перебирать узлы a, b и c, чтобы потом обнаружить, что a и -a являются детьми одного узла, а значит их можно сократить вместе с узлом. Но, понятное дело, что перебор деревьев это не такая простая задача, так как куда проще производить все действия в цикле сразу со всеми детьми, как это представлено на рисунке справа. Кстати, такой перебор нам позволяют сделать математические свойства ассоциативности и коммутативности.


    Во время процесса симплификации возникает задача сравнения двух узлов, которые, в свою очередь, могут иметь один из четырех типов. Это нужно, например, для того, чтобы сократить такие выражения как sin(x + y) и -sin(x + y). Понятно, что можно по узлам сравнить сами узлы и все их потомки. Но проблема заключается в том, что данный метод не справится с ситуацией, когда слагаемые или множители переставлены местами, например sin(x + y) и -sin(y + x). Для разрешения данной проблемы используется предварительная сортировка слагаемых или множителей, для которых выполняется свойство коммутативности (т.е. для сложения и вычитания). Сравниваются узлы так, как показано на картинке ниже, т.е. значения меньше констант, константы меньше переменных и т.д. Для функций все немного сложнее, поскольку нужно сравнить не только их названия, но и количество аргументов и сами аргументы.


    Таким образом, после всех вышеописанных преобразований и переборов, исходное выражение неплохо упрощается.

    Обработка рациональных дробей


    В найденных мною реализациях рациональных не было конвертации обычного строкового типа, например 0.666666, в тип дроби с определенным числителем и знаменателем, т.е. в 2 / 3 с определенной точностью. Тогда я решил все-таки написать свою реализацию с большим количеством опций. Функции ниже могут, например, определить, является ли какое-то число чисто иррациональным, или оно имеет период или конечную дробную часть и может быть преобразовано в рациональное, например sin(pi) в 0 при определенной точности. А вообще другие подробности смотрите в моем ответе на stackoverflow.com. Краткое графическое пояснения метода представлено на рисунке ниже, а код — в листинге ниже. Стоит отметить, что все-таки точности стандартных математических функций и типа double не хватает для нормального распознавания рациональных и вещественных чисел, однако теоретически все работает.


    Преобразование decimal в дробь
    /// <summary>
    /// Convert decimal to fraction
    /// </summary>
    /// <param name="value">decimal value to convert</param>
    /// <param name="result">result fraction if conversation is succsess</param>
    /// <param name="decimalPlaces">precision of considereation frac part of value</param>
    /// <param name="trimZeroes">trim zeroes on the right part of the value or not</param>
    /// <param name="minPeriodRepeat">minimum period repeating</param>
    /// <param name="digitsForReal">precision for determination value to real if period has not been founded</param>
    /// <returns></returns>
    public static bool FromDecimal(decimal value, out Rational<T> result, 
        int decimalPlaces = 28, bool trimZeroes = false, decimal minPeriodRepeat = 2, int digitsForReal = 9)
    {
        var valueStr = value.ToString("0.0000000000000000000000000000", CultureInfo.InvariantCulture);
        var strs = valueStr.Split('.');
    
        long intPart = long.Parse(strs[0]);
        string fracPartTrimEnd = strs[1].TrimEnd(new char[] { '0' });
        string fracPart;
    
        if (trimZeroes)
        {
            fracPart = fracPartTrimEnd;
            decimalPlaces = Math.Min(decimalPlaces, fracPart.Length);
        }
        else
            fracPart = strs[1];
    
        result = new Rational<T>();
        try
        {
            string periodPart;
            bool periodFound = false;
    
            int i;
            for (i = 0; i < fracPart.Length; i++)
            {
                if (fracPart[i] == '0' && i != 0)
                    continue;
    
                for (int j = i + 1; j < fracPart.Length; j++)
                {
                    periodPart = fracPart.Substring(i, j - i);
                    periodFound = true;
                    decimal periodRepeat = 1;
                    decimal periodStep = 1.0m / periodPart.Length;
                    var upperBound = Math.Min(fracPart.Length, decimalPlaces);
                    int k;
                    for (k = i + periodPart.Length; k < upperBound; k += 1)
                    {
                        if (periodPart[(k - i) % periodPart.Length] != fracPart[k])
                        {
                            periodFound = false;
                            break;
                        }
                        periodRepeat += periodStep;
                    }
    
                    if (!periodFound && upperBound - k <= periodPart.Length && periodPart[(upperBound - i) % periodPart.Length] > '5')
                    {
                        var ind = (k - i) % periodPart.Length;
                        var regroupedPeriod = (periodPart.Substring(ind) + periodPart.Remove(ind)).Substring(0, upperBound - k);
                        ulong periodTailPlusOne = ulong.Parse(regroupedPeriod) + 1;
                        ulong fracTail = ulong.Parse(fracPart.Substring(k, regroupedPeriod.Length));
                        if (periodTailPlusOne == fracTail)
                            periodFound = true;
                    }
    
                    if (periodFound && periodRepeat >= minPeriodRepeat)
                    {
                        result = FromDecimal(strs[0], fracPart.Substring(0, i), periodPart);
                        break;
                    }
                    else
                        periodFound = false;
                }
    
                if (periodFound)
                    break;
            }
    
            if (!periodFound)
            {
                if (fracPartTrimEnd.Length >= digitsForReal)
                    return false;
                else
                {
                    result = new Rational<T>(long.Parse(strs[0]), 1, false);
                    if (fracPartTrimEnd.Length != 0)
                        result = new Rational<T>(ulong.Parse(fracPartTrimEnd), TenInPower(fracPartTrimEnd.Length));
                    return true;
                }
            }
    
            return true;
        }
        catch
        {
            return false;
        }
    }
    
    public static Rational<T> FromDecimal(string intPart, string fracPart, string periodPart)
    {
        Rational<T> firstFracPart;
        if (fracPart != null && fracPart.Length != 0)
        {
            ulong denominator = TenInPower(fracPart.Length);
            firstFracPart = new Rational<T>(ulong.Parse(fracPart), denominator);
        }
        else
            firstFracPart = new Rational<T>(0, 1, false);
    
        Rational<T> secondFracPart;
        if (periodPart != null && periodPart.Length != 0)
            secondFracPart =
                new Rational<T>(ulong.Parse(periodPart), TenInPower(fracPart.Length)) *
                new Rational<T>(1, Nines((ulong)periodPart.Length), false);
        else
            secondFracPart = new Rational<T>(0, 1, false);
    
        var result = firstFracPart + secondFracPart;
        if (intPart != null && intPart.Length != 0)
        {
            long intPartLong = long.Parse(intPart);
            result = new Rational<T>(intPartLong, 1, false) + (intPartLong == 0 ? 1 : Math.Sign(intPartLong)) * result;
        }
    
        return result;
    }
    
    private static ulong TenInPower(int power)
    {
        ulong result = 1;
        for (int l = 0; l < power; l++)
            result *= 10;
        return result;
    }
    
    private static decimal TenInNegPower(int power)
    {
        decimal result = 1;
        for (int l = 0; l > power; l--)
            result /= 10.0m;
        return result;
    }
    
    private static ulong Nines(ulong power)
    {
        ulong result = 9;
        if (power >= 0)
            for (ulong l = 0; l < power - 1; l++)
                result = result * 10 + 9;
        return result;
    }
    


    Компиляция выражения


    Для преобразования полученного семантического дерева в код IL после этапа симплификации, я использовал Mono.Cecil.

    В начале процесса создается сборка, класс и метод, в который будут записываться команды. Потом для каждого FuncNode рассчитывается, сколько он раз встречается в программе. Например, если есть функция sin(x^2) * cos(x^2), то в ней функция возведения x в степень 2 встречается два раза, а функции sin и cos — по одному. В дальнейшем, данная информация о повторах вычислений функций используется следующим образом (т.е. таким образом второго раза вычисления одной и той же функции не происходит):

    if (!func.Calculated)
    {
    	EmitFunc(funcNode);
    	func.Calculated = true;
    }
    else
    	IlInstructions.Add(new OpCodeArg(OpCodes.Ldloc, funcNode.Number));
    


    После генерации всего этого кода, возможны еще некоторые оптимизации IL кода, представленные на рисунке ниже.

    На данном рисунке:
    1. Ldarg — Операция загрузки определенного аргумента функции в стек.
    2. Ldc_R8 — Загрузка определенного значения double в стек.
    3. Stloc.n — Извлечение из стека последнего значения и его сохранение в локальную переменную n.
    4. Ldloc.n — Помещение локальной переменной n в стек.

    Инструкции в бежевых прямоугольниках можно удалить, если выполнены определенные условия. Например, случай на левом верхнем изображении описывается так: Если текущая инструкция — это загрузка аргумента из функции или загрузка константы, и следующая инструкция — ее сохранение в локальную переменную n, то данный блок инструкций можно удалить, предварительно заменив инструкции загрузки локальной переменной n на загрузку аргумента функции или загрузки константы. Процесс замены продолжать до первой инструкции сохранения в локальную переменную n. Аналогичным образом объясняются остальные три случая. например последовательность Ldloc.n; Stloc.n сразу может быть удалена.

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

    Быстрое возведение в степень

    Я думаю, что почти все знают про алгоритм быстрого возведения числа в степень. Но ниже представлен данный алгоритм на уровне компиляции:
    Алгоритм быстрого возведения в степень, реализованный в IL коде
    if (power <= 3)
    {
    	IlInstructions.Add(new OpCodeArg(OpCodes.Stloc, funcNode.Number));
    	IlInstructions.Add(new OpCodeArg(OpCodes.Ldloc, funcNode.Number));
    	for (int i = 1; i < power; i++)
    	{
    		IlInstructions.Add(new OpCodeArg(OpCodes.Ldloc, funcNode.Number));
    		IlInstructions.Add(new OpCodeArg(OpCodes.Mul));
    	}
    }
    else if (power == 4)
    {
    	IlInstructions.Add(new OpCodeArg(OpCodes.Stloc, funcNode.Number));
    	IlInstructions.Add(new OpCodeArg(OpCodes.Ldloc, funcNode.Number));
    	IlInstructions.Add(new OpCodeArg(OpCodes.Ldloc, funcNode.Number));
    	IlInstructions.Add(new OpCodeArg(OpCodes.Mul));
    	IlInstructions.Add(new OpCodeArg(OpCodes.Stloc, funcNode.Number));
    	IlInstructions.Add(new OpCodeArg(OpCodes.Ldloc, funcNode.Number));
    	IlInstructions.Add(new OpCodeArg(OpCodes.Ldloc, funcNode.Number));
    	IlInstructions.Add(new OpCodeArg(OpCodes.Mul));
    }
    else
    {
    	// result: funcNode.Number
    	// x: funcNode.Number + 1
    
    	//int result = x;
    	IlInstructions.Add(new OpCodeArg(OpCodes.Stloc, funcNode.Number + 1));
    	IlInstructions.Add(new OpCodeArg(OpCodes.Ldloc, funcNode.Number + 1));
    
    	power--;
    	do
    	{
    		if ((power & 1) == 1)
    		{
    			IlInstructions.Add(new OpCodeArg(OpCodes.Ldloc, funcNode.Number + 1));
    			IlInstructions.Add(new OpCodeArg(OpCodes.Mul));
    		}
    
    		if (power <= 1)
    			break;
    
    		//x = x * x;
    		IlInstructions.Add(new OpCodeArg(OpCodes.Ldloc, funcNode.Number + 1));
    		IlInstructions.Add(new OpCodeArg(OpCodes.Ldloc, funcNode.Number + 1));
    		IlInstructions.Add(new OpCodeArg(OpCodes.Mul));
    		IlInstructions.Add(new OpCodeArg(OpCodes.Stloc, funcNode.Number + 1));
    
    		power = power >> 1;
    	}
    	while (power != 0);
    }
    


    Стоит отметить, что алгоритм быстрого возведения в степень не является самым оптимальным алгоритмом возведения в степень. Например, ниже представлено перемножение одинаковых переменных пять раз двумя способами. Оптимизации типа x^4 + x^3 + x^2 + x = x*(x*(x*(x + 1) + 1) + 1) тоже у меня не реализованы.

    • var t = a ^ 2; a*a*a*a*a*a = t ^ 2 * t — обычный «быстрый» алгоритм.
    • a*a*a*a*a*a = (a ^ 3) ^ 2 — оптимальный «быстрый» алгоритм.


    Кстати, еще стоит упомянуть, что для обычных фиксированных чисел double, результат перемножения для упомянутых выше перемножений в разном порядке, будет различный (т.е. (a ^ 2) ^ 2 * a ^ 2 != (a ^ 3) ^ 2). И из-за этого некоторые компиляторы не оптимизируют многие подобные выражения. По этому поводу есть интересные Q&A в stackoverlofw: Why doesn't GCC optimize a*a*a*a*a*a to (a*a*a)*(a*a*a)? и Why Math.Pow(x, 2) not optimized to x * x neither compiler nor JIT?.

    Оптимизация локальных переменных

    Как известно из предыдущих этапов, локальные переменные используются для хранения результата всех функций, которые встречаются в исходном выражении более 1 раза. Для работы с локальными переменными используются всего две простые инструкции: stloc и ldloc, которые используют один аргумент, отвечающий за номер данной локальной переменной. Но если номер локальной переменной инкрементировать каждый раз при появлении повторяющегося результата вычислений (для ее создания), то локальных переменных может быть очень много. Для нивелирования данной проблемы, был реализован алгоритм сжатия жизненных циклов локальных переменных, процесс которого можно наглядно увидеть на рисунке ниже. Как видим, вместо 5 локальных переменных в выражении можно использовать всего 3. Однако этот «жадный» алгоритм не является самым оптимальным перестановки, однако он вполне подходит к реализуемой задаче.

    Компиляция неопределенных функций и их производных

    В разработанной библиотеке можно использовать не только простые функции от одного переменного вида f(x), но и другие, такие как f(x,a,b(x)), где a — неизвестная константа, а b(x) — неизвестная функция, передаваемая в виде делегата. Как известно, определение производной от любой функции выглядит следующим образом: b(x)' = (b(x + dx) — b(x)) / dx. И когда модуль компиляции встречает неопределенную функцию, он генерирует следующий код:
    ldarg.1
    ldarg.0
    callvirt TResult System.Func`2<System.Double,System.Double>::Invoke(T)
    ret
    

    Код вычисления производной неопределенной функции (dx = 0.000001)
    ldarg.1
    ldarg.0
    ldc.r8 1E-06
    add
    callvirt TResult System.Func`2<System.Double,System.Double>::Invoke(T)
    ldarg.1
    ldarg.0
    callvirt TResult System.Func`2<System.Double,System.Double>::Invoke(T)
    sub
    ldc.r8 0.000001
    div
    ret
    


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



    Тестирование


    Тестирование проводилось для нескольких этапов обработки, которые можно посмотреть в проекте MathFunction.Tests. Из интересных моментов можно упомянуть использование WolframAlpha.NET для тестирования аналитических производных и загрузку и выгрузку сборок с помощью доменов.

    WolframAlpha

    WolframAlpha.NET — это API оболочка известного математического сервиса wolframalpha. В моем проекте она пригодилась для того, чтобы сверять заведомо правильную производную, полученную на этом сервисе, с производной, полученной с помощью моей библиотеки. Использовалось следующим образом:

    public static bool CheckDerivative(string expression, string derivative)
    {
    	return CheckEquality("diff(" + expression + ")", derivative);
    }
    
    public static bool CheckEquality(string expression1, string expression2)
    {
    	WolframAlpha wolfram = new WolframAlpha(ConfigurationManager.AppSettings["WolframAlphaAppId"]);
    	string query = "(" + expression1.Replace(" ", "") + ")-(" + expression2.Replace(" ", "") + ")";
    	QueryResult result = wolfram.Query(query);
    	result.RecalculateResults();
    	try
    	{
    		double d;
    		return double.TryParse(result.GetPrimaryPod().SubPods[0].Plaintext, out d) && d == 0.0;
    	}
    	catch
    	{
    		return false;
    	}
    }
    


    Загрузка и выгрузка сборок с помощью доменов

    Создание домена, сборки из файла, создание экземпляра определенного типа и получение объектов с информацией о методах (MethodInfo):

    Domain = AppDomain.CreateDomain("MathFuncDomain");
    MathFuncObj = _domain.CreateInstanceFromAndUnwrap(FileName, NamespaceName + "." + ClassName);
    Type mathFuncObjType = _mathFuncObj.GetType();
    Func = mathFuncObjType.GetMethod(FuncName);
    FuncDerivative = mathFuncObjType.GetMethod(FuncDerivativeName);
    


    Вычисление результата скомпилированной функции:
    return (double)Func.Invoke(_mathFuncObj, new object[] { x })
    


    Выгрузка домена и удаление файла со сборкой:
    if (_domain != null)
        AppDomain.Unload(Domain);
    File.Delete(FileName);
    


    Сравнение сгенерированного IL кода

    Финальный IL код является более оптимальным, по сравнению с кодом, сгенерированным стандартным C# компилятором csc.exe в Release режиме, достаточно всего лишь взглянуть на сравнение следующих двух листингов например для функции
    x ^ 3 + sin(3 * ln(x * 1)) + x ^ ln(2 * sin(3 * ln(x))) - 2 * x ^ 3

    csc.exe .NET 4.5.1 MathExpressions.NET
    IL_0000: ldarg.0
    IL_0001: ldarg.0
    IL_0002: mul
    IL_0003: ldarg.0
    IL_0004: mul
    IL_0005: ldc.r8 3
    IL_000e: ldarg.0
    IL_000f: ldc.r8 1
    IL_0018: mul
    IL_0019: call float64 [mscorlib]System.Math::Log(float64)
    IL_001e: mul
    IL_001f: call float64 [mscorlib]System.Math::Sin(float64)
    IL_0024: add
    IL_0025: ldarg.0
    IL_0026: ldc.r8 2
    IL_002f: ldc.r8 3
    IL_0038: ldarg.0
    IL_0039: call float64 [mscorlib]System.Math::Log(float64)
    IL_003e: mul
    IL_003f: call float64 [mscorlib]System.Math::Sin(float64)
    IL_0044: mul
    IL_0045: call float64 [mscorlib]System.Math::Log(float64)
    IL_004a: call float64 [mscorlib]System.Math::Pow(float64, float64)
    IL_004f: add
    IL_0050: ldc.r8 2
    IL_0059: ldarg.0
    IL_005a: mul
    IL_005b: ldarg.0
    IL_005c: mul
    IL_005d: ldarg.0
    IL_005e: mul
    IL_005f: sub
    IL_0060: ret
    

    IL_0000: ldc.r8 3
    IL_0009: ldarg.0
    IL_000a: call float64 [mscorlib]System.Math::Log(float64)
    IL_000f: mul
    IL_0010: call float64 [mscorlib]System.Math::Sin(float64)
    IL_0015: stloc.0
    IL_0016: ldloc.0
    IL_0017: ldarg.0
    IL_0018: ldarg.0
    IL_0019: mul
    IL_001a: ldarg.0
    IL_001b: mul
    IL_001c: sub
    IL_001d: ldarg.0
    IL_001e: ldc.r8 2
    IL_0027: ldloc.0
    IL_0028: mul
    IL_0029: call float64 [mscorlib]System.Math::Log(float64)
    IL_002e: call float64 [mscorlib]System.Math::Pow(float64, float64)
    IL_0033: add
    IL_0034: ret
    


    Любопытно при этом посмотреть, какой при этом получается дизассемблированный C# код в ILSpy для функции, являющейся производной к вышеупомянутой функции, т.е. к
    (ln(2 * sin(3 * ln(x))) * x ^ -1 + 3 * ln(x) * cos(3 * ln(x)) * sin(3 * ln(x)) ^ -1 * x ^ -1) * x ^ ln(2 * sin(3 * ln(x))) + 3 * cos(3 * ln(x)) * x ^ -1 + -(3 * x ^ 2)

    double arg_24_0 = 2.0;
    double arg_1A_0 = 3.0;
    double num = Math.Log(x);
    double num2 = arg_1A_0 * num;
    double num3 = Math.Sin(num2);
    double num4 = Math.Log(arg_24_0 * num3);
    double arg_3B_0 = num4;
    double num5 = 1.0 / x;
    double arg_54_0 = arg_3B_0 * num5;
    double arg_4F_0 = 3.0 * num;
    num = Math.Cos(num2);
    return (arg_54_0 + arg_4F_0 * num / num3 * num5) * Math.Pow(x, num4) + num * num5 * 3.0 - x * x * 3.0;
    

    Как видим, создалось большое количество локальных переменных для хранения уже рассчитанных результатов функций. Хотя в реальности локальных переменных меньше, потому что, например, для double arg_24_0 = 2.0; не создается локальная переменная, это просто константа.

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

    Заключение


    Во время реализации всего вышеперечисленного на C#, я понял, что функциональный язык больше подходит для таких задач, где требуется произвести что-то связанное с символьными вычислениями и подстановкой. Даже, например, OpenSource система компьютерной алгебры maxima написана на lisp. Кстати, на F# есть реализация символьных вычислений на codeproject, и кода там явно меньше, чем у меня в проекте. Однако реализация данного проекта позволила мне получить больше опыта в разборе выражений, оптимизации, генерации низкоуровневого кода, работе с числами с плавающей точкой и просто в математике.

    Для всех заинтересованных я выложил ее на Github: source. Также и доступна сама программа: MathExpressions.NET. При острой необходимости могу ее отрефакторить. Пулл-реквесты приветствуются.

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

    UPDATE: Так как никто не отгадал, что зашифровано в сердечке, то выкладываю результат. Данная формула является частным случаем многочлена Матиясевича, множество положительных значений которых при неотрицательных значениях переменных совпадает с множеством простых чисел. На википедии также есть информация.
    Share post
    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More
    Ads

    Comments 34

      +4
      Основательный подход, великолепная статья, и вдвойне приятно, что не перевод!

      P.S. отдельное спасибо за грамматики — обновил знания.
        +5
          0
          1. В .NET нет длинной арифметики. В Java плохая по эффективности реализация (плата за стандарт). Самая хорошая здесь (всегда заморочка при сборке под Windows)
          gmplib.org/

          2. Lisp как язык программирования лучше подходит для универсальных (есть специализированные. например, Singular, Magma) систем компьютерной алгебры. Так Maxima (произошла от Macsyma) и Reduce значительно быстрее написанных на C и наиболее распространенных Maple, Mathematica. Классический коэффициент, после компиляции Lisp в 5-6 раз медленнее C, но математикам удобнее писать на Lisp свои алгоритмы, чем объяснять их программистам для встраивания их на C.

          3. В Reduce встроен следующий оператор
          decompose(xˆ8-88*xˆ7+2924*xˆ6-43912*xˆ5+263431*xˆ4-218900*xˆ3+65690*xˆ2-7700*x+234)
          -> {U^2 + 35*U + 234, U=V^2 + 10*V, V=X^2 — 22*X}

          decompose(uˆ2+vˆ2+2u*v+1)
          -> {W^2 + 1, W=U + V}

          На Reduce также есть очень неплохой пакет scope. Пример его использования
          z:=a^2*b^2+10*a^2*m^6+a^2*m^2+2*a*b*m^4+2*b^2*m^6+b^2*m^2;
          Результат
          g1 := b*a
          g5 := m*m
          g2 := g5*b*b
          g3 := g5*a*a
          g4 := g5*g5
          z := g2 + g3 + g1*(2*g4 + g1) + g4*(2*g2 + 10*g3)

            +4
            > 1. В .NET нет длинной арифметики…

            System.Numerics.BigInteger, на хабре была статья.
              0
              Согласен исправить фразу.
              1. В .NET нет эффективной реализации длинной арифметики.…

              А такую по Вашей ссылке можно написать на С за 2-3 дня если знаешь алгоритмы. Если мучится с ассемблером под разные процессоры, то просто долго. Если на GMP, что я дал ссылку несколько лет и небольшой коллектив с хорошим математическим образованием и алгоритмическим мышлением. Maple несколько лет назад тоже перешел GMP со своей реализации. Вроде в Magma хорошая реализация длинной арифметики, но там все коды закрыты и об этом можно судить косвенно по скорости выполнения символьных вычислений.
                +4
                В длинной арифметике ничего сложного нет. Сами в школе писали, как говорится. А C# очень часто проигрывает чистому Си. Например, BigInteger реализованы как immutable, что очень негативно сказывается на производительности. Преимущество C# в том что программист может быстро сделать, а не в том что код будет идеально быстро работать.
                  0
                  1. Наивный алгоритм умножения O(n^2) проигрывает методу Карацубы O(n^log_2(3)=1.58..) при правильной реализации при размере чисел в 7-8 машинных слов. Умножения Шёнхаге — Штрассена O(N·logN·log logN) использует быстрое преобразование Фурье. Примерно 100-400 машинных слов (сильно зависит от архитектуры). Алгоритм Фюрера 2007 год O(N·logN). Есть еще обобщение алгоритма Карацубы алгоритм Тоома-Кука.
                  2. Деление на число с машинное слово, на число в два машинных слова. Для больших чисел целочисленное деление заменяют умножением (в зависимости от размера соответствующий быстрый алгоритм умножения) на обратный элемент, который получают используя метод Ньютона.
                  3. Вычисление НОД (gcd) постоянно необходимо при реализации рациональных чисел. Тут бинарный алгоритм.
                  4. Про факторизацию и дискретное логарифмирование просто молчу.
                  5. Эффективная реализация китайской теоремы об остатках очень нетривиальная задача. Необходима для поддержки модуляроной арифметики.

                  Ну и где тут школа?
                    0
                    Это не про меня.
                    ALU не проектирую, криптографией не занимаюсь, теорию чисел знаю не намного лучше продвинутого школьника. Хватило спортивного программирования 12 лет назад.
                    Речь о том, что когда нужен калькулятор с длинной арифметикой, то BigInteger будет достаточно. А для числодробилки нужны не только алгоритмы, но может еще и FPGA сгодится.
                    Полюбопытствую, а Вам зачем такие подробности длинной арифметики знать?
                      0
                      FPGA приделать достаточно сложно и вряд ли имеет смысл (тут все плохо параллелится), ну разве модулярная арифметика.
                      Одно из моих увлечений разработка алгоритмов в области компьютерной алгебры. До версии GMP < 2.0 имел собственную библиотеку длинной арифметики, которая была быстрее. Но у меня другое направление, хотя воспоминания остались очень приятные. Вдобавок это был мой первый эксперимент по массовому применению метода утверждений (Тройки Хоара) + ООП на C++.
                      С тех пор постоянно использую два подхода. Криптографией тоже увлекался и придумал пару схем открытого ключа. Но на все время не хватает. Например нет хорошей открытой реализации на C/C++ вычисления НОД (gcd) для полиномов от нескольких переменных. Хорошая есть в Maple и Magma. Тут лучший алгоритм
                      www.cecm.sfu.ca/CAG/papers/brown.ps
                      Мало кто знает, что на заре компьютерной техники в ЭВМ при проектировании сразу закладывались возможности эффективной реализации длинной арифметики.
                  0
                  В .NET нет эффективной реализации длинной арифметики.…

                  В чистом .NET — нет. Однако есть библиотека под .NET, которая называется IntXLib (разработчик, кстати, из Украины), в которой используется дискретное преобразование Хартли для эффективной реализации операций с большими числами.

                  Понятное дело, что GMP будет быстрее .NET версии, потому что в ней, наверное, непосредственно используются специальные команды ассемблера, векторизация, более агрессивная оптимизация (C/C++). В проекте IntXLib есть сравнение производительности с GMP.
              0
              Интересный проект. Есть мысль попытаться добавить поддержку кастомных функций.
                +1
                Ну кастомные функции есть, и даже производные от них берутся в аналитическом виде на этапе дифференцирования (добавляется штрих) и численном виде на этапе компиляции (по определению производной с заданной delta) следующим образом:

                using (var mathAssembly = new MathAssembly("b(x) + 10 * x * a", "x"))
                {
                  var b = new Func<double, double>(x => x * x);
                  var funcResult = mathAssembly.Func(5, 2, b); // x = 5; a = 2; b = x ^ 2
                  // funcResult == 5 ^ 2 + 10 * 5 * 2 = 125
                  var funcDerResult = mathAssembly.FuncDerivative(5, 2, b); // x = 5; a = 2; b = x ^ 2
                  // funcDerResult == (b(x + dx) - b(x)) / dx + 10 * a = 30
                }
                

                Или вы что-то другое имеете в виду? :)
                0
                Очень хорошая задача для глубокого понимания многих вопросов CS.

                В реальности, для максимального использования проверенного кода, можно сделать проще:
                1. Парсить выражения можно с помощью Roslyn (но только стандартные выражения C# ?)
                2. Можно парсить выражение в Expression<Func<double[], double>> и затем вызывать .Compile();

                А вот в аналитических преобразованиях я дальше производных тоже не двинулся. Интересно было бы посмотреть LGPL код подобной библиотеки.
                  0
                  Кроме того, можно еще парсить и с помощью NRefactory, который к тому же еще и OpenSource. И я, кстати, использую его для еще одного «теоретического» проекта для минификации C# кода: CSharp-Minifier. А он в свою очередь используется для минификации кода в квайнах. Об этом, впрочем, я собираюсь написать в будущем.

                  Но использование Roslyn и NRefactory для математических выражений — большой оверхед, и они не позволяют работать с IL кодом напрямую для лучшей оптимизации.

                  А вот в аналитических преобразованиях я дальше производных тоже не двинулся. Интересно было бы посмотреть LGPL код подобной библиотеки.

                  Из наиболее известных и мощных есть система Axiom с документацией и исходниками.

                  Однако стоит отметить, что даже упрощение выражений (т.е. приведение выражению к виду, содержащему наименьшее количество операций), является довольно сложной задачей, где нужно так или иначе производить перебор (с эвристиками или без). А если говорить об аналитических интегралах и другой символьной математике, то это уже чуть ли не творчество :)
                  +1
                  В сердечке трудно различить L и 1. Видны какие-то формулы из теоремы косинусов, трёхгранные углы… но некоторые формулы непонятны вообще. И главное, непонятно, зачем сумму квадратов (или 4-х степеней?) вычитать из 1 :) Забавно, что слагаемые можно упорядочить так, что в каждом следующем появляется хотя бы одна новая переменная. В целом, выглядит как какая-то сложная конфигурация точек и углов на сфере.
                    0
                    Хм, хорошая попытка, но не угадали совсем. Вероятно, отгадать сможет только тот, кто видел ее хотя бы раз. Я скажу о правильном результате позже, когда кто-нибудь отгадает или пройдет еще какое-то время :)
                      0
                      Можно ещё предположить, что это что-то, связанное с квантовыми состояниями атома гелия. Но там я совсем не разбираюсь.
                        0
                        Тоже нет, это чистая математика. Добавил правильный ответ в конец статьи.
                    +2
                    Финальный IL код является более оптимальным, по сравнению с кодом, сгенерированным стандартным C# компилятором csc.exe в Release режиме, достаточно всего лишь взглянуть на сравнение следующих двух листингов например для функции

                    Стоит отметить, что csc практически ничего не оптимизирует, т.к. все оптимизации возложены на JIT. Он умеет совершать только совсем небольшие оптимизации. Стоит посмотреть на результирующий ассемблер (и скорость его выполнения, конечно же). Прилагаю ассемблер для случая x64 Release:
                    00000000  push        ebp 
                    00000001  mov         ebp,esp 
                    00000003  sub         esp,10h 
                    00000006  fld         qword ptr [ebp+8] 
                    00000009  fld         st(0) 
                    0000000b  fld1 
                    0000000d  fmulp       st(1),st 
                    0000000f  sub         esp,8 
                    00000012  fstp        qword ptr [esp] 
                    00000015  fstp        qword ptr [ebp-8] 
                    00000018  call        723800FB 
                    0000001d  fld         qword ptr [ebp-8] 
                    00000020  fxch        st(1) 
                    00000022  fmul        dword ptr ds:[013A3608h] 
                    00000028  fsin 
                    0000002a  fld         st(1) 
                    0000002c  fmul        st,st(2) 
                    0000002e  fmul        st,st(2) 
                    00000030  faddp       st(1),st 
                    00000032  fld         st(1) 
                    00000034  sub         esp,8 
                    00000037  fstp        qword ptr [esp] 
                    0000003a  fld         st(1) 
                    0000003c  sub         esp,8 
                    0000003f  fstp        qword ptr [esp] 
                    00000042  fstp        qword ptr [ebp-8] 
                    00000045  fstp        qword ptr [ebp-10h] 
                    00000048  call        723800FB 
                    0000004d  fld         qword ptr [ebp-10h] 
                    00000050  fxch        st(1) 
                    00000052  fmul        dword ptr ds:[013A3610h] 
                    00000058  fsin 
                    0000005a  fmul        dword ptr ds:[013A3618h] 
                    00000060  sub         esp,8 
                    00000063  fstp        qword ptr [esp] 
                    00000066  fstp        qword ptr [ebp-10h] 
                    00000069  call        723800FB 
                    0000006e  fld         qword ptr [ebp-10h] 
                    00000071  sub         esp,8 
                    00000074  fxch        st(1) 
                    00000076  fstp        qword ptr [esp] 
                    00000079  fstp        qword ptr [ebp-10h] 
                    0000007c  call        724D0466 
                    00000081  fld         qword ptr [ebp-10h] 
                    00000084  fld         qword ptr [ebp-8] 
                    00000087  faddp       st(2),st 
                    00000089  fld         st(0) 
                    0000008b  fmul        dword ptr ds:[013A3620h] 
                    00000091  fmul        st,st(1) 
                    00000093  fmulp       st(1),st 
                    00000095  fsubp       st(1),st 
                    00000097  mov         esp,ebp 
                    00000099  pop         ebp 
                    0000009a  ret         8 
                    
                      +2
                      Хотя забавно, что debug на моей машине выдает более быстрый код (юзает SSE) :)
                      00000038  vmovsd      xmm0,qword ptr [rbp+000000B0h] 
                      00000041  vmovsd      xmm1,qword ptr [rbp+000000B0h] 
                      0000004a  vmulsd      xmm0,xmm0,xmm1 
                      0000004f  vmovsd      xmm1,qword ptr [rbp+000000B0h] 
                      00000058  vmulsd      xmm0,xmm0,xmm1 
                      0000005d  vmovsd      qword ptr [rbp+78h],xmm0 
                      00000063  vmovsd      xmm0,qword ptr [000001B8h] 
                      0000006c  vmovsd      qword ptr [rbp+70h],xmm0 
                      00000072  vmovsd      xmm0,qword ptr [rbp+000000B0h] 
                      0000007b  vmulsd      xmm0,xmm0,mmword ptr [000001C0h] 
                      00000084  call        000000005FCD3BA0 
                      00000089  vmovsd      qword ptr [rbp+68h],xmm0 
                      0000008f  vmovsd      xmm0,qword ptr [rbp+70h] 
                      00000095  vmovsd      xmm1,qword ptr [rbp+68h] 
                      0000009b  vmulsd      xmm0,xmm0,xmm1 
                      000000a0  call        000000005F65C6F0 
                      000000a5  vmovsd      qword ptr [rbp+60h],xmm0 
                      000000ab  vmovsd      xmm0,qword ptr [rbp+78h] 
                      000000b1  vmovsd      xmm1,qword ptr [rbp+60h] 
                      000000b7  vaddsd      xmm0,xmm0,xmm1 
                      000000bc  vmovsd      qword ptr [rbp+58h],xmm0 
                      000000c2  vmovsd      xmm0,qword ptr [rbp+000000B0h] 
                      000000cb  vmovsd      qword ptr [rbp+50h],xmm0 
                      000000d1  vmovsd      xmm0,qword ptr [000001C8h] 
                      000000da  vmovsd      qword ptr [rbp+48h],xmm0 
                      000000e0  vmovsd      xmm0,qword ptr [000001D0h] 
                      000000e9  vmovsd      qword ptr [rbp+40h],xmm0 
                      000000ef  vmovsd      xmm0,qword ptr [rbp+000000B0h] 
                      000000f8  call        000000005FCD3BA0 
                      000000fd  vmovsd      qword ptr [rbp+38h],xmm0 
                      00000103  vmovsd      xmm0,qword ptr [rbp+40h] 
                      00000109  vmovsd      xmm1,qword ptr [rbp+38h] 
                      0000010f  vmulsd      xmm0,xmm0,xmm1 
                      00000114  call        000000005F65C6F0 
                      00000119  vmovsd      qword ptr [rbp+30h],xmm0 
                      0000011f  vmovsd      xmm0,qword ptr [rbp+48h] 
                      00000125  vmovsd      xmm1,qword ptr [rbp+30h] 
                      0000012b  vmulsd      xmm0,xmm0,xmm1 
                      00000130  call        000000005FCD3BA0 
                      00000135  vmovsd      qword ptr [rbp+28h],xmm0 
                      0000013b  vmovsd      xmm0,qword ptr [rbp+50h] 
                      00000141  vmovsd      xmm1,qword ptr [rbp+28h] 
                      00000147  call        000000005FCD3BC8 
                      0000014c  vmovsd      qword ptr [rbp+20h],xmm0 
                      00000152  vmovsd      xmm0,qword ptr [rbp+58h] 
                      00000158  vmovsd      xmm1,qword ptr [rbp+20h] 
                      0000015e  vaddsd      xmm0,xmm0,xmm1 
                      00000163  vmovsd      xmm1,qword ptr [rbp+000000B0h] 
                      0000016c  vmulsd      xmm1,xmm1,mmword ptr [000001D8h] 
                      00000175  vmovsd      xmm2,qword ptr [rbp+000000B0h] 
                      0000017e  vmulsd      xmm1,xmm1,xmm2 
                      00000183  vmovsd      xmm2,qword ptr [rbp+000000B0h] 
                      0000018c  vmulsd      xmm1,xmm1,xmm2 
                      00000191  vsubsd      xmm0,xmm0,xmm1 
                      00000196  vmovsd      qword ptr [rbp+00000080h],xmm0 
                      0000019f  nop 
                      000001a0  jmp         00000000000001A2 
                      000001a2  vmovsd      xmm0,qword ptr [rbp+00000080h] 
                      000001ab  lea         rsp,[rbp+00000090h] 
                      000001b2  pop         rsi 
                      000001b3  pop         rdi 
                      000001b4  pop         rbp 
                      000001b5  ret 
                      
                        0
                        Да, он юзает SSE, но при этом вызовов call всего 4 (как и в моей оптимизированной версии), в отличие от оптимизированной версии, где их 6. Так что быстрее или нет это еще вопрос.

                        Но, в любом случае, спасибо. Протестирую и добавлю листинги чистого ассемблера в статью.
                        +1
                        А как вы получили этот ассемблер? Если из студии, то там он тоже неоптимизированный. Нужно как-то через WinDbg или через еще что-то.
                          0
                          Debug -> Options

                          1) Снимаем галку Enable Just My Code
                          2) Снимаем галку Suppress JIT Optimisations

                          Можно еще всяких галок понаставить, но эти основные
                            0
                            Нет, это не то. Функции, как минимум, не инлайнятся, это объяснимо почему: не будут работать точки останова внутри них. Может сейчас, конечно, как-то по-другому, но когда я тестил, все было так, и без дебаггера работало быстрее (Ctrl + F5).
                              +2
                              Это то.
                                      static void Main(string[] args)
                                      {
                                          Console.WriteLine(Get4());
                                          Console.ReadKey();
                                      }
                              
                                      private static int Get4()
                                      {
                                          return 4;
                                      }
                              

                              Выполните сказанные мной манипуляции, поставьте сборку release и бряк на return 4. Он не сработает, сразу выведется 4
                        0
                        habrahabr.ru/post/270443 — я ссылаюсь на Вас в своей статье. Было бы интересно если Вы предложите дополнительные идеи.
                        Кстати — там может даже будет нужно ухудшать код программы.
                        Обработка приватных данных на публичных вычислительных сетях
                        С удовольствием парочку интересных идей украл бы
                          0
                          Хм, спасибо, интересно.
                          Буду изучать ваши материалы.
                            0
                            абстрактный вычислитель Маркова — замена в сроке символов шаблона на последовательность.

                            абстрактный вычислитель Маркова, абстрактный вычислитель Тьюринга, машина Поста — всё они эквивалентны.
                            Соответственно, это модель когда на ленте абстрактного вычислителя записан текст программы.
                            Шаги замен переводят в эквивалентны тексты программ.
                            Соответственно можно искать оптимизированную программ в каком либо смысле — например — меньше команд или быстрее исполнение.
                            Глянте в сторону и языка РЕФАЛ с его супер компилятором.
                            Только надо делать оптимизированную программ в каком либо смысле — а можно чередовать шаги — сперва делаем меньше команд — потом — быстрее. Чередуем несколько раз — и вот алгоритм для обфукации кода.
                            Может старап забацать чтобы пару баксов срубить?
                            Сделать habrahabr.ru/post/270443 Коммерческая тулза для обработки приватных данных на публичных сетях — например в облаке MS Azure?
                              0
                              Легко рассуждать, сложно сделать. Да и обфускаторов сейчас достаточно, причем бесплатных.
                                0
                                буду пытаться для консольного приложения.
                                пока буду собирать обфукаторы и оптимизаторы разные
                                считаю что надо реализовывать для платформы .net
                                шаг 1. слить все dll в один выполняемый утилитой ILMerge
                                шаг 2. применить обфукаторы разные
                                шаг 3. передать данные и программу в Azure на исполнение
                                так что знакомые с подводными камнями нужны советы их
                          0
                          например, алгоритм можно записать как последовательное применение операндов — соответственно интересно находить коммутативные элементы, чтобы переупорядочить эту последовательность операндов — чтобы по содержимому данных и кода программы невозможно было бы понять какие исходные данные поступили на обработку и какие данные были возвращены
                            +1
                            Стоит отметить, что все-таки точности стандартных математических функций и типа double не хватает для нормального распознавания рациональных и вещественных чисел, однако теоретически все работает.

                            Можно использовать double-double арифметику, я при реализации вот этим документом вдохновлялся.
                              +1
                              Для преобразования чисел в виде рациональной дроби цепные дроби не рассматривали? Есть и готовые реализации.
                                0
                                Не рассматривал, но при желании можно реализовать. В репозитории, кстати, сменил GOLD парсер на ANTLR. Но остаются другие рефакторинги, чтобы шаги преобразования были как можно более иммутабельными.

                                Возможно позже опубликую обновленную статью только уже на английском хабре.

                              Only users with full accounts can post comments. Log in, please.