В наше время, а дальше больше, знание умножения и сложения становится стратегическим. Вот то самое умножение в столбик, знакомое с начальной школы и сложение с переносом, знакомое еще раньше умножения теперь стратегический ресурс.
И причина простая как 2х2. Весь обмен информацией теперь шифруется и подписывается. А это, как ни крути, умножение и сложение больших чисел
Любая передача информации, денег, команд и знаний нынче требует стандартных шагов: зашифровать-расшифровать и проверить. А это, в основном, именно те самые умножения и сложения больших чисел в большом количестве.
Ну и поскольку самый наш любимый размер это 256 битовые числа, то попробуем кое-что тут посмотреть интересного.
В это статье будем рассматривать умножение и сложение только буквально “в столбик”, именно как в школе учили.
Алгоритм великого А.А.Карацубы рассматривать не будем сейчас.
Начнем. Обычный алгоритм умножения числа на большое число в сложение с другим большим числом. Это обычная и очень популярная операция. Вот простое умножение двух больших чисел по модулю другого большого числа требует несколько таких операций.
В качестве железа для тестов возьмем видеокарту NVIDIA и все вычисления будем делать на С и PTX ассемблере. Просто такие у автора есть, было бы другое железо, то пример и основа были бы другими.
Вот обычная типичная реализация на ассемблере PTX. (это псевдо ассемблер, который еще раз переводит ассемблерный код в реально ассемблерный код). На этом PTX пишут все эти ИИ библиотеки, если хотят что-бы было быстро.
__device__ inline void cu_32_mul_add_to_288_b(uint32_t *res, const uint32_t *a, const uint32_t b) {
uint64_t mul2 = 0;
uint32_t l32 = 0, h32 = 0, buf[9];
asm volatile(
"mul.wide.u32 %9, %12, %20;nt" // %9 = a[0]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32}
"add.cc.u32 %0, 0, %10;nt" // buf[0] += %lo32
"addc.cc.u32 %1, 0, %11;nt" // buf[1] += %hi32
"addc.u32 %2, 0, 0;nt"
"mul.wide.u32 %9, %13, %20;nt" // %mul1 = a[1]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32}
"add.cc.u32 %1, %1, %10;nt" // buf[1] += %lo32
"addc.cc.u32 %2, %2, %11;nt" // buf[2] += %hi32
"addc.u32 %3, 0, 0;nt"
"mul.wide.u32 %9, %14, %20;nt" // %mul1 = a[2]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32}
"addc.cc.u32 %2, %2, %10;nt" // buf[2] += %lo32
"addc.cc.u32 %3, %3, %11;nt" // buf[3] += %hi32
"addc.u32 %4, 0, 0;nt"
"mul.wide.u32 %9, %15, %20;nt" // %mul1 = a[3]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32}
"addc.cc.u32 %3, %3, %10;nt" // buf[3] += %lo32
"addc.cc.u32 %4, %4, %11;nt" // buf[4] += %hi32
"addc.u32 %5, 0, 0;nt"
"mul.wide.u32 %9, %16, %20;nt" // %mul1 = a[4]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32}
"addc.cc.u32 %4, %4, %10;nt" // buf[4] += %lo32
"addc.cc.u32 %5, %5, %11;nt" // buf[5] += %hi32
"addc.u32 %6, 0, 0;nt"
"mul.wide.u32 %9, %17, %20;nt" // %mul1 = a[5]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32}
"addc.cc.u32 %5, %5, %10;nt" // buf[5] += %lo32
"addc.cc.u32 %6, %6, %11;nt" // buf[6] += %hi32
"addc.u32 %7, 0, 0;nt"
"mul.wide.u32 %9, %18, %20;nt" // %mul1 = a[6]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32}
"addc.cc.u32 %6, %6, %10;nt" // buf[6] += %lo32
"addc.cc.u32 %7, %7, %11;nt" // buf[7] += %hi32
"addc.u32 %8, 0, 0;nt" //
"mul.wide.u32 %9, %19, %20;nt" // %mul1 = a[7]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32c}
"addc.cc.u32 %7, %7, %10;nt" // buf[7] += %lo32
"addc.cc.u32 %8, %8, %11;nt" // buf[8] += %hi32
: "+r"(buf[0]), "+r"(buf[1]), "+r"(buf[2]), "+r"(buf[3]),
"+r"(buf[4]), "+r"(buf[5]), "+r"(buf[6]), "+r"(buf[7]), "+r"(buf[8]),
"=l"(mul2), "+r"(l32), "+r"(h32)
: "r"(a[0]), "r"(a[1]), "r"(a[2]), "r"(a[3]),
"r"(a[4]), "r"(a[5]), "r"(a[6]), "r"(a[7]),
"r"(b)
);
asm volatile(
"add.cc.u32 %0, %0, %9;nt" // res[0] += %lo32
"addc.cc.u32 %1, %1, %10;nt" // res[0] += %lo32
"addc.cc.u32 %2, %2, %11;nt" // res[0] += %lo32
"addc.cc.u32 %3, %3, %12;nt" // res[0] += %lo32
"addc.cc.u32 %4, %4, %13;nt" // res[0] += %lo32
"addc.cc.u32 %5, %5, %14;nt" // res[0] += %lo32
"addc.cc.u32 %6, %6, %15;nt" // res[0] += %lo32
"addc.cc.u32 %7, %7, %16;nt" // res[0] += %lo32
"addc.cc.u32 %8, %8, 0;nt" // res[0] += %lo32
: "+r"(res[0]), "+r"(res[1]), "+r"(res[2]), "+r"(res[3]),
"+r"(res[4]), "+r"(res[5]), "+r"(res[6]), "+r"(res[7]), "+r"(res[8])
: "r"(buf[0]), "r"(buf[1]), "r"(buf[2]), "r"(buf[3]),
"r"(buf[4]), "r"(buf[5]), "r"(buf[6]), "r"(buf[7])
);
}
Код обычный, популярный, его напишет на несколько микросекунд любой ИИ.
Ну казалось бы что тут улучшить?! Повторю опять, что алгоритм А.А,Карацубы в данной ситуации не рассматриваем, просто ищем улучшение простого алгоритма умножение-с-накоплением реализованного в столбик.
Вы удивитесь, но такое улучшение есть, весьма изящное и элегантное. Поиск по авторам алгоритмов автор статьи не проводил, первый или не первый раз публикуется не проверял. Но тем не менее придумал этот вариант сам, без какого-то ни было ИИ.
__device__ inline void cu_32_mul_add_to_512a(uint32_t *res, const uint32_t *a, const uint32_t b) {
uint64_t mul2 = 0;
uint32_t l32 = 0, h32 = 0;
asm volatile(
// even part
"mul.wide.u32 %9, %12, %20;nt" // %9 = a[0]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32}
"add.cc.u32 %0, %0, %10;nt" // res[0] += %lo32
"addc.cc.u32 %1, %1, %11;nt" // res[1] += %hi32
"mul.wide.u32 %9, %14, %20;nt" // %mul1 = a[2]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32}
"addc.cc.u32 %2, %2, %10;nt" // res[2] += %lo32
"addc.cc.u32 %3, %3, %11;nt" // res[3] += %hi32
"mul.wide.u32 %9, %16, %20;nt" // %mul1 = a[4]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32}
"addc.cc.u32 %4, %4, %10;nt" // res[4] += %lo32
"addc.cc.u32 %5, %5, %11;nt" // res[5] += %hi32
"mul.wide.u32 %9, %18, %20;nt" // %mul1 = a[6]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32}
"addc.cc.u32 %6, %6, %10;nt" // res[6] += %lo32
"addc.cc.u32 %7, %7, %11;nt" // res[7] += %hi32
"addc.u32 %8, %8, 0;nt" // res[8] += carry
//odd part
"mul.wide.u32 %9, %13, %20;nt" // %mul1 = a[1]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32}
"add.cc.u32 %1, %1, %10;nt" // res[1] += %lo32
"addc.cc.u32 %2, %2, %11;nt" // res[2] += %hi32
"mul.wide.u32 %9, %15, %20;nt" // %mul1 = a[3]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32}
"addc.cc.u32 %3, %3, %10;nt" // res[3] += %lo32
"addc.cc.u32 %4, %4, %11;nt" // res[4] += %hi32
"mul.wide.u32 %9, %17, %20;nt" // %mul1 = a[5]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32}
"addc.cc.u32 %5, %5, %10;nt" // res[5] += %lo32
"addc.cc.u32 %6, %6, %11;nt" // res[6] += %hi32
"mul.wide.u32 %9, %19, %20;nt" // %mul1 = a[7]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32c}
"addc.cc.u32 %7, %7, %10;nt" // res[7] += %lo32
"addc.u32 %8, %8, %11;nt" // res[7] += %hi32
: "+r"(res[0]), "+r"(res[1]), "+r"(res[2]), "+r"(res[3]),
"+r"(res[4]), "+r"(res[5]), "+r"(res[6]), "+r"(res[7]), "+r"(res[8]),
"=l"(mul2), "+r"(l32), "+r"(h32)
: "r"(a[0]), "r"(a[1]), "r"(a[2]), "r"(a[3]),
"r"(a[4]), "r"(a[5]), "r"(a[6]), "r"(a[7]),
"r"(b)
);
}Решение изящное, элегантное, экономит операции и регистры. Умножаются и прибавляются сначала четные слова, далее нечетные и перенос распространяется естественно.
Те. кому оно понравилось, могут вписать моё имя маленькими буквами рядом с именем Анатолия Александровича. Совсем маленькими буквами, но изящно и красиво.
Решение это хорошее, но на данном железе абсолютно бестолковое (об этом далее), увы. На другом не проверял, т.к. нет его у меня.
Проверим теперь быстродействие на реальном железе.
Проверялось на вот таком
Устройство 1
Compute Capability: 8.9
Имя: NVIDIA RTX 500 Ada Generation Laptop GPU
Общая глобальная память: -347209728
Разделяемая память на блок: 49152
Количество регистров на блок: 65536
Общая константная память: 65536
Максимальное количество нитей на блоке: 1024
Максимальное количество нитей в сетке: 1024 1024 64
Максимальный размер сетки: 2147483647 65535 65535
Размер wrap-a: 32
[Vector addition of 50000 elements]
CUDA kernel launch with 196 blocks of 256 threads
Test PASSEDДопишем умножение на скаляр до умножения двух больших 256 бит чисел по нашему любимому модулю P = “fffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f”, это модуль биткоин, весьма популярный для вычислений. Алгоритм быстрого подсчета модуля обычный и расписан везде, где можно, не буду его объяснять.
Далее будет много кода разного. Умножение по модулю 256 бит чисел в базе 32 и 64 и два алгоритма. Итого четыре варианта программ. Еще программа тестирования их и замера скорости. В коде всё будет видно.
вариант Б, классический, с базой 32 бит
#include "habr_256.h"
__device__ inline void cu_u32_set(uint32_t *res, const uint32_t *a) {
#pragma unroll
for (int i = 0; i < 8; i++) {
res[i] = a[i];
}
}
__device__ inline void cu_32_mul_add_to_512b(uint32_t *res, const uint32_t *a, const uint32_t b) {
uint64_t mul2 = 0;
uint32_t l32 = 0, h32 = 0, buf[9];
asm volatile(
"//mul to 288;nt"
"mul.wide.u32 %9, %12, %20;nt" // %9 = a[0]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32}
"add.cc.u32 %0, 0, %10;nt" // buf[0] += %lo32
"addc.cc.u32 %1, 0, %11;nt" // buf[1] += %hi32
"addc.u32 %2, 0, 0;nt"
"mul.wide.u32 %9, %13, %20;nt" // %mul1 = a[1]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32}
"add.cc.u32 %1, %1, %10;nt" // buf[1] += %lo32
"addc.cc.u32 %2, %2, %11;nt" // buf[2] += %hi32
"addc.u32 %3, 0, 0;nt"
"mul.wide.u32 %9, %14, %20;nt" // %mul1 = a[2]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32}
"addc.cc.u32 %2, %2, %10;nt" // buf[2] += %lo32
"addc.cc.u32 %3, %3, %11;nt" // buf[3] += %hi32
"addc.u32 %4, 0, 0;nt"
"mul.wide.u32 %9, %15, %20;nt" // %mul1 = a[3]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32}
"addc.cc.u32 %3, %3, %10;nt" // buf[3] += %lo32
"addc.cc.u32 %4, %4, %11;nt" // buf[4] += %hi32
"addc.u32 %5, 0, 0;nt"
"mul.wide.u32 %9, %16, %20;nt" // %mul1 = a[4]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32}
"addc.cc.u32 %4, %4, %10;nt" // buf[4] += %lo32
"addc.cc.u32 %5, %5, %11;nt" // buf[5] += %hi32
"addc.u32 %6, 0, 0;nt"
"mul.wide.u32 %9, %17, %20;nt" // %mul1 = a[5]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32}
"addc.cc.u32 %5, %5, %10;nt" // buf[5] += %lo32
"addc.cc.u32 %6, %6, %11;nt" // buf[6] += %hi32
"addc.u32 %7, 0, 0;nt"
"mul.wide.u32 %9, %18, %20;nt" // %mul1 = a[6]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32}
"addc.cc.u32 %6, %6, %10;nt" // buf[6] += %lo32
"addc.cc.u32 %7, %7, %11;nt" // buf[7] += %hi32
"addc.u32 %8, 0, 0;nt" //
"mul.wide.u32 %9, %19, %20;nt" // %mul1 = a[7]*b
"mov.b64 {%10, %11}, %9;nt" // %mul1 -> {%lo32, %hi32c}
"addc.cc.u32 %7, %7, %10;nt" // buf[7] += %lo32
"addc.cc.u32 %8, %8, %11;nt" // buf[8] += %hi32
: "+r"(buf[0]), "+r"(buf[1]), "+r"(buf[2]), "+r"(buf[3]),
"+r"(buf[4]), "+r"(buf[5]), "+r"(buf[6]), "+r"(buf[7]), "+r"(buf[8]),
"=l"(mul2), "+r"(l32), "+r"(h32)
: "r"(a[0]), "r"(a[1]), "r"(a[2]), "r"(a[3]),
"r"(a[4]), "r"(a[5]), "r"(a[6]), "r"(a[7]),
"r"(b)
);
asm volatile(
"add.cc.u32 %0, %0, %9;nt" // res[0] += buf
"addc.cc.u32 %1, %1, %10;nt" // res[1] += buf
"addc.cc.u32 %2, %2, %11;nt" // res[2] += buf
"addc.cc.u32 %3, %3, %12;nt" // res[3] += buf
"addc.cc.u32 %4, %4, %13;nt" // res[4] += buf
"addc.cc.u32 %5, %5, %14;nt" // res[5] += buf
"addc.cc.u32 %6, %6, %15;nt" // res[6] += buf
"addc.cc.u32 %7, %7, %16;nt" // res[7] += buf
"addc.cc.u32 %8, %8, 0;nt" // res[8] += buf
: "+r"(res[0]), "+r"(res[1]), "+r"(res[2]), "+r"(res[3]),
"+r"(res[4]), "+r"(res[5]), "+r"(res[6]), "+r"(res[7]), "+r"(res[8])
: "r"(buf[0]), "r"(buf[1]), "r"(buf[2]), "r"(buf[3]),
"r"(buf[4]), "r"(buf[5]), "r"(buf[6]), "r"(buf[7])
);
}
__device__ inline void cu_32_mul_mod_pb(uint32_t *res, const uint32_t *a, const uint32_t *b) {
uint32_t tmp[18] = {};
cu_32_mul_add_to_512b(&tmp[0], a, b[0]);
cu_32_mul_add_to_512b(&tmp[1], a, b[1]);
cu_32_mul_add_to_512b(&tmp[2], a, b[2]);
cu_32_mul_add_to_512b(&tmp[3], a, b[3]);
cu_32_mul_add_to_512b(&tmp[4], a, b[4]);
cu_32_mul_add_to_512b(&tmp[5], a, b[5]);
cu_32_mul_add_to_512b(&tmp[6], a, b[6]);
cu_32_mul_add_to_512b(&tmp[7], a, b[7]);
uint64_t mul = 0;
uint32_t l32 = 0, h32 = 0;
uint32_t t8=0, t9=0, t10=0, t0=0, t1=0, t2=0, t3=0;
// tmp[8-18] * 3d1 + tmp[0] -> tmp[0]
asm volatile(
"mul.wide.u32 %16, %8, 0x03d1;nt"
"mov.b64 {%17, %18}, %16;nt"
"add.cc.u32 %0, %0, %17;nt"
"addc.cc.u32 %1, %1, %18;nt"
"mul.wide.u32 %16, %10, 0x03d1;nt"
"mov.b64 {%17, %18}, %16;nt"
"addc.cc.u32 %2, %2, %17;nt"
"addc.cc.u32 %3, %3, %18;nt"
"mul.wide.u32 %16, %12, 0x03d1;nt"
"mov.b64 {%17, %18}, %16;nt"
"addc.cc.u32 %4, %4, %17;nt"
"addc.cc.u32 %5, %5, %18;nt"
"mul.wide.u32 %16, %14, 0x03d1;nt"
"mov.b64 {%17, %18}, %16;nt"
"addc.cc.u32 %6, %6, %17;nt"
"addc.cc.u32 %7, %7, %18;nt"
"addc.u32 %19, 0, 0;nt"
// -------------------------------------------------------
"mul.wide.u32 %16, %9, 0x03d1;nt"
"mov.b64 {%17, %18}, %16;nt"
"add.cc.u32 %1, %1, %17;nt"
"addc.cc.u32 %2, %2, %18;nt"
"mul.wide.u32 %16, %11, 0x03d1;nt"
"mov.b64 {%17, %18}, %16;nt"
"addc.cc.u32 %3, %3, %17;nt"
"addc.cc.u32 %4, %4, %18;nt"
"mul.wide.u32 %16, %13, 0x03d1;nt"
"mov.b64 {%17, %18}, %16;nt"
"addc.cc.u32 %5, %5, %17;nt"
"addc.cc.u32 %6, %6, %18;nt"
"mul.wide.u32 %16, %15, 0x03d1;nt"
"mov.b64 {%17, %18}, %16;nt"
"addc.cc.u32 %7, %7, %17;nt"
"addc.cc.u32 %19, %19, %18;nt"
"addc.u32 %20, 0, 0;nt"
: "+r"(tmp[0]), "+r"(tmp[1]), "+r"(tmp[2]), "+r"(tmp[3]),
"+r"(tmp[4]), "+r"(tmp[5]), "+r"(tmp[6]), "+r"(tmp[7]),
"+r"(tmp[8]), "+r"(tmp[9]), "+r"(tmp[10]), "+r"(tmp[11]),
"+r"(tmp[12]), "+r"(tmp[13]), "+r"(tmp[14]), "+r"(tmp[15]),
"+l"(mul), "+r"(l32), "+r"(h32),
"+r"(t8), "+r"(t9), "+r"(t10)
);
// *+ hi part of 1000003d1
// tmp[8-18] * 3d1 + tmp[0] -> tmp[0]
uint32_t carry;
asm volatile(
"add.cc.u32 %1, %1, %8;nt"
"addc.cc.u32 %2, %2, %9;nt"
"addc.cc.u32 %3, %3, %10;nt"
"addc.cc.u32 %4, %4, %11;nt"
"addc.cc.u32 %5, %5, %12;nt"
"addc.cc.u32 %6, %6, %13;nt"
"addc.cc.u32 %7, %7, %14;nt"
"addc.cc.u32 %19, %19, %15;nt"
"addc.u32 %20, 0, 0;nt"
: "+r"(tmp[0]), "+r"(tmp[1]), "+r"(tmp[2]), "+r"(tmp[3]),
"+r"(tmp[4]), "+r"(tmp[5]), "+r"(tmp[6]), "+r"(tmp[7]),
"+r"(tmp[8]), "+r"(tmp[9]), "+r"(tmp[10]), "+r"(tmp[11]),
"+r"(tmp[12]), "+r"(tmp[13]), "+r"(tmp[14]), "+r"(tmp[15]),
"+l"(mul), "+r"(l32), "+r"(h32),
"+r"(t8), "+r"(t9), "+r"(t10)
);
asm volatile(
// t8,t9,10 (buf[4]) нужно умножить на 100000000 + 000003d1 и добавить в tmp
"mov.u32 %1, %7;nt"
"mov.u32 %2, %8;nt"
"mov.u32 %3, %9;nt"
"mul.wide.u32 %4, %7, 0x000003d1;nt" // mul = %h32 * 3d1
"mov.b64 {%5, %6}, %4;nt"
"add.cc.u32 %0, %0, %5;nt"
"addc.cc.u32 %1, %1, %6;nt"
"mul.wide.u32 %4, %9, 0x000003d1;nt" // mul = %h32 * 3d1
"mov.b64 {%5, %6}, %4;nt"
"addc.cc.u32 %2, %2, %5;nt"
"addc.u32 %3, %3, %6;nt"
"mul.wide.u32 %4, %8, 0x000003d1;nt" // mul = %h32 * 3d1
"mov.b64 {%5, %6}, %4;nt"
"add.cc.u32 %1, %1, %5;nt"
"addc.cc.u32 %2, %2, %6;nt"
"addc.u32 %3, 0, 0;nt"
: "+r"(t0), "+r"(t1), "=r"(t2), "=r"(t3),
"+l"(mul), "+r"(l32), "+r"(h32)
: "r"(t8), "r"(t9), "r"(t10)
);
asm volatile( //tmp += t0t1t2t3 (t8t9t10*100003d1)
"add.cc.u32 %0, %0, %9;nt"
"addc.cc.u32 %1, %1, %10;nt"
"addc.cc.u32 %2, %2, %11;nt"
"addc.cc.u32 %3, %3, %12;nt"
"addc.cc.u32 %4, %4, 0;nt"
"addc.cc.u32 %5, %5, 0;nt"
"addc.cc.u32 %6, %6, 0;nt"
"addc.cc.u32 %7, %7, 0;nt"
"addc.u32 %8, 0, 0;nt"
: "+r"(tmp[0]), "+r"(tmp[1]), "+r"(tmp[2]), "+r"(tmp[3]),
"+r"(tmp[4]), "+r"(tmp[5]), "+r"(tmp[6]), "+r"(tmp[7]),
"=r"(carry)
: "r"(t0), "r"(t1), "r"(t2), "r"(t3)
);
if (carry != 0) {
asm volatile(
"add.cc.u32 %0, %9, 0x000003d1;nt"
"addc.cc.u32 %1, %10, 0x1;nt"
"addc.cc.u32 %2, %11, 0;nt"
"addc.cc.u32 %3, %12, 0;nt"
"addc.cc.u32 %4, %13, 0;nt"
"addc.cc.u32 %5, %14, 0;nt"
"addc.cc.u32 %6, %15, 0;nt"
"addc.cc.u32 %7, %16, 0;nt"
"addc.u32 %8, 0, 0;nt" // Add 0 + 0 + carry flag, write result to %1
: "=r"(tmp[0]), "=r"(tmp[1]), "=r"(tmp[2]), "=r"(tmp[3]),
"=r"(tmp[4]), "=r"(tmp[5]), "=r"(tmp[6]), "=r"(tmp[7]),
"=r"(carry)
: "r"(res[0]), "r"(res[1]), "r"(res[2]), "r"(res[3]),
"r"(res[4]), "r"(res[5]), "r"(res[6]), "r"(res[7])
);
}
asm volatile(
"add.cc.u32 %0, %9, 0x000003d1;nt"
"addc.cc.u32 %1, %10, 0x1;nt"
"addc.cc.u32 %2, %11, 0;nt"
"addc.cc.u32 %3, %12, 0;nt"
"addc.cc.u32 %4, %13, 0;nt"
"addc.cc.u32 %5, %14, 0;nt"
"addc.cc.u32 %6, %15, 0;nt"
"addc.cc.u32 %7, %16, 0;nt"
"addc.u32 %8, 0, 0;nt" // Add 0 + 0 + carry flag, write result to %1
: "=r"(tmp[8]), "=r"(tmp[9]), "=r"(tmp[10]), "=r"(tmp[11]),
"=r"(tmp[12]), "=r"(tmp[13]), "=r"(tmp[14]), "=r"(tmp[15]),
"=r"(carry)
: "r"(tmp[0]), "r"(tmp[1]), "r"(tmp[2]), "r"(tmp[3]),
"r"(tmp[4]), "r"(tmp[5]), "r"(tmp[6]), "r"(tmp[7])
);
if (carry == 0) {
#pragma unroll
for (int i = 0; i < 8; ++i) {
res[i] = tmp[i];
}
} else {
#pragma unroll
for (int i = 0; i < 8; ++i) {
res[i] = tmp[i+8];
}
}
}
__global__ void test_kernel_mul_mod_32b(Data * dev_r, Data * dev_a, Data * dev_b) {
uint64_t tid = threadIdx.x + blockIdx.x * blockDim.x;
if (tid >= CUDA_NUM_THREADS*CUDA_NUM_BLOCKS) return;
uint32_t a[8],b[8],r[8];
for (uint i = 0; i < CUDA_STRIDE; i++) {
uint64_t idx = tid*CUDA_STRIDE+i;
cu_u32_set(a, dev_a[idx].w32);
cu_u32_set(b, dev_b[idx].w32);
for (uint j = 0; j < T_SIZE; j++) {
cu_32_mul_mod_pb(r, a, b);
a[0] = a[0] ^ r[0] ;
}
cu_u32_set(dev_r[idx].w32, r);
}
}
Комментарии не нужны, тут и так ясно, что вызывается cu_32_mul_mod_pb на всех точках и на каждой выполняется T_SIZE раз. Сама cu_32_mul_mod_pb просто умножает два 256 бит числа a и b по модулю P и прибавляет их к результату r. Data тут простой union для того, чтобы спокойно далее запустить перемножение в 64 бит базе тех же самых 256 бит чисел.
Строка a[0] = a[0] ^ r[0] ; нужна для того, что иногда оптимизатор выполняет весь код в цикле по j всего один раз и наверно для этого есть основания, т.к. данные не меняются.
Хедер
#ifndef HABR_256_HABR_256_H
#define HABR_256_HABR_256_H
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <sys/random.h>
#include <errno.h>
#include <string.h>
#include <cuda_runtime.h>
#include <gmp.h>
#include <chrono>
#define CUDA_NUM_THREADS 32
#define CUDA_NUM_BLOCKS 32
#define CUDA_STRIDE 16*256
#define N_SIZE CUDA_NUM_THREADS*CUDA_NUM_BLOCKS*CUDA_STRIDE
#define T_SIZE 1
#define CUDA_CHECK(ans) { gpuAssert((ans), __FILE__, __LINE__); }
union Data {
uint64_t w64[4];
uint32_t w32[8];
};
void mpz_64_to_uint256(const mpz_t x, uint64_t * res);
void mpz_32_to_uint256(const mpz_t x, uint32_t * res);
__global__ void test_kernel_add_32(Data * dev_r, Data * dev_a, Data * dev_b);
__global__ void test_kernel_add_64(Data * dev_r, Data * dev_a, Data * dev_b);
__global__ void test_kernel_mul_32(Data * dev_r, Data * dev_a, uint32_t * dev_b);
__global__ void test_kernel_mul_64(Data * dev_r, Data * dev_a, const uint64_t * dev_b);
__global__ void test_kernel_mul_mod_64_b(Data * dev_r, Data * dev_a, const Data * dev_b);
__global__ void test_kernel_mul_mod_64_a(Data * dev_r, Data * dev_a, const Data * dev_b);
__global__ void test_kernel_mul_mod_32a(Data * dev_r, Data * dev_a, Data * dev_b);
__global__ void test_kernel_mul_mod_32b(Data * dev_r, Data * dev_a, Data * dev_b);
__device__ void u64_512_printf(const char *txt, uint64_t *a);
__device__ inline void cu_u64_set(uint64_t *res, const uint64_t *a);
__device__ inline void cu_u32_set(uint32_t *res, const uint32_t *a);
__device__ inline void cu_u64_mul_mod_p(uint64_t *res, const uint64_t *a, const uint64_t *b);
__device__ inline void cu_u64_mul_mod_pa(uint64_t *res, const uint64_t *a, const uint64_t *b);
__device__ void cu_u64_printf(const char *txt, const uint64_t *a);
__device__ void cu_u32_printf(const char *txt, const uint32_t *a);
__device__ void u64_512_printf(char *txt, uint64_t *a);
__device__ void cu_u32_512_printf(char *txt, uint32_t *a);
__device__ inline void cu_32_mul_mod_pa(uint32_t *res, const uint32_t *a, const uint32_t *b);
__device__ inline uint64_t cu_u64_add_3d1_with_carry(uint64_t *res, const uint64_t *x);
void u64_printf(const char *txt, const uint64_t *a);
void u32_printf(const char *txt, const uint32_t *a);
void u64_to_mpz(const uint64_t * x, mpz_t y);
void u32_to_mpz(const uint32_t * x, mpz_t y);
void gpuAssert(cudaError_t code, const char *file, int line);
#endif //HABR_256_HABR_256_H
И все оставшиеся три варианта, два для 64 битовых слов и один для 32 бит
вариант А с базой 32 бит
#include "habr_256.h"
__device__ inline void cu_32_mul_add_to_512a(uint32_t *res, const uint32_t *a, const uint32_t b) {
uint64_t mul2 = 0;
uint32_t l32 = 0, h32 = 0;
asm volatile(
: "+r"(res[0]), "+r"(res[1]), "+r"(res[2]), "+r"(res[3]),
"+r"(res[4]), "+r"(res[5]), "+r"(res[6]), "+r"(res[7]), "+r"(res[8]),
"=l"(mul2), "+r"(l32), "+r"(h32)
: "r"(a[0]), "r"(a[1]), "r"(a[2]), "r"(a[3]),
"r"(a[4]), "r"(a[5]), "r"(a[6]), "r"(a[7]),
"r"(b)
);
}
__device__ inline void cu_32_mul_mod_pa(uint32_t *res, const uint32_t *a, const uint32_t *b) {
uint32_t tmp[18] = {};
cu_32_mul_add_to_512a(&tmp[0], a, b[0]);
cu_32_mul_add_to_512a(&tmp[1], a, b[1]);
cu_32_mul_add_to_512a(&tmp[2], a, b[2]);
cu_32_mul_add_to_512a(&tmp[3], a, b[3]);
cu_32_mul_add_to_512a(&tmp[4], a, b[4]);
cu_32_mul_add_to_512a(&tmp[5], a, b[5]);
cu_32_mul_add_to_512a(&tmp[6], a, b[6]);
cu_32_mul_add_to_512a(&tmp[7], a, b[7]);
uint64_t mul = 0;
uint32_t l32 = 0, h32 = 0;
uint32_t t8=0, t9=0, t10=0, t0=0, t1=0, t2=0, t3=0;
// tmp[8-18] * 3d1 + tmp[0] -> tmp[0]
asm volatile(
: "+r"(tmp[0]), "+r"(tmp[1]), "+r"(tmp[2]), "+r"(tmp[3]),
"+r"(tmp[4]), "+r"(tmp[5]), "+r"(tmp[6]), "+r"(tmp[7]),
"+r"(tmp[8]), "+r"(tmp[9]), "+r"(tmp[10]), "+r"(tmp[11]),
"+r"(tmp[12]), "+r"(tmp[13]), "+r"(tmp[14]), "+r"(tmp[15]),
"+l"(mul), "+r"(l32), "+r"(h32),
"+r"(t8), "+r"(t9), "+r"(t10)
);
// *+ hi part of 1000003d1
// tmp[8-18] * 3d1 + tmp[0] -> tmp[0]
uint32_t carry;
asm volatile(
"add.cc.u32 %1, %1, %8;nt"
"addc.cc.u32 %2, %2, %9;nt"
"addc.cc.u32 %3, %3, %10;nt"
"addc.cc.u32 %4, %4, %11;nt"
"addc.cc.u32 %5, %5, %12;nt"
"addc.cc.u32 %6, %6, %13;nt"
"addc.cc.u32 %7, %7, %14;nt"
"addc.cc.u32 %19, %19, %15;nt"
"addc.u32 %20, 0, 0;nt"
: "+r"(tmp[0]), "+r"(tmp[1]), "+r"(tmp[2]), "+r"(tmp[3]),
"+r"(tmp[4]), "+r"(tmp[5]), "+r"(tmp[6]), "+r"(tmp[7]),
"+r"(tmp[8]), "+r"(tmp[9]), "+r"(tmp[10]), "+r"(tmp[11]),
"+r"(tmp[12]), "+r"(tmp[13]), "+r"(tmp[14]), "+r"(tmp[15]),
"+l"(mul), "+r"(l32), "+r"(h32),
"+r"(t8), "+r"(t9), "+r"(t10)
);
asm volatile(
// t8,t9,10 (buf[4]) нужно умножить на 100000000 + 000003d1 и добавить в tmp
"mov.u32 %1, %7;nt"
"mov.u32 %2, %8;nt"
"mov.u32 %3, %9;nt"
"mul.wide.u32 %4, %7, 0x000003d1;nt"
"mov.b64 {%5, %6}, %4;nt"
"add.cc.u32 %0, %0, %5;nt"
"addc.cc.u32 %1, %1, %6;nt"
"mul.wide.u32 %4, %9, 0x000003d1;nt"
"mov.b64 {%5, %6}, %4;nt"
"addc.cc.u32 %2, %2, %5;nt"
"addc.u32 %3, %3, %6;nt"
"mul.wide.u32 %4, %8, 0x000003d1;nt"
"mov.b64 {%5, %6}, %4;nt"
"add.cc.u32 %1, %1, %5;nt"
"addc.cc.u32 %2, %2, %6;nt"
"addc.u32 %3, 0, 0;nt"
: "+r"(t0), "+r"(t1), "=r"(t2), "=r"(t3),
"+l"(mul), "+r"(l32), "+r"(h32)
: "r"(t8), "r"(t9), "r"(t10)
);
asm volatile( //tmp += t0t1t2t3 (t8t9t10*100003d1)
"add.cc.u32 %0, %0, %9;nt"
"addc.cc.u32 %1, %1, %10;nt"
"addc.cc.u32 %2, %2, %11;nt"
"addc.cc.u32 %3, %3, %12;nt"
"addc.cc.u32 %4, %4, 0;nt"
"addc.cc.u32 %5, %5, 0;nt"
"addc.cc.u32 %6, %6, 0;nt"
"addc.cc.u32 %7, %7, 0;nt"
"addc.u32 %8, 0, 0;nt"
: "+r"(tmp[0]), "+r"(tmp[1]), "+r"(tmp[2]), "+r"(tmp[3]),
// 0, 1, 2, 3,
"+r"(tmp[4]), "+r"(tmp[5]), "+r"(tmp[6]), "+r"(tmp[7]),
// 4, 5, 6, 7,
"=r"(carry)
// 8,
: "r"(t0), "r"(t1), "r"(t2), "r"(t3)
// 9, 10, 11, 12
);
if (carry != 0) {
asm volatile(
"add.cc.u32 %0, %9, 0x000003d1;nt"
"addc.cc.u32 %1, %10, 0x1;nt"
"addc.cc.u32 %2, %11, 0;nt"
"addc.cc.u32 %3, %12, 0;nt"
"addc.cc.u32 %4, %13, 0;nt"
"addc.cc.u32 %5, %14, 0;nt"
"addc.cc.u32 %6, %15, 0;nt"
"addc.cc.u32 %7, %16, 0;nt"
"addc.u32 %8, 0, 0;nt"
: "=r"(tmp[0]), "=r"(tmp[1]), "=r"(tmp[2]), "=r"(tmp[3]),
"=r"(tmp[4]), "=r"(tmp[5]), "=r"(tmp[6]), "=r"(tmp[7]),
"=r"(carry)
: "r"(res[0]), "r"(res[1]), "r"(res[2]), "r"(res[3]),
"r"(res[4]), "r"(res[5]), "r"(res[6]), "r"(res[7])
);
}
asm volatile(
"add.cc.u32 %0, %9, 0x000003d1;nt"
"addc.cc.u32 %1, %10, 0x1;nt"
"addc.cc.u32 %2, %11, 0;nt"
"addc.cc.u32 %3, %12, 0;nt"
"addc.cc.u32 %4, %13, 0;nt"
"addc.cc.u32 %5, %14, 0;nt"
"addc.cc.u32 %6, %15, 0;nt"
"addc.cc.u32 %7, %16, 0;nt"
"addc.u32 %8, 0, 0;nt" // Add 0 + 0 + carry flag, write result to %1
: "=r"(tmp[8]), "=r"(tmp[9]), "=r"(tmp[10]), "=r"(tmp[11]),
"=r"(tmp[12]), "=r"(tmp[13]), "=r"(tmp[14]), "=r"(tmp[15]),
"=r"(carry)
: "r"(tmp[0]), "r"(tmp[1]), "r"(tmp[2]), "r"(tmp[3]),
"r"(tmp[4]), "r"(tmp[5]), "r"(tmp[6]), "r"(tmp[7])
);
if (carry == 0) {
#pragma unroll
for (int i = 0; i < 8; ++i) {
res[i] = tmp[i];
}
} else {
#pragma unroll
for (int i = 0; i < 8; ++i) {
res[i] = tmp[i+8];
}
}
}
__global__ void test_kernel_mul_mod_32a(Data * dev_r, Data * dev_a, Data * dev_b) {
uint64_t tid = threadIdx.x + blockIdx.x * blockDim.x;
if (tid >= CUDA_NUM_THREADS*CUDA_NUM_BLOCKS) return;
uint32_t a[8],b[8],r[8];
for (uint i = 0; i < CUDA_STRIDE; i++) {
uint64_t idx = tid*CUDA_STRIDE+i;
cu_u32_set(a, dev_a[idx].w32);
cu_u32_set(b, dev_b[idx].w32);
for (uint j = 0; j < T_SIZE; j++) {
cu_32_mul_mod_pa(r, a, b);
a[0] = a[0] ^ r[0];
}
cu_u32_set(dev_r[idx].w32, r);
}
}
Для полноты еще есть варианты А и В для базы 64 бит. Команды ассемблерные умножения другие, но алгоритм тот же
вариант Б, классический, с базой 64 бит
#include "habr_256.h"
__device__ inline void cu_u64_set(uint64_t *res, const uint64_t *a) {
#pragma unroll
for (int i = 0; i < 4; i++) {
res[i] = a[i];
}
}
__device__ inline uint64_t cu_u64_add_3d1_with_carry(uint64_t *res, const uint64_t *x) {
uint64_t carry;
asm volatile(
"add.cc.u64 %0, %5, 0x1000003d1;nt"
"addc.cc.u64 %1, %6, 0;nt"
"addc.cc.u64 %2, %7, 0;nt"
"addc.cc.u64 %3, %8, 0;nt"
"addc.u64 %4, 0, 0;nt"
: "=l"(res[0]), "=l"(res[1]), "=l"(res[2]), "=l"(res[3]),
"=l"(carry)
: "l"(x[0]), "l"(x[1]), "l"(x[2]), "l"(x[3])
);
return carry;
}
__device__ inline void cu_u64_mul_add_to_512(uint64_t *res, const uint64_t *a, const uint64_t b) {
uint64_t buf[5] = {};
asm volatile(
"mul.lo.u64 %0, %5, %9;nt"
"mul.hi.u64 %%rd0, %5, %9;nt"
"// Накопление со сдвигом по цепочкеnt"
"mad.lo.cc.u64 %1, %6, %9, %%rd0;nt"
"madc.hi.cc.u64 %%rd0, %6, %9, 0;nt"
"madc.lo.cc.u64 %2, %7, %9, %%rd0;nt"
"madc.hi.cc.u64 %%rd0, %7, %9, 0;nt"
"madc.lo.cc.u64 %3, %8, %9, %%rd0;nt"
"madc.hi.cc.u64 %4, %8, %9, 0;nt"
: "=l"(buf[0]), "=l"(buf[1]), "=l"(buf[2]), "=l"(buf[3]), "=l"(buf[4])
: "l"(a[0]), "l"(a[1]), "l"(a[2]), "l"(a[3]),
"l"(b)
);
asm volatile(
"add.cc.u64 %0, %0, %5;nt"
"addc.cc.u64 %1, %1, %6;nt"
"addc.cc.u64 %2, %2, %7;nt"
"addc.cc.u64 %3, %3, %8;nt"
"addc.cc.u64 %4, %4, %9;nt"
: "+l"(res[0]), "+l"(res[1]), "+l"(res[2]), "+l"(res[3]), "+l"(res[4])
: "l"(buf[0]), "l"(buf[1]), "l"(buf[2]), "l"(buf[3]), "l"(buf[4])
);
}
__device__ inline void cu_u64_mul_mod_p(uint64_t *res, const uint64_t *a, const uint64_t *b) {
uint64_t tmp[9] = {};
uint64_t buf[5] = {};
uint64_t carry = 0;
cu_u64_mul_add_to_512(&tmp[0], a, b[0]);
cu_u64_mul_add_to_512(&tmp[1], a, b[1]);
cu_u64_mul_add_to_512(&tmp[2], a, b[2]);
cu_u64_mul_add_to_512(&tmp[3], a, b[3]);
asm volatile(
"// tmp[4-7] * 10..03d1 -> buf[0-4]nt"
"mul.lo.u64 %0, %5, 0x1000003d1;nt"
"mul.hi.u64 %%rd0, %5, 0x1000003d1;nt"
"mad.lo.cc.u64 %1, %6, 0x1000003d1, %%rd0;nt"
"madc.hi.cc.u64 %%rd0, %6, 0x1000003d1, 0;nt"
"madc.lo.cc.u64 %2, %7, 0x1000003d1, %%rd0;nt"
"madc.hi.cc.u64 %%rd0, %7, 0x1000003d1, 0;nt"
"madc.lo.cc.u64 %3, %8, 0x1000003d1, %%rd0;nt"
"madc.hi.cc.u64 %4, %8, 0x1000003d1, 0;nt"
: "=l"(buf[0]), "=l"(buf[1]), "=l"(buf[2]), "=l"(buf[3]), "=l"(buf[4])
: "l"(tmp[4]), "l"(tmp[5]), "l"(tmp[6]), "l"(tmp[7])
);
asm volatile(
"// buf[0-4] + tmp[0-3] -> buf[0-4] nt"
"add.cc.u64 %0, %0, %5;nt"
"addc.cc.u64 %1, %1, %6;nt"
"addc.cc.u64 %2, %2, %7;nt"
"addc.cc.u64 %3, %3, %8;nt"
"addc.u64 %4, %4, 0;nt"
: "+l"(buf[0]), "+l"(buf[1]), "+l"(buf[2]), "+l"(buf[3]), "+l"(buf[4])
: "l"(tmp[0]), "l"(tmp[1]), "l"(tmp[2]), "l"(tmp[3]), "l"(tmp[4])
);
asm volatile(
"mad.lo.cc.u64 %0, %4, 0x1000003d1, %0;nt"
"madc.hi.cc.u64 %%rd0, %4, 0x1000003d1, 0;nt"
"addc.cc.u64 %1, %1, %%rd0;nt"
"addc.cc.u64 %2, %2, 0;nt"
"addc.cc.u64 %3, %3, 0;nt"
"addc.u64 %4, 0, 0;nt"
: "+l"(buf[0]), "+l"(buf[1]), "+l"(buf[2]), "+l"(buf[3]), "+l"(buf[4])
);
if (buf[4] != 0)
carry = cu_u64_add_3d1_with_carry(res, buf);
carry = cu_u64_add_3d1_with_carry(res, buf);
if (carry == 1) return;
cu_u64_set(res, buf);
}
__global__ void test_kernel_mul_mod_64_b(Data * dev_r, Data * dev_a, const Data * dev_b) {
uint_fast64_t tid = threadIdx.x + blockIdx.x * blockDim.x;
if (tid >= CUDA_NUM_THREADS*CUDA_NUM_BLOCKS) return;
uint64_t a[4],b[4],r[4];
for (uint i = 0; i < CUDA_STRIDE; i++) {
uint64_t idx = tid*CUDA_STRIDE+i;
cu_u64_set(a, dev_a[idx].w64);
cu_u64_set(b, dev_b[idx].w64);
for (uint j = 0; j < T_SIZE; j++) {
cu_u64_mul_mod_p(r, a, b);
a[0] = a[0] ^ (r[0] & 0x00000000ffffffffULL);
}
cu_u64_set(dev_r[idx].w64, r);
}
}
еще немного кода
Вариант А, модерновый, с базой 64 бита
#include "habr_256.h"
__device__ inline void cu_u64_set(uint64_t *res, const uint64_t *a) {
#pragma unroll
for (int i = 0; i < 4; i++) {
res[i] = a[i];
}
}
__device__ inline uint64_t cu_u64_add_3d1_inplace(uint64_t *res) {
uint64_t carry;
asm volatile(
"add.cc.u64 %0, %0, 0x1000003d1;nt"
"addc.cc.u64 %1, %1, 0;nt"
"addc.cc.u64 %2, %2, 0;nt"
"addc.cc.u64 %3, %3, 0;nt"
"addc.u64 %4, 0, 0;nt" // Add 0 + 0 + carry flag, write result to %1
: "+l"(res[0]), "+l"(res[1]), "+l"(res[2]), "+l"(res[3]),
"=l"(carry)
);
return carry;
}
__device__ inline uint64_t cu_u64_add_3d1_with_carry(uint64_t *res, const uint64_t *x) {
uint64_t carry;
asm volatile(
"add.cc.u64 %0, %5, 0x1000003d1;nt"
"addc.cc.u64 %1, %6, 0;nt"
"addc.cc.u64 %2, %7, 0;nt"
"addc.cc.u64 %3, %8, 0;nt"
"addc.u64 %4, 0, 0;nt" // Add 0 + 0 + carry flag, write result to %1
: "=l"(res[0]), "=l"(res[1]), "=l"(res[2]), "=l"(res[3]),
"=l"(carry)
: "l"(x[0]), "l"(x[1]), "l"(x[2]), "l"(x[3])
);
return carry;
}
__device__ __forceinline__ void cu_u64_mul_add_to_512a(uint64_t *res, const uint64_t *a, const uint64_t b) {
asm volatile(
"mad.lo.cc.u64 %0, %5, %9, %0;nt"
"madc.hi.cc.u64 %1, %5, %9, %1;nt"
"madc.lo.cc.u64 %2, %7, %9, %2;nt"
"madc.hi.cc.u64 %3, %7, %9, %3;nt"
"addc.u64 %4, 0,0;nt"
"mad.lo.cc.u64 %1, %6, %9, %1;nt"
"madc.hi.cc.u64 %2, %6, %9, %2;nt"
"madc.lo.cc.u64 %3, %8, %9, %3;nt"
"madc.hi.cc.u64 %4, %8, %9, %4;nt"
: "+l"(res[0]), "+l"(res[1]), "+l"(res[2]), "+l"(res[3]), "+l"(res[4])
: "l"(a[0]), "l"(a[1]), "l"(a[2]), "l"(a[3]),
"l"(b)
);
}
__device__ inline void cu_u64_mul_mod_pa(uint64_t *res, const uint64_t *a, const uint64_t *b) {
uint64_t tmp[9] = {0};
uint64_t carry = 0;
cu_u64_mul_add_to_512a(&tmp[0], a, b[0]);
cu_u64_mul_add_to_512a(&tmp[1], a, b[1]);
cu_u64_mul_add_to_512a(&tmp[2], a, b[2]);
cu_u64_mul_add_to_512a(&tmp[3], a, b[3]);
asm volatile(
"// tmp[4-7] * 10..03d1 + tmp[0-4] -> tmp[0-4]nt"
"mad.lo.cc.u64 %0, %4, 0x1000003d1, %0;nt"
"madc.hi.cc.u64 %1, %4, 0x1000003d1, %1;nt"
"// Накопление со сдвигом по цепочкеnt"
"madc.lo.cc.u64 %2, %6, 0x1000003d1, %2;nt"
"madc.hi.cc.u64 %3, %6, 0x1000003d1, %3;nt"
"madc.lo.cc.u64 %1, %5, 0x1000003d1, %1;nt"
"madc.hi.cc.u64 %2, %5, 0x1000003d1, %2;nt"
"madc.lo.cc.u64 %3, %7, 0x1000003d1, %3;nt"
"madc.hi.cc.u64 %4, %7, 0x1000003d1, 0;nt"
: "+l"(tmp[0]), "+l"(tmp[1]), "+l"(tmp[2]), "+l"(tmp[3]), "+l"(tmp[4])
: "l"(tmp[5]), "l"(tmp[6]), "l"(tmp[7])
);
asm volatile(
"// tmp[4] нужно умножить на 1000003d1 и добавить в tmpnt"
"mad.lo.cc.u64 %0, %4, 0x1000003d1, %0;nt"
"madc.hi.cc.u64 %%rd0, %4, 0x1000003d1, 0;nt"
"addc.cc.u64 %1, %1, %%rd0;nt"
"addc.cc.u64 %2, %2, 0;nt"
"addc.cc.u64 %3, %3, 0;nt"
"addc.u64 %4, 0, 0;nt"
: "+l"(tmp[0]), "+l"(tmp[1]), "+l"(tmp[2]), "+l"(tmp[3]), "+l"(tmp[4])
);
if (tmp[4] != 0)
carry = cu_u64_add_3d1_inplace(tmp);
carry = cu_u64_add_3d1_with_carry(res, tmp);
if (carry == 1) {
return;
}
cu_u64_set(res, tmp);
}
__global__ void test_kernel_mul_mod_64_a(Data * dev_r, Data * dev_a, const Data * dev_b) {
uint_fast64_t tid = threadIdx.x + blockIdx.x * blockDim.x;
if (tid >= CUDA_NUM_THREADS*CUDA_NUM_BLOCKS) return;
uint64_t a[4],b[4],r[4];
for (uint i = 0; i < CUDA_STRIDE; i++) {
uint64_t idx = tid*CUDA_STRIDE+i;
cu_u64_set(a, dev_a[idx].w64);
cu_u64_set(b, dev_b[idx].w64);
for (uint j = 0; j < T_SIZE; j++) {
cu_u64_mul_mod_pa(r, a, b);
a[0] = a[0] ^ (r[0] & 0x00000000ffffffffULL);
}
cu_u64_set(dev_r[idx].w64, r);
}
}
Теперь еще для полноты картины main.cu и хедер
main.cu
void test_mul_mod() {
timespec start{}, end{};
cudaError_t code;
long long start_ns;
long long end_ns;
long long diff_ns;
gmp_randstate_t state;
gmp_randinit_default(state);
const unsigned long seed = std::chrono::system_clock::now().time_since_epoch().count();
gmp_randseed_ui(state, seed);
Data *a = (Data *)malloc(sizeof(Data) * N_SIZE);
Data *b = (Data *)malloc(sizeof(Data) * N_SIZE);
Data *r32a = (Data *)malloc(sizeof(Data) * N_SIZE);
Data *r64a = (Data *)malloc(sizeof(Data) * N_SIZE);
Data *r32b = (Data *)malloc(sizeof(Data) * N_SIZE);
Data *r64b = (Data *)malloc(sizeof(Data) * N_SIZE);
mpz_t gmp_a, gmp_b, gmp_r, gmp_t, P;
mpz_inits(gmp_a, gmp_b, P, gmp_r, gmp_t, nullptr);
for (int j = 0; j < N_SIZE; j++) {
mpz_urandomm(gmp_a, state, P);
mpz_urandomm(gmp_b, state, P);
mpz_32_to_uint256(gmp_a, a[j].w32);
mpz_32_to_uint256(gmp_b, b[j].w32);
}
// Выделяем память на устройстве
Data *dev_a, *dev_r32a, *dev_r64a = nullptr, *dev_r32b, *dev_r64b= nullptr;
Data *dev_b;
CUDA_CHECK(cudaMalloc(&dev_a, N_SIZE * sizeof(Data)));
CUDA_CHECK(cudaMalloc(&dev_b, N_SIZE * sizeof(Data)));
CUDA_CHECK(cudaMalloc(&dev_r32a, N_SIZE * sizeof(Data)));
CUDA_CHECK(cudaMalloc(&dev_r64a, N_SIZE * sizeof(Data)));
CUDA_CHECK(cudaMalloc(&dev_r32b, N_SIZE * sizeof(Data)));
CUDA_CHECK(cudaMalloc(&dev_r64b, N_SIZE * sizeof(Data)));
for (int j = 0; j < N_SIZE; j++) {
r32a[j].w64[0] = 0;
r32a[j].w64[1] = 0;
r32a[j].w64[2] = 0;
r32a[j].w64[3] = 0;
r32b[j].w64[0] = 0;
r32b[j].w64[1] = 0;
r32b[j].w64[2] = 0;
r32b[j].w64[3] = 0;
r64b[j].w64[0] = 0;
r64a[j].w64[1] = 0;
r64a[j].w64[2] = 0;
r64a[j].w64[3] = 0;
r64b[j].w64[0] = 0;
r64b[j].w64[1] = 0;
r64b[j].w64[2] = 0;
r64b[j].w64[3] = 0;
}
CUDA_CHECK(cudaMemcpy(dev_a, a, N_SIZE * sizeof(Data), cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(dev_b, b, N_SIZE * sizeof(Data), cudaMemcpyHostToDevice));
// Запуск ядра 32a
clock_gettime(CLOCK_MONOTONIC, &start);
test_kernel_mul_mod_32a<<<CUDA_NUM_THREADS, CUDA_NUM_BLOCKS>>>(dev_r32a, dev_a, dev_b);
code = cudaGetLastError();
if (code != cudaSuccess) {
fprintf(stderr, "CUDA error: %s %s %dn", cudaGetErrorString(code), __FILE__, __LINE__);
exit(code);
}
CUDA_CHECK(cudaDeviceSynchronize()); // ждём и проверяем ошибку ядра
clock_gettime(CLOCK_MONOTONIC, &end);
start_ns = start.tv_sec * 1000000000LL + start.tv_nsec;
end_ns = end.tv_sec * 1000000000LL + end.tv_nsec;
diff_ns = end_ns - start_ns;
printf("mul_mod 32a bit время GPU: %lld нс (%.6f мс) точек %dn", diff_ns, (diff_ns / 1e6), N_SIZE);
// Получение результата
CUDA_CHECK(cudaMemcpy(r32a, dev_r32a, N_SIZE * sizeof(Data), cudaMemcpyDeviceToHost));
// ~~~~~~~~~~~~~~ //
// Запуск ядра 32b
clock_gettime(CLOCK_MONOTONIC, &start);
test_kernel_mul_mod_32b<<<CUDA_NUM_THREADS, CUDA_NUM_BLOCKS>>>(dev_r32b, dev_a, dev_b);
code = cudaGetLastError();
if (code != cudaSuccess) {
fprintf(stderr, "CUDA error: %s %s %dn", cudaGetErrorString(code), __FILE__, __LINE__);
exit(code);
}
CUDA_CHECK(cudaDeviceSynchronize()); // ждём и проверяем ошибку ядра
clock_gettime(CLOCK_MONOTONIC, &end);
start_ns = start.tv_sec * 1000000000LL + start.tv_nsec;
end_ns = end.tv_sec * 1000000000LL + end.tv_nsec;
diff_ns = end_ns - start_ns;
printf("mul_mod 32b bit время GPU: %lld нс (%.6f мс) точек %dn", diff_ns, (diff_ns / 1e6), N_SIZE);
// Получение результата
CUDA_CHECK(cudaMemcpy(r32b, dev_r32b, N_SIZE * sizeof(Data), cudaMemcpyDeviceToHost));
// ~~~~~~~~~~~~~~ //
// Запуск ядра 64a
clock_gettime(CLOCK_MONOTONIC, &start);
test_kernel_mul_mod_64_a<<<CUDA_NUM_THREADS, CUDA_NUM_BLOCKS>>>(dev_r64a, dev_a, dev_b);
code = cudaGetLastError();
if (code != cudaSuccess) {
fprintf(stderr, "CUDA error: %s %s %dn", cudaGetErrorString(code), __FILE__, __LINE__);
exit(code);
}
CUDA_CHECK(cudaDeviceSynchronize()); // ждём и проверяем ошибку ядра
clock_gettime(CLOCK_MONOTONIC, &end);
start_ns = start.tv_sec * 1000000000LL + start.tv_nsec;
end_ns = end.tv_sec * 1000000000LL + end.tv_nsec;
diff_ns = end_ns - start_ns;
printf("mul_mod 64a bit время GPU: %lld нс (%.6f мс) точек %dn", diff_ns, (diff_ns / 1e6), N_SIZE);
// Получение результата
CUDA_CHECK(cudaMemcpy(r64a, dev_r64a, N_SIZE * sizeof(Data), cudaMemcpyDeviceToHost));
// ~~~~~~~~~~~~~~ //
// Запуск ядра 64b
clock_gettime(CLOCK_MONOTONIC, &start);
test_kernel_mul_mod_64_b<<<CUDA_NUM_THREADS, CUDA_NUM_BLOCKS>>>(dev_r64b, dev_a, dev_b);
code = cudaGetLastError();
if (code != cudaSuccess) {
fprintf(stderr, "CUDA error: %s %s %dn", cudaGetErrorString(code), __FILE__, __LINE__);
exit(code);
}
CUDA_CHECK(cudaDeviceSynchronize()); // ждём и проверяем ошибку ядра
clock_gettime(CLOCK_MONOTONIC, &end);
start_ns = start.tv_sec * 1000000000LL + start.tv_nsec;
end_ns = end.tv_sec * 1000000000LL + end.tv_nsec;
diff_ns = end_ns - start_ns;
printf("mul_mod 64b bit время GPU: %lld нс (%.6f мс) точек %dn", diff_ns, (diff_ns / 1e6), N_SIZE);
// Получение результата
CUDA_CHECK(cudaMemcpy(r64b, dev_r64b, N_SIZE * sizeof(Data), cudaMemcpyDeviceToHost));
// ~~~~~~~~~~~~~~ //
// Запуск ядра 32a
clock_gettime(CLOCK_MONOTONIC, &start);
test_kernel_mul_mod_32a<<<CUDA_NUM_THREADS, CUDA_NUM_BLOCKS>>>(dev_r32a, dev_a, dev_b);
code = cudaGetLastError();
if (code != cudaSuccess) {
fprintf(stderr, "CUDA error: %s %s %dn", cudaGetErrorString(code), __FILE__, __LINE__);
exit(code);
}
CUDA_CHECK(cudaDeviceSynchronize()); // ждём и проверяем ошибку ядра
clock_gettime(CLOCK_MONOTONIC, &end);
start_ns = start.tv_sec * 1000000000LL + start.tv_nsec;
end_ns = end.tv_sec * 1000000000LL + end.tv_nsec;
diff_ns = end_ns - start_ns;
printf("mul_mod 32a bit время GPU: %lld нс (%.6f мс) точек %dn", diff_ns, (diff_ns / 1e6), N_SIZE);
// Получение результата
CUDA_CHECK(cudaMemcpy(r32a, dev_r32a, N_SIZE * sizeof(Data), cudaMemcpyDeviceToHost));
// ~~~~~~~~~~~~~~ //
// Запуск ядра 32b
clock_gettime(CLOCK_MONOTONIC, &start);
test_kernel_mul_mod_32b<<<CUDA_NUM_THREADS, CUDA_NUM_BLOCKS>>>(dev_r32b, dev_a, dev_b);
code = cudaGetLastError();
if (code != cudaSuccess) {
fprintf(stderr, "CUDA error: %s %s %dn", cudaGetErrorString(code), __FILE__, __LINE__);
exit(code);
}
CUDA_CHECK(cudaDeviceSynchronize()); // ждём и проверяем ошибку ядра
clock_gettime(CLOCK_MONOTONIC, &end);
start_ns = start.tv_sec * 1000000000LL + start.tv_nsec;
end_ns = end.tv_sec * 1000000000LL + end.tv_nsec;
diff_ns = end_ns - start_ns;
printf("mul_mod 32b bit время GPU: %lld нс (%.6f мс) точек %dn", diff_ns, (diff_ns / 1e6), N_SIZE);
// Получение результата
CUDA_CHECK(cudaMemcpy(r32b, dev_r32b, N_SIZE * sizeof(Data), cudaMemcpyDeviceToHost));
// ~~~~~~~~~~~~~~ //
// Запуск ядра 64a
clock_gettime(CLOCK_MONOTONIC, &start);
test_kernel_mul_mod_64_a<<<CUDA_NUM_THREADS, CUDA_NUM_BLOCKS>>>(dev_r64a, dev_a, dev_b);
code = cudaGetLastError();
if (code != cudaSuccess) {
fprintf(stderr, "CUDA error: %s %s %dn", cudaGetErrorString(code), __FILE__, __LINE__);
exit(code);
}
CUDA_CHECK(cudaDeviceSynchronize()); // ждём и проверяем ошибку ядра
clock_gettime(CLOCK_MONOTONIC, &end);
start_ns = start.tv_sec * 1000000000LL + start.tv_nsec;
end_ns = end.tv_sec * 1000000000LL + end.tv_nsec;
diff_ns = end_ns - start_ns;
printf("mul_mod 64a bit время GPU: %lld нс (%.6f мс) точек %dn", diff_ns, (diff_ns / 1e6), N_SIZE);
// Получение результата
CUDA_CHECK(cudaMemcpy(r64a, dev_r64a, N_SIZE * sizeof(Data), cudaMemcpyDeviceToHost));
// ~~~~~~~~~~~~~~ //
// Запуск ядра 64b
clock_gettime(CLOCK_MONOTONIC, &start);
test_kernel_mul_mod_64_b<<<CUDA_NUM_THREADS, CUDA_NUM_BLOCKS>>>(dev_r64b, dev_a, dev_b);
code = cudaGetLastError();
if (code != cudaSuccess) {
fprintf(stderr, "CUDA error: %s %s %dn", cudaGetErrorString(code), __FILE__, __LINE__);
exit(code);
}
CUDA_CHECK(cudaDeviceSynchronize()); // ждём и проверяем ошибку ядра
clock_gettime(CLOCK_MONOTONIC, &end);
start_ns = start.tv_sec * 1000000000LL + start.tv_nsec;
end_ns = end.tv_sec * 1000000000LL + end.tv_nsec;
diff_ns = end_ns - start_ns;
printf("mul_mod 64b bit время GPU: %lld нс (%.6f мс) точек %dn", diff_ns, (diff_ns / 1e6), N_SIZE);
// Получение результата
CUDA_CHECK(cudaMemcpy(r64b, dev_r64b, N_SIZE * sizeof(Data), cudaMemcpyDeviceToHost));
// ~~~~~~~~~~~~~~ //
for (uint i = 0; i < N_SIZE; i++) {
if ((
r64a[i].w64[3] ^ r64b[i].w64[3] |
r64a[i].w64[2] ^ r64b[i].w64[2] |
r64a[i].w64[1] ^ r64b[i].w64[1] |
r64a[i].w64[0] ^ r64b[i].w64[0]) == 0
) {
continue;
}
printf("Error n");
}
for (uint i = 0; i < min(4, N_SIZE); i++) {
u32_printf("r32a ", r32a[i].w32);
u32_printf("r32b ", r32a[i].w32);
u64_printf("r64a ", r64a[i].w64);
u64_printf("r64b ", r64a[i].w64);
}
CUDA_CHECK(cudaFree(dev_a));
CUDA_CHECK(cudaFree(dev_b));
CUDA_CHECK(cudaFree(dev_r64a));
CUDA_CHECK(cudaFree(dev_r32a));
CUDA_CHECK(cudaFree(dev_r64b));
CUDA_CHECK(cudaFree(dev_r32b));
}
int main() {
test_mul_mod();
return (0);
}
Результат получился совсем не тот, что ожидал. Дефолтный алгоритм с базой 64 бита оказался быстрее всех. Последние строки для визуально контроля результатов разных программ. Множители одинаковые и случайные для каждого экземпляра программы и грузятся один раз.
mul_mod 32a bit время GPU: 1784910771 нс (1784.910771 мс) точек 4194304
mul_mod 32b bit время GPU: 2278021284 нс (2278.021284 мс) точек 4194304
mul_mod 64a bit время GPU: 2635622032 нс (2635.622032 мс) точек 4194304
mul_mod 64b bit время GPU: 1536215011 нс (1536.215011 мс) точек 4194304
mul_mod 32a bit время GPU: 1739873265 нс (1739.873265 мс) точек 4194304
mul_mod 32b bit время GPU: 2277945774 нс (2277.945774 мс) точек 4194304
mul_mod 64a bit время GPU: 2635468058 нс (2635.468058 мс) точек 4194304
mul_mod 64b bit время GPU: 1541266850 нс (1541.266850 мс) точек 4194304
r32a 718620095f369811cbbbb4593c894bfc2adca86e5a50abb4a0dae1677b73aa60
r32b 718620095f369811cbbbb4593c894bfc2adca86e5a50abb4a0dae1677b73aa60
r64a 718620095f369811cbbbb4593c894bfc2adca86e5a50abb4a0dae1677b73aa60
r64b 718620095f369811cbbbb4593c894bfc2adca86e5a50abb4a0dae1677b73aa60
r32a 9580d7384bdd722df8d7614eb65fca9d51acb13aea79f1cbff8da83d19a5716c
r32b 9580d7384bdd722df8d7614eb65fca9d51acb13aea79f1cbff8da83d19a5716c
r64a 9580d7384bdd722df8d7614eb65fca9d51acb13aea79f1cbff8da83d19a5716c
r64b 9580d7384bdd722df8d7614eb65fca9d51acb13aea79f1cbff8da83d19a5716c
r32a ed452e82729b459144326bc75124b3088caaeca4ccfd806a1900c094a83c5613
r32b ed452e82729b459144326bc75124b3088caaeca4ccfd806a1900c094a83c5613
r64a ed452e82729b459144326bc75124b3088caaeca4ccfd806a1900c094a83c5613
r64b ed452e82729b459144326bc75124b3088caaeca4ccfd806a1900c094a83c5613
r32a bb4ee241c7154a1ec0e34e9ca29f05c0b394a800934ca85a8e58e782ccc5884f
r32b bb4ee241c7154a1ec0e34e9ca29f05c0b394a800934ca85a8e58e782ccc5884f
r64a bb4ee241c7154a1ec0e34e9ca29f05c0b394a800934ca85a8e58e782ccc5884f
r64b bb4ee241c7154a1ec0e34e9ca29f05c0b394a800934ca85a8e58e782ccc5884fУдивительные тонкости мира CUDA.
Поведение NVIDIA RTX 500 Ada Generation Laptop GPU оказалось весьма хитрым и, наверно, чип обладает своим каким-то неведомым искусственным интеллектом. Если закомментировать строку // a[0] = a[0] ^ (r[0] & 0x00000000ffffffffULL); во всех вариантах, то с теми же размерностями результат такой
mul_mod 32a bit время GPU: 534502115 нс (534.502115 мс) точек 4194304
mul_mod 32b bit время GPU: 179423263 нс (179.423263 мс) точек 4194304
mul_mod 64a bit время GPU: 35610833 нс (35.610833 мс) точек 4194304
mul_mod 64b bit время GPU: 8107617 нс (8.107617 мс) точек 4194304
mul_mod 32a bit время GPU: 472195051 нс (472.195051 мс) точек 4194304
mul_mod 32b bit время GPU: 179256781 нс (179.256781 мс) точек 4194304
mul_mod 64a bit время GPU: 35555436 нс (35.555436 мс) точек 4194304
mul_mod 64b bit время GPU: 8095543 нс (8.095543 мс) точек 4194304Строка эта не может использовать столько времени, это точно результат оптимизатора. Данные в том же формате, операций сложения/умножения столько же.
Обращений в файлах .ptx к глобальной памяти командой ld.global.u64 в 64 бит вариантах всего 8, в 32 бит всего 16. Т.е. изменение времени не из-за обращений к медленной глобальной памятию
Возможно у автора тут есть ошибки, и буду весьма признателен если кто найдет. Но проверял и прогонял тесты многократно, даже заново всё переписал как-то.
Если сделать T_SIZE = 1, т.е. убрать цикл и повторять один раз, то картина совсем другая, они все выполняются приблизительно одинаково и вариант “а” быстрее немного.
mul_mod 32a bit время GPU: 8847056 нс (8.847056 мс) точек 4194304
mul_mod 32b bit время GPU: 8656183 нс (8.656183 мс) точек 4194304
mul_mod 64a bit время GPU: 8132406 нс (8.132406 мс) точек 4194304
mul_mod 64b bit время GPU: 8190300 нс (8.190300 мс) точек 4194304
mul_mod 32a bit время GPU: 8304138 нс (8.304138 мс) точек 4194304
mul_mod 32b bit время GPU: 8612442 нс (8.612442 мс) точек 4194304
mul_mod 64a bit время GPU: 8081219 нс (8.081219 мс) точек 4194304
mul_mod 64b bit время GPU: 8102045 нс (8.102045 мс) точек 4194304Конечно же измерения сделаны простовато, нужно провести расчеты для кэшей и других ухищрений, но всё таки разница в вычислении простого умножения по модулю или незначительна или существенна и зависит от совсем ну постороннего параметра.
Вот такие тонкости нашего большого дела.
Автор: ChePeter
