В наше время, а дальше больше, знание умножения и сложения становится стратегическим. Вот то самое умножение в столбик, знакомое с начальной школы и сложение с переносом, знакомое еще раньше умножения теперь стратегический ресурс.
И причина простая как 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;\n\t" // %9 = a[0]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32} "add.cc.u32 %0, 0, %10;\n\t" // buf[0] += %lo32 "addc.cc.u32 %1, 0, %11;\n\t" // buf[1] += %hi32 "addc.u32 %2, 0, 0;\n\t" "mul.wide.u32 %9, %13, %20;\n\t" // %mul1 = a[1]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32} "add.cc.u32 %1, %1, %10;\n\t" // buf[1] += %lo32 "addc.cc.u32 %2, %2, %11;\n\t" // buf[2] += %hi32 "addc.u32 %3, 0, 0;\n\t" "mul.wide.u32 %9, %14, %20;\n\t" // %mul1 = a[2]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32} "addc.cc.u32 %2, %2, %10;\n\t" // buf[2] += %lo32 "addc.cc.u32 %3, %3, %11;\n\t" // buf[3] += %hi32 "addc.u32 %4, 0, 0;\n\t" "mul.wide.u32 %9, %15, %20;\n\t" // %mul1 = a[3]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32} "addc.cc.u32 %3, %3, %10;\n\t" // buf[3] += %lo32 "addc.cc.u32 %4, %4, %11;\n\t" // buf[4] += %hi32 "addc.u32 %5, 0, 0;\n\t" "mul.wide.u32 %9, %16, %20;\n\t" // %mul1 = a[4]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32} "addc.cc.u32 %4, %4, %10;\n\t" // buf[4] += %lo32 "addc.cc.u32 %5, %5, %11;\n\t" // buf[5] += %hi32 "addc.u32 %6, 0, 0;\n\t" "mul.wide.u32 %9, %17, %20;\n\t" // %mul1 = a[5]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32} "addc.cc.u32 %5, %5, %10;\n\t" // buf[5] += %lo32 "addc.cc.u32 %6, %6, %11;\n\t" // buf[6] += %hi32 "addc.u32 %7, 0, 0;\n\t" "mul.wide.u32 %9, %18, %20;\n\t" // %mul1 = a[6]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32} "addc.cc.u32 %6, %6, %10;\n\t" // buf[6] += %lo32 "addc.cc.u32 %7, %7, %11;\n\t" // buf[7] += %hi32 "addc.u32 %8, 0, 0;\n\t" // "mul.wide.u32 %9, %19, %20;\n\t" // %mul1 = a[7]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32c} "addc.cc.u32 %7, %7, %10;\n\t" // buf[7] += %lo32 "addc.cc.u32 %8, %8, %11;\n\t" // 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;\n\t" // res[0] += %lo32 "addc.cc.u32 %1, %1, %10;\n\t" // res[0] += %lo32 "addc.cc.u32 %2, %2, %11;\n\t" // res[0] += %lo32 "addc.cc.u32 %3, %3, %12;\n\t" // res[0] += %lo32 "addc.cc.u32 %4, %4, %13;\n\t" // res[0] += %lo32 "addc.cc.u32 %5, %5, %14;\n\t" // res[0] += %lo32 "addc.cc.u32 %6, %6, %15;\n\t" // res[0] += %lo32 "addc.cc.u32 %7, %7, %16;\n\t" // res[0] += %lo32 "addc.cc.u32 %8, %8, 0;\n\t" // 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;\n\t" // %9 = a[0]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32} "add.cc.u32 %0, %0, %10;\n\t" // res[0] += %lo32 "addc.cc.u32 %1, %1, %11;\n\t" // res[1] += %hi32 "mul.wide.u32 %9, %14, %20;\n\t" // %mul1 = a[2]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32} "addc.cc.u32 %2, %2, %10;\n\t" // res[2] += %lo32 "addc.cc.u32 %3, %3, %11;\n\t" // res[3] += %hi32 "mul.wide.u32 %9, %16, %20;\n\t" // %mul1 = a[4]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32} "addc.cc.u32 %4, %4, %10;\n\t" // res[4] += %lo32 "addc.cc.u32 %5, %5, %11;\n\t" // res[5] += %hi32 "mul.wide.u32 %9, %18, %20;\n\t" // %mul1 = a[6]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32} "addc.cc.u32 %6, %6, %10;\n\t" // res[6] += %lo32 "addc.cc.u32 %7, %7, %11;\n\t" // res[7] += %hi32 "addc.u32 %8, %8, 0;\n\t" // res[8] += carry //odd part "mul.wide.u32 %9, %13, %20;\n\t" // %mul1 = a[1]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32} "add.cc.u32 %1, %1, %10;\n\t" // res[1] += %lo32 "addc.cc.u32 %2, %2, %11;\n\t" // res[2] += %hi32 "mul.wide.u32 %9, %15, %20;\n\t" // %mul1 = a[3]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32} "addc.cc.u32 %3, %3, %10;\n\t" // res[3] += %lo32 "addc.cc.u32 %4, %4, %11;\n\t" // res[4] += %hi32 "mul.wide.u32 %9, %17, %20;\n\t" // %mul1 = a[5]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32} "addc.cc.u32 %5, %5, %10;\n\t" // res[5] += %lo32 "addc.cc.u32 %6, %6, %11;\n\t" // res[6] += %hi32 "mul.wide.u32 %9, %19, %20;\n\t" // %mul1 = a[7]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32c} "addc.cc.u32 %7, %7, %10;\n\t" // res[7] += %lo32 "addc.u32 %8, %8, %11;\n\t" // 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;\n\t" "mul.wide.u32 %9, %12, %20;\n\t" // %9 = a[0]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32} "add.cc.u32 %0, 0, %10;\n\t" // buf[0] += %lo32 "addc.cc.u32 %1, 0, %11;\n\t" // buf[1] += %hi32 "addc.u32 %2, 0, 0;\n\t" "mul.wide.u32 %9, %13, %20;\n\t" // %mul1 = a[1]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32} "add.cc.u32 %1, %1, %10;\n\t" // buf[1] += %lo32 "addc.cc.u32 %2, %2, %11;\n\t" // buf[2] += %hi32 "addc.u32 %3, 0, 0;\n\t" "mul.wide.u32 %9, %14, %20;\n\t" // %mul1 = a[2]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32} "addc.cc.u32 %2, %2, %10;\n\t" // buf[2] += %lo32 "addc.cc.u32 %3, %3, %11;\n\t" // buf[3] += %hi32 "addc.u32 %4, 0, 0;\n\t" "mul.wide.u32 %9, %15, %20;\n\t" // %mul1 = a[3]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32} "addc.cc.u32 %3, %3, %10;\n\t" // buf[3] += %lo32 "addc.cc.u32 %4, %4, %11;\n\t" // buf[4] += %hi32 "addc.u32 %5, 0, 0;\n\t" "mul.wide.u32 %9, %16, %20;\n\t" // %mul1 = a[4]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32} "addc.cc.u32 %4, %4, %10;\n\t" // buf[4] += %lo32 "addc.cc.u32 %5, %5, %11;\n\t" // buf[5] += %hi32 "addc.u32 %6, 0, 0;\n\t" "mul.wide.u32 %9, %17, %20;\n\t" // %mul1 = a[5]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32} "addc.cc.u32 %5, %5, %10;\n\t" // buf[5] += %lo32 "addc.cc.u32 %6, %6, %11;\n\t" // buf[6] += %hi32 "addc.u32 %7, 0, 0;\n\t" "mul.wide.u32 %9, %18, %20;\n\t" // %mul1 = a[6]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32} "addc.cc.u32 %6, %6, %10;\n\t" // buf[6] += %lo32 "addc.cc.u32 %7, %7, %11;\n\t" // buf[7] += %hi32 "addc.u32 %8, 0, 0;\n\t" // "mul.wide.u32 %9, %19, %20;\n\t" // %mul1 = a[7]*b "mov.b64 {%10, %11}, %9;\n\t" // %mul1 -> {%lo32, %hi32c} "addc.cc.u32 %7, %7, %10;\n\t" // buf[7] += %lo32 "addc.cc.u32 %8, %8, %11;\n\t" // 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;\n\t" // res[0] += buf "addc.cc.u32 %1, %1, %10;\n\t" // res[1] += buf "addc.cc.u32 %2, %2, %11;\n\t" // res[2] += buf "addc.cc.u32 %3, %3, %12;\n\t" // res[3] += buf "addc.cc.u32 %4, %4, %13;\n\t" // res[4] += buf "addc.cc.u32 %5, %5, %14;\n\t" // res[5] += buf "addc.cc.u32 %6, %6, %15;\n\t" // res[6] += buf "addc.cc.u32 %7, %7, %16;\n\t" // res[7] += buf "addc.cc.u32 %8, %8, 0;\n\t" // 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;\n\t" "mov.b64 {%17, %18}, %16;\n\t" "add.cc.u32 %0, %0, %17;\n\t" "addc.cc.u32 %1, %1, %18;\n\t" "mul.wide.u32 %16, %10, 0x03d1;\n\t" "mov.b64 {%17, %18}, %16;\n\t" "addc.cc.u32 %2, %2, %17;\n\t" "addc.cc.u32 %3, %3, %18;\n\t" "mul.wide.u32 %16, %12, 0x03d1;\n\t" "mov.b64 {%17, %18}, %16;\n\t" "addc.cc.u32 %4, %4, %17;\n\t" "addc.cc.u32 %5, %5, %18;\n\t" "mul.wide.u32 %16, %14, 0x03d1;\n\t" "mov.b64 {%17, %18}, %16;\n\t" "addc.cc.u32 %6, %6, %17;\n\t" "addc.cc.u32 %7, %7, %18;\n\t" "addc.u32 %19, 0, 0;\n\t" // ------------------------------------------------------- "mul.wide.u32 %16, %9, 0x03d1;\n\t" "mov.b64 {%17, %18}, %16;\n\t" "add.cc.u32 %1, %1, %17;\n\t" "addc.cc.u32 %2, %2, %18;\n\t" "mul.wide.u32 %16, %11, 0x03d1;\n\t" "mov.b64 {%17, %18}, %16;\n\t" "addc.cc.u32 %3, %3, %17;\n\t" "addc.cc.u32 %4, %4, %18;\n\t" "mul.wide.u32 %16, %13, 0x03d1;\n\t" "mov.b64 {%17, %18}, %16;\n\t" "addc.cc.u32 %5, %5, %17;\n\t" "addc.cc.u32 %6, %6, %18;\n\t" "mul.wide.u32 %16, %15, 0x03d1;\n\t" "mov.b64 {%17, %18}, %16;\n\t" "addc.cc.u32 %7, %7, %17;\n\t" "addc.cc.u32 %19, %19, %18;\n\t" "addc.u32 %20, 0, 0;\n\t" : "+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;\n\t" "addc.cc.u32 %2, %2, %9;\n\t" "addc.cc.u32 %3, %3, %10;\n\t" "addc.cc.u32 %4, %4, %11;\n\t" "addc.cc.u32 %5, %5, %12;\n\t" "addc.cc.u32 %6, %6, %13;\n\t" "addc.cc.u32 %7, %7, %14;\n\t" "addc.cc.u32 %19, %19, %15;\n\t" "addc.u32 %20, 0, 0;\n\t" : "+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;\n\t" "mov.u32 %2, %8;\n\t" "mov.u32 %3, %9;\n\t" "mul.wide.u32 %4, %7, 0x000003d1;\n\t" // mul = %h32 * 3d1 "mov.b64 {%5, %6}, %4;\n\t" "add.cc.u32 %0, %0, %5;\n\t" "addc.cc.u32 %1, %1, %6;\n\t" "mul.wide.u32 %4, %9, 0x000003d1;\n\t" // mul = %h32 * 3d1 "mov.b64 {%5, %6}, %4;\n\t" "addc.cc.u32 %2, %2, %5;\n\t" "addc.u32 %3, %3, %6;\n\t" "mul.wide.u32 %4, %8, 0x000003d1;\n\t" // mul = %h32 * 3d1 "mov.b64 {%5, %6}, %4;\n\t" "add.cc.u32 %1, %1, %5;\n\t" "addc.cc.u32 %2, %2, %6;\n\t" "addc.u32 %3, 0, 0;\n\t" : "+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;\n\t" "addc.cc.u32 %1, %1, %10;\n\t" "addc.cc.u32 %2, %2, %11;\n\t" "addc.cc.u32 %3, %3, %12;\n\t" "addc.cc.u32 %4, %4, 0;\n\t" "addc.cc.u32 %5, %5, 0;\n\t" "addc.cc.u32 %6, %6, 0;\n\t" "addc.cc.u32 %7, %7, 0;\n\t" "addc.u32 %8, 0, 0;\n\t" : "+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;\n\t" "addc.cc.u32 %1, %10, 0x1;\n\t" "addc.cc.u32 %2, %11, 0;\n\t" "addc.cc.u32 %3, %12, 0;\n\t" "addc.cc.u32 %4, %13, 0;\n\t" "addc.cc.u32 %5, %14, 0;\n\t" "addc.cc.u32 %6, %15, 0;\n\t" "addc.cc.u32 %7, %16, 0;\n\t" "addc.u32 %8, 0, 0;\n\t" // 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;\n\t" "addc.cc.u32 %1, %10, 0x1;\n\t" "addc.cc.u32 %2, %11, 0;\n\t" "addc.cc.u32 %3, %12, 0;\n\t" "addc.cc.u32 %4, %13, 0;\n\t" "addc.cc.u32 %5, %14, 0;\n\t" "addc.cc.u32 %6, %15, 0;\n\t" "addc.cc.u32 %7, %16, 0;\n\t" "addc.u32 %8, 0, 0;\n\t" // 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;\n\t" "addc.cc.u32 %2, %2, %9;\n\t" "addc.cc.u32 %3, %3, %10;\n\t" "addc.cc.u32 %4, %4, %11;\n\t" "addc.cc.u32 %5, %5, %12;\n\t" "addc.cc.u32 %6, %6, %13;\n\t" "addc.cc.u32 %7, %7, %14;\n\t" "addc.cc.u32 %19, %19, %15;\n\t" "addc.u32 %20, 0, 0;\n\t" : "+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;\n\t" "mov.u32 %2, %8;\n\t" "mov.u32 %3, %9;\n\t" "mul.wide.u32 %4, %7, 0x000003d1;\n\t" "mov.b64 {%5, %6}, %4;\n\t" "add.cc.u32 %0, %0, %5;\n\t" "addc.cc.u32 %1, %1, %6;\n\t" "mul.wide.u32 %4, %9, 0x000003d1;\n\t" "mov.b64 {%5, %6}, %4;\n\t" "addc.cc.u32 %2, %2, %5;\n\t" "addc.u32 %3, %3, %6;\n\t" "mul.wide.u32 %4, %8, 0x000003d1;\n\t" "mov.b64 {%5, %6}, %4;\n\t" "add.cc.u32 %1, %1, %5;\n\t" "addc.cc.u32 %2, %2, %6;\n\t" "addc.u32 %3, 0, 0;\n\t" : "+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;\n\t" "addc.cc.u32 %1, %1, %10;\n\t" "addc.cc.u32 %2, %2, %11;\n\t" "addc.cc.u32 %3, %3, %12;\n\t" "addc.cc.u32 %4, %4, 0;\n\t" "addc.cc.u32 %5, %5, 0;\n\t" "addc.cc.u32 %6, %6, 0;\n\t" "addc.cc.u32 %7, %7, 0;\n\t" "addc.u32 %8, 0, 0;\n\t" : "+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;\n\t" "addc.cc.u32 %1, %10, 0x1;\n\t" "addc.cc.u32 %2, %11, 0;\n\t" "addc.cc.u32 %3, %12, 0;\n\t" "addc.cc.u32 %4, %13, 0;\n\t" "addc.cc.u32 %5, %14, 0;\n\t" "addc.cc.u32 %6, %15, 0;\n\t" "addc.cc.u32 %7, %16, 0;\n\t" "addc.u32 %8, 0, 0;\n\t" : "=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;\n\t" "addc.cc.u32 %1, %10, 0x1;\n\t" "addc.cc.u32 %2, %11, 0;\n\t" "addc.cc.u32 %3, %12, 0;\n\t" "addc.cc.u32 %4, %13, 0;\n\t" "addc.cc.u32 %5, %14, 0;\n\t" "addc.cc.u32 %6, %15, 0;\n\t" "addc.cc.u32 %7, %16, 0;\n\t" "addc.u32 %8, 0, 0;\n\t" // 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;\n\t" "addc.cc.u64 %1, %6, 0;\n\t" "addc.cc.u64 %2, %7, 0;\n\t" "addc.cc.u64 %3, %8, 0;\n\t" "addc.u64 %4, 0, 0;\n\t" : "=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;\n\t" "mul.hi.u64 %%rd0, %5, %9;\n\t" "// Накопление со сдвигом по цепочке\n\t" "mad.lo.cc.u64 %1, %6, %9, %%rd0;\n\t" "madc.hi.cc.u64 %%rd0, %6, %9, 0;\n\t" "madc.lo.cc.u64 %2, %7, %9, %%rd0;\n\t" "madc.hi.cc.u64 %%rd0, %7, %9, 0;\n\t" "madc.lo.cc.u64 %3, %8, %9, %%rd0;\n\t" "madc.hi.cc.u64 %4, %8, %9, 0;\n\t" : "=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;\n\t" "addc.cc.u64 %1, %1, %6;\n\t" "addc.cc.u64 %2, %2, %7;\n\t" "addc.cc.u64 %3, %3, %8;\n\t" "addc.cc.u64 %4, %4, %9;\n\t" : "+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]\n\t" "mul.lo.u64 %0, %5, 0x1000003d1;\n\t" "mul.hi.u64 %%rd0, %5, 0x1000003d1;\n\t" "mad.lo.cc.u64 %1, %6, 0x1000003d1, %%rd0;\n\t" "madc.hi.cc.u64 %%rd0, %6, 0x1000003d1, 0;\n\t" "madc.lo.cc.u64 %2, %7, 0x1000003d1, %%rd0;\n\t" "madc.hi.cc.u64 %%rd0, %7, 0x1000003d1, 0;\n\t" "madc.lo.cc.u64 %3, %8, 0x1000003d1, %%rd0;\n\t" "madc.hi.cc.u64 %4, %8, 0x1000003d1, 0;\n\t" : "=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] \n\t" "add.cc.u64 %0, %0, %5;\n\t" "addc.cc.u64 %1, %1, %6;\n\t" "addc.cc.u64 %2, %2, %7;\n\t" "addc.cc.u64 %3, %3, %8;\n\t" "addc.u64 %4, %4, 0;\n\t" : "+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;\n\t" "madc.hi.cc.u64 %%rd0, %4, 0x1000003d1, 0;\n\t" "addc.cc.u64 %1, %1, %%rd0;\n\t" "addc.cc.u64 %2, %2, 0;\n\t" "addc.cc.u64 %3, %3, 0;\n\t" "addc.u64 %4, 0, 0;\n\t" : "+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;\n\t" "addc.cc.u64 %1, %1, 0;\n\t" "addc.cc.u64 %2, %2, 0;\n\t" "addc.cc.u64 %3, %3, 0;\n\t" "addc.u64 %4, 0, 0;\n\t" // 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;\n\t" "addc.cc.u64 %1, %6, 0;\n\t" "addc.cc.u64 %2, %7, 0;\n\t" "addc.cc.u64 %3, %8, 0;\n\t" "addc.u64 %4, 0, 0;\n\t" // 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;\n\t" "madc.hi.cc.u64 %1, %5, %9, %1;\n\t" "madc.lo.cc.u64 %2, %7, %9, %2;\n\t" "madc.hi.cc.u64 %3, %7, %9, %3;\n\t" "addc.u64 %4, 0,0;\n\t" "mad.lo.cc.u64 %1, %6, %9, %1;\n\t" "madc.hi.cc.u64 %2, %6, %9, %2;\n\t" "madc.lo.cc.u64 %3, %8, %9, %3;\n\t" "madc.hi.cc.u64 %4, %8, %9, %4;\n\t" : "+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]\n\t" "mad.lo.cc.u64 %0, %4, 0x1000003d1, %0;\n\t" "madc.hi.cc.u64 %1, %4, 0x1000003d1, %1;\n\t" "// Накопление со сдвигом по цепочке\n\t" "madc.lo.cc.u64 %2, %6, 0x1000003d1, %2;\n\t" "madc.hi.cc.u64 %3, %6, 0x1000003d1, %3;\n\t" "madc.lo.cc.u64 %1, %5, 0x1000003d1, %1;\n\t" "madc.hi.cc.u64 %2, %5, 0x1000003d1, %2;\n\t" "madc.lo.cc.u64 %3, %7, 0x1000003d1, %3;\n\t" "madc.hi.cc.u64 %4, %7, 0x1000003d1, 0;\n\t" : "+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 и добавить в tmp\n\t" "mad.lo.cc.u64 %0, %4, 0x1000003d1, %0;\n\t" "madc.hi.cc.u64 %%rd0, %4, 0x1000003d1, 0;\n\t" "addc.cc.u64 %1, %1, %%rd0;\n\t" "addc.cc.u64 %2, %2, 0;\n\t" "addc.cc.u64 %3, %3, 0;\n\t" "addc.u64 %4, 0, 0;\n\t" : "+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 %d\n", 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 мс) точек %d\n", 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 %d\n", 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 мс) точек %d\n", 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 %d\n", 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 мс) точек %d\n", 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 %d\n", 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 мс) точек %d\n", 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 %d\n", 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 мс) точек %d\n", 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 %d\n", 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 мс) точек %d\n", 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 %d\n", 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 мс) точек %d\n", 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 %d\n", 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 мс) точек %d\n", 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
Конечно же измерения сделаны простовато, нужно провести расчеты для кэшей и других ухищрений, но всё таки разница в вычислении простого умножения по модулю или незначительна или существенна и зависит от совсем ну постороннего параметра.
Вот такие тонкости нашего большого дела.
