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;