就地整数乘法

Question

我正在编写一个用C语言编写的程序，旨在尽可能短的时间计算大数的幂。我将数字表示为数字向量，因此所有操作都必须手动编写。

如果没有所有中间结果的分配和释放，程序速度将快得多。有没有任何算法可以“原地”执行整数乘法？例如以下函数：

void BigInt_Times(BigInt *a, const BigInt *b);

要在不使用中间值的情况下，将 a 和 b 的乘积结果储存在 a 中。

- Paul Manta

---

那么问题是什么？你到目前为止尝试了什么？ - Macmade

@Macmade：有没有一种在原地进行整数乘法的算法？ - zw324

@Macmade 这个问题是：“是否有一种算法可以进行整数乘法，且是原地操作的？”；到目前为止，我已经编写了一个执行整数乘法的函数，但不是原地操作的；还有一个计算幂的函数（具有log(f)复杂度，其中f是乘法函数的复杂度）。 - Paul Manta

你写这篇文章的原因是什么？为什么不使用已经为C语言编写的众多现成的高性能、经过测试、有文档、可立即工作的大数库呢？ - Dour High Arch

@DourHighArch 学习项目。 - Paul Manta

这是一道相对容易回答的问题，但是我需要一个白板才能解答它... - Mysticial

Answer 1

标准算法包括将“a”的每个数字（单词）与“b”的每个数字相乘，并将它们总和到结果的适当位置。 因此，“a”的第i位数字进入结果的第i到i+n个数字中的每个数字。 因此，为了在原地执行此操作，您需要从最高有效位向下计算输出数字。 这比从最低位到最高位进行计算要棘手一些，但并不是很难...

- Chris Dodd

Answer 2

这里的muln()是无符号整数的2n（实际上是n）乘以n = 2n原地乘法。您可以将其调整为使用32位或64位“数字”而不是8位进行操作。模运算符保留以便清晰明了。

muln2()是n乘以n = n的原地乘法（如此处所示），也在8位“数字”上进行操作。

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <limits.h>

typedef unsigned char uint8;
typedef unsigned short uint16;
#if UINT_MAX >= 0xFFFFFFFF
typedef unsigned uint32;
#else
typedef unsigned long uint32;
#endif
typedef unsigned uint;

void muln(uint8* dst/* n bytes + n extra bytes for product */,
          const uint8* src/* n bytes */,
          uint n)
{
  uint c1, c2;

  memset(dst + n, 0, n);

  for (c1 = 0; c1 < n; c1++)
  {
    uint8 carry = 0;

    for (c2 = 0; c2 < n; c2++)
    {
      uint16 p = dst[c1] * src[c2] + carry + dst[(c1 + n + c2) % (2 * n)];
      dst[(c1 + n + c2) % (2 * n)] = (uint8)(p & 0xFF);
      carry = (uint8)(p >> 8);
    }

    dst[c1] = carry;
  }

  for (c1 = 0; c1 < n; c1++)
  {
    uint8 t = dst[c1];
    dst[c1] = dst[n + c1];
    dst[n + c1] = t;
  }
}

void muln2(uint8* dst/* n bytes */,
           const uint8* src/* n bytes */,
           uint n)
{
  uint c1, c2;

  if (n >= 0xFFFF) abort();

  for (c1 = n - 1; c1 != ~0u; c1--)
  {
    uint16 s = 0;
    uint32 p = 0; // p must be able to store ceil(log2(n))+2*8 bits

    for (c2 = c1; c2 != ~0u; c2--)
    {
      p += dst[c2] * src[c1 - c2];
    }

    dst[c1] = (uint8)(p & 0xFF);

    for (c2 = c1 + 1; c2 < n; c2++)
    {
      p >>= 8;
      s += dst[c2] + (uint8)(p & 0xFF);
      dst[c2] = (uint8)(s & 0xFF);
      s >>= 8;
    }
  }
}

int main(void)
{
  uint8 a[4] = { 0xFF, 0xFF, 0x00, 0x00 };
  uint8 b[2] = { 0xFF, 0xFF };

  printf("0x%02X%02X * 0x%02X%02X = ", a[1], a[0], b[1], b[0]);
  muln(a, b, 2);
  printf("0x%02X%02X%02X%02X\n", a[3], a[2], a[1], a[0]);

  a[0] = -2; a[1] = -1;
  b[0] = -3; b[1] = -1;
  printf("0x%02X%02X * 0x%02X%02X = ", a[1], a[0], b[1], b[0]);
  muln2(a, b, 2);
  printf("0x%02X%02X\n", a[1], a[0]);

  return 0;
}

输出：

0xFFFF * 0xFFFF = 0xFFFE0001
0xFFFE * 0xFFFD = 0x0006

我认为这是我们能够就地完成的最佳方案。我不喜欢muln2()的一件事情是它必须累加更大的中间乘积，然后传播更大的进位。

- Alexey Frunze

Answer 3

听起来你并不需要一个算法，相反，你需要更好地利用语言的特性。

为什么不直接创建你在答案中提到的那个函数呢？使用它并享受吧！（该函数很可能最终会返回对 a 的引用作为其结果。）

- John Fisher

Answer 4

通常，大整数的表示根据所表示的值而变化长度；一般来说，结果会比任何一个操作数都要长。 特别是，对于乘法，结果表示的大小大致等于两个参数的大小之和。

如果您确定内存管理真正成为您特定平台的瓶颈，您可以考虑实施一个更新第三个值的乘法函数。就您上面的C风格函数原型而言：

void BigInt_Times_Update(const BigInt* a, const BigInt* b, BigInt* target);

这样，您可以像C++ std::vector<>容器一样处理内存管理：只有在现有大小太小时，您的更新目标才需要重新分配其堆数据。

- comingstorm

---

我已经有类似于你建议的函数（除了它返回值）。此外，我有一个限制，只能使用C语言特性。 - Paul Manta

我并不是想要“使用”std::vector<>。我的意思是，如果你有更新BigInt的函数（而不是返回一个新的BigInt），只要你不需要更多的空间，就可以保留堆分配的内存，而不是重新分配它。这就是std::vector<>进行动态分配的方式；我建议你遵循它的例子... - comingstorm

我已经编辑了我的回答，请看一下是否更清晰。 - comingstorm