Вычисление длины 4D-вектора

Точки 4D характеризуются 4-мя координатами. Например, первая точка имеет координаты (x1, y1, z1, w1), а вторая - (x2, y2, z2, w2). В этом случае расстояние между двумя этими точками находится по формуле:

d = sqrt( (y2-y1)2 + (x2-x1)2 + (z2-z1)2+ (w2-w1)2)

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

// Функция вычисления расстояния между двумя точками 4D 
// Для определения координат используется тип single (single precision floating-point)
// Для вычисления применяется сопроцессор NEON
//
// Входные параметры:
// X0 - указатель на 8 чисел с плавающей точкой
// в виде (x1, x2, x3, x4), (y1, y2, y3, y4)
// Результат:
// W0 - вычисленное расстояние (также применяется тип single)

.global distance
distance:
    LDP Q2, Q3, [X0]                // загружаем 8 32-разрядных чисел, по 4 в каждый регистр
    FSUB V1.4S, V2.4S, V3.4S        // V1 = V2 - V3
    FMUL V1.4S, V1.4S, V1.4S        // V1 = V1 * V1 = (xi-yi)^2
    FADDP V0.4S, V1.4S, V1.4S       // S0 = S0 + S1 + S2 + S3
    FADDP V0.4S, V0.4S, V0.4S
    FSQRT S4, S0                    // sqrt(S0)
    FMOV W0, S4                     // результат помещаем в W0
    RET

На входе в функцию distance мы ожидаем, что через регистр X0 будет передаваться адрес 4 точек в виде (x1, x2, x3, x4) и (y1, y2, y3, y4). Условно числа объединены в 2 вектора. Каждый вектор состоит из 4-х 32-разрядных чисел с плавающей точкой, поэтому один вектор можно поместить в один 128-разрядный регистр V, который будет обрабатывать их одновременно как 4 дорожки. И вначале мы загружаем один вектор из 4 точек в регистр V2 и второй вектор в регистр V3.

LDP Q2, Q3, [X0]         // загружаем 4 числа

С помощью инструкции FSUB производим вычитание всех четырех компонентов одновременно

FSUB V1.4S, V2.4S, V3.4S        // V1 = V2 - V3

Далее также одновременно для всех четырех чисел производится возведение в квадрат с помощью инструкции FMUL

FMUL V1.4S, V1.4S, V1.4S     // V1 = V1 * V1 = (xi-yi)^2

В итоге в регистре V1 в четырех дорожках будет по одному 32-разрядному числу, каждый из которых преставляет квадрат разности. Теперь нам надо сложить все эти числа. Разом все четыре числа мы сложим не можем, но можем это сделать по этапно. Сложение всех вычисленных квадратов производится в регистре V0:

FADDP V0.4S, V1.4S, V1.4S

здесь сначала складываем каждую первую пару и вторую пару 32-разрядных чисел с плавающей точкой. Опять же их две операции сложения будут производиться параллельно. То есть в итоге получим два 32-разрядных чисел с плавающей точкой, которые будут занимать две дорожки регистра V0.

С помощью второй инструкции сложения складываем выше полученных два числа в одно

FADDP V0.4S, V0.4S, V0.4S

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

FSQRT S4, S0     // sqrt(S0)

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

Пусть основной файл программы называется main.s. Используем в нем выше определенную функцию distance:

.global main
main:
    STP X19, X20, [SP, #-16]!
    STR LR, [SP, #-16]!
    LDR X20, =points        // загружаем указатель на точки
    MOV W19, #3             // для итерации по 3 пар точкек
loop: 
    MOV X0, X20             // указатель на первое число
    BL distance             // вызываем функцию distance
    FMOV S2, W0             // для преобразования single передаем в fpu
    FCVT D0, S2             // преобразуем single в double
    FMOV X1, D0             // передаем double для форматированного вывода в printf
    LDR X0, =prtstr         // строка форматирования для вывода на консоль
    BL printf               // выводим на консоль расстояние
    ADD X20, X20, #(8*4)    // 2 точки - 8 чисел, каждое по 4 байта
    SUBS W19, W19, #1       // уменьшаем счетчик цикла
    B.NE loop               // повторяем цикл и переходим к следующей паре точек
    
    MOV X0, #0 // код возврата
    LDR LR, [SP], #16
    LDP X19, X20, [SP], #16
    RET
.data
    points: 
        .single 0.0, 0.0, 0.0, 0.0, 17.0, 4.0, 2.0, 1.0
        .single 1.3, 5.4, 3.1, -1.5, -2.4, 0.323, 3.4, -0.232
        .single 1.323e10, -1.2e-4, 34.55, 5454.234, 10.9, -3.6, 4.2, 1.3
    prtstr: .asciz "Distance = %f\n"

Для тестирования в секции .data определяется 6 точек под меткой points - каждая точка имеет 4 компоненты. И в программе проходим в цикле по всем этим точкам и передает поочередно каждую пару точек в функцию distance и вычисляем расстояние между ними. То есть мы ожидаем, что функция вычислит 3 расстояния.

Скомпилируем приложение командой

aarch64-none-linux-gnu-gcc main.s distance.s -o main -static

и запустим на выполнение, и мы должны получить следующий вывод:

Distance = 17.606817
Distance = 6.415898
Distance = 13230000128.000000