Пишем «калькулятор». Часть II. Решаем уравнения, рендерим в LaTeX, ускоряем функции до сверхсветовой

    Привет!

    Итак, в первой части мы уже неплохо поработали, сделав производную, упрощение, и компиляцию. Так, а что еще должен уметь наш простенький калькулятор? Ну хотя бы уравнения вида

    $(x - b)(\tan(\sin(x))^2 - 3\tan(\sin(x)) + c) = 0$

    точно должен решать. А еще красиво нарисовать это дело в латехе, и будет прямо хорошо! Погнали!





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


    Ускорение компиляции


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

    $\sin(x^2) + \cos(x^2) + x^2 + \sin(x^2)$


    На окончание первой части скорость этой функции была такой:
    Метод Время (в наносекундах)
    Substitute 6800
    Наша скомпилированная функция 650
    Эта жа функция, написанная прямо в коде 430

    Что?
    Substitute — классический способ подсчета значения выражения, то есть например
    var x = MathS.Var("x");
    var expr = x + 3 * x;
    Console.WriteLine(expr.Substitute(x, 5).Eval());
    >>> 20
    

    Наша скомпилированная функция — это когда мы делаем так же, но предварительно скомпилировав ее
    var x = MathS.Var("x");
    var expr = x + 3 * x;
    var func = expr.Compile(x);
    Console.WriteLine(func.Substitute(5));
    

    Функция, написанная прям в коде, это когда мы делаем
    static Complex MyFunc(Complex x)
                => x + 3 * x;
    


    Как мы могли заметить, в этой функции есть повторяющиеся части, например, $x^2$, и неплохо было бы их закешировать.

    Для этого введем еще две инструкции PULLCACHE и TOCACHE. Первая будет класть в стек число, лежащее в кеше по адресу, который мы ей передаем. Вторая будет копировать(stack.Peek()) последнее число из стека в кеш, тоже по определенному адресу.

    Остается сделать таблицу, в которую при компиляции мы будем записывать функции для кеширования. Что мы не будем кешировать? Ну во-первых то, что встречается один раз. Лишняя инструкция обращения к кешу — нехорошо. Во-вторых, слишком простые операции тоже смысла кешировать не имеет. Ну например обращение к переменной или числу.

    При интерпретации списка инструкций у нас будет заранее созданный массив для кеширования. Теперь инструкции для этой функции выглядят как

    PUSHCONST (2, 0)
    PUSHVAR 0
    CALL powf
    TOCACHE 0       #здесь мы посчитали квадрат, и так как он нам нужен более одного раза, закешируем-ка мы его.
    CALL sinf
    TOCACHE 1         #синус от квадрата нам тоже понадобится
    PULLCACHE 0     # Оп, встретили использование квадрата. Тащим из кеша
    PULLCACHE 0
    CALL cosf
    PULLCACHE 1
    CALL sumf
    CALL sumf
    CALL sumf
    


    После этого мы получаем уже явно более хороший результат:
    Метод Время (в наносекундах)
    Subtitute 6800
    Наша скомпилированная функция 330 (было 650)
    Эта жа функция, написанная прямо в коде 430

    В репе и компиляция, и интерпретация инструкций реализованы тут.

    LaTeX


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

    Имея дерево выражений сделать рендеринг в латех очень просто. Как сделать это с точки зрения логики? Итак, у нас есть вершина дерева. Если это число или переменная, тут все просто. Если эта вершина, например, оператор деления, нам больше хочется $\frac{a}{b}$, чем $a/b$, поэтому для деления мы пишем что-то типа

    public static class Div
    {
        internal static string Latex(List<Entity> args)
        => @"\frac{" + args[0].Latexise() + "}{" + args[1].Latexise() + "}";
    }
    


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

    $\left(\left(\left(\left(x + 3\right) + \left(a + b\right)\right)\right) + c\right)$


    Однако внимательный человек сразу (а не как я) заметит, что если их убрать совсем, то при построения выражения вида $c * (a + b)$, у нас будет печататься

    $c * a + b$



    Решается просто, вводим приоритеты операторов а-ля
    args[0].Latexise(args[0].Priority < Const.PRIOR_MUL) + "*" + args[1].Latexise(args[1].Priority < Const.PRIOR_MUL);
    

    И вуаля, вы прекрасны.

    Латексиризация сделана тут. И слова latexise не существует, я его сам придумал, не надо пинать меня за это.

    Решение уравнений


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

    Численное решение. Метод Ньютона


    Он предельно простой, зная функцию $f(x)$ мы будем искать корень по итеративной формуле

    $x_{n+1} = x_n - \frac{f(x)}{f(x)'}$


    Так как мы все решаем в комплексной плоскости, можно в принципе написать двумерный цикл, который будет искать решения, а затем вернуть уникальные. При этом производную от функции мы теперь можем найти аналитически, а затем обе функции $f(x)$ и $f(x)'$ скомпилировать.

    Ньютон сидит тут.

    Аналитическое решение


    Первые мысли довольно очевидны. Так, выражение, корни уравнения $a(x) * b(x)$ равны совокупности корней $a(x)$ и $b(x)$, аналогично для деления:

    internal static void Solve(Entity expr, VariableEntity x, EntitySet dst)
    {
        if (expr is OperatorEntity)
        {
            switch (expr.Name)
            {
                case "mul":
                    Solve(expr.Children[0], x, dst);
                    Solve(expr.Children[1], x, dst);
                    return;
                case "div":
                    Solve(expr.Children[0], x, dst);
                    return;
            }
        }
    ...
    


    Для синуса это будет выглядить чуть по-другому:
    case "sinf":
        Solve(expr.Children[0] + "n" * MathS.pi, x, dst);
        return;
    

    Ведь мы хотим найти все корни, а не только те, что равны 0.

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

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

    $\sin(x)^2 + \sin(x) - 0.4 = 0$


    не существует никакого шаблона, а вот к паре

    $ \begin{cases} t^2 + t - 0.4 = 0\\ t=\sin(x) \end{cases} $


    еще как существует.

    Поэтому мы делаем функцию что-то типа GetMinimumSubtree, которая вернула бы самую эффективную замену переменной. Что такое эффективная замена? Эта такая замена, в которой мы
    1. Максимизируем количество употреблений этой замены
    2. Максимизируем глубину дерева (чтобы в уравнении $\sin(\sin(x))^2 + 3$ у нас была замена $\sin(\sin(x))$)
    3. Убеждаемся, что этой заменой мы заменили все упоминания переменной, например если в уравнении $\sin(x)^2 + \sin(x) + x$ мы сделаем замену $t=\sin(x)$, то все станет плохо. Поэтому в этом уравнении лучшая замена по $x$ это $x$ (то есть нет хорошей замены), а вот например в $\sin(x)^2 + \sin(x) + c$ мы можем смело делать замену $t=\sin(x)$.

    После замены уравнение выглядит много проще.

    Полином


    Итак, первое, что мы делаем, мы делаем expr.Expand() — раскроем все скобочки, чтобы избавиться от гадости вида $x(x + 2)$, превратив в $x^2 + 2x$. Теперь, после раскрытия у нас получится что-то типа

    $c + x^3 + 3x^2 - a * x^3 + x$


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

    Что у нас есть о одночлене? Одночлен — это произведение множителей, один из которых — переменная, либо опетатор степени от переменной и целого числа. Поэтому мы заведем две переменных, в одной будет степень, в другой — коэффициент. Далее просто идем по множителям, и каждый раз убеждаемся, что либо там $x$ в целой степени, либо без степени вовсе. Если встречаем бяку — возвращаемся с null-ом.
    бяка
    Если мы вдруг встретили какой-нибудь $x^{3.2}$, $log(x, 2)$ и прочие штуки, в которых упоминается $x$, но не подходит под паттерн одночлена — это плохо.


    Окей, мы собрали словарик, в котором ключ — это степень (целое число), а значение — коэффициент одночлена. Так выглядит он для предыдущего примера:
    0 => c
    1 => 1
    2 => 3
    3 => 1 - a
    


    Вот, к примеру, реализация решения квадратного уравнения
    if (powers[powers.Count - 1] == 2)
    {
        var a = GetMonomialByPower(2);
        var b = GetMonomialByPower(1);
        var c = GetMonomialByPower(0);
        var D = MathS.Sqr(b) - 4 * a * c;
        res.Add((-b - MathS.Sqrt(D)) / (2 * a));
        res.Add((-b + MathS.Sqrt(D)) / (2 * a));
        return res;
    }
    

    Этот момент еще не доделан (например, нужно правильно сделать для кубических полиномов), но вообще происходит это тут.

    Обратная развертка замены


    Итак, вот мы сделали какую-нибудь замену типа $t=\arcsin(x^3 + a^2)$, и теперь нам хочется получить оттуда x. Здесь у нас просто пошаговое развертывание функций и операторов, типа

    $t = \arcsin(x^3 + a^2)$


    $\sin(t) = x^3 + a^2$


    $\sin(t) - a^2 = x^3$


    $(\sin(t) - a^2)^{\frac{1}{3}} = x$



    Отрывок кода выглядит примерно так:
    switch (func.Name)
    {
        case "sinf":
            // sin(x) = value => x = arcsin(value)
            return FindInvertExpression(a, MathS.Arcsin(value), x);
        case "cosf":
            // cos(x) = value => x = arccos(value)
            return FindInvertExpression(a, MathS.Arccos(value), x);
        case "tanf":
            // tan(x) = value => x = arctan(value)
            return FindInvertExpression(a, MathS.Arctan(value), x);
    ...
    

    Код этих функций здесь.

    Все, конечное решение уравнений (пока!) выглядит примерно так
    1. Если знаем ноль функции или оператора, отправляем Solve туда (например, если $a * b = 0$, то запустим Solve(a) и Solve(b))
    2. Найти замену
    3. Решить как полином, если это возможно
    4. Для всех решений развернуть замену, получив конечное решение
    5. Если как полином оно не решилось и у нас только одна переменная, решаем методом Ньютона


    Итак, за сегодня мы:
    • Ускорили работу скомпилированной функции благодаря кешу
    • Отрендерили в LaTeX
    • Решили небольшую часть случаев уравнений


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

    В следующей части...
    Если еще хоть кому-то интересно, попробуем усложнить решение уравнений, добавив кубическую и четвертую степени, а так же более умную сворачивалку. Может быть будем подбирать паттерны, ведь любой школьник решит

    $1 - \sin(x)^2 + \cos(x) + 0.5$


    Но только не наш алгоритм. Более того, может быть будет неопределенный интеграл (начало). Расширение числа в кватернионы или матрицы пока смысла не имеет, но как-нибудь тоже можно будет. Если есть конкретная идея для части III, напишите в личные сообщения или комментарии.


    О проекте
    Те, кто посмотрели код, могли заметить, насколько же он сырой и не отрефакторенный. И, разумеется, я буду рефакторить, поэтому основная идея этой статьи — донести информацию теоретически, и уж потом подглядывая в код. А еще приветствуются pull-реквесты, хотя наверное до wolfram.alpha нам все равно не добраться. Проект открытый.


    А вот парочка примеров того, что мы сегодня сделали
    var x = MathS.Var("x");
    var y = MathS.Var("y");
    var expr = x.Pow(y) + MathS.Sqrt(x + y / 4) * (6 / x);
    Console.WriteLine(expr.Latexise());
    >>> {x}^{y}+\sqrt{x+\frac{y}{4}}*\frac{6}{x}
    

    ${x}^{y}+\sqrt{x+\frac{y}{4}}*\frac{6}{x}$

    var x = MathS.Var("x");
    var expr = MathS.Sin(x) + MathS.Sqrt(x) / (MathS.Sqrt(x) + MathS.Cos(x)) + MathS.Pow(x, 3);
    var func = expr.Compile(x);
    Console.WriteLine(func.Substitute(3));
    >>> 29.4752368584034
    


    Entity expr = "(sin(x)2 - sin(x) + a)(b - x)((-3) * x + 2 + 3 * x ^ 2 + (x + (-3)) * x ^ 3)";
    foreach (var root in expr.Solve("x"))
        Console.WriteLine(root);
    >>> arcsin((1 - sqrt(1 + (-4) * a)) / 2)
    >>> arcsin((1 + sqrt(1 + (-4) * a)) / 2)
    >>> b
    >>> 1
    >>> 2
    >>> i
    >>> -i
    


    Спасибо за внимание!

    Только зарегистрированные пользователи могут участвовать в опросе. Войдите, пожалуйста.

    Делать следующую часть?

    • 83,6%Да (с парсингом)56
    • 7,5%Да (без парсинга)5
    • 9,0%Нет6
    Реклама
    AdBlock похитил этот баннер, но баннеры не зубы — отрастут

    Подробнее

    Комментарии 10

      +4

      В Латехе.

        0

        Латех, латекс… да какая разница?)

          +3

          Как писал Кнут в TeXBook, "TEX (pronounced tecks) is the admirable Text EXecutive processor developed by Honeywell Information Systems", т.е. совсем другая штука.


          Разумеется, сейчас мало кто перепутает, но всё-таки стоит учитывать рекомендации создателя TeX, особенно когда этому посвящена вся первая глава.


          Небольшой отрывок

          English words like ‘technology’ stem from a Greek root beginning with the letters τεχ...; and this same Greek word means art as well as technology. Hence the name TeX, which is an uppercase form of τεχ.


          Insiders pronounce the χ of TeX as a Greek chi, not as an ‘x’, so that TeX rhymes with the word blecchhh. It’s the ‘ch’ sound in Scottish words like loch or German words like ach; it’s a Spanish ‘j’ and a Russian ‘kh’. When you say it correctly to your computer, the terminal may become slightly moist.

            +1

            Окей, поправил. Хотя все-таки странное докапывание, если честно

              +1

              Это не докапывание. Просто указали на небольшую ошибку, а после объяснили, что разница всё-таки есть.


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

        0
        Ещё б добавили как графики строить из c#, включая 3д
          +1

          Не знаю, поможет ли автору мой комментарий, но поскольку сам был вдохновлен подобной идеей когда-то (до практики руки не дошли), накопились некоторые соображения по теме. Я имел неплохой опыт работы с Wolfram Language и большинство замечаний будут так или иначе связаны с ним. Также спешу заметить, что в репозиторий с кодом я так и не заглянул, потому если что-то уже именно так и реализовано, но в статье не отражено (или я не понял) — оно и к лучшему (надеюсь).


          Во-первых. Pattern matching — прекрасная концепция функциональных языков. И она же по сути является краеугольным камнем любых символьных манипуляторов. (Может быть дело пойдет проще в F#?) В примере из предыдущей статьи с "уплощением" дерева для операции + (т.е. перевод x + (y + ...) в x + y + ...) автор вынужден был использовать слишком общие паттерны вроде any, const, Num(1) и т.д. Т.е. никак нельзя взять и сопоставить целое поддерево. Было бы куда проще, если бы можно было задавать правила не вида


          if expr == plus(const1, any1) && any1 == plus(const2, any2)

          а прямо-таки


          if expr == plus(any1, plus(any2, any3))

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


          MatchQ[plus[1, plus[2, plus[3, 4]]], plus[x_, plus[y_, z_]]] == True

          Несложно и заменить все вхождения данного паттерна в дереве:


          plus[1, plus[2, plus[3, 4]]] //. plus[x___, plus[y__], z___] :> plus[x, y, z]
          (* prints plus[1, 2, 3, 4] *)

          Плюс, поскольку такие общие (x___ в частности) паттерны все-таки довольно накладны, математика из коробки предлагает т.н. атрибуты. Например, атрибут Flat, будучи присвоен функции plus, лишает нас необходимости самим ее уплощать.


          Если бы можно было писать произвольные паттерны и доставать из них все упоминающиеся переменные (plus(a_, b_) => {a = ..., b = ...}), все бы существенно упростилось.


          Во-вторых. Что касается упрощений, тут дело гораздо сложнее. Расширенный pattern matching, конечно, поможет, но шаблоны могут быть довольно сложные. Опять же, математика решает этот вопрос в Simplify и FullSimplify путем перебора всех возможных комбинаций подвыражений данного выражения: к каждой комбинации применяется та или иная упрощающая тактика, оценивается качество упрощения. Про качество тоже отдельный разговор. Например, можно минимизировать количество листов дерева выражения. Или длину максимальной ветви дерева. Играя с этим, можно делать интересные вещи.


          В Simplify применяются некоторые эвристики, так что не все комбинации выражений действительно берутся в расчет. Вот пример вывода шагов упрощения выражения x^2 + 2x + a + 1 до a + (1 + x)^2 (вывод FullSimplify будет куда длиннее; к сожалению, я не знаю способа вывести это в формате "до-после", тут только "до"):


          1 + a + 2 x + x^2
          2 x
          x^2
          1 + 2 x + x^2
          (1 + x)^2
          1 + x
          a + (1 + x)^2
          ...

          Код получения такого вывода:


          Simplify[x^2 + 2 x + a + 1, TransformationFunctions -> {Echo, Automatic}]

          Как выражать правили замены? Думаю, также, как и сейчас — паттернами. Создается список элементарных упростителей. Далее перебором комбинаций подвыражений (поддеревьев т.е.) производятся попытки упрощения каждым упростителем. Из-за закрытых исходников математики мой интерес в этой области остался неудовлетворенным.


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


          В-четвертых. Раньше уже упоминали, что можно компилировать выражения прямо в IL-код или в какие-нибудь нативные представления AST (Expressions). Математика тоже так может, плюс она может и в машинный код. Где-то на хабре была статья о компиляции математических функций в IL, но, боюсь, ту, о которой помню, не найду. Вот другая.


          В-пятых. Вообще я был вдохновлен подходом интерактивных proof assistant'ов (Coq, The KeY Project, Z3 и т.д.). И прямо-таки горел желанием что-то свое написать. Но не срослось, поскольку я оценил объем работ неподъемным. Идея была в том, чтобы на каждом из этапов упрощения позволять пользователю выбирать ту тактику упрощения, которая бы ему была более приемлема. Из вычислителя-упростителя система тогда превращается в настоящую систему компьютерного доказательства. Мы выдвигаем гипотезу, что \forall x . sin(x)^2 + cos(x)^2 = 1, доказываем ее как-либо и оформляем в виде теоремы основн_триг_тождество. Далее уже можно на нее ссылаться, чтобы упростить sqrt(1 - sin^2(x)), заменив, например, 1. С группировкой тактик и теорем по категориям и с добавлением изрядной доли эвристик можно даже пользователя привлекать минимально к процессу доказательства. Проблем здесь много — чего только стоит контекст (например, теорема верна, если только x \in \Reals или x \in \Complex), потому я не буду сильно много здесь размышлять.


          В-шестых. Хорошо бы на каком-либо этапе более детально осмотреться. Я имею в виду, например, подсмотреть у Mathematica и у готовых open source продуктов (у Axiom есть прямо-таки теоретические материалы, правда на поверку весьма скудные, где, кроме прочего, есть и аннотированные исходники). К сожалению, по вопросам систем компьютерной алгебры (СКА, CAS) крайне мало подробной информации (во всяком случае мне попадалось).

            0

            if expr == plus(any1, plus(any2, any3))
            У меня нет такового ужаса, и все ваши комментарии в сторону того, что нельзя хардкодить паттерны не очень уместны, так как я их и не хардкодил. Вот их список: https://github.com/Angourisoft/MathS/blob/master/AngouriMath/Functions/Evaluation/Patterns/Patterns.cs
            И вообще, на данный момент simplify более менее нормально работает, ваш пример проработает в том числе.
            Простенький язык я придумывать не буду так как смысла в ран-тайме добавлять функционал не вижу. Хотя кто знает, может и понадобится когда-нибудь. Но тогда это нужно очень узкой аудитории.


            Насчет компиляции в IL — ок, сделаем.


            Пятое: да, действительно умная упрощалка, которая будет работать не за экспоненту — это сложно. Хотя может и в эту сторону будем развиваться.


            Чекну Axiom. У mathematica я подсмотреть не смогу по очевидным причинам

              0

              Этот "ужас" — просто псевдокод, чуть-чуть приближенный к математике по синтаксису. Демонстрация подхода. Про хардкод скорее относится к фрагментам кода конкретно из этой части статьи, а именно к аналитическим солверам. Вы сами демонстрируете фрагменты с switch (func.Name) { case "sinf" ... }, собственно, вот. Здесь есть над чем поработать.


              У Mathematica прекрасная документация. Так что многие вещи можно вынести оттуда. Есть достаточно много открытых проектов и на ней. В т.ч. всякие пакеты символьных вычислений над тензорами — вот где нужен хороший pattern matching! Я в свое время тоже с этим игрался в математике (маленький самописный тензорный движок) и был почти доволен.

            0

            НЛО прилетело и опубликовало эту надпись здесь

            Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

            Самое читаемое