多精度算术：一个娱乐项目

Andrea Simonassi

5.00/5 (2投票s)

2022年6月11日

CPOL

11分钟阅读

6521

121

展示本系列文章中使用的算法和代码

下载源代码 - 73.9 KB

本文是系列文章的一部分

多精度算术：一个有趣的爱好项目

引言

在本文中，我将回顾本系列文章随附的代码。

该仓库（上面的链接，代码快照已附在文章中以备将来参考，因为链接容易失效）包含两个实现：大的是“C”语言版本 + 一些x64汇编代码，小的是一个JavaScript原型，旨在尽可能简单，我甚至将其用作调试工具来验证大实现的正确性。

C代码并非旨在成为可用的API，它只包含核心算法，它是“学习”代码，是娱乐编程的极致，当然可以进行打磨并用作一个功能贫乏但可用的API，但我认为创建一个新的API没有任何意义，因为已经有很多高性能、功能丰富的MP API，例如GMP或类似的，它们是由备受推崇的专业数学家开发的。

这个程序更像是一个“编译脚本”而不是一个真正的程序：编译它并运行！它的输出只取决于执行计算机的速度：每次运行在同一台计算机上都会产生相同的输出。一旦您收集了用于分析的输出数据，您就可以丢弃可执行文件。

要在同一台机器上更改输出：编辑源代码，重新编译并运行……从这个意义上说，它更像是一个脚本而不是一个程序。

这一切的意义在于能够编辑代码并重新编译它，以找到新的输出，从而丰富代码和程序员的个人文化。如果有人想在现实生活中使用部分代码呢？为什么不呢！

程序将诊断消息输出到 stderr，结果输出到 stdout，以制表符分隔文件。

这是R-Studio加载的制表符分隔文件输出

Excel加载的同一文件

程序用法

C 版本

编译代码后，运行程序并将标准输出重定向到txt文件

几分钟（或更长时间……测试并非并行执行）后，您可以使用r-studio、Excel或任何其他工具打开txt文件进行分析。

示例：使用数据透视条形图比较按 operand2_size 进行除法的实现速度

JavaScript 版本

JS实现比“大程序”更“真实”，因为它能接受输入并根据输入产生输出，但当然，性能无法比较。

要打开JS实现，请查看 sources 目录并打开 BigIntSimple/wwwroot/default.html，该页面使用“简朴”一词来形容恰如其分。由于采用UX可用性最佳实践（LOL）创建，其使用方式应不言自明。

这是一个基于10的乘法示例，有趣的是它显示了中间步骤，这在调试大型实现时可能很有用，特别是如果设置为使用 base16，但老实说，用这个工具来乘2K字大小的数字相当困难。

这是一个基于16的除法示例

可移植性问题

程序的核心使用可移植的STD “C” 编写，但有些部分不可移植。

汇编语言

项目中的“C”代码具有一定程度的可移植性，但汇编代码的可移植性较差，如果您没有汇编器或您的汇编器不支持x64 Intel语法，或者您的架构不同，请删除ASM代码：链接器会抱怨缺少代码，您应该找出要删除的代码以使C代码编译。汇编部分是一个“附加组件”，可以安全删除，只是您可能需要重新定义纯“C”代码中的一些小部分。

不可移植的 C 代码

即使某些C代码也不可移植，我添加了一些保护措施，并在可能的情况下添加了一些回退可移植例程，例如，高精度时钟例程就是这样实现的

test.h

#if defined(_WIN32) || defined(WIN32)
   #include <Windows.h>
   #define CLOCK_T unsigned long long
   #undef FAILED
   #undef BitScanReverse
#else
   #warning maybe you want to define a better CLOCK_T for your platform, 
   #also checkout test_common.c
   #define CLOCK_T clock_t
#endif

test_common.c

#if defined(_WIN32) || defined(WIN32)

   CLOCK_T precise_clock() {
       unsigned long long o;
       if (QueryPerformanceCounter((LARGE_INTEGER*) &o) != 0){
             return o;
       }     
       return clock_zero();
   }

   double seconds_from_clock(CLOCK_T clock){
       unsigned long long o;
       if (QueryPerformanceFrequency(
             (LARGE_INTEGER*)&o
       ) == 0)return 0.0;

       double d = (double)(clock) /(double)o;
       assert(d > 0);
       return d;
    }

    CLOCK_T clock_zero() {     
       return 0ULL;
    }

#else

#warning maybe you want to define a better clock function

    CLOCK_T precise_clock() {
       return clock();
    }

    double seconds_from_clock(CLOCK_T clock){
       return clock / CLOCKS_PER_SEC;
    }

    CLOCK_T clock_zero(){
       return 0L;
    }

#endif

编译器应该会发出警告，让您有机会编写不可移植的部分。

项目结构

我使用了MS VC++编译器，所以仓库有一个VC++项目文件。如果您不使用MSVC++，将所有 *.c 文件和 *.asm 文件编译成一个可执行文件应该没问题。

将 Typedef reg_t 对应于目标平台的原生无符号 WORD

在每个平台上，您应该检查一些 typedef，查看 BigIntSimple.h 和 test.h 中的一些平台特定部分，请参见 BigIntSimple.h 的以下摘录。

/* define _R such that appends L to int literal if necessary to keep portability,
but should not be necessary since compiler automatically promotes int to longs  */

#define _R(a) ((reg_t)(a ## L))
typedef uint32_t numsize_t;
typedef uint64_t reg_t;
typedef int _result_t;

/* define that if the target compiler does not have a dword unsigned integer type */
#define NO_DWORD_INTS
typedef uint32_t multiply_small;
typedef uint64_t multiply_big;

以下是 test.h 中一些不可移植的示例部分

#define NOMEM L"NO MEMORY"
#define _fprintf fwprintf 
#define _sprintf  swprintf
typedef wchar_t _char_t;
#define STR(x) L ## x     
#define EXPAND_(a) STR(a)
#define WFILE EXPAND_(__FILE__)  
#define LOG_INFO(x, ...) _fprintf(stderr, L"[INFO] " x L"\n", ##__VA_ARGS__)
#define MY_ASSERT(c, x, ...) if (!(c)) 
{ _fprintf(stderr, L"assertion failed: (%s %d) - " x L"\n",  
  WFILE, __LINE__, ##__VA_ARGS__); abort(); }
#define LOG_ERROR(x, ...) _fprintf(stderr, L"[ERROR] " x L"\n", ##__VA_ARGS__)

其他

抱歉，我没有时间尝试在任何Linux Box上编译这个，这是我的错，我一定会这样做，但我投入到这个项目的时间非常少。

Using the Code

算术挑战

现在我向你挑战，创建你的“C”语言实现，实现四个基本运算中的一个或多个，并尝试让其运行时间快于 ref 实现，我知道做得比我好得多是可能的。

奖品什么都没有，尽情享受吧。

我鼓励强大的志愿者使用任何非便携式/异构硬件特定、并行、GPU、向量、任何微优化技巧来击败 ref 实现。如果能看到内联函数的使用，并比较同一机器上不同编译器的输出，那将是件好事。

我希望看到不同的算法在工作，我也有意创建一个不同的算法并将其添加进来（你可能会从我的提交中注意到，我很少在这方面工作），比如算法的一些并行版本，更好的乘法算法，十进制字符串输出，针对单个数字优化的除法等等。

当然，我们将通过渐近复杂度类别来比较算法。

在添加新的算术算法之前需要了解的重要事项

使用的位置系统

选择数字表示方式极大地影响算法的工作方式和所使用的算术，例如，可以选择使用负基数，或者带有符号数字的正基数，但这对我来说仍然太复杂了。我决定使用正基数和正数字（只实现无符号算术）。只需很少的努力，代码就可以扩展以表示负数。当然，如果底层计算机没有处理寄存器大小无符号数字的指令，这可能是一个问题，但据我所知，所有主要架构都支持无符号算术。

在这个项目中，数字被表示为无符号数字的 string，每个数字的基数为 2^{(8*sizeof(reg_t))}（为了最佳性能，请为您的平台正确考虑 typedef reg_t），最低有效数字位于位置 0。string 的大小总是作为单独的无符号整数值，类型为 numsize_t 提供，最好定义为

typedef uint32_t numsize_t;

但如果您觉得 uint32_t 太小，可以随意使用更大的值。

Zero

零通过传递零大小值来表示，无论A数组的内容如何（在这种情况下A可以为 NULL）。此外，一个大小为K且K位数字都为0的数字等于零，但这需要扫描该数字，因此任何封装此数字系统的API代码都应该规范化零，以便只有一个可能的零表示。

如果您需要围绕此代码构建一个完整的算术API，您可以将数字定义为如下结构体

struct BigUnsignedInteger
{
       numsize_t magnitude;
       reg_t* data;
};
struct BigInteger
{
       char sign;
       numsize_t magnitude;
       reg_t* data;
};
struct BigRational
{
       char sign;
       struct BigUnsignedInteger p;
       struct BigUnsignedInteger q;
};

等等。

在继续实现之前还需要了解的其他重要事项

在源头文件上，有两种操作签名和不同的调用约定。

求和、减法和乘法签名

EXTERN numsize_t YourOperation(reg_t * A, numsize_t ASize, 
       reg_t * B, numsize_t BSize, reg_t* R);

在所有三种情况下，A是一个数字，ASize 是A的大小；B是第二个数字，BSize 是它的大小；R是预分配（由调用者）用于存储结果的空间（数字按照上述“使用的位置系统”段落表示）。

加法/减法调用约定

约定是，对于加法和减法，调用者将R数组作为已分配的数组传递，其容量至少为 MAX(Asize, BSize)+1，当您接收R数组作为输入时，必须将其值视为未定义，在您的实现返回之前，您应该始终将其设置为结果，如果结果为零，则您的实现返回零，并且您可以将R的内容保留为未定义，否则您必须返回有效数字的数量（基数为 2^{(8*sizeof(reg_t))}）并相应地设置R值。

在减法中，如果B大于A，则结果是未定义的，您不必检查A>=B，那是调用者的责任，您需要在这里尽可能快：这是此函数可以包装到其他函数中以构建最终API的原因之一。

实现并非必需，但某些加法和减法算法可能能够在原地工作，因此例如，R可以指向与A或B相同的数组；如果您打算将您的加法/减法例程用作乘法或除法的子例程以节省额外的缓冲区，这可能会很有用，但这并非必需。

乘法调用约定

在乘法中，调用者必须传递已分配的 R 数组，其容量为 ASize + Bsize。当您接收 R 作为输入时，您必须将其值视为未定义，如果您的实现需要将 R 填充为零，则您的实现必须完成此操作。

除法签名

_div_result_t LongDivision(reg_t *A, numsize_t m, reg_t *B, numsize_t n, reg_t * Q, 
                           numsize_t * q, reg_t * R, numsize_t * r);

这个签名稍微复杂一些，因为我们需要标记 DIVIDE_BY_0 或其他可能的错误情况：加法、减法和乘法不需要分配任何额外的空间（嗯，乘法有时需要分配，但可以回退到其他不需要额外缓冲区空间的较慢算法）。

所以我们有A及其大小m，B及其大小n，Q是用于存储结果的空间和q，用于存储输出q大小的位置，R是用于存储余数的空间，r是用于存储余数大小的位置。

调用者必须保证A的容量为 m+1，B的容量为 n，这样你的实现就可以自由地将B左移一些位，直到它的第一个左位变为1，因为归一化B可以让你的算法更好地工作，你的实现必须尽可能快，所以你不想分配东西或使你的算法复杂化，最好把这个责任交给调用者，当然，一个假设的API应该以一种对用户透明的方式封装数字。

调用者保证 Q 足够大以容纳结果，即至少 m-n+2 位数字的空间（2 是因为可能的移位）。

调用者还保证R足够大以包含余数，因此它必须与n位数字一样大。

创建您自己的算术算法

以示例为例，让我们创建一个新的乘法算法示例：同样的过程适用于加法、减法和除法。

编辑 BigIntSimple.h 文件，并为您超快的乘法算法添加签名。

numsize_t MySuperFastMultiplication
    (reg_t* A, numsize_t m, reg_t* B, numsize_t n, reg_t* R);

然后向项目中添加一个新 .c 文件，您将在其中实现乘法，将其添加到项目中，然后添加以下代码

#include "BigIntSimple.h"

numsize_t MySuperFastMultiplication
    (reg_t* A, numsize_t m, reg_t* B, numsize_t n, reg_t* R)
{
       /*your superfast algorithm here, we will implement it later, 
         we can compile it as is, of course tests are going to fail.. */
       return 0;
}

编辑 test_config.c 文件，并将您的算法插入乘法算法列表（搜索乘法初始化部分）

static _operationdescriptor multiplications [] =
{
       {
             STR("C ref"),
             OP, /* op for sum, mul and sub, DIVOP for div */
             LongMultiplicationV2,
             {0,0,NULL} /* statistics (0 initial size, 0 count, NULL pointer) */
       },
       {
             STR("My super fast multiplication"), /*INSERT YOUR ALGORITHM AFTER 
                                                    THE FIRST REFERENCE ALGORITHM 
                                                    IN ANY POSITION*/
             OP, /* op for sum, mul and sub, DIVOP for div */
             MySuperFastMultiplication,
             {0,0,NULL} /* statistics (0 initial size, 0 count, NULL pointer) */
       },
       /* etcetera */
 };

就是这样！

现在编译并运行程序。

以下是 STDERR 输出的摘录，当然，测试失败了，因为我们总是对我们的超快速乘法返回0（但它通过了A*0必须等于0的测试，并且它还具有交换律和结合律。😊）

这是一个Excel图表的样本

哇，“我的超快速乘法”真的很快，对于32K字的运算数，它只用了0秒，而其他算法用了8秒……但问题是“我的超快速算法”并没有工作……它还需要实现。

调试您的实现

您可以随心所欲地调试它，但我建议这样做：在 test_config.c 中注释掉所有其他乘法算法，只保留您自己的和参考算法；在 main 函数中，注释掉所有测试用例，只保留您正在处理的操作的测试用例。

int main(){
       setlocale(LC_ALL, "");
       initTest();        

       /*testParse();
       testBSR();
       testCompare();
       testSum();
       testSub();*/
       testMul(); 
       /*testDiv();*/
       write_summary();
       return 0;
}

添加新单元测试

如果您需要额外的单元测试，您可以添加任意数量的测试，请遵循git仓库上的手册获取建议，或者随心所欲地做：https://github.com/andreasimonassi/bigint/blob/master/BigIntSimple/creating_unit_tests_readme.md

历史

2020年5月：首次发布
2022 年 6 月 11 日：在线发布

多精度算术： 一个娱乐项目