This is the solution to tutorial in m4-tutorial-croatia2023 .
brew install qemu
Please make sure that the ‘helloworld’ example successfully outputs
Hello Croatia 2023!
when run in qemu prior to attending the tutorial. Installation instructions are available above.
You can build and run the example using the following steps:
cd helloworld
make
make run-qemu
- Armv7E-M architecture, released in 2010,
- stages (fetch, decode, execute), branching breaks pipeline (consider also branch prediction and speculative execution)
- r0-r15, r13 (sp, stack pointer), r14 (lr, link register), r15 (pc, program counter), r0-r12 are general purpose, and r14 can also be freely used after being saved to the stack
- Instr Rd, Rn(, Rm), many instructions have a variant that sets flags by appending s
- barrel shifter
ARM 独有的特征,是一种常见的优化方法,实现一条指令实现两条指令的效果,仅花费一个周期。
例如
mov r0 , #42 mov r1 , #37 ror r1 , r1 , #1 orr r2 , r0 , r1 lsl r2 , r2 , #1 eor r0 , r2
又可以以这种优化的方式实现
mov r0, #42 mov r1, #37 orr r2, r0, r1, ror #1 eor r0, r2, r2, lsl #1
barrel shifter 不会更新 Rm 中的值,只是会得到临时的结果用于 Instr
- conditional branchs cmp r0, r1/immediate bxx label, jump with flag corresponding to xx (\==, !=, >, >=, <, <=)
- conditional execution if-then-else(ITE), up to 4 instructions
- stack Used when data does not fit in registers, push {r0, r1}, pop {r0, r2} (assigned register by user)
- memory
- memory pipeline
- a single ldr take 2 cycles (when not stalled), N consecutive ldr take N+1 cycles
- str takes 1 cycle, does not pipeline
- ldrd/strd/ldm/stm do not pipeline, takes N+1 cycles when load N words
- For more details look at https://developer.arm.com/documentation/ddi0439/b/Programmers-Model/Instruction-set-summary/Load-store-timings
- lr keep track of ‘return address’
- arguments of function 通过 r0-r3 用于存储实参,如果多于 4 个参数,则先将第五个参数存在 r0 中再将 r0 压入栈中。也就是在最右边的会先被压入栈。 long 和 double 类型的参数需要使用两个寄存器存储
- AAPCS (ARM Architecture Procedure Call Standard)
当运行参考实现时, qemu 上的执行时间大致在这附近波动
====== START ====== crypto_stream_chacha20: 19300 cycles (note that these are meaningless on qemu) OK chacha20 ====== END ======
实现这部分初步体会到了 berrel shift 的妙用,对于一块操作,如
*a = *a + *b;
*d = *d ^ *a;
*d = rotate(*d, 16);
我最开始的设想是,
mov r8, r7, lsr #16
mov r7, r7, lsl #16
eor r7, r7, r8
我觉得已经充分使用了,但是随后又发现其实还可以减少一个周期,变为
mov r8, r7, lsr #16
eor r7, r8, r7, lsl #16
总的来说,是把原本一个 5 cycles 的操作变为了 2 cycles.
经过这部分优化之后,执行时间降低到 16000 附近
====== START ======
crypto_stream_chacha20: 17375 cycles (note that these are meaningless on qemu)
OK chacha20
====== END ======
初步的想法就是将 16 个32比特的状态作为参数传入,然后根据实参来决定执行前4个后后嗣个 quarterround。只要指定好每次传入的实参,那其实就是同样的操作执行4次。不过多余的参数已经被压入栈了,所以每操作完四个参数,再从栈中 pop 出随后要处理的4个参数即可。
但是在执行时出现了 HardFault_Handler, 检查了一下 mps2.ld, 栈大小是 0x400 也就是 1kB, 所以不是栈的问题。
应该是 pop 的问题。是这样的,调用函数时先把参数从又开始压至栈中,再执行 push {r4-r11, lr}, 所以参数的地址实际上在 sp + 9*4 的位置,我通过 gdb 进行 debug 发现的确如此。经过调整后结果正确,不再有问题。
====== START ======
crypto_stream_chacha20: 16800 cycles (note that these are meaningless on qemu)
OK chacha20
====== END ======
也就是把 20 轮循环打包成单个函数,再加一个 flag 和循环跳转。不过这次不再能手动指定寄存器的位置,所以需要一些额外的寄存器来存储地址。因为每一行的第一个值都是 x0-x3, 所以就按照 x0-x15 的顺序传参,这样一来, x0-x3 在一开始存在寄存器中,每次只需要到栈中取一次地址,然后存储在 r10, r11, r14 三个寄存器中。遵循前面的习惯,以 r4-r7 存储四个传入的值,r8 存储计算过程中的中间值。为了使用循环,还需要使用一个寄存器 r9 存储计数器。这一步也比较简单。不过在执行时间上未体现出明显的优势。按理说这里少了许多的 push 和 pop, 应该能节约很多时间。
这里是 32 位的蒙哥马利模乘,输入是两个参数,参数1是普通域上的,参数2是蒙哥马利域上的,都是 32 位的有符号数。根据 Dilithium 的参考实现,其使用的蒙哥马利模约减如下:
int32_t PQCLEAN_DILITHIUM2_CLEAN_montgomery_reduce(int64_t a) {
int32_t t;
t = (int32_t)((uint64_t)a * (uint64_t)QINV);
t = (a - (int64_t)t * Q) >> 32;
return t;
}
即首先将模乘的两项相乘,得到一个 64 位的中间结果,随后乘以 qinv 并取低 32 位得到 t。最后计算 a - t*q 的高 32 位作为结果。蒙哥马利方面, $R = 232 \mod q = -4186625$,
其复杂的点在于需要使用两个寄存器来保存中间状态。比如 a * qinv 这里,其实只需要低 32 位,而结果中的低32位一定是由两个乘数的低32位计算得来的,所以使用普通的 mul 指令即可,将结果保存在一个普通的 32 位寄存器中也就得到了低32位的结果。随后的 t = a-t*q 则可以通过一条指令 smlal 实现,不过这里只只是加法不支持乘法,所以我们实际上计算的 t 是 a * (-qinv), 这样一来这里的加法就变成了乘法,最后取高 32 位即位移位的操作。
还有两个汇编的小 tips:
- 就是返回的值默认放在寄存器 r0 中,之前 chacha20 的部分没有涉及到返回值
- 当把一个 32 比特的立即数放入寄存器中时,需要先 movw lowpart, 再 movt highpart
butterfly 实现单个蝶型运算,输入的参数第一个是一个数组,也就是一个指针,包含两个 32 位的有符号数,第二个是本次蝶型运算使用的 twiddle factor, 同样是一个 32 比特的有符号数。计算就是普通的蝶型运算,首先使用一次模乘,随后是一个加法和一个减法。使用 bl 跳转到有返回值的函数 modmul, 需要注意的是, a[0] 作为随后会用到的值,不要存放到 r0-r3 中,会被覆盖掉。然后 modmul 的参数应该被按照 ABI 放置在 r0 和 r1 中。返回值,即 a1 * tw 则会存放在 r0 中留待 butterfly 使用。
执行第一层的 ntt 操作,使用一个固定的蝶型因子。复杂的地方在于需要使用一个 for 循环。输入是一个指针 a, 该地址之后有 256 个有符号的 32 位整形。每次循环是 a + i 和 a + 128 + i 两个数的蝶型运算。所以需要做两个累加操作,一个是 counter 用于处理循环结束,二是 a0 的地址,需要每次加4以得到下一个系数,此外 a1 的地址即为 a0 的地址加上 128*4 . 还有一个就是我们的模约减是蒙哥马利模乘,所以 tw 应该是蒙哥马利域上的,不要使用 4808194 而是 25847.
ntt 是最复杂的部分,需要实现一个三重循环,只要细心处理,问题倒是不大; 此外,还涉及到访问一个全局变量的问题,可以使用 =varname 的方式取得其地址
kyber 使用的是 ches’23 上提出的定制版本的 plantard reduction, 相较于使用蒙哥马利模乘,所需的乘法减少了一个。如果不考虑 lazy reduction 的话,其实 NTT 实现和蒙哥马利的实现方式区别就仅仅是
assignment 中给出了实现的指令,直接实现即可,将两个常量赋值到对应的寄存器中,最后将结果的高 16 位存储在 r0 中。
注意是半字的操作数,此外 ntt 只有 7 层。