Как писать на ассемблере в 2021 году

    Несмотря на наличие множества языков различной степени высокоуровневости, сегодня ассемблер не потерял своей актуальности и в индексе TIOBE находится на почётном 10-ом месте (на февраль 2021), обогнав такие модные языки как Go и Rust. Одна из причин его привлекательности – в простоте и максимальной близости к железу; с другой стороны, программирование на ассемблере всё ещё может рассматриваться как искусство и даёт совершенно особые эмоции.



    Однако писать программы целиком и полностью на ассемблере — не просто долго, муторно и сложно — но ещё и несколько глупо — ведь высокоуровневые абстракции для того и были придуманы, чтобы сократить время разработки и упростить процесс программирования. Поэтому чаще всего на ассемблере пишут отдельно взятые хорошо оптимизированные функции, которые затем вызываются из языков более высокого уровня, таких как с++ и c#.

    Исходя из этого, наиболее удобной средой для программирования будет Visual Studio, в состав которой уже входит MASM. Подключить его к проекту на с/c++ можно через контекстное меню проекта Build Dependencies – Build Customizations…, поставив галочку напротив masm, а сами программы на ассемблере будут располагаться в файлах с расширением .asm (в свойствах которого Item Type должно иметь значение Microsoft Macro Assembler). Это позволит не просто компилировать и вызывать программы на ассемблере без лишних телодвижений – но и осуществлять сквозную отладку, «проваливаясь» в ассемблерный исходник непосредственно из c++ или c# (в том числе и по точке останова внутри ассемблерного листинга), а также отслеживать состояния регистров наряду с обычными переменными в окне Watch.

    Подсветка синтаксиса


    B Visual Studio нет встроенной подсветки синтаксиса для ассемблера и прочих достижений современного IDE-строения; но её можно обеспечить с помощью сторонних расширений.

    AsmHighlighter — исторически первый с минимальным функционалом и неполным набором команд — отсутсвуют не только AVX, но и некоторые из стандартных, в частности fsqrt. Именно этот факт побудил к написанию собственного расширения —

    ASM Advanced Editor. В нём, помимо подсветки и сворачивания участков кода (с использованием комментариев ";[", ";[+" и ";]") реализована привязка подсказок к регистрам, всплывающих по наведению курсора ниже по коду (также через комментарии). Выглядит это так:

    ;rdx=ссылка на текущий элемент входного массива

    или так:

    mov rcx, 8;=счётчик цикла

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

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

    Asm Dude — обнаружился чуть позже. В нём автор пошёл другим путём и сфокусировал силы на встроенном справочнике команд и автодополнении, в том числе и с отслеживанием меток. Сворачивание кода там также присутствует (по «#region / #end region»), но привязки комментариев к регистрам вроде ещё нет.

    32 vs. 64


    С тех пор, как появилась 64-битная платформа, стало нормой писать по 2 варианта приложений. Пора завязывать с этим! Сколько можно тянуть легаси. Это же касается и расширений — найти процессор без SSE2 можно разве что в музее – к тому же, без SSE2 64-битные приложения и не заработают. Никакого удовольствия от программирования не будет, если писать по 4 варианта оптимизированных функций для каждой платформы. Только 64 бит/AVX, только хардкор! Хотя может быть и прямо противоположный взгляд — новые процессоры и так работают быстро, и оптимизацию стоит делать под старые. В общем, всё зависит от конкретной задачи.

    Преимущество 64-битной платформе вовсе не в «широких» регистрах – а в том, что этих самых регистров стало в 2 раза больше – по 16 штук как общего назначения, так и XMM/YMM. Это не только упрощает программирование, но и позволяет значительно сократить обращения к памяти.

    FPU


    Если ранее без FPU было никуда, т.к. функции с вещественными числами оставляли результат на вершине стека, то на 64-битной платформе обмен проходит уже без его участия с использованием регистров xmm расширения SSE2. Intel в своих руководствах также активно рекомендует отказаться от FPU в пользу SSE2. Однако есть нюанс: FPU позволяет производить вычисления с 80-битной точностью — что в некоторых случаях может оказаться критически важным. Поэтому поддержка FPU никуда не делаcь, и рассматривать её как устаревшую технологию совершенно точно не стоит. Например, вычисление гипотенузы можно делать «в лоб» без опасения переполнения,

    а именно
    fld x
    fmul st(0), st(0)
    fld y
    fmul st(0), st(0)
    faddp st(1), st(0)
    fsqrt
    fstp hypot

    Основная сложность при программировании FPU — это его стековая организация. Для упрощения была написана небольшая утилитка, автоматически генерирующей комментарии с текущим состоянием стека (планировалось добавить подобную функциональность непосредственно в основное расширение для подсветки синтаксиса — но до этого руки так и не дошли)

    Пример оптимизации: преобразование Хартли


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

    Итак, в качестве задачи возьмём алгоритм (медленного) преобразования Хартли:

    код
    static void ht_csharp(double[] data, double[] result)
    {
        int n = data.Length;
        double phi = 2.0 * Math.PI / n;
        for (int i = 0; i < n; ++i)
        {
            double sum = 0.0;
            for (int j = 0; j < n; ++j)
            {
                double w = phi * i * j;
                sum += data[j] * (Math.Cos(w) + Math.Sin(w));
            }
            result[i] = sum / Math.Sqrt(n);
        }
    }

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

    код (комментарии удалены)
    ht_asm PROC
    local sqrtn:REAL10
    local _2pin:REAL10
    local k:DWORD
    local n:DWORD
    
    and r8, 0ffffffffh
    mov n, r8d
    mov r11, rcx
    xor rcx, rcx
    
    mov r9, r8
    dec r9
    shr r9, 1
    
    mov r10, r8
    sub r10, 2
    shl r10, 3
    
    finit
    fld _2pi
    fild n
    fdivp st(1), st
    fstp _2pin
    
    fld1
    fild n
    fsqrt
    fdivp st(1), st
    ;
    mov rax, r11
    mov rcx, r8
    
    fldz
    @loop0:
    	fadd QWORD PTR [rax]
    	add rax, 8
    	loop @loop0
    fmul st, st(1)
    fstp QWORD PTR [rdx]
    fstp sqrtn
    add rdx, 8
    
    mov k, 1
    @loop1:
    	mov rax, r11
    	fld QWORD PTR [rax]
    	fld st(0)
    	add rax, 8
    	fld _2pin
    
    	fild k
    	fmulp st(1),st
    	fsincos
    
    	fld1;=u
    	fldz;=v
    
    	mov rcx, r8
    	dec rcx
    	@loop2:
    		fld st(1)
    		fmul st(0),st(4)
    		fld st(1)
    		fmul st,st(4)
    		faddp st(1),st
    		fxch st(1)
    		fmul st, st(4)
    		fxch st(2)
    		fmul st,st(3)
    		fsubrp st(2),st
    		
    		fld st(0)
    		fadd st, st(2)
    		fmul QWORD PTR [rax]
    		faddp st(5), st
    		
    		fld st(0)
    		fsubr st, st(2)
    		fmul QWORD PTR [rax]
    		faddp st(6), st
    		
    		add rax, 8
    		loop @loop2
    	
    	fcompp
    	fcompp
    	fld sqrtn
    	fmul st(1), st
    	fxch st(1)
    	fstp QWORD PTR [rdx]
    	fmulp st(1), st
    	fstp QWORD PTR [rdx+r10]
    	
    	add rdx,8
    	sub r10, 16
    	inc k
    	dec r9
    	jnz @loop1
    
    test r10, r10
    jnz @exit
    
    mov rax, r11
    fldz
    mov rcx, r8
    shr rcx, 1
    
    @loop3:;[
    	fadd QWORD PTR [rax]
    	fsub QWORD PTR [rax+8]
    	add rax, 16
    	loop @loop3;]
    fld sqrtn
    fmulp st(1), st
    fstp QWORD PTR [rdx]
    
    @exit:
    ret
    ht_asm ENDP

    Обратите внимание: тут нет ни разворачивания цикла, ни SSE/AVX, ни косинусных таблиц, ни снижения сложности за счёт «быстрого» алгоритма преобразования. Единственная явная оптимизация — это итеративное вычисление синуса/косинуса во внутреннем цикле алгоритма непосредственно в регистрах FPU.

    Поскольку речь идёт об интегральном преобразовании, то помимо скорости, нас ещё интересует точность вычисления и уровень накопленных погрешностей. В данном случае посчитать её очень просто — делая два преобразования подряд, мы должны получить (в теории) исходные данные. На практике они будут слегка отличаться, и можно будет посчитать ошибку через среднеквадратическое отклонения полученного результата от аналитического.

    Результаты авто-оптимизации программы на c++ также могут сильно зависеть от настроек параметров компилятора и выбора допустимого расширенного набора инструкций (SSE/AVX/etc). При этом есть два нюанса:

    1. Современные компиляторы склонны всё возможное вычислять на этапе компиляции – поэтому вполне возможно в откомпилированном коде вместо алгоритма увидеть заранее посчитанное значение, что при замере производительности даст преимущество компилятору в 100500 раз. Во избежание этого в моих замерах используется внешняя функция zero(), добавляющая неопределённости входным параметрам.
    2. Если компилятору указать «не использовать AVX» — это ещё не значит, что в полученном коде будет отсутствовать AVX. Внешние библиотеки сами проверяют доступный набор команд на текущей платформе и выбирают соответствующую реализацию. Поэтому единственно надёжный способ сравнения производительности в таком случае – испытывать код на платформе, где AVX отсутствует в принципе.

    Самый интересный параметр оптимизации – это Floating Point Model, принимающий значения Precise|Strict|Fast. В случае с Fast компилятору разрешается делать любые математические преобразования на своё усмотрение (в том числе и итеративные вычисления) – собственно, только в этом режиме и происходит автоматическая векторизация.

    Итак, компилятор Visual Studio 2019, целевая платформа AVX2, Floating Point Model=Precise. Чтобы было ещё интереснее — будет измерять из проекта на c# на массиве из 10000 элементов:



    C# ожидаемо оказался медленнее с++, а функция на ассемблере оказалась быстрее в 9 раз! Однако ещё рано радоваться — установим Floating Point Model=Fast:



    Как видно, это помогло значительно ускорить код и отставание от ручной оптимизации составило всего лишь в 1.8 раз. Но вот что не изменилось – так это погрешность. Что тот, что другой вариант дал ошибку в 4 значащих цифры – а это немаловажно при математических вычислениях.

    В данном случае наш вариант оказался и быстрее, и точнее. Но так бывает не всегда – и выбирая FPU для хранения результатов мы неизбежно будем терять в возможности оптимизации векторизацией. Также никто не запрещает комбинировать FPU и SSE2 в тех случаях, когда это имеет смысл (в частности, такой подход я использовал в реализации double-double арифметики, получив 10-кратное ускорение при умножении).

    Дальнейшая оптимизация преобразования Хартли лежит уже в другой плоскости и (для произвольного размера) требует алгоритма Блюстейна, который также критичен к точности промежуточных вычислений. Ну а этот проект можно скачать на GitHub, и в качестве бонуса там также можно найти реализацию функций для суммирования/масштабирования массивов на FPU/SSE2/AVX.

    Что почитать


    Литературы по ассемблеру навалом. Но можно выделить несколько ключевых источников:
    1. Официальная документация от Intel. Ничего лишнего, вероятность опечаток минимальна (кои в печатной литературе встречаются повсеместно).
    2. Официальная документация от Microsoft.
    3. Онлайн справочник, спарсенный из официальной документации.
    4. Сайт Агнера Фога, признанного эксперта по оптимизации. Также содержит образцы оптимизированного кода на C++ с использованием интринсиков.
    5. SIMPLY FPU.
    6. 40 Basic Practices in Assembly Language Programming.
    7. Все, что нужно знать, чтобы начать программировать для 64-разрядных версий Windows.

    Appendix: Почему бы просто не использовать интринсики (Intrinsics)?


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

    Длинный ответ:

    • Не вся оптимизация делается векторизацией. Если во времена DOS оптимизация делалась за счёт экономии тактов и количества обращений к памяти, то сейчас основной инструмент оптимизации – это организация оптимальной работы с памятью во избежание промахов кеша.
    • Интринсики к новым инструкциям появляются с некоторым запаздыванием. Недостающие интринсики нельзя добавить простым включением заголовочного файла – их поддержка должна быть реализована на уровне компилятора.
    • Не на все инструкции возможно сделать интринсики. Когда оптимизация касается точности математических выражений, можно добиться значительных улучшений, программируя непосредственно на математическом сопроцессоре (FPU). Для команд сопроцессора же интринсики отсутствуют. Они также отсутствуют для команд, на результат которых влияют флаги переноса/переполнения/etc.
    • Интринсики не представляют из себя какой-либо высокоуровневый фреймворк – по факту, это не более чем лишь ещё один уровень мнемоник над ассемблерными инструкциями, мало чем отличающийся от ассемблерных вставок. Что автоматически означает, что код вовсе не становится более читаемым, а при его написании так или иначе нужно придерживаться ассемблерного стиля мышления.
    • Также, из C++ кода вовсе не очевидно, что существуют какие-то ограничения на количество одновременно используемых SIMD-регистров. В 32-х битном ассемблере вы не можете использовать XMM8 наравне с XMM0 / XMM7 – код просто не откомпилируется. А интринсиками вы можете определить хоть тысячу одновременно используемых регистров — и компилятор вынужден будет сам решать, что хранить в регистрах, что в памяти, и как между собой их оптимально комбинировать. В таком случае применение интринсиков не даёт никаких преимуществ в плане увеличения производительности – компилятор работает с ними так же, как и с обычными структурами на C++.
    • Программировать на ассемблере не так уж и сложно, особенно если не пытаться делать на нём хитрые высокоуровневые конструкции. Многое можно почерпнуть, изучая ассемблерный листинг простых алгоритмов в режиме отладки и с различными уровнями оптимизации – компилятор от Microsoft производит достаточно лёгкий для понимания и в то же время эффективный код, комментируя его исходным кодом на C++.

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

      +2
      Честно говоря пару раз пробовал переписать сишную функцию на ассемблере. Как я не старался прироста в производительности не получил.
        +6
        Так то Intel, там и компилятор нормальный. А я вот как не напишу функцию на ассемблере для своего VLIW DSP — все раз в 10 быстрее чем на Си.
          0

          Напишите проверку что аргумент ф-ции делится на 10 без остатка на чистом асме и сравните с Си.

          0
          Так язык не увеличивает производительность, он расширяет возможности. Перед переписыванием нужно понимать что ты сможешь сделать лучше компилятора.
            0
            Ну и откуда взять это понимание?
            Я помню времена, когда просто посчитав количество инструкций в функции можно было бы примерно оценить ее производительность. Но сейчас такое не работает. Глядя на количество строк кода невозможно в уме прикинуть быстр этот код или нет.
            А вот компилятор скорее всего все это знает.
              +1
              Очевидно, из архитектуры процессора и эмпирической оценки времени выполнения команд. Как-то смотришь на результат компиляции и находишь команды, которые можно объединить, переставить или убрать исходя из общего смысла кода.
              По крайне мере, у меня получалось векторизация в циклах лучше, чем gcc 4.8.x c -O3.
              А целые программы ассемблере писать можно только с большого горя — когда стека нет, например.
                0
                Вот, например, отреверсенная программа после С-компилятора, архитектура со стеком и всеми делами:
                Извиняюсь за форматирование, у меня не работает кнопка исходный код, а хабр не умеет выравнивать ни табами, ни пробелами
                LCD_WC_DB0
                ;FSR0 to Arg1
                movf FSR2L, W, ACCESS
                addlw 0xFE
                movwf FSR0L, ACCESS
                movlw 0xFF
                addwfc FSR2H, W, ACCESS
                movwf FSR0H, ACCESS
                ;W = Arg1
                movf INDF0, W, ACCESS
                andlw b'00000001'
                bnz LCD_WC_DB0ON
                LCD_WC_DB0OFF
                bcf LCD_DB0 ;DB0 = 0
                bra LCD_WC_DB1
                LCD_WC_DB0ON
                bsf LCD_DB0 ;DB0 = 1
                LCD_WC_DB1
                ;FSR0 to Arg1
                movf FSR2L, W, ACCESS
                addlw 0xFE
                movwf FSR0L, ACCESS
                movlw 0xFF
                addwfc FSR2H, W, ACCESS
                movwf FSR0H, ACCESS
                ;W = Arg1
                movf INDF0, W, ACCESS
                andlw b'00000010'
                bnz LCD_WC_DB1ON
                LCD_WC_DB1OFF
                bcf LCD_DB1 ;DB1 = 0
                bra LCD_WC_DB2
                LCD_WC_DB1ON
                bsf LCD_DB1 ;DB1 = 1

                Штуки КАПСОМ на всю строку без команд — это метки.

                Так вот, задача функцыи простая: взять свой агрумент из стека в виде байта, в соответсвии с каждым битом из байта устарновить лог «0» или «1» на выходе микроконтроллера.
                Имеем после компилятора: используется косвенная адресацыя, установили указатель на аргумент (6 команд, т.к. архитектура 8 бит, а указатели — слово), поместили в аккумулятор, наложыли маску в соответствии с целевым битом (2 команды), в зависимости от значения установили лог «0» или «1» (4 команды). И так 8 раз — для каждого бита, итого 12х8=96 команд (192 байта памяти).
                Как бы это выглядело на ассемблере: установили указатель на агрумент (6 команд, но можно упростить до 3х, если учесть что старшая часть адреса неизменна), поместили аргумент в память (2 команды), потом 8 раз: сдвиг на 1 бит (1 команда) и установка «0»/«1» по результатам (4 команды, но можно и в 3 соптимизировать). Итого 6+8*(1+4) = 46 команд или даже 39 после простой оптимизацыи. Или 92 и 78 байт памяти соответственно.
                В итоге в отреверсенном коде из 11 кБ половина — это всякие прыжки указателями вокруг стека в т.ч проверка в каждой функцыи по завершению, не вышли ли указатели за границы стека и их коррекцыя. И таким образом программу можно уменьшыть по размеру в 2 раза!
                Конечно, есть подозрения, что это сделано бесплатной версией компилятора, не умеющей оптимизацыю, но в итоге этот код в продакшене.
                  0
                  И что из этого следует? Если нужна оптимизация, то надо для начала смотреть код и настройки компилятора.
            0

            Попробуйте посчитать количество вхождений символа в строку. Если использовать интринсики, то становится сильно быстрее, чем в случае очевидного цикла, даже с -O3 -march=native.

              0
              Просто переписать недостаточно. Нужно сделать это лучше компилятора. Нужно иметь хитрый план по оптимизации. Нужно знать то, что не знает компилятор, и учесть это (например, если размер циклического буфера кратен степени двойки — можно использовать and вместо остатка для деления, если выровнен по 32 байта — можно использовать vmovapd вместо vmovupd, если размер массива ограничен сверху небольшим значением — можно развернуть его полностью и использовать динамически сформированные адреса для перехода, и т.д. и т.п.). Нужно решить оптимизационную задачу — как минимизировать количество обращений к памяти и промахи кеша, используя регистры по максимуму. В частности, используя интрисики AVX в с++ коде ускорения в 4 раза обычно не получается только потому, что компилятор один фиг хранит их (ymm регистры) в памяти и полностью все 16 штук сразу не будет использовать никогда.

              Конечно, многие оптимизационные стратегии можно применить непосредственно на си. Но опять же, есть нюанс — для этого нужно уметь мыслить на ассемблере. А чтобы уметь мыслить на ассемблере, нужно уметь программировать на ассемблере и иметь соответствующий опыт.
              0
              Вы сравниваете теплое с мягким — очевидно, что тривиальная программа при прочих равных медленнее/менее точна чем программа с итеративным вычислением тригонометрии. Так что непонятно, на что списывать разницу в результатах — на ассемблер или на алгоритм. Нужно реализовать итеративную тригонометрию на C# и сравнить ассемблер с ней.
              Ну и странно, что fsincos вызывается n раз — можно посчитать функции только для phi, а дальше все по формулам.
                +1
                Сравнение конечно же несправедливое. Компилятор использует новейший набор инструкций AVX2, а я — устаревшую технологию, медленную по определению. Никто же этим фактом не возмутился. Сравниваются же подходы к оптимизации в первую очередь — машинное против человеческого. Итеративное вычисление ко/синусов компилятору вполне доступно (как минимум в теории) — ведь это же по сути просто комплексное умножение.

                fsincos вызывается n раз во внешнем цикле — потому что регистров сопроцессора не хватило бы, чтобы их все там держать. Пришлось бы сохранять их во внешних переменных — а это слегка бы нарушило чистоту эксперимента. А так мы получаем повышение точности вычислений без явного определений переменных расширенной точности.
                  0
                  Вот сравнение с исправленным вариантом от teology (cpp2 и csharp2) также на 10000 элементах:



                  Результат стал получше, но не в разы — причём с си-шарпом практически сравнялись. А вот падение точности — катастрофическое, 3 значащих цифры только осталось.
                    0
                    1. kovserg увидел ошибку в моем коде, надо ввести временную переменную tmp, в которой сохранить предыдущее значение cos_w и cos_dw и использовать tmp при вычислении sin_w, sin_dw.


                    2. Оптимизатор настраиваете при компиляции на максимум -O3?


                      0
                      1. Да, tmp я и добавил. 2. Оптимизация — Maximum Optimization (Favor Speed) (/O2), /O3 в вариантах выбора нету. Есть /Ox, но по результату принципиально не отличается — и то в худшую сторону.
                        0

                        Может long double?)

                          0
                          long double в компиляторе от Microsoft от double ничем не отличается. Вообще, основная задача 80-битной точности — это хранить промежуточные данные и соответственно их использование должно быть проблемой компилятора, а не программиста.
                      0
                      Странные у Вас тайминги. Я не смог запустить Ваш ассемблер, но сварганил сравнение изначального C++ с оптимизированным — так у меня оптимизированный в 9 быстрее. Но точность — да, проседает, причем я даже не знаю почему, вроде бы вычисление тригонометрии должно быть менее точным чем сложения/умножения.
                      Вот код:
                      Тест
                      #include <stdio.h>
                      #include <stdlib.h>
                      #define _USE_MATH_DEFINES
                      #include <math.h>
                      #include <chrono>
                      
                      void ht_cpp0(double* data, double* result, int n)
                      {
                          double phi = 2.0 * M_PI / n;
                          for (int i = 0; i < n; ++i)
                          {
                              double sum = 0.0;
                              for (int j = 0; j < n; ++j)
                              {
                                  double w = phi * i * j;
                                  sum += data[j] * (cos(w) + sin(w));
                              }
                              result[i] = sum / sqrt(n);
                          }
                      }
                      
                      void ht_cpp1(double* data, double* result, int n)
                      {
                          double phi = 2.0 * M_PI / n;
                          double cos_phi = cos(phi);
                          double sin_phi = sin(phi);
                          double sqrtn = 1 / sqrt(n);
                          double sin_dw = 0.0;
                          double cos_dw = 1.0;
                          for (int i = 0; i < n; ++i)
                          {
                              double sin_w = 0.0;
                              double cos_w = 1.0;
                              double sum = 0.0;
                      
                              for (int j = 0; j < n; ++j)
                              {
                                  sum += data[j] * (cos_w + sin_w);
                                  double new_cos_w = cos_w * cos_dw - sin_w * sin_dw;
                                  double new_sin_w = sin_w * cos_dw + cos_w * sin_dw;
                                  cos_w = new_cos_w;
                                  sin_w = new_sin_w;
                              }
                      
                              result[i] = sum * sqrtn;
                              double new_cos_dw = cos_dw * cos_phi - sin_dw * sin_phi;
                              double new_sin_dw = sin_dw * cos_phi + cos_dw * sin_phi;
                              cos_dw = new_cos_dw;
                              sin_dw = new_sin_dw;
                          }
                      }
                      
                      double MeanSquaredError(double* data1, double* data2, int n)
                      {
                          double sum = 0;
                          for (int i = 0; i < n; i++)
                          {
                              double df = data1[i] - data2[i];
                              sum += df * df;
                          }
                          return sqrt(sum);
                      }
                      
                      int main()
                      {
                          enum : int { data_length = 10000 };
                          double data_src[data_length];
                          for (int i = 0; i < data_length; i++)
                              data_src[i] = i + 1;
                          double data_in[data_length];
                          double data_out[data_length];
                      
                          {
                              auto t0 = std::chrono::high_resolution_clock::now();
                              ht_cpp0(data_src, data_in, data_length);
                              ht_cpp0(data_in, data_out, data_length);
                              auto t1 = std::chrono::high_resolution_clock::now();
                              std::chrono::duration<double> ms_double = t1 - t0;
                              double err = MeanSquaredError(data_src, data_out, data_length);
                              printf("cpp0 : %lf\terr : %1.10lf\n", ms_double, err);
                          }
                      
                          {
                              auto t0 = std::chrono::high_resolution_clock::now();
                              ht_cpp1(data_src, data_in, data_length);
                              ht_cpp1(data_in, data_out, data_length);
                              auto t1 = std::chrono::high_resolution_clock::now();
                              std::chrono::duration<double> ms_double = t1 - t0;
                              double err = MeanSquaredError(data_src, data_out, data_length);
                              printf("cpp1 : %lf\terr : %1.10lf\n", ms_double, err);
                          }
                      
                          system("pause");
                      }
                        0
                        Тайминги в секундах. Мерил и из шарпа, и из плюсов — почти совпадают. А где, кстати, ваши? У вас же, я так понимаю, другой компилятор и возможно даже другая платформа.
                          0
                          Windows/VS 2019. Но я понял, в чем было дело — я не включил поддержку AVX2. С ней происходит векторизация внутреннего цикла изначального варианта, вызываются функции "__vdecl_sin4"/"__vdecl_cos4", которые видимо обрабатывают по 4 аргумента за раз, и тогда соотношение времени выполнения двух алгоритмов у меня примерно совпадает с Вашим. Но тут стоит заметить, что для оптимизированного алгоритма компилятор векторизацию не сделал (хотя unroll на 10 итераций сделать смог), если доделать ее — результат должен улучшиться.
                          0
                          Но точность — да, проседает, причем я даже не знаю почему, вроде бы вычисление тригонометрии должно быть менее точным чем сложения/умножения.


                          Все просто объяснимо. Синус-косинус считается каждый раз заново и погрешность там возникает при вычислении 10-20 слагаемых и эта ошибка более-менее постоянна в течение всего алгоритма. А алгоритм с рекурсивной формулой накапливает ошибку при n^2 пересчетах. Проделайте вычисления самой первой итерации большого цикла вручную и поймете.
                      –3
                      Шарп то тут при чем? Давайте бульдозером картофель копать… ;) Кроме того, если уж такая проблема и возникает, ничего не мешает шарпом подгрузить dll.
                        0
                        Тут и так шарпом подгружается dll. И расширения для Visual Studio тоже на шарпе написаны (работают на 2012+).
                        +1
                        Попробуйте:
                        1. Вместо умножений использовать инкремент.
                        2. Вместо вычислений sin/cos углов вычислять sin/cos суммы углов (текущий угол = предыдущий угол + инкремент).
                        3. На всякий случай закэшировать вычисление 1/sqrt(n), вдруг компилятор слабо оптимизирует.

                        static void ht_csharp(double[] data, double[] result)
                        {
                            int n = data.Length;
                            double phi = 2.0 * Math.PI / n;
                            double cos_phi = Math.Cos(phi); // кэшируем
                            double sin_phi = Math.Sin(phi); // кэшируем
                            double sqrtn = 1 / Math.Sqrt(n); // кэшируем
                            double dw = 0.0; // угол приращения
                            double sin_dw = 0.0;
                            double cos_dw = 1.0;
                            for (int i = 0; i < n; ++i)
                            {
                                double w = 0.0;
                                double sin_w = 0.0; // начальный синус
                                double cos_w = 1.0; // начальный косинус
                                double sum = 0.0;
                                for (int j = 0; j < n; ++j)
                                {
                                    sum += data[j] * (cos_w + sin_w);
                                    w += dw;
                                    cos_w = cos_w * cos_dw - sin_w * sin_dw; // пересчитываем косинус по формуле косинуса суммы
                                    sin_w = sin_w * cos_dw + cos_w * sin_dw; // пересчитываем синус по формуле синуса суммы
                                }
                                result[i] = sum * sqrtn;
                                dw += phi;
                                cos_dw = cos_dw * cos_phi - sin_dw * sin_phi; // пересчитываем косинус по формуле косинуса суммы
                                sin_dw = sin_dw * cos_phi + cos_dw * sin_phi; // пересчитываем синус по формуле синуса суммы
                            }
                        }
                        

                        Можете еще выбросить вычисление углов w и dw, но компилятор должен это сам сделать.
                        Итого: в циклах нет ни одного вычисления синуса и косинуса. Ускорение должно составить 40 крат?
                          0
                          В ассемблере так и сделано, но только для внутреннего цикла, во внешнем считаются sin/cos каждую итерацию.
                            +2
                            cos_w = cos_w * cos_dw - sin_w * sin_dw;
                            sin_w = sin_w * cos_dw + cos_w * sin_dw;
                            ...
                            cos_dw = cos_dw * cos_phi - sin_dw * sin_phi;
                            sin_dw = sin_dw * cos_phi + cos_dw * sin_phi;

                            1. Вам понадобиться временная переменная иначе, результат будет немного странный.
                            2. При больших n можно огрести проблемы с численной неустойчивостью.
                            3. Зачем w и dw если не используются

                            ps: БПФ ускорит гораздо лучше особенно при больших n
                              0

                              Надо придумать преобразование без хранения во временной переменной.)

                              0
                              Ваш код даёт некорректный результат, в нём где-то ошибка. Даже если на это закрыть глаза — то он работает всего лишь в 2 раза быстрее (с компилятором c#) наивной реализации.
                                0
                                UPD: тормоза в вашем коде были связаны с тем, что там появлялся NaN. Исправленный вариант работает нормально (по скорости), результаты сравнения выше.
                              0
                              Тогда этот код C# может оказаться быстрее ассемблера.
                              P.S. Я проверил, код эквивалентно работает, ошибок нет.
                                +1
                                По мне лучше интринсиками пользоваться, чем асм вставками. Или весь проект писать на асме, ну или dll-ку на асме. Тем более писать на высокоуровневом ассемблере не слишком сложно.
                                Вот пример.
                                UASM
                                ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
                                ;;	Quick sort		v1.01
                                ;;	http://rosettacode.org/wiki/Sorting_algorithms/
                                ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
                                .386
                                .model flat, stdcall
                                option casemap:none
                                
                                include msvcrt.inc
                                include macros.asm
                                
                                .code
                                align_proc
                                quicksort proc uses esi edi ebx pBegin:ptr sdword, pEnd:ptr sdword
                                	ASSUME	edi:ptr sdword, esi:ptr sdword
                                	mov		edi, pBegin		;указатель на первый элемент массива
                                	mov		esi, pEnd		;указатель на последний элемент массива
                                	.if (edi!=esi)
                                		;pivot = pBegin[(pEnd - pBegin + 1)/2];//ecx
                                		lea		edx, [esi+1*4]
                                		.for (edx-=edi, edx>>=3, ecx=[edi+edx*4]: : edi+=4, esi-=4)
                                			.for (: sdword ptr [edi] < ecx: edi+=4)
                                			.endfor
                                			.for (: sdword ptr [esi] > ecx: esi-=4)
                                			.endfor
                                			.break .if (edi >= esi)
                                			swap	[edi], [esi]
                                		.endfor
                                		quicksort(pBegin, &[edi-1*4])
                                		quicksort(edi, pEnd)
                                	.endif
                                	ASSUME	edi:nothing, esi:nothing
                                	ret
                                quicksort endp
                                
                                createSpisok(spisok, sdword, 4, 65, 2, -31, 0, 99, 2, 83, 782, 1)
                                
                                main proc C argc:sdword, argv:ptr ptr, envp:ptr 
                                	printf("Quick sort.\n")
                                	.for (esi = spisok.pFirst: esi <= spisok.pLast: esi+=4)
                                		printf("%d ", dword ptr [esi])
                                	.endfor
                                	printf("\n")
                                	quicksort(spisok.pFirst, spisok.pLast)
                                	.for (esi = spisok.pFirst: esi <= spisok.pLast: esi+=4)
                                		printf("%d ", dword ptr [esi])
                                	.endfor
                                	printf("\n")
                                	xor		eax, eax
                                	ret
                                main endp
                                
                                main_startup3_END

                                Библиотеки макросы свои.

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

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