Arb 任意精度计算库

贡献者： addis

预备知识 1　C++ 基础

未完成：根据 arb 的主页，arb 在 2023 年已经融入到 flint 的开发中，官网和文档不再更新。

　　 Arb 是一个支持任意精度计算的 C/C++ 库（也提供 python 接口），支持对计算进行严谨的误差估计，即每个任意精度浮点数 $z$ 都会伴随一个误差半径 $r$，使得精确结果必定落在复平面上以 $z$ 为圆心半径为 $r$ 的圆盘中。通过增加浮点数的精度，就可以用数值方法无限逼近理论值。Arb 还提供了许多特殊函数的计算，例如 $\Gamma$ 函数，超几何函数等，以及常用的线性代数功能和离散傅里叶变换等。官方主页 arblib.org 包含详细的文档。

1. 安装

　　以 Ubuntu 为例，最方便的安装方式就是使用 apt 安装。首先安装 dependency

sudo apt install libmpfr-dev libflint-dev

其中 MPFR 和 flint 两个包分别用于任意精度浮点数以及数论。然后安装 Arb

sudo apt install libflint-arb-dev libflint-arb2

　　但目前这并不是 Arb 的最新版本（例如没有实现库仑函数）。要获得最新版本，可以直接从 GitHub 下载源码编译即可（默认使用 gcc 编译器）。目前笔者使用的版本是 release 2.19.0（ubuntu 20.04 和 22.04 均可直接用 apt 安装更高版本）。

　　可以用 ./configure --help 查看编译选项，若所有的包都安装在默认目录则不需要编译选项。

./configure [编译选项];
make -j4;
sudo make install;

其中 -j4 是使用 4 线程进行编译，也可以改成其他数字。

　　要在 windows 上安装，参考这个。这里面的 dll 是可以直接拿来用的，无需编译。

2. 编译

　　在 Ubuntu 中如果你使用 apt 安装，在编译程序是需要加 -lflint-arb 选项。如果你直接从源码编译，则需要加 -larb 选项。对一些编译器（例如 intel 的 icpc），可能还需要加上 -lflint, -lmpfr, -lgmp，如果 link 阶段遇到问题可以试试加上。

　　一个 C++ 例程：该程序用 80 bit 初始精度计算复参数的超几何函数 $_1F_1$ 并于 Mathematica 的结果比较。如果 Arb 估计的结果精度小于 16 位有效数字则显示警告。初始精度越高，结果的有效数字也越高，具体取决于 $_1F_1$ 的参数。读者可以尝试在程序内加入一个循环，若结果有效数字不够，则把初始精度加倍再次计算，直到达到精度为止。

代码 1：test.cpp

#include "arb_hypgeom.h"
#include "acb_hypgeom.h"
#include <complex>
#include <iostream>

using namespace std;

typedef double Doub, Doub_O;
typedef const double Doub_I;
typedef complex<double> Comp, Comp_O;
typedef complex<double> Comp_I;
typedef int Int;
typedef const int Int_I;
typedef int & Int_O, Int_IO;

// 1F1 hypdergeometric function with complex arguments
Comp hypergeom1F1(Comp_I a, Comp_I b, Comp_I z)
{
	slong prec = 80; // set precision bit (log10/log2 = 3.322)
	Comp res;
	arb_t temp1;
	arb_init(temp1);
	acb_t a1, b1, z1, res1;
	acb_init(a1); acb_init(b1); acb_init(z1); acb_init(res1);
	// error range is set to 0
	acb_set_d_d(a1, real(a), imag(a));
	acb_set_d_d(b1, real(b), imag(b));
	acb_set_d_d(z1, real(z), imag(z));
	acb_hypgeom_1f1(res1, a1, b1, z1, 0, prec);
	// acb_printn(res1, 100, 0); printf("\n"); // print result
	int digits = acb_rel_accuracy_bits(res1)/3.321928;
	if (digits < 16) {
		cout << "warning: hypergeom1F1 error too large" << endl;
	}
	acb_get_real(temp1, res1);
	res.real(arf_get_d(arb_midref(temp1), ARF_RND_NEAR));
	acb_get_imag(temp1, res1);
	res.imag(arf_get_d(arb_midref(temp1), ARF_RND_NEAR));
	acb_clear(a1); acb_clear(b1); acb_clear(z1);
    acb_clear(res1); arb_clear(temp1);
	return res;
}

int main()
{
	cout << "hypergeom1F1(1.23+1.23I, 1.23+1.23I, 1.23+1.23I) = " << endl;
	cout << hypergeom1F1(Comp(1.23,1.23),Comp(1.23,1.23),Comp(1.23,1.23))
    cout << endl;
	printf("Mathematica: 1.143503984180676 + 3.224470526790991i\n");
}

编译：

g++ -o test.x test.cpp -larb -lflint

　　确认版本：arb.h 头文件中定义了版本宏

#define __ARB_VERSION 2
#define __ARB_VERSION_MINOR 23

常用函数和常数

　　数学常数（prec 是结果的二进制精度，跟上面的 hypergeom1F1 不同）：

void arb_const_pi(arb_t z, slong prec) 圆周率
void arb_const_sqrt_pi(arb_t z, slong prec) 根号圆周率
arb_const_log2
void arb_const_e(arb_t z, slong prec)

3. 内部实现

预备知识 2　FLINT 库笔记

要配合 gdb 学习可以自己编译一个 debug 版本的 libarb 库，见 Automake 笔记。
ARF_PREC_EXACT 的定义是 std::numeric_limits<slong>::max()

arf_t 的数据结构

arf_t 是任意精度浮点数，文档。
arf_t 的数据结构，看懂下面的代码就了解了。基本就是小数部分 .d（例如 0.100101110，小数点后面用一个 flint 大整数表示），指数部分 .exp（一个 flint 大整数），和 .size 部分（最后一 bit 是符号位，1 表示复数，.size >> 1 是用到的 limb 的个数，这可能和 alloc 的 limb 的个数是不一样的）。
当 .size == 0 时，arf_t 是某个 special 值，即：.exp == ARF_EXP_ZERO 时，arf 为零；.exp == ARF_EXP_POS_INF 时，arf 为正无穷；.exp == ARF_EXP_NEG_INF 时，arf 为负无穷；.exp == ARF_EXP_NAN 是，arf 为 nan；如果不是 special，就叫做 normal。
相关函数：arf_is_special，arf_is_zero，arf_is_pos_inf，arf_is_neg_inf，arf_is_nan，arf_is_normal，arf_is_finite，arf_is_inf，arf_is_one

#define ARF_NOPTR_LIMBS 2
#define ARF_NOPTR_D(x) ((x)->d.noptr.d) // 小数部分的 limb 指针（非 alloc）
#define ARF_PTR_D(x) ((x)->d.ptr.d) // 小数部分的 limb 指针（alloc）

typedef struct
{ // mp_limb_t 是 GMP 整数的一个 limb
    mp_limb_t d[ARF_NOPTR_LIMBS];
}
mantissa_noptr_struct; // arf 的小数部分至少有两个 limb

typedef struct
{
    mp_size_t alloc; // GMP limb 的 alloc 的个数， 并不一定全部用到。
    mp_ptr d; // mp_limb_t （GMP 整数的一个 limb）的指针
}
mantissa_ptr_struct;

typedef union
{
    mantissa_noptr_struct noptr;
    mantissa_ptr_struct ptr;
}
mantissa_struct;

typedef struct
{
    fmpz exp; // 指数（fmpz 是 flint 的任意精度整数）
    mp_size_t size; // 最后 1bit 是符号， >>1 是小数部分 limb 的数量
    mantissa_struct d; // 小数部分
	// 当 size <= ARF_NOPTR_LIMBS 时使用 d.noptr， 否则使用 d.ptr
}
arf_struct;

typedef arf_struct arf_t[1];
typedef arf_struct * arf_ptr;
typedef const arf_struct * arf_srcptr;

/* Raw size field (encodes both limb size and sign). */
#define ARF_XSIZE(x) ((x)->size)

// 获取用到的 limb 的个数
#define ARF_SIZE(x) (ARF_XSIZE(x) >> 1)
// 获取符号位， 1 为负数
#define ARF_SGNBIT(x) (ARF_XSIZE(x) & 1)

// x 是 arf_t， 获取其小数部分的 limb 指针和 limb 个数
// 当 limb 个数小于 ARF_NOPTR_LIMBS = 2 时， 不存在额外 alloc 的空间。
#define ARF_GET_MPN_READONLY(xptr, xn, x)   \
    do {                                    \
        xn = ARF_SIZE(x);                   \
        if (xn <= ARF_NOPTR_LIMBS)          \
            xptr = ARF_NOPTR_D(x);          \
        else                                \
            xptr = ARF_PTR_D(x);            \
    } while (0)


// get a quad precision number from arf_t type
// similar to arf_get_d()
inline void arf_get_q(Qdoub_O v, const arf_t x, arf_rnd_t rnd)
{
	arf_t t;
	// mp_limb_t is the type of a limb, with FLINT_BITS bits
	// typedef const mp_limb_t *mp_srcptr;
	mp_srcptr tp; // pointer to least significant limb
	mp_size_t tn; // # of limbs

	arf_init(t);
	arf_set_round(t, x, 113, rnd);
	ARF_GET_MPN_READONLY(tp, tn, t);
	if (tn == 1)
		v = (Qdoub)(tp[0]);
	else if (tn == 2)
	    // FLINT_BITS 是每个 limb 的 bit 数， 即 sizeof(mp_limb_t)*8
		v = (Qdoub)(tp[1]) + (Qdoub)(tp[0]) * ldexpq(1,-FLINT_BITS);
	else if (tn == 3)
		v = (Qdoub)(tp[2]) + (Qdoub)(tp[1]) * ldexpq(1,-FLINT_BITS)
		    + (Qdoub)(tp[0]) * ldexpq(1,-2*FLINT_BITS);
	else if (tn == 4)
		v = (Qdoub)(tp[3]) + (Qdoub)(tp[2]) * ldexpq(1,-FLINT_BITS)
		    + (Qdoub)(tp[1]) * ldexpq(1,-2*FLINT_BITS)
			+ (Qdoub)(tp[0]) * ldexpq(1,-3*FLINT_BITS);
	else
		SLS_ERR("not implemented!");

	v = ldexpq(v, ARF_EXP(t) - FLINT_BITS);

	if (ARF_SGNBIT(t)) // 1 for negative
		v = -v;
	arf_clear(t);
}

mag_t 的数据结构

typedef struct
{
    fmpz exp; // 指数部分， flint 整数
    mp_limb_t man; // 小数部分， 单个 limb
}
mag_struct; // 用于表示误差半径， 非负浮点数

typedef mag_struct mag_t[1];
typedef mag_struct * mag_ptr;
typedef const mag_struct * mag_srcptr;

MAG_INLINE int
mag_is_inf(const mag_t x)
{
    return (MAG_MAN(x) == 0) && (MAG_EXP(x) != 0);
}

　　它表示的值为 x->man * 2^ (x->exp - MAG_BITS) 其中 #define MAG_BITS 30。

arb_t 的数据结构

　　文档。

typedef struct
{
    arf_struct mid; // 区间中点
    mag_struct rad; // 误差半径
}
arb_struct;

typedef arb_struct arb_t[1];
typedef arb_struct * arb_ptr;
typedef const arb_struct * arb_srcptr;

acb_t 数据结构

　　就是两个 arb 类型

typedef struct
{
    arb_struct real;
    arb_struct imag;
}
acb_struct;

typedef acb_struct acb_t[1];
typedef acb_struct * acb_ptr;
typedef const acb_struct * acb_srcptr;

致读者：小时百科一直以来坚持所有内容免费无广告，这导致我们处于严重的亏损状态。长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。因此，我们请求广大读者热心打赏 ，使网站得以健康发展。如果看到这条信息的每位读者能慷慨打赏 20 元，我们一周就能脱离亏损，并在接下来的一年里向所有读者继续免费提供优质内容。但遗憾的是只有不到 1% 的读者愿意捐款，他们的付出帮助了 99% 的读者免费获取知识，我们在此表示感谢。