2
\$\begingroup\$

This post is a second part of my original post, Add two huge base-10 numbers, which deals with adding two huge base-10 numbers. However, in this case, I'm multiplying two non-negative whole numbers. In both cases, the input and output are represented as null-terminated ASCII strings, with digits in big-endian order. I use my own custom strlen function which only modifies rcx. I've tried to incorporate most of the feedback I've received and managed to shorten the code from 170 to 130 lines. Is there anything else I can do which will make it run faster?

Here's the code:

extern add_whole
extern strlen

section .text
global _multiply_whole
_multiply_whole:
  ; Input:
  ;   - char *a    -> rdi
  ;   - char *b    -> rsi
  ;   - char *res  -> rdx  (allocate strlen(a) + strlen(b) + 1 (for null terminator) bytes)
  ;   - char *buf1 -> rcx  (temporary buffer for program to use)
  ;   - char *buf2 -> r8   (both buf1 and buf2 have the same length as res)

  ; Registers used:
  ;   - rax
  ;   - rbx
  ;   - r9
  ;   - r10
  ;   - r11
  ;   - r12
  ;   - r13
  ;   - r14
  ;   - r15
  
  push   rbx                  ; Push used callee-saved registers.
  push   r12
  push   r13
  push   r14
  push   r15
  mov    r9, rdi               ; Move char* a to r9.
  mov    r10, rcx              ; Save rcx since strlen() doesn't preserve it.
  call   strlen                ; rax = strlen(a)
  mov    rbx, rax              ; Save rax in order to prevent it from being overwritten.
  mov    rdi, rsi              ; strlen's first argument is now b.
  call   strlen                ; Call strlen with b.
  dec    rax
  xchg   rax, rbx              ; rbx = strlen(b) - 1
  mov    rdi, r9               ; Restore rdi.
  mov    rcx, r10              ; Restore rcx.
  xor    r10d, r10d            ; r10b is the carry.
  mov    byte [rdx], '0'       ; Set the default value of the result string to '0'.
  xor    r14d, r14d            ; The number of leading zeroes to add.
.loop_1:
  lea    r13, [rax+1]          ; r13 = strlen(a) + 1
  add    rcx, rax              ; rcx = &a[strlen(a)]
  mov    r9, rax               ; r9 = strlen(a)
  mov    r15d, 10              ; The number we're dividing by.
  push   rax                   ; Save rax since we use it for mul and div later.
.loop_2:
  movzx  r11d, byte [rsi+rbx]  ; r11b is a byte of char *b
  sub    r11b, '0'
  movzx  eax, byte [rdi+r9-1]  ; al is a byte of char *a
  sub    al, '0'
  mul    r11b
  mov    r11b, al              ; r11b is the multiplication of char *a and char *b.
  add    r11b, r10b            ; Add carry.
  push   rdx
  mov    al, r11b
  div    r15b
  mov    r10b, al              ; carry = r11b / 10
  mov    al, ah
  mov    byte [rcx], al
  add    byte [rcx], '0'       ; Add the remainder to the result string.
  pop    rdx
  dec    rcx
  dec    r9
  jnz    .loop_2               ; Loop r9 times.
  mov    byte [rcx], r10b      ; Add the final carry (might be 0)
  add    byte [rcx], '0'       ; to the string.
  xor    r10b, r10b            ; Reset carry.
  xor    r11d, r11d            ; r11 = 0
  test   r14, r14
  jle    .after_loop_3         ; if r14 <= 0 then skip loop_3
.loop_3:                       ; Adds r14 trailing zeroes to buf1, r14 is the amount of times loop_1 has run.
  mov    byte [rcx+r13], '0'   ; buf1[r13] = '0'
  inc    r13                   ; r13++
  inc    r11                   ; r11++
  cmp    r11, r14
  js     .loop_3               ; Loop r14 times.
.after_loop_3:
  mov    byte [rcx+r13], 0     ; Zero-terminate the string to pass it to add_whole().
  push   rcx                   ; Save used caller-saved registers.
  push   rdi
  push   rsi
  push   rdx
  mov    rdi, rdx              ; param1 = rdx (char *res)
  mov    rsi, rcx              ; param2 = rcx (char *buf1)
  mov    rdx, r8               ; param3 = r8 (char *buf2)
  mov    r15, r8               ; Saving the r8 register.
  call   add_whole             ; buf2 = add_whole(res, buf1) (this doesn't actually return anything)
  mov    r8, r15               ; Restoring the r8 register.
  pop    rdx
  pop    rsi
  xor    r11d, r11d
  mov    rdi, r8
  call   strlen                ; r9 = strlen(buf2) (length of the addition result)
  mov    r9, rax
  pop    rdi
  pop    rcx
  pop    rax
.loop_4:
  mov    r15b, byte [r8+r11]
  mov    byte [rdx+r11], r15b
  inc    r11
  cmp    r11, r9
  js     .loop_4
  mov    byte [rdx+r9], 0
  xor    r11d, r11d
.loop_5:
  mov    byte [rcx+r11], 0
  mov    byte [r8+r11], 0
  inc    r11
  cmp    r11, r13
  js     .loop_5
  inc    r14
  mov    r15, rbx
  dec    rbx
  test   r15, r15
  jg     .loop_1
  mov    rdi, rdx
  call   strlen
  mov    byte [rdx+rax], 0
  pop    r15
  pop    r14
  pop    r13
  pop    r12
  pop    rbx
  ret

This is basically the pen-and-paper multiplication algorithm implemented in assembly. I've used movzx in some cases and even avoided pushing and popping registers (instead, I've used two movs as suggested by https://stackoverflow.com/questions/73996728/is-moving-into-another-register-faster-or-slower-than-push-and-pop.)

Currently, this code can calculate 1000! 27 times per second. I used this code for that benchmark.

#include "multiply.h"
#include <chrono>
#include <iostream>

std::string factorial(int x) {
  std::string answer = "1";
  for (auto i = 2; i <= x; i++) {
    std::string toMultiply = std::to_string(i);
    char *partAns =
        (char *)calloc(answer.length() + toMultiply.length() + 1, 1);
    multiply(answer.c_str(), toMultiply.c_str(), partAns);
    answer = partAns;
    free(partAns);
  }

  return answer;
}

int main() {
  double time = 0;
  int x = 0;
  while (time < 1) {
    auto start = std::chrono::high_resolution_clock::now();
    std::string s = factorial(1000);
    auto end = std::chrono::high_resolution_clock::now();
    time += std::chrono::duration_cast<std::chrono::microseconds>(end - start)
                .count() *
            1e-6;
    x++;
  }

  std::cout << "1000!'s per second: " << x << '\n';
}

multiply.h:

/// @brief Multiplies the first rational argument with the second rational
/// argument and stores the result in the third argument.
/// @param a The first non-negative rational number as a decimal.
/// @param b The first non-negative rational number as a decimal.
/// @param res Where a * b will be stored.
void multiply(const char *a, const char *b, char *res) {
  extern size_t strlen(const char *str);
  size_t a_length = strlen(a);
  size_t b_length = strlen(b);
  char *a_copy = (char *)calloc(a_length + 1, 1);
  char *b_copy = (char *)calloc(b_length + 1, 1);
  size_t ptr = 0;
  int displacement = 0;
  unsigned char flag = 0;
  for (size_t i = 0; i < a_length; i++) {
    if (flag == 1)
      displacement++;
    if (a[i] == '.')
      flag = 1;
    else {
      a_copy[ptr] = a[i];
      ptr++;
    }
  }
  ptr = 0;
  flag = 0;
  for (size_t i = 0; i < b_length; i++) {
    if (flag == 1)
      displacement++;
    if (b[i] == '.')
      flag = 1;
    else {
      b_copy[ptr] = b[i];
      ptr++;
    }
  }

  size_t bufsize = strlen(a) + strlen(b) + 1;
  char *buf1 = (char *)calloc(bufsize, 1);
  char *buf2 = (char *)calloc(bufsize, 1);
  _multiply_whole(a_copy, b_copy, res, buf1, buf2);

  size_t reslength = strlen(res);
  for (size_t i = 0; i < displacement; i++)
    res[reslength - i] = res[reslength - i - 1];
  if (displacement)
    res[reslength - displacement] = '.';

  free(buf1);
  free(buf2);
  free(a_copy);
  free(b_copy);
}

Any way to improve this?

\$\endgroup\$
2
  • \$\begingroup\$ Will your fourth question be about transcendental functions base ten? \$\endgroup\$
    – greybeard
    Oct 15, 2022 at 5:39
  • \$\begingroup\$ @greybeard I don't know; currently I haven't made any transcendental functions, but if I do in the future, then I might post a question on it. \$\endgroup\$
    – avighnac
    Oct 15, 2022 at 5:42

1 Answer 1

3
\$\begingroup\$

Unfortunate compare/branch pair

A cmp/js pair cannot macro-fuse on various Intel processors. Since it's meant to be a less-than comparison, usually it would be done with cmp/jl or cmp/jb (depending on signedness) both of which can macro-fuse. Also, cmp/js in general can go wrong in case of overflow. That case doesn't look very relevant here in this code, I just mention it so you can avoid it when it does matter.


avoided pushing and popping registers

Not entirely, but I'll leave that up to you.


div and 8-bit registers

Note that the dividend for div r15b is ax, not al. ah is zero at the time of the division, but seemingly by accident, and in such a way that it creates a situation where the high-byte register and the low-byte register may have gotten renamed separately and need to be "unified" with an implicit extra µop. That would only happen on microarchitectures that rename those registers separately. On others, writing to byte registers may have an implicit dependency on the corresponding 64-bit register, to "merge" the new value into it. movzx has been used in some cases, but there are still other operations that risk triggering the curse of partial register writes.

The curse of partial register writes also applies to xor r10b, r10b for example, which (unlike xor r10d, r10d) is not a special zeroing idiom (so it is a real xor instruction, with associated execution cost, and not a dependency-breaker). Unless there is some overriding reason not to, you should zero the corresponding 32-bit register. Using xor r11d, r11d to set r11 to zero is good, keep doing that.

By the way the div can be avoided, which is probably still a good idea, despite significant improvements to the performance of integer division in recent x86 processors. There are various ways to approach this.

For example for a number in range 0..91 held in eax, division by 10 may be implemented like this:

imul edx, eax, 0CDh   ; approximately 1/10 * 2^11
shr edx, 11

which works for that entire range of numbers and actually more, the first time it goes wrong is when eax = 0x0405.

Then the remainder can be determined by subtracting 10 times the quotient from the dividend,

imul ecx, edx, -10
add eax, ecx          ; this add may be merged with the conversion to ASCII into a 3-part lea

Or two leas could be used to implement a multiplication by 10, but then there would be a sub which cannot be merged with the conversion to ASCII.

I used some relatively arbitrary registers for these examples, I do not mean to imply that these snippets can be plugged into the code as-is. Please take them as descriptions of techniques rather than as concrete suggestions.


mov byte [rcx], al
add byte [rcx], '0'

A read-modify-write operation is relatively "heavy" in various ways. Not super heavy, but heavy enough to think about avoiding it. It could be avoided by using lea to do a non-destructive addition, or even by adding '0' to eax (note that I'm choosing eax instead of al on purpose, to avoid the wrath of partial register writes), storing al, and then subtracting '0' again. That add and sub would be part of a loop-carried dependency chain, but that does not look significant to me compared to the rest of the loop body.


Sweeping changes

One of the biggest things you could do for performance, is multiplying several digits at once. The basic multiplication algorithm does a quadratic number of multiplications after all. Switching from base 10 to base 109 (so that each limb fits in 32 bits, and each product fits in 64 bits) would not reduce multiplications by a factor of 9, but rather a factor of 81[1]. I chose a power of ten for the base, since that most easily converts to and from the base ten input and output. Even with that, it is not so simple to do efficiently, or even to do it at all. The temporary results into which the additions take place, should remain in base 109 for as long as possible, since converting that to base ten is (while simpler than general base conversions, since each limb would map nicely to 9 output digits) actually not cheap.

You could take this one step further and take advantage of 64bit-to-128bit multiplication to do it all in base 1018 for another 4x reduction in multiplications, but the resulting 128bit partial products are more complicated to deal with. I don't think you should do it.

There are various more complex algorithms for big-integer multiplication that are worth using when numbers have about the size that you're already dealing with (over a thousand decimal digits), but they apply only when both inputs are sufficiently big, which is not the case when calculating a factorial.


note 1: but that factor is lower when calculating a factorial, since one of the multiplicands is always small in that case.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.