numth.h 文件参考

#include <stdbool.h>
#include "lmmp.h"

numth.h 的引用(Include)关系图:

此图展示该文件直接或间接的被哪些文件引用了:

浏览源代码.

类型定义
typedef int8_t	schar
typedef int8_t *	scharp
typedef int32_t	sint
typedef int32_t *	sintp
typedef int64_t	slong
typedef int64_t *	slongp
typedef int16_t	sshort
typedef int16_t *	sshortp
typedef uint8_t	uchar
	Copyright (C) 2026 HJimmyK(Jericho Knox)
typedef uint8_t *	ucharp
typedef uint32_t	uint
typedef uint32_t *	uintp
typedef uint64_t	ulong
typedef uint64_t *	ulongp
typedef uint16_t	ushort
typedef uint16_t *	ushortp

函数
mp_size_t	lmmp_2factorial_ (mp_ptr dst, mp_bitcnt_t bits, mp_size_t rn, uint n)
	计算 n!! 双阶乘
mp_size_t	lmmp_2factorial_size_ (uint n, mp_bitcnt_t *bits)
	计算 n!! 双阶乘的 limb 缓冲区长度
mp_size_t	lmmp_arith_seqprod_ (mp_ptr dst, mp_size_t rn, uint x, uint n, uint m)
	计算等差数列乘积 x(x+m)...(x+n*m)
mp_size_t	lmmp_arith_seqprod_size_ (uint x, uint n, uint m)
	计算等差数列乘积 x(x+m)...(x+n*m)的 limb 缓冲区长度
void	lmmp_binvert_ (mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_size_t n)
	计算 [numa,na] 在 B^n 下的逆元
void	lmmp_binvert_2_ (mp_ptr dst, mp_srcptr numa)
	计算 [numa,2] 在 B^2 下的逆元
void	lmmp_binvert_3_ (mp_ptr dst, mp_srcptr numa)
	计算 [numa,3] 在 B^3 下的逆元
void	lmmp_binvert_4_ (mp_ptr dst, mp_srcptr numa)
	计算 [numa,4] 在 B^4 下的逆元
void	lmmp_binvert_n_dc_ (mp_ptr dst, mp_srcptr numa, mp_size_t n, mp_ptr tp)
	计算 [numa,n] 在 B^n 下的逆元
uint	lmmp_binvert_uint_ (uint a)
	计算 a 在2^32下的逆元
ulong	lmmp_binvert_ulong_ (ulong a)
	计算 a 在2^64下的逆元
void	lmmp_binvert_unbalanced_ (mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_size_t n, mp_ptr tp)
	计算 [numa,na] 在 B^n 下的逆元
void	lmmp_binvert_unbalanced_1_ (mp_ptr dst, mp_limb_t a, mp_size_t n)
	计算 a 在 B^n 下的逆元
void	lmmp_binvert_unbalanced_2_ (mp_ptr dst, mp_srcptr numa, mp_size_t n)
	计算 [numa,2] 在 B^n 下的逆元
mp_limb_t	lmmp_cbrt_3_ (mp_limb_t a0, mp_limb_t a1, mp_limb_t a2)
	计算算数立方根 floor(cbrt(a0+a1*B+a2*B^2))
void	lmmp_cbrt_6_ (mp_ptr dst, mp_srcptr numa, mp_size_t na)
	计算算数立方根 floor(cbrt([numa,na]))
ulong	lmmp_cbrt_chebyshev_ (ulong n)
	计算算数立方根 floor(cbrt(n))
ulong	lmmp_cbrt_ulong_ (ulong n)
	计算算数立方根 floor(cbrt(n))
mp_limb_t	lmmp_cbrtapprox_3_ (mp_limb_t a0, mp_limb_t a1, mp_limb_t a2)
	计算近似立方根 floor(cbrt(a0+a1*B+a2*B^2))-[0|1]
void	lmmp_cbrtapprox_6_ (mp_ptr dst, mp_srcptr numa, mp_size_t na)
	计算近似立方根 floor(cbrt([numa,na]))-[0|1]
void	lmmp_divexact_ (mp_ptr dst, mp_srcptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn)
	精确除法（[dst,nn]=[np,nn]/[dp,dn]，且余数必须为0）
void	lmmp_divexact_1_ (mp_ptr dst, mp_srcptr np, mp_size_t nn, mp_limb_t d, mp_limb_t dinv)
	精确除法（[dst,nn]=[np,nn]/d，且余数必须为0）
void	lmmp_divexact_2_ (mp_ptr dst, mp_srcptr np, mp_size_t nn, mp_srcptr dp, mp_srcptr dinv)
	精确除法（[dst,nn]=[np,nn]/[dp,2]，且余数必须为0）
void	lmmp_divexact_basecase_ (mp_ptr dst, mp_ptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn, mp_limb_t dinv)
	精确除法（[dst,nn]=[np,nn]/[dp,dn]，且余数必须为0），朴素算法
void	lmmp_divexact_divide_ (mp_ptr dst, mp_srcptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn)
	精确除法（[dst,nn]=[np,nn]/[dp,dn]，且余数必须为0），分治算法
void	lmmp_divexact_unbalanced_ (mp_ptr dst, mp_srcptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn, mp_ptr dinv)
	精确除法（[dst,nn]=[np,nn]/[dp,dn]，且余数必须为0）
mp_size_t	lmmp_factorial_ (mp_ptr dst, mp_bitcnt_t bits, mp_size_t rn, uint n)
	计算 n! 阶乘
mp_size_t	lmmp_factorial_size_ (uint n, mp_bitcnt_t *bits)
	计算 n! 阶乘的 limb 缓冲区长度
mp_limb_t	lmmp_gcd_11_ (mp_limb_t u, mp_limb_t v)
	计算两个无符号整数的最大公约数
mp_limb_t	lmmp_gcd_1_ (mp_srcptr up, mp_size_t un, mp_limb_t vlimb)
	计算两个无符号整数的最大公约数
mp_size_t	lmmp_gcd_22_ (mp_ptr dst, mp_srcptr up, mp_srcptr vp)
	计算两个无符号整数的最大公约数
mp_size_t	lmmp_gcd_2_ (mp_ptr dst, mp_srcptr up, mp_size_t un, mp_srcptr vp)
	计算两个无符号整数的最大公约数
mp_size_t	lmmp_gcd_basecase_ (mp_ptr dst, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn)
	计算两个无符号整数的最大公约数（不建议使用此算法，更高版本可能被彻底弃用）
mp_size_t	lmmp_gcd_lehmer_ (mp_ptr dst, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn)
	计算两个无符号整数的最大公约数（Lehmer算法）
mp_size_t	lmmp_hyperfac_ (mp_ptr dst, mp_bitcnt_t bits, mp_size_t rn, ushort n)
	计算hyper阶乘（k^k累乘至n）
mp_size_t	lmmp_hyperfac_size_ (ushort n, mp_bitcnt_t *bits)
	计算hyper阶乘的 limb 缓冲区长度
bool	lmmp_is_prime_notrial_ (ulong n)
	判断素数（无试除法）
bool	lmmp_is_prime_uint_ (uint n)
	判断素数
bool	lmmp_is_prime_ulong_ (ulong n)
	判断素数
ulong	lmmp_mulmod_ulong_ (ulong a, ulong b, ulong mod, ulongp q)
	计算两个无符号整数的乘积，对mod取模，商放入 q 中
mp_size_t	lmmp_multinomial_ (mp_ptr dst, mp_size_t rn, uint n, const uintp r, uint m)
	计算多项式系数
mp_size_t	lmmp_multinomial_size_ (const uintp r, uint m, ulong *n)
	计算多项式系数的 limb 缓冲区长度
mp_size_t	lmmp_nCr_ (mp_ptr dst, mp_bitcnt_t bits, mp_size_t rn, uint n, uint r)
	计算 nCr 组合数 ( nCr = n! / (r!(n-r)!) )
mp_size_t	lmmp_nCr_size_ (uint n, uint r, mp_bitcnt_t *bits)
	计算 nCr 组合数的 limb 缓冲区长度
ulong	lmmp_next_prime_ulong_ (ulong n)
	大于n的下一个素数
mp_size_t	lmmp_nPr_ (mp_ptr dst, mp_bitcnt_t bits, mp_size_t rn, ulong n, ulong r)
	计算 nPr 排列数 ( nPr = n! / (n-r)! )
mp_size_t	lmmp_nPr_size_ (ulong n, ulong r, mp_bitcnt_t *bits)
	计算 nPr 排列数的 limb 缓冲区长度
ulong	lmmp_nthroot_ulong_ (ulong n, ulong root)
	计算 floor(n^(1/root))
mp_size_t	lmmp_odd_factorial_uint_ (mp_ptr dst, mp_size_t rn, uint n)
	计算 n! 阶乘的奇数部分
mp_size_t	lmmp_odd_nCr_uint_ (mp_ptr dst, mp_size_t rn, uint n, uint r)
	计算 nCr 组合数的奇数部分
mp_size_t	lmmp_odd_nCr_ushort_ (mp_ptr dst, mp_size_t rn, uint n, uint r)
	计算 nCr 组合数的奇数部分
mp_size_t	lmmp_odd_nPr_uint_ (mp_ptr dst, mp_size_t rn, ulong n, ulong r)
	计算 nPr 排列数的奇数部分
mp_size_t	lmmp_odd_nPr_ulong_ (mp_ptr dst, mp_size_t rn, ulong n, ulong r)
	计算 nPr 排列数的奇数部分
mp_size_t	lmmp_odd_nPr_ushort_ (mp_ptr dst, mp_size_t rn, ulong n, ulong r)
	计算 nPr 排列数的奇数部分
mp_size_t	lmmp_pow_ (mp_ptr dst, mp_size_t rn, mp_srcptr base, mp_size_t n, ulong exp)
	计算大整数幂 [dst,rn] = [base,n] ^ exp
mp_size_t	lmmp_pow_1_ (mp_ptr dst, mp_size_t rn, mp_limb_t base, ulong exp)
	计算幂次方 [dst,rn] = [base,1] ^ exp
mp_size_t	lmmp_pow_1_size_ (mp_limb_t base, ulong exp)
	计算幂次方需要的limb缓冲区长度 base ^ exp
mp_size_t	lmmp_pow_basecase_ (mp_ptr dst, mp_size_t rn, mp_srcptr base, mp_size_t n, ulong exp)
	计算奇数次幂算法 [dst,rn] = [base,n] ^ exp
mp_size_t	lmmp_pow_size_ (mp_srcptr base, mp_size_t n, ulong exp)
	计算幂次方需要的limb缓冲区长度 [base,n] ^ exp
mp_size_t	lmmp_pow_win2_ (mp_ptr dst, mp_size_t rn, mp_srcptr base, mp_size_t n, ulong exp)
	计算幂次方2比特窗口快速幂算法 [dst,rn] = [base,n] ^ exp
uint	lmmp_powmod_uint_ (uint base, ulong exp, uint mod)
	计算 base^exp 对 mod 取模
ulong	lmmp_powmod_ulong_ (ulong base, ulong exp, ulong mod)
	计算 base^exp 对 mod 取模
ulong	lmmp_prev_prime_ulong_ (ulong n)
	小于等于n的上一个素数
mp_size_t	lmmp_primefac_ (mp_ptr dst, mp_size_t rn, uint n)
	计算质数阶乘（不超过n的质数累乘）
mp_size_t	lmmp_primefac_size_ (uint n)
	计算质数阶乘的 limb 缓冲区长度
mp_size_t	lmmp_remove_ (mp_ptr np, mp_size_t *nn, mp_srcptr dp, mp_size_t dn)
	除去[np,nn]中的[dp,dn]的因子
ulong	lmmp_sqrt_ulong_ (ulong a)
	计算算术平方根 floor(sqrt(a))
mp_size_t	lmmp_superfac_ (mp_ptr dst, mp_bitcnt_t bits, mp_size_t rn, ushort n)
	计算super阶乘（k!累乘至n）
mp_size_t	lmmp_superfac_size_ (ushort n, mp_bitcnt_t *bits)
	计算super阶乘的 limb 缓冲区长度
ushortp	lmmp_trialdiv_ (mp_srcptr np, mp_size_t nn, ushort N, ushort *rn)
	试除法
mp_size_t	lmmp_u16_pow_1_ (mp_ptr dst, mp_size_t rn, ulong base, ulong exp)
	计算幂次方 [dst,rn] = [base,1] ^ exp
mp_size_t	lmmp_u32_pow_1_ (mp_ptr dst, mp_size_t rn, ulong base, ulong exp)
	计算幂次方 [dst,rn] = [base,1] ^ exp
mp_size_t	lmmp_u4_pow_1_ (mp_ptr dst, mp_size_t rn, ulong base, ulong exp)
	计算幂次方 [dst,rn] = [base,1] ^ exp
mp_size_t	lmmp_u64_pow_1_ (mp_ptr dst, mp_size_t rn, ulong base, ulong exp)
	计算幂次方 [dst,rn] = [base,1] ^ exp
mp_size_t	lmmp_u8_pow_1_ (mp_ptr dst, mp_size_t rn, ulong base, ulong exp)
	计算幂次方 [dst,rn] = [base,1] ^ exp

类型定义说明

schar

typedef int8_t schar

在文件 numth.h 第 28 行定义.

scharp

typedef int8_t* scharp

在文件 numth.h 第 36 行定义.

sint

typedef int32_t sint

在文件 numth.h 第 33 行定义.

sintp

typedef int32_t* sintp

在文件 numth.h 第 40 行定义.

slong

typedef int64_t slong

在文件 numth.h 第 34 行定义.

slongp

typedef int64_t* slongp

在文件 numth.h 第 42 行定义.

sshort

typedef int16_t sshort

在文件 numth.h 第 30 行定义.

sshortp

typedef int16_t* sshortp

在文件 numth.h 第 38 行定义.

uchar

typedef uint8_t uchar

Copyright (C) 2026 HJimmyK(Jericho Knox)

This file is part of LAMMP.

LAMMP is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License (LGPL) as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.

This program is distributed WITHOUT ANY WARRANTY.

See https://www.gnu.org/licenses/.

在文件 numth.h 第 27 行定义.

ucharp

typedef uint8_t* ucharp

在文件 numth.h 第 35 行定义.

uint

typedef uint32_t uint

在文件 numth.h 第 31 行定义.

uintp

typedef uint32_t* uintp

在文件 numth.h 第 39 行定义.

ulong

typedef uint64_t ulong

在文件 numth.h 第 32 行定义.

ulongp

typedef uint64_t* ulongp

在文件 numth.h 第 41 行定义.

ushort

typedef uint16_t ushort

在文件 numth.h 第 29 行定义.

ushortp

typedef uint16_t* ushortp

在文件 numth.h 第 37 行定义.

函数说明

lmmp_2factorial_()

mp_size_t lmmp_2factorial_	(	mp_ptr	dst,
		mp_bitcnt_t	bits,
		mp_size_t	rn,
		uint	n
	)

计算 n!! 双阶乘

参数

dst	结果指针
bits	n!! 的2的因子数
rn	结果指针的 limb 长度
n	双阶乘的阶数

警告

dst!=NULL, rn>0

注解

0的双阶乘为1，n为偶数时，n!!=2*4*...*n，n为奇数时，n!!=1*3*...*n

返回

返回 dst 的实际 limb 长度

lmmp_2factorial_size_()

mp_size_t lmmp_2factorial_size_	(	uint	n,
		mp_bitcnt_t *	bits
	)

计算 n!! 双阶乘的 limb 缓冲区长度

参数

n	双阶乘的阶数
bits	被修改为 n!! 的2的因子数

警告

bits!=NULL

返回

n!! 双阶乘的 limb 缓冲区长度

lmmp_arith_seqprod_()

mp_size_t lmmp_arith_seqprod_	(	mp_ptr	dst,
		mp_size_t	rn,
		uint	x,
		uint	n,
		uint	m
	)

计算等差数列乘积 x(x+m)...(x+n*m)

参数

dst	结果指针
rn	结果指针的 limb 长度
x	首项
n	等差数列共n+1项
m	公差（大于1）

警告

x>0, m>1, n>0, dst!=NULL, rn>0, x+n*m<=0xffffffff

返回

结果的实际的 limb 缓冲区长度

lmmp_arith_seqprod_size_()

mp_size_t lmmp_arith_seqprod_size_	(	uint	x,
		uint	n,
		uint	m
	)

计算等差数列乘积 x(x+m)...(x+n*m)的 limb 缓冲区长度

参数

x	首项
n	等差数列共n+1项
m	公差

警告

x>0, m>1, n>0, x+n*m<=0xffffffff

返回

等差数列乘积的 limb 缓冲区长度（比实际长度多 1-2 个 limb）

计算等差数列乘积 x(x+m)...(x+n*m)的 limb 缓冲区长度

This file is part of LAMMP.

LAMMP is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License (LGPL) as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.

This program is distributed WITHOUT ANY WARRANTY.

See https://www.gnu.org/licenses/.

在文件 arith_seqprod.c 第 24 行定义.

                                                           {
    lmmp_param_assert(x > 0);
    lmmp_param_assert(n > 0);
    lmmp_param_assert(m > 1);
    // x(x+m)(x+2m)...(x+nm) <= (x+m*n/2)^(n+1)
    ulong t = (x + ((ulong)m * n + 1) / 2);
    mp_size_t rn1 = lmmp_pow_1_size_(t, n + 1);
    // 共有n+1个数，每个数的位数最多为uint，一个limb最多可以容纳两个uint乘积
    mp_size_t rn2 = (n + 1) / 2 + 1;
    return LMMP_MIN(rn1, rn2);
}

引用了 LMMP_MIN, lmmp_param_assert, lmmp_pow_1_size_(), n , 以及 t.

函数调用图:

lmmp_binvert_()

void lmmp_binvert_	(	mp_ptr	dst,
		mp_srcptr	numa,
		mp_size_t	na,
		mp_size_t	n
	)

计算 [numa,na] 在 B^n 下的逆元

参数

dst	结果指针（长度为 n 个limb）
numa	待求逆元指针（长度为 na 个limb）
na	待求逆元的 limb 长度
n	结果的 limb 长度

警告

n>=na>0, numa!=NULL, dst!=NULL, numa[0]%2==1, sep(dst,numa)

lmmp_binvert_2_()

void lmmp_binvert_2_	(	mp_ptr	dst,
		mp_srcptr	numa
	)

计算 [numa,2] 在 B^2 下的逆元

参数

numa	待求逆元指针（长度为 2 个limb）
dst	结果指针（长度为 2 个limb）

警告

numa!=NULL, dst!=NULL, numa[0]%2==1, eqsep(dst,numa)

在文件 binvert_1.c 第 56 行定义.

                                                 {
    lmmp_param_assert(numa[0] % 2 == 1);
    mp_limb_t k, t;
    mp_limb_t a1 = numa[1];
    mp_limb_t a0 = numa[0];
    mp_limb_t xn = lmmp_binvert_ulong_(a0);
    mp_limb_t z;
    /*
     xn * a0 == 1 + k * B
     yn := xn * (2 - a * xn) mod B^2
        := xn * (2 - a0 * xn - a1 * xn * B) mod B^2
        := xn * (2 - 1 - k*B - a1 * xn * B) mod B^2
        := xn * (1 - k*B - a1 * xn * B) mod B^2
        := (xn - xn*k * B - a1 * xn^2 * B) mod B^2
    */
    _umul64to128_(a0, xn, &t, &k);
    z = xn * k;
    z += a1 * xn * xn;
    dst[0] = xn;
    dst[1] = -z;
}

引用了 _umul64to128_(), a0, a1, k, lmmp_binvert_ulong_(), lmmp_param_assert, n, t , 以及 xn.

被这些函数引用 lmmp_binvert_3_(), lmmp_binvert_4_(), lmmp_binvert_n_dc_(), lmmp_binvert_unbalanced_2_() , 以及 lmmp_divexact_().

函数调用图:

这是这个函数的调用关系图:

lmmp_binvert_3_()

void lmmp_binvert_3_	(	mp_ptr	dst,
		mp_srcptr	numa
	)

计算 [numa,3] 在 B^3 下的逆元

参数

numa	待求逆元指针（长度为 3 个limb）
dst	结果指针（长度为 3 个limb）

警告

numa!=NULL, dst!=NULL, numa[0]%2==1, sep(dst,numa)

被这些函数引用 lmmp_binvert_n_dc_().

这是这个函数的调用关系图:

lmmp_binvert_4_()

void lmmp_binvert_4_	(	mp_ptr	dst,
		mp_srcptr	numa
	)

计算 [numa,4] 在 B^4 下的逆元

参数

numa	待求逆元指针（长度为 4 个limb）
dst	结果指针（长度为 4 个limb）

警告

numa!=NULL, dst!=NULL, numa[0]%2==1, sep(dst,numa)

被这些函数引用 lmmp_binvert_n_dc_().

这是这个函数的调用关系图:

lmmp_binvert_n_dc_()

void lmmp_binvert_n_dc_	(	mp_ptr	dst,
		mp_srcptr	numa,
		mp_size_t	n,
		mp_ptr	tp
	)

计算 [numa,n] 在 B^n 下的逆元

参数

numa	待求逆元指针（长度为 n 个limb）
dst	结果指针（长度为 n 个limb）
n	结果的 limb 长度
tp	临时工作区指针（长度为 5*(n+1)/2 个limb）

警告

numa!=NULL, dst!=NULL, numa[0]%2==1, sep(dst,numa,tp)

被这些函数引用 lmmp_divexact_divide_() , 以及 lmmp_divexact_unbalanced_().

这是这个函数的调用关系图:

lmmp_binvert_uint_()

uint lmmp_binvert_uint_

(

uint

a

)

计算 a 在2^32下的逆元

参数

a

待求逆元

警告

a%2==1

返回

逆元

在文件 binvert_1.c 第 30 行定义.

                                {
    lmmp_param_assert(a % 2 == 1);
    uint r, y;
 
    r = binv_tab[(a / 2) & 0x7F]; /* 8 bits */
    y = 1 - a * r;
    r = r * (1 + y); /* 16 bits */
    y *= y;
    r = r * (1 + y); /* 32 bits */
    return r;
}

引用了 binv_tab, lmmp_param_assert , 以及 n.

lmmp_binvert_ulong_()

ulong lmmp_binvert_ulong_

(

ulong

a

)

计算 a 在2^64下的逆元

参数

a

待求逆元

警告

a%2==1

返回

逆元

在文件 binvert_1.c 第 42 行定义.

                                   {
    lmmp_param_assert(a % 2 == 1);
    ulong r, y;
 
    r = binv_tab[(a / 2) & 0x7F]; /* 8 bits */
    y = 1 - a * r;
    r = r * (1 + y); /* 16 bits */
    y *= y;
    r = r * (1 + y); /* 32 bits */
    y *= y;
    r = r * (1 + y); /* 64 bits */
    return r;
}

引用了 binv_tab, lmmp_param_assert , 以及 n.

被这些函数引用 lmmp_binvert_2_(), lmmp_binvert_n_dc_(), lmmp_binvert_unbalanced_1_(), lmmp_divexact_(), lmmp_is_prime_notrial_() , 以及 lmmp_powmod_ulong_().

这是这个函数的调用关系图:

lmmp_binvert_unbalanced_()

void lmmp_binvert_unbalanced_	(	mp_ptr	dst,
		mp_srcptr	numa,
		mp_size_t	na,
		mp_size_t	n,
		mp_ptr	tp
	)

计算 [numa,na] 在 B^n 下的逆元

参数

dst	结果指针（长度为 n 个limb）
numa	待求逆元指针（长度为 na 个limb）
na	待求逆元的 limb 长度
n	结果的 limb 长度
tp	临时工作区指针（长度为 (9*na+5)/2 个limb）

警告

numa[0]%2==1, n>na, dst!=NULL, numa!=NULL, tp!=NULL, sep(dst,numa,tp)

lmmp_binvert_unbalanced_1_()

void lmmp_binvert_unbalanced_1_	(	mp_ptr	dst,
		mp_limb_t	a,
		mp_size_t	n
	)

计算 a 在 B^n 下的逆元

参数

dst	结果指针（长度为 n 个limb）
a	待求逆元
n	结果的 limb 长度

警告

a%2==1, n>1, dst!=NULL

被这些函数引用 lmmp_binvert_().

这是这个函数的调用关系图:

lmmp_binvert_unbalanced_2_()

void lmmp_binvert_unbalanced_2_	(	mp_ptr	dst,
		mp_srcptr	numa,
		mp_size_t	n
	)

计算 [numa,2] 在 B^n 下的逆元

参数

dst	结果指针（长度为 n 个limb）
numa	待求逆元指针（长度为 2 个limb）
n	结果的 limb 长度

警告

numa[0]%2==1, n>2, dst!=NULL, numa!=NULL, sep(dst,numa)

被这些函数引用 lmmp_binvert_().

这是这个函数的调用关系图:

lmmp_cbrt_3_()

mp_limb_t lmmp_cbrt_3_	(	mp_limb_t	a0,
		mp_limb_t	a1,
		mp_limb_t	a2
	)

计算算数立方根 floor(cbrt(a0+a1*B+a2*B^2))

参数

a0	低位 limb
a1	中位 limb
a2	高位 limb

警告

a1>0

注解

a2可以为0，但a1需要大于0，即这个数至少应有65个bit

返回

floor(cbrt(a0+a1*B+a2*B^2))

在文件 cbrt.c 第 74 行定义.

                                                                 {
    lmmp_param_assert(a1 > 0);
 
    mp_limb_t r = lmmp_cbrtapprox_3_(a0, a1, a2);
    if (r == LIMB_MAX)
        return LIMB_MAX;
    mp_limb_t t[3], a[3] = {a0, a1, a2};
    lmmp_cube_3_(t, r + 1);
    int cmp = lmmp_cmp_(t, a, 3);
    // approx的结果至多只会低估1
    if (cmp <= 0)
        return r + 1;
    else
        return r;
}

引用了 a0, a1, a2, LIMB_MAX, lmmp_cbrtapprox_3_(), lmmp_cmp_(), lmmp_cube_3_(), lmmp_param_assert, n , 以及 t.

函数调用图:

lmmp_cbrt_6_()

void lmmp_cbrt_6_	(	mp_ptr	dst,
		mp_srcptr	numa,
		mp_size_t	na
	)

计算算数立方根 floor(cbrt([numa,na]))

参数

dst	结果指针（长度为 2 个limb）
numa	被开方数指针
na	被开方数的 limb 长度

警告

dst!=NULL, numa!=NULL, 3<na<=6, numa[na-1]!=0, eqsep(dst,numa)

在文件 cbrt.c 第 152 行定义.

                                                            {
    mp_limb_t ret[2];
    lmmp_cbrtapprox_6_(ret, numa, na);
 
    if (ret[1] == LIMB_MAX && ret[0] == LIMB_MAX) {
        dst[0] = LIMB_MAX;
        dst[1] = LIMB_MAX;
    } else {
        mp_limb_t r[2];
        r[0] = ret[0] + 1;
        r[1] = ret[1] + (r[0] == 0 ? 1 : 0);
        mp_limb_t t[10];
        mp_size_t tn = lmmp_cube_(t, r, 2, t + 6);
        if (tn < na) {
            dst[0] = r[0];
            dst[1] = r[1];
        } else if (tn > na) {
            dst[0] = ret[0];
            dst[1] = ret[1];
        } else if (tn == na) {
            int cmp = lmmp_cmp_(t, numa, tn);
            // approx的结果至多只会低估1
            if (cmp <= 0) {
                dst[0] = r[0];
                dst[1] = r[1];
            } else {
                dst[0] = ret[0];
                dst[1] = ret[1];
            }
        }
    }
}

引用了 LIMB_MAX, lmmp_cbrtapprox_6_(), lmmp_cmp_(), lmmp_cube_(), n , 以及 t.

函数调用图:

lmmp_cbrt_chebyshev_()

ulong lmmp_cbrt_chebyshev_

(

ulong

n

)

计算算数立方根 floor(cbrt(n))

参数

n

被开方数

返回

floor(cbrt(n))

警告

n>0

注解

使用Chebyshev估计，当n较大时，此算法更占优势

在文件 cbrt_1.c 第 40 行定义.

                                    {
    lmmp_param_assert(n > 0);
    typedef union {
        ulong uword_val;
        double double_val;
    } uni;
 
    int rem, mul;
    double factor, root, dec, dec2;
    ulong ret, expo, table_index;
    uni alias;
 
    /* upper_limit is the max cube root possible for one word */
 
    const ulong upper_limit = 2642245;              /* 2642245 < (2^64)^(1/3) */
    const ulong expo_mask = 0x7FF0000000000000;     /* exponent bits in double */
    const ulong mantissa_mask = 0x000FFFFFFFFFFFFF; /* mantissa bits in float */
    const ulong table_mask = 0x000F000000000000;    /* first 4 bits of mantissa */
    const uint mantissa_bits = 52;
    const ulong bias_hex = 0x3FE0000000000000;
    const uint bias = 1022;
    alias.double_val = (double)n;
 
    expo = alias.uword_val & expo_mask; /* extracting exponent */
    expo >>= mantissa_bits;
    expo -= bias; /* Subtracting bias */
 
    /* extracting first 4 bits of mantissa, this will help select correct poly */
    /* note mantissa of 0.5 is 0x0000000000000 not 0x1000000000000 */
 
    table_index = alias.uword_val & table_mask;
    table_index >>= (mantissa_bits - 4);
 
    /* extracting decimal part, 0.5 <= dec <= 1 */
    ret = alias.uword_val & mantissa_mask;
    ret |= bias_hex;
    alias.uword_val = ret;
    dec = alias.double_val;
 
    rem = expo % 3;
    expo /= 3;                  /* cube root of 2^expo */
    factor = factor_table[rem]; /* select factor */
 
    /* Calculating cube root of dec using chebyshev approximation polynomial */
    /* Evaluating approx polynomial at (dec) by Estrin's scheme */
 
    dec2 = dec * dec;
    root = (coeff[table_index][0] + coeff[table_index][1] * dec);
    root += (coeff[table_index][2] * dec2);
 
    mul = (ulong)1 << expo; /* mul = 2^expo */
    root *= mul;            /* dec^(1/3) * 2^(expo/3) */
    root *= factor;         /* root*= (expo%3)^(1/3) */
    ret = root;
 
    /* In case ret^3 or (ret+1)^3 will cause overflow */
 
    if (ret >= upper_limit) {
        if (n >= upper_limit * upper_limit * upper_limit)
            return upper_limit;
        ret = upper_limit - 1;
    }
    while (ret * ret * ret <= n) {
        (ret) += 1;
        if (ret == upper_limit)
            break;
    }
    while (ret * ret * ret > n) (ret) -= 1;
 
    return ret;
}

引用了 coeff, factor_table, lmmp_param_assert , 以及 n.

被这些函数引用 lmmp_cbrt_ulong_().

这是这个函数的调用关系图:

lmmp_cbrt_ulong_()

ulong lmmp_cbrt_ulong_

(

ulong

n

)

计算算数立方根 floor(cbrt(n))

参数

n

被开方数

返回

floor(cbrt(n))

注解

会依据n的大小，选择合适的算法

在文件 cbrt_1.c 第 133 行定义.

                                {
    double val, x, xcub, num, den;
    ulong ret, upper_limit;
 
    if (n < 125)
        return (n >= 1) + (n >= 8) + (n >= 27) + (n >= 64);
    if (n < 1331)
        return 5 + (n >= 216) + (n >= 343) + (n >= 512) + (n >= 729) + (n >= 1000);
    if (n < 4913)
        return 11 + (n >= 1728) + (n >= 2197) + (n >= 2744) + (n >= 3375) + (n >= 4096);
 
    val = (double)n;
 
    if (n >= 1ULL << 46)
        return lmmp_cbrt_chebyshev_(n);
 
    upper_limit = 2642245; /* 2642245 = floor((2^64)^(1/3)) */
 
    x = lmmp_cbrt_estimate((double)n);
 
    /* Kahan's iterations to get cube root */
 
    xcub = x * x * x;
    num = (xcub - val) * x;
    den = (xcub + xcub + val);
    x -= (num / den);
    ret = x;
 
    if (ret >= upper_limit) {
        if (n >= upper_limit * upper_limit * upper_limit)
            return upper_limit;
        ret = upper_limit - 1;
    }
    while (ret * ret * ret <= n) {
        (ret) += 1;
        if (ret == upper_limit)
            break;
    }
    while (ret * ret * ret > n) (ret) -= 1;
 
    return ret;
}

引用了 lmmp_cbrt_chebyshev_(), lmmp_cbrt_estimate() , 以及 n.

被这些函数引用 lmmp_nthroot_ulong_().

函数调用图:

这是这个函数的调用关系图:

lmmp_cbrtapprox_3_()

mp_limb_t lmmp_cbrtapprox_3_	(	mp_limb_t	a0,
		mp_limb_t	a1,
		mp_limb_t	a2
	)

计算近似立方根 floor(cbrt(a0+a1*B+a2*B^2))-[0|1]

参数

a0	低位 limb
a1	中位 limb
a2	高位 limb

警告

a1>0

注解

a2可以为0，但a1需要大于0，即这个数至少应有65个bit

返回

floor(cbrt(a0+a1*B+a2*B^2))-[0|1]

在文件 cbrt.c 第 32 行定义.

                                                                       {
    lmmp_param_assert(a1 > 0);
    mp_limb_t x[2];
    /* exact high 65 bits */
    mp_limb_t a_hi;
    mp_bitcnt_t bits;
    if (a2 == 0) {
        mp_bitcnt_t a1_bits = lmmp_limb_bits_(a1);
        bits = LIMB_BITS + a1_bits;
        a1_bits--;
        if (a1_bits == 0)
            a_hi = a0;
        else
            a_hi = (a1 << (LIMB_BITS - a1_bits)) | (a0 >> a1_bits);
    } else {
        mp_bitcnt_t a2_bits = lmmp_limb_bits_(a2);
        bits = LIMB_BITS * 2 + a2_bits;
        a2_bits--;
        if (a2_bits == 0)
            a_hi = a1;
        else
            a_hi = (a2 << (LIMB_BITS - a2_bits)) | (a1 >> a2_bits);
    }
    lmmp_debug_assert(bits >= 65);
 
    x[1] = bits - 1;
    x[0] = log2_fixed_64(a_hi);
 
    mp_limb_t r = lmmp_div_1_(x, x, 2, 3);
    if (2 * r >= 3) // round
        lmmp_inc(x);
 
    mp_bitcnt_t shift = x[1];
    x[0] = exp2_fixed_64(x[0]);
 
    lmmp_debug_assert(shift <= 64);
    if (shift == 64)
        return LIMB_MAX;
    else
        return (x[0] >> (64 - shift)) | (1ULL << shift);
}

引用了 a0, a1, a2, exp2_fixed_64(), LIMB_BITS, LIMB_MAX, lmmp_debug_assert, lmmp_div_1_(), lmmp_inc, lmmp_limb_bits_, lmmp_param_assert, log2_fixed_64() , 以及 n.

被这些函数引用 lmmp_cbrt_3_().

函数调用图:

这是这个函数的调用关系图:

lmmp_cbrtapprox_6_()

void lmmp_cbrtapprox_6_	(	mp_ptr	dst,
		mp_srcptr	numa,
		mp_size_t	na
	)

计算近似立方根 floor(cbrt([numa,na]))-[0|1]

参数

dst	结果指针（长度为 2 个limb）
numa	被开方数指针
na	被开方数的 limb 长度

警告

dst!=NULL, numa!=NULL, 3<na<=6, numa[na-1]!=0, eqsep(dst,numa)

在文件 cbrt.c 第 106 行定义.

                                                                  {
    lmmp_param_assert(na > 3 && na <= 6);
    lmmp_param_assert(dst != NULL && numa != NULL);
    lmmp_param_assert(numa[na - 1] != 0);
    /* extract the first 129 bits */
    int bits = lmmp_limb_bits_(numa[na - 1]);
    mp_bitcnt_t n = bits - 1;
    mp_limb_t high, low;
    if (bits == 1) {
        high = numa[na - 2];
        low = numa[na - 3];
    } else {
        bits--;
        high = (numa[na - 1] << (64 - bits)) | (numa[na - 2] >> bits);
        low = (numa[na - 2] << (64 - bits)) | (numa[na - 3] >> bits);
    }
 
    n += LIMB_BITS * (na - 1);
    mp_limb_t x[3] = {0, 0, n};
 
    log2_fixed_128(x, high, low);
    mp_limb_t r = lmmp_div_1_(x, x, 3, 3);
    if (2 * r >= 3) // round
       lmmp_inc(x);
 
    n = x[2];
    high = x[1];
    low = x[0];
 
    exp2_fixed_128(x, high, low);
 
    lmmp_debug_assert(n >= 64 && n <= 128);
    if (n == 64) {
        dst[0] = x[1];
        dst[1] = 1;
    } else if (n < 128) {
        n -= 64;
        mp_limb_t t = 1ULL << n;
        dst[1] = (x[1] >> (64 - n)) | t;
        dst[0] = (x[1] << n) | (x[0] >> (64 - n));
    } else {
        dst[1] = LIMB_MAX;
        dst[0] = LIMB_MAX;
    }
}

引用了 exp2_fixed_128(), LIMB_BITS, LIMB_MAX, lmmp_debug_assert, lmmp_div_1_(), lmmp_inc, lmmp_limb_bits_, lmmp_param_assert, log2_fixed_128(), n , 以及 t.

被这些函数引用 lmmp_cbrt_6_().

函数调用图:

这是这个函数的调用关系图:

lmmp_divexact_()

void lmmp_divexact_	(	mp_ptr	dst,
		mp_srcptr	np,
		mp_size_t	nn,
		mp_srcptr	dp,
		mp_size_t	dn
	)

精确除法（[dst,nn]=[np,nn]/[dp,dn]，且余数必须为0）

参数

dst	结果指针（长度为 nn-dn+1 个limb）
np	被除数指针（长度为 nn 个limb）
nn	被除数的 limb 长度
dp	除数指针（长度为 dn 个limb）
dn	结除数的 limb 长度

警告

dp[0]%2==1, nn>=dn>0, dst!=NULL, np!=NULL, dp!=NULL, eqsep(dst,np), sep([dst|np],dp)

被这些函数引用 lmmp_odd_nCr_div_().

这是这个函数的调用关系图:

lmmp_divexact_1_()

void lmmp_divexact_1_	(	mp_ptr	dst,
		mp_srcptr	np,
		mp_size_t	nn,
		mp_limb_t	d,
		mp_limb_t	dinv
	)

精确除法（[dst,nn]=[np,nn]/d，且余数必须为0）

参数

dst	结果指针（长度为 nn 个limb）
np	被除数指针（长度为 nn 个limb）
nn	被除数的 limb 长度
d	除数
dinv	除数的逆元（d*dinv==1 mod 2^64）

警告

d%2==1, d*dinv==1 mod 2^64, nn>0, dst!=NULL, np!=NULL, eqsep(dst,np)

精确除法（[dst,nn]=[np,nn]/d，且余数必须为0）

This file is part of LAMMP.

LAMMP is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License (LGPL) as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.

This program is distributed WITHOUT ANY WARRANTY.

See https://www.gnu.org/licenses/.

在文件 divexact.c 第 24 行定义.

                                                                                           {
    lmmp_param_assert(d % 2 == 1);
    lmmp_param_assert(dinv * d == 1);
    lmmp_param_assert(nn > 0);
    lmmp_debug_assert(dst != NULL && np != NULL);
    mp_limb_t c = 0;
    mp_limb_t l, s, lo, hi, q;
    mp_size_t i;
    for (i = 0; i < nn - 1; i++) {
        s = np[i];
        l = s - c;
        c = l > s;
        q = l * dinv;
        dst[i] = q;
        _umul64to128_(q, d, &lo, &hi);
        c += hi;
    }
    s = np[i];
    l = s - c;
    c = l > s;
    q = l * dinv;
    dst[i] = q;
    lmmp_debug_assert(c == 0);
}

引用了 _umul64to128_(), c, l, lmmp_debug_assert, lmmp_param_assert, lo, n , 以及 s.

被这些函数引用 lmmp_divexact_().

函数调用图:

这是这个函数的调用关系图:

lmmp_divexact_2_()

void lmmp_divexact_2_	(	mp_ptr	dst,
		mp_srcptr	np,
		mp_size_t	nn,
		mp_srcptr	dp,
		mp_srcptr	dinv
	)

精确除法（[dst,nn]=[np,nn]/[dp,2]，且余数必须为0）

参数

dst	结果指针（长度为 nn-1 个limb）
np	被除数指针（长度为 nn 个limb）
nn	被除数的 limb 长度
dp	除数指针（长度为 2 个limb）
dinv	除数的逆元指针（长度为 2 个limb）

警告

dp[0]%2==1, dp*dinv==1 mod 2^128, nn>1, dst!=NULL, np!=NULL, eqsep(dst,np), sep(dp,dinv,[dst|np])

lmmp_divexact_basecase_()

void lmmp_divexact_basecase_	(	mp_ptr	dst,
		mp_ptr	np,
		mp_size_t	nn,
		mp_srcptr	dp,
		mp_size_t	dn,
		mp_limb_t	dinv
	)

精确除法（[dst,nn]=[np,nn]/[dp,dn]，且余数必须为0），朴素算法

参数

dst	结果指针（长度为 nn-dn+1 个limb）
np	被除数指针（长度为 nn 个limb），将会被覆写为全零
nn	被除数的 limb 长度
dp	除数指针（长度为 dn 个limb）
dn	结除数的 limb 长度
dinv	除数低位dp[0]关于B的逆元

警告

dp[0]%2==1, d[0]*dinv==1 mod 2^64, nn>=dn>0, dst!=NULL, np!=NULL, dp!=NULL, eqsep(dst,np), sep(dp,[dst|np])

注解

若dst==np，将会覆写 [dst,nn] 区域，其中 [dst,nn-dn+1] 为商，[dst+nn-dn+1,dn-1] 将会被覆写为0

lmmp_divexact_divide_()

void lmmp_divexact_divide_	(	mp_ptr	dst,
		mp_srcptr	np,
		mp_size_t	nn,
		mp_srcptr	dp,
		mp_size_t	dn
	)

精确除法（[dst,nn]=[np,nn]/[dp,dn]，且余数必须为0），分治算法

参数

dst	结果指针（长度为 nn-dn+1 个limb）
np	被除数指针（长度为 nn 个limb）
nn	被除数的 limb 长度
dp	除数指针（长度为 dn 个limb）
dn	结除数的 limb 长度

警告

dp[0]%2==1, nn>=dn>0, dn>=nn-dn+1, dst!=NULL, np!=NULL, dp!=NULL, sep(dst,np,dp)

lmmp_divexact_unbalanced_()

void lmmp_divexact_unbalanced_	(	mp_ptr	dst,
		mp_srcptr	np,
		mp_size_t	nn,
		mp_srcptr	dp,
		mp_size_t	dn,
		mp_ptr	dinv
	)

精确除法（[dst,nn]=[np,nn]/[dp,dn]，且余数必须为0）

参数

dst	结果指针（长度为 nn-dn+1 个limb）
np	被除数指针（长度为 nn 个limb），将会被修改为全零
nn	被除数的 limb 长度
dp	除数指针（长度为 dn 个limb）
dn	结除数的 limb 长度
dinv	除数[dp,dn]关于B^dn的逆元，若为NULL，则自动计算

警告

dp[0]%2==1, nn>=dn>0, dst!=NULL, np!=NULL, dp!=NULL, eqsep(dst,np), sep(dp,dinv,[dst|np])

注解

若dst==np，只会覆写 [dst,nn-dn+1] 区域

lmmp_factorial_()

mp_size_t lmmp_factorial_	(	mp_ptr	dst,
		mp_bitcnt_t	bits,
		mp_size_t	rn,
		uint	n
	)

计算 n! 阶乘

参数

dst	结果指针
bits	n! 的2的因子数
rn	结果指针的 limb 长度
n	阶乘的阶数

警告

dst!=NULL, rn>0

返回

返回 dst 的实际 limb 长度

lmmp_factorial_size_()

mp_size_t lmmp_factorial_size_	(	uint	n,
		mp_bitcnt_t *	bits
	)

计算 n! 阶乘的 limb 缓冲区长度

参数

n	阶乘的阶数
bits	被修改为 n! 的2的因子数

警告

bits!=NULL

返回

n! 阶乘的 limb 缓冲区长度（比实际长度多）

被这些函数引用 lmmp_odd_nCr_uint_().

这是这个函数的调用关系图:

lmmp_gcd_11_()

mp_limb_t lmmp_gcd_11_	(	mp_limb_t	u,
		mp_limb_t	v
	)

计算两个无符号整数的最大公约数

参数

u	第一个无符号整数
v	第二个无符号整数

返回

最大公约数

警告

u!=0, v!=0

计算两个无符号整数的最大公约数

This file is part of LAMMP.

LAMMP is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License (LGPL) as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.

This program is distributed WITHOUT ANY WARRANTY.

See https://www.gnu.org/licenses/.

在文件 gcd_1.c 第 24 行定义.

                                                 {
    lmmp_param_assert(u > 0 && v > 0);
    int count, k;
    k = lmmp_tailing_zeros_(u | v);
    u >>= lmmp_tailing_zeros_(u);
    v >>= lmmp_tailing_zeros_(v);
    while (u != v) {
        if (u > v) {
            u -= v;
            count = lmmp_tailing_zeros_(u); 
            u >>= count;
        } else {
            v -= u;
            count = lmmp_tailing_zeros_(v);
            v >>= count;
        }
    }
    return u << k;
}

引用了 k, lmmp_param_assert, lmmp_tailing_zeros_ , 以及 n.

被这些函数引用 lmmp_gcd_1_(), lmmp_gcd_22_() , 以及 lmmp_gcd_lehmer_().

这是这个函数的调用关系图:

lmmp_gcd_1_()

mp_limb_t lmmp_gcd_1_	(	mp_srcptr	up,
		mp_size_t	un,
		mp_limb_t	vlimb
	)

计算两个无符号整数的最大公约数

参数

up	第一个无符号整数指针
un	第一个无符号整数的 limb 长度
v	第二个无符号整数

警告

v!=0, up!=NULL, un>0

返回

最大公约数

在文件 gcd_1.c 第 44 行定义.

                                                                   {
    lmmp_param_assert(un > 0);
    lmmp_param_assert(vlimb > 0);
    mp_limb_t ulimb;
    if (un == 1) {
        ulimb = up[0];
    } else {
        ulimb = lmmp_mod_1_(up, un, vlimb);
    }
    if (ulimb == 0)
        return vlimb;
    else
        return lmmp_gcd_11_(ulimb, vlimb);
}

引用了 lmmp_gcd_11_(), lmmp_mod_1_(), lmmp_param_assert , 以及 n.

被这些函数引用 lmmp_gcd_lehmer_().

函数调用图:

这是这个函数的调用关系图:

lmmp_gcd_22_()

mp_size_t lmmp_gcd_22_	(	mp_ptr	dst,
		mp_srcptr	up,
		mp_srcptr	vp
	)

计算两个无符号整数的最大公约数

参数

up	第一个无符号整数指针，长度为 2
vp	第二个无符号整数指针，长度为 2
dst	结果指针（长度为 2，两个 limb 都会进行写入，即使最高位可能为0）

警告

up!=NULL, vp!=NULL, [up,2]!=0, [vp,2]!=0, dst!=NULL, eqsep(dst,[up|vp])

注解

我们不要求 up 和 vp 的高位不为 0，但要求两个数均不可以高低位全为 0

返回

dst 的实际 limb 长度

计算两个无符号整数的最大公约数

This file is part of LAMMP.

LAMMP is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License (LGPL) as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.

This program is distributed WITHOUT ANY WARRANTY.

See https://www.gnu.org/licenses/.

在文件 gcd_2.c 第 27 行定义.

                                                               {
    lmmp_param_assert(dst != NULL);
    lmmp_param_assert(up != NULL);
    lmmp_param_assert(vp != NULL);
    lmmp_param_assert(!(up[1] == 0 && up[0] == 0));
    lmmp_param_assert(!(vp[1] == 0 && vp[0] == 0));
    mp_limb_t u[2] = { up[0], up[1] };
    mp_limb_t v[2] = { vp[0], vp[1] };
    int k, cnt;
 
    if (u[1] == 0 && v[1] == 0) {
        dst[0] = lmmp_gcd_11_(u[0], v[0]);
        dst[1] = 0;
        return 1;
    } else if (u[1] == 0) {
        cnt = lmmp_tailing_zeros_(u[0] | v[0]);
        u[0] = u[0] >> lmmp_tailing_zeros_(u[0]);
        goto gcd_1_2;
    } else if (v[1] == 0) {
        cnt = lmmp_tailing_zeros_(u[0] | v[0]);
        v[0] = v[0] >> lmmp_tailing_zeros_(v[0]);
        goto gcd_2_1;
    }
 
    if (u[0] == 0 && v[0] == 0) {
        dst[0] = lmmp_gcd_11_(u[1], v[1]); 
        dst[1] = 0;
        return 1;
    } else if (u[0] == 0) {
        u[0] = u[1] >> lmmp_tailing_zeros_(u[1]);
        cnt = lmmp_tailing_zeros_(v[0]);
        goto gcd_1_2;
    } else if (v[0] == 0) {
        v[0] = v[1] >> lmmp_tailing_zeros_(v[1]);
        cnt = lmmp_tailing_zeros_(u[0]);
        goto gcd_2_1;
    }
    cnt = lmmp_tailing_zeros_(u[0] | v[0]);
    k = lmmp_tailing_zeros_(u[0]);
 
    if (k > 0)
        _u128lshr(u, u, k);
    k = lmmp_tailing_zeros_(v[0]);
    if (k > 0)
        _u128lshr(v, v, k);
    while (!(u[0] == v[0] && u[1] == v[1])) {
        if (u[1] == 0 && v[1] != 0) goto gcd_1_2;
        if (v[1] == 0 && u[1] != 0) goto gcd_2_1;
        if (u[1] == 0 && v[1] == 0) goto gcd_1_1;
 
        if (_u128cmp(u, v)) {
            _u128sub(v, v, u);
            if (v[0] == 0) {
                v[0] = v[1] >> lmmp_tailing_zeros_(v[1]);
                goto gcd_2_1;
            } else if (v[1] == 0) {
                v[0] = v[0] >> lmmp_tailing_zeros_(v[0]);
                goto gcd_2_1;
            }
            k = lmmp_tailing_zeros_(v[0]);
            // k > 0
            _u128lshr(v, v, k);
        } else {
            _u128sub(u, u, v);
            if (u[0] == 0) {
                u[0] = u[1] >> lmmp_tailing_zeros_(u[1]);
                goto gcd_1_2;
            } else if (u[1] == 0) {
                u[0] = u[0] >> lmmp_tailing_zeros_(u[0]);
                goto gcd_1_2;
            }
            k = lmmp_tailing_zeros_(u[0]);
            // k > 0
            _u128lshr(u, u, k);
        }
    }
    dst[0] = u[0];
    dst[1] = u[1];
    if (cnt > 0)
        _u128lshl(dst, dst, cnt);
    return 2;
 
    gcd_1_2:   // [u,1]  ,  [v,2]
        k = lmmp_tailing_zeros_(v[0]);
        if (k > 0)
            _u128lshr(v, v, k);
        while (v[1] != 0) {
            _u128sub64(v, v, u[0]);
            if (v[1] == 0) 
                goto gcd_1_1;
            if (v[0] == 0) {
                v[0] = v[1] >> lmmp_tailing_zeros_(v[1]);
                goto gcd_1_1;
            }
            k = lmmp_tailing_zeros_(v[0]);
            // k > 0
            _u128lshr(v, v, k);
        }
        goto gcd_1_1;
 
    gcd_2_1:   // [u,2]  ,  [v,1]
        k = lmmp_tailing_zeros_(u[0]);
        if (k > 0)
            _u128lshr(u, u, k);
        while (u[1] != 0) {
            _u128sub64(u, u, v[0]);
            if (u[1] == 0)
                goto gcd_1_1;
            if (u[0] == 0) {
                u[0] = u[1] >> lmmp_tailing_zeros_(u[1]);
                goto gcd_1_1;
            }
            k = lmmp_tailing_zeros_(u[0]);
            // k > 0
            _u128lshr(u, u, k);
        }
        goto gcd_1_1;
 
    gcd_1_1:   // [u,1]  ,  [v,1]
        while (u[0] != v[0]) {
            if (u[0] > v[0]) {
                u[0] -= v[0];
                u[0] >>= lmmp_tailing_zeros_(u[0]);
            } else {
                v[0] -= u[0];
                v[0] >>= lmmp_tailing_zeros_(v[0]);
            }
        }
        dst[0] = u[0];
        dst[1] = 0;
        if (cnt > 0)
            _u128lshl(dst, dst, cnt);
        return dst[1] == 0 ?  1  :  2;
}

引用了 _u128cmp, _u128lshl, _u128lshr, _u128sub, _u128sub64, k, lmmp_gcd_11_(), lmmp_param_assert, lmmp_tailing_zeros_ , 以及 n.

被这些函数引用 lmmp_gcd_2_().

函数调用图:

这是这个函数的调用关系图:

lmmp_gcd_2_()

mp_size_t lmmp_gcd_2_	(	mp_ptr	dst,
		mp_srcptr	up,
		mp_size_t	un,
		mp_srcptr	vp
	)

计算两个无符号整数的最大公约数

参数

up	第一个无符号整数指针
un	第一个无符号整数的 limb 长度
vp	第二个无符号整数指针，长度为 2
dst	结果指针（长度至少为 2，两个 limb 都会进行写入，即使最高位可能为0）

警告

up!=NULL, un>2, vp!=NULL, vp[1]!=0, dst!=NULL, eqsep(dst,[up|vp])

返回

dst 的实际 limb 长度

在文件 gcd_2.c 第 162 行定义.

                                                                            {
    lmmp_param_assert(dst != NULL);
    lmmp_param_assert(up != NULL);
    lmmp_param_assert(vp != NULL);
    lmmp_param_assert(un > 2);
    lmmp_param_assert(vp[1] != 0);
    mp_limb_t u[2] = {vp[0], vp[1]};
    lmmp_mod_2_(up, un, u);
    if (u[1] == 0 && u[0] == 0) {
        dst[0] = vp[0];
        dst[1] = vp[1];
        return 2;
    } else {
        return lmmp_gcd_22_(dst, vp, u);
    }
}

引用了 lmmp_gcd_22_(), lmmp_mod_2_(), lmmp_param_assert , 以及 n.

函数调用图:

lmmp_gcd_basecase_()

mp_size_t lmmp_gcd_basecase_	(	mp_ptr	dst,
		mp_srcptr	up,
		mp_size_t	un,
		mp_srcptr	vp,
		mp_size_t	vn
	)

计算两个无符号整数的最大公约数（不建议使用此算法，更高版本可能被彻底弃用）

参数

dst	结果指针（长度至少为 min(un,vn)）
up	第一个无符号整数指针
un	第一个无符号整数的 limb 长度
vp	第二个无符号整数指针
vn	第二个无符号整数的 limb 长度

警告

up!=NULL, un>0, vp!=NULL, vn>0, eqsep(dst,[up|vp]), dst!=NULL

注解

朴素的辗转相除法，与Lehmer算法具有相似的渐进时间复杂度，但Lehmer算法绝大多数场合更加优秀

返回

dst 的实际 limb 长度

计算两个无符号整数的最大公约数（不建议使用此算法，更高版本可能被彻底弃用）

This file is part of LAMMP.

LAMMP is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License (LGPL) as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.

This program is distributed WITHOUT ANY WARRANTY.

See https://www.gnu.org/licenses/.

在文件 gcd_basecase.c 第 20 行定义.

                                                                                                 {
    lmmp_param_assert(un > 0 && vn > 0);
    if (un < vn) {
        LMMP_SWAP(up, vp, mp_srcptr);
        LMMP_SWAP(un, vn, mp_size_t);
    } else if (un == vn) {
        int cmp = lmmp_cmp_(up, vp, un);
        if (cmp < 0) {
            LMMP_SWAP(up, vp, mp_srcptr);
        } else if (cmp == 0) {
            lmmp_copy(dst, up, un);
            return un;
        }
    }
    TEMP_DECL;
#define an un
#define bn vn
    mp_ptr a = TALLOC_TYPE(an, mp_limb_t);
    mp_ptr b = TALLOC_TYPE(bn, mp_limb_t);
    lmmp_copy(a, up, an);
    lmmp_copy(b, vp, bn);
    while (bn > 0 || (bn == 1 && b[0] == 0)) {
        // dst = a % b;
        lmmp_div_(NULL, dst, a, an, b, bn);
        lmmp_copy(a, b, bn);
        an = bn;
        while (dst[bn - 1] == 0 && bn > 0) {
            --bn;
        }
        lmmp_copy(b, dst, bn);
    }
    lmmp_copy(dst, a, an);
    TEMP_FREE;
    return an;
#undef an
#undef bn
}

引用了 an, bn, lmmp_cmp_(), lmmp_copy, lmmp_div_(), lmmp_param_assert, LMMP_SWAP, n, TALLOC_TYPE, TEMP_DECL , 以及 TEMP_FREE.

函数调用图:

lmmp_gcd_lehmer_()

mp_size_t lmmp_gcd_lehmer_	(	mp_ptr	dst,
		mp_srcptr	up,
		mp_size_t	un,
		mp_srcptr	vp,
		mp_size_t	vn
	)

计算两个无符号整数的最大公约数（Lehmer算法）

参数

dst	结果指针（长度至少为 min(un,vn)）
up	第一个无符号整数指针
un	第一个无符号整数的 limb 长度
vp	第二个无符号整数指针
vn	第二个无符号整数的 limb 长度

警告

up!=NULL, un>0, vp!=NULL, vn>0, eqsep(dst,[up|vp]), dst!=NULL

返回

dst 的实际 limb 长度

在文件 gcd_lehmer.c 第 242 行定义.

                                                                                               {
    lmmp_param_assert(un > 0 && vn > 0);
    lmmp_param_assert(up != NULL && vp != NULL);
    lmmp_param_assert(dst != NULL);
 
    if (un < vn) {
        LMMP_SWAP(up, vp, mp_srcptr);
        LMMP_SWAP(un, vn, mp_size_t);
    } else if (un == vn) {
        int cmp = lmmp_cmp_(up, vp, un);
        if (cmp == 0) {
            lmmp_copy(dst, up, un);
            return un;
        } else if (cmp < 0) {
            LMMP_SWAP(up, vp, mp_srcptr);
        }
    }
    // u > v
 
    mp_gcd_lehmer_t M;
    slong x = 0, y = 0;
 
#define an un
#define bn vn
    TEMP_B_DECL;
    // [a,an+1] [b,bn+1]
    // A * a_old may overlow
    mp_ptr a = BALLOC_TYPE(an + 1, mp_limb_t);
    mp_ptr b = BALLOC_TYPE(bn + 1, mp_limb_t);
    lehmer_stack_t ms;
    mp_ptr temp = BALLOC_TYPE((an + 1) * 3, mp_limb_t);
    ms.tp = temp;
    ms.mp = temp + (an + 1);
    ms.np = temp + (an + 1) * 2;
 
    lmmp_copy(a, up, an);
    lmmp_copy(b, vp, bn);
 
    bool bzero = false;
    while (bzero == false) {
        if (an > 1 && bn == 1) {
            dst[0] = lmmp_gcd_1_(a, an, b[0]);
            return 1;
        } else if (an == 1 && bn == 1) {
            dst[0] = lmmp_gcd_11_(a[0], b[0]);
            return 1;
        } 
        // a > b
        lmmp_lehmer_extract_(a, an, b, bn, &x, &y);
        lmmp_gcd_lehmer_step_(x, y, &M);
 
        if (M.m21 == 0) {
            lmmp_div_(NULL, dst, a, an, b, bn);
            lmmp_copy(a, b, bn);
            an = bn;
            while (dst[bn - 1] == 0 && bn > 0) {
                --bn;
            }
            if (bn == 0)
                bzero = true;
            else 
                lmmp_copy(b, dst, bn);
        } else {
            bzero = lmmp_lehmer_mul_(a, &an, b, &bn, &M, &ms);
            if ((an < bn) || (an == bn && lmmp_cmp_(a, b, an) < 0)) {
                LMMP_SWAP(a, b, mp_ptr);
                LMMP_SWAP(an, bn, mp_size_t);
            }
        }
    }
    lmmp_copy(dst, a, an);
    TEMP_B_FREE;
    return an;
#undef an
#undef bn
}

引用了 an, BALLOC_TYPE, bn, lmmp_cmp_(), lmmp_copy, lmmp_div_(), lmmp_gcd_11_(), lmmp_gcd_1_(), lmmp_gcd_lehmer_step_(), lmmp_lehmer_extract_(), lmmp_lehmer_mul_(), lmmp_param_assert, LMMP_SWAP, n, TEMP_B_DECL , 以及 TEMP_B_FREE.

函数调用图:

lmmp_hyperfac_()

mp_size_t lmmp_hyperfac_	(	mp_ptr	dst,
		mp_bitcnt_t	bits,
		mp_size_t	rn,
		ushort	n
	)

计算hyper阶乘（k^k累乘至n）

参数

dst	结果指针
bits	hyper阶乘的2的因子数
rn	结果指针的 limb 长度
n	hyper阶乘的阶数

警告

dst!=NULL, rn>0

返回

dst 的实际 limb 长度

lmmp_hyperfac_size_()

mp_size_t lmmp_hyperfac_size_	(	ushort	n,
		mp_bitcnt_t *	bits
	)

计算hyper阶乘的 limb 缓冲区长度

参数

n	hyper阶乘的阶数
bits	被修改为 hyper阶乘的2的因子数

警告

bits!=NULL

返回

hyper阶乘的 limb 缓冲区长度

lmmp_is_prime_notrial_()

bool lmmp_is_prime_notrial_

(

ulong

n

)

判断素数（无试除法）

参数

n	待判断的数（建议为极有可能为素数的数）

注解

不进行试除法过滤，适用于判断已被小素数试除法过滤的数或强伪素数

警告

n>=2

返回

若 n 为素数，返回 true，否则返回 false

在文件 is_prime_ulong.c 第 293 行定义.

                                     {
    lmmp_param_assert(n > 1);
    if (n < 684630005672341) {
        ushort bases[2];
        bases[0] = 2;
        bases[1] = dj_base49[((0x3AC69A35UL * n) & 0xFFFFFFFFUL) >> 21] + 3;
        if (n % bases[0] == 0)
            return false;
        if (n % bases[1] == 0)
            return false;
 
        ulong one = 1;
        ulong m_1 = n - 1;
        ulong m_inv = lmmp_binvert_ulong_(n);
        m_inv = -m_inv;
        ulong R2 = mont63_R2(n);
        one = to_mont63(one, R2, n, m_inv);
        m_1 = to_mont63(m_1, R2, n, m_inv);
 
        ulong u = n - 1, t = 0;
        while (u % 2 == 0) u /= 2, ++t;
 
        if (miller_rabin_63(bases[0], t, u, n, m_inv, one, m_1))
            if (miller_rabin_63(bases[1], t, u, n, m_inv, one, m_1))
                return true;
            else
                return false;
        else
            return false;
    } else {
        ushort bases[3];
        ulong bbmask = dj_base64[((0x3AC69A35UL * n) & 0xFFFFFFFFUL) >> 18];
        bases[0] = 2;
        bases[1] = (bbmask & 0x8000) ? 26460 : 9375;
        bases[2] = (bbmask & 0x7FFF) + 3;
 
        if (n % bases[0] == 0)
            return false;
        if (n % bases[1] == 0)
            return false;
        if (n % bases[2] == 0)
            return false;
 
        ulong one = 1;
        ulong m_1 = n - 1;
        ulong m_inv = lmmp_binvert_ulong_(n);
        m_inv = -m_inv;
 
        if (n <= MONT63_MAX) {
            ulong R2 = mont63_R2(n);
            one = to_mont63(one, R2, n, m_inv);
            m_1 = to_mont63(m_1, R2, n, m_inv);
 
            ulong u = n - 1, t = 0;
            while (u % 2 == 0) u /= 2, ++t;
 
            if (miller_rabin_63(bases[0], t, u, n, m_inv, one, m_1))
                if (miller_rabin_63(bases[1], t, u, n, m_inv, one, m_1))
                    if (miller_rabin_63(bases[2], t, u, n, m_inv, one, m_1))
                        return true;
                    else
                        return false;
                else
                    return false;
            else
                return false;
        } else {
            ulong R2 = mont64_R2(n);
            one = to_mont64(one, R2, n, m_inv);
            m_1 = to_mont64(m_1, R2, n, m_inv);
 
            ulong u = n - 1, t = 0;
            while (u % 2 == 0) u /= 2, ++t;
 
            if (miller_rabin_64(bases[0], t, u, n, m_inv, one, m_1))
                if (miller_rabin_64(bases[1], t, u, n, m_inv, one, m_1))
                    if (miller_rabin_64(bases[2], t, u, n, m_inv, one, m_1))
                        return true;
                    else
                        return false;
                else
                    return false;
            else
                return false;
        }
    }
}

引用了 dj_base49, dj_base64, lmmp_binvert_ulong_(), lmmp_param_assert, miller_rabin_63(), miller_rabin_64(), MONT63_MAX, mont63_R2(), mont64_R2(), n, t, to_mont63() , 以及 to_mont64().

被这些函数引用 lmmp_is_prime_ulong_(), lmmp_next_prime_ulong_() , 以及 lmmp_prev_prime_ulong_().

函数调用图:

这是这个函数的调用关系图:

lmmp_is_prime_uint_()

bool lmmp_is_prime_uint_

(

uint

n

)

判断素数

参数

n

待判断的数

返回

若 n 为素数，返回 true，否则返回 false

判断素数

Hashed 2-bases for n < 684630005672341 (slightly more than 2^49) Hashed 3-bases for n < 2^64

Based on Steve Worley's 2^32 example: http://www.mersenneforum.org/showthread.php?t=12209 With a 3-base encoding idea from Bradley Berg.

Copyright 2014, Dana Jacobsen dana@.nosp@m.acm..nosp@m.org

在文件 is_prime_ulong.c 第 259 行定义.

                                 {
    if (n % 2 == 0 || n <= 1)
        return false;
    int judge = lmmp_is_prime_table_(n);
    if (judge == 0) {
        return false;
    } else if (judge == 1) {
        return true;
    }
 
    if (trial_div35711(n))
        return false;
 
    ushort bases[2];
    bases[0] = 2;
    bases[1] = dj_base49[((0x3AC69A35UL * n) & 0xFFFFFFFFUL) >> 21] + 3;
    if (n % bases[0] == 0)
        return false;
    if (n % bases[1] == 0)
        return false;
 
    ulong u = n - 1, t = 0;
    while (u % 2 == 0) u /= 2, ++t;
 
    _udiv64_t binv = _udiv64_gen(n);
    if (miller_rabin_32(bases[0], t, u, n, &binv))
        if (miller_rabin_32(bases[1], t, u, n, &binv))
            return true;
        else
            return false;
    else
        return false;
}

引用了 _udiv64_gen(), dj_base49, lmmp_is_prime_table_(), miller_rabin_32(), n, t , 以及 trial_div35711().

函数调用图:

lmmp_is_prime_ulong_()

bool lmmp_is_prime_ulong_

(

ulong

n

)

判断素数

参数

n

待判断的数

注解

如果 n 的实际值小于2^32，此函数不会调用 lmmp_is_prime_uint_， 如果你可以保证 n 的实际值小于2^32，使用 lmmp_is_prime_uint_ 将会更快

返回

若 n 为素数，返回 true，否则返回 false

在文件 is_prime_ulong.c 第 381 行定义.

                                   {
    if (n % 2 == 0 || n <= 1)
        return false;
    int judge = lmmp_is_prime_table_(n);
    if (judge == 0) {
        return false;
    } else if (judge == 1) {
        return true;
    }
    if (trial_div35711(n))
        return false;
    return lmmp_is_prime_notrial_(n);
}

引用了 lmmp_is_prime_notrial_(), lmmp_is_prime_table_(), n , 以及 trial_div35711().

函数调用图:

lmmp_mulmod_ulong_()

ulong lmmp_mulmod_ulong_	(	ulong	a,
		ulong	b,
		ulong	mod,
		ulongp	q
	)

计算两个无符号整数的乘积，对mod取模，商放入 q 中

参数

a	第一个无符号整数
b	第二个无符号整数
q	商的结果指针
mod	取模数

警告

a < mod, b < mod, q!=NULL

返回

余数

lmmp_multinomial_()

mp_size_t lmmp_multinomial_	(	mp_ptr	dst,
		mp_size_t	rn,
		uint	n,
		const uintp	r,
		uint	m
	)

计算多项式系数

参数

dst	结果指针
rn	结果指针的 limb 长度
n	r[i] 的总和
r	需要计算的系数的数组
m	系数的个数

警告

m>1, n>0

注解

多项式系数为 ( r1+r2+...+rm )! / ( r1! * r2! * ... * rm!)

返回

返回 dst 的实际 limb 长度

lmmp_multinomial_size_()

mp_size_t lmmp_multinomial_size_	(	const uintp	r,
		uint	m,
		ulong *	n
	)

计算多项式系数的 limb 缓冲区长度

参数

r	需要计算的系数的数组
m	系数的个数
n	输出变量，将会被修改为 r[i] 的总和，即r1+r2+...+rm

返回

多项式系数的 limb 缓冲区长度（比实际长度多 1-2 个 limb）

注解

多项式系数为 ( r1+r2+...+rm )! / ( r1! * r2! * ... * rm!)

警告

我们使用 ulong* n 来同时计算 r[i] 的总和，因为 n 可能超过 0xffffffff。 我们预计算 n，这不仅可以作为后续多项式系数函数的参数传入。 同时也请调用者注意判断 n 是否超过了 0xffffffff 这是 lmmp_multinomial_ 函数的限制。

lmmp_nCr_()

mp_size_t lmmp_nCr_	(	mp_ptr	dst,
		mp_bitcnt_t	bits,
		mp_size_t	rn,
		uint	n,
		uint	r
	)

计算 nCr 组合数 ( nCr = n! / (r!(n-r)!) )

参数

dst	结果指针
bits	nCr 的2的因子数
rn	结果指针的 limb 长度
n	组合数的总数
r	组合数的选择数

返回

返回 dst 的实际 limb 长度

警告

r<=n/2, n<=0xffffffff, dst!=NULL, rn>0

lmmp_nCr_size_()

mp_size_t lmmp_nCr_size_	(	uint	n,
		uint	r,
		mp_bitcnt_t *	bits
	)

计算 nCr 组合数的 limb 缓冲区长度

参数

n	组合数的总数
r	组合数的选择数
bits	被修改为 nCr 的2的因子数

警告

r<=n/2, bits!=NULL

返回

nCr 组合数的 limb 缓冲区长度（比实际长度多 1-2 个 limb）

lmmp_next_prime_ulong_()

ulong lmmp_next_prime_ulong_

(

ulong

n

)

大于n的下一个素数

参数

n

起始点（不含）

警告

如果 n 大于等于ulong可表示最大的质数，则返回ulong_max

返回

大于n的下一个素数

在文件 is_prime_ulong.c 第 455 行定义.

                                      {
    if (n < prime_short_table[PRIME_SHORT_TABLE_SIZE - 1]) {
        ushort idx = lmmp_prime_cnt16_(n);
        return prime_short_table[idx];
    } else if (n >= ULONG_PRIME_MAX) {
        return MP_ULONG_MAX;
    } else {
        n += (n % 2 == 0) ? 1 : 2;
        while (1) {
            if (trial_div35711(n)
             || trial_div13(n)
             || trial_div17(n)
             || trial_div19(n)
             || trial_div23(n)
             || trial_div29(n)
             || trial_div31(n)
             || trial_div37(n)
             || trial_div41(n)) {
                n += 2; 
            } else {
                if (lmmp_is_prime_notrial_(n)) {
                    return n;
                } else {
                    n += 2;
                }
            }
        }
    }
}

引用了 lmmp_is_prime_notrial_(), lmmp_prime_cnt16_(), MP_ULONG_MAX, n, prime_short_table, PRIME_SHORT_TABLE_SIZE, trial_div13(), trial_div17(), trial_div19(), trial_div23(), trial_div29(), trial_div31(), trial_div35711(), trial_div37(), trial_div41() , 以及 ULONG_PRIME_MAX.

函数调用图:

lmmp_nPr_()

mp_size_t lmmp_nPr_	(	mp_ptr	dst,
		mp_bitcnt_t	bits,
		mp_size_t	rn,
		ulong	n,
		ulong	r
	)

计算 nPr 排列数 ( nPr = n! / (n-r)! )

参数

dst	结果指针
bits	nPr 的2的因子数
rn	结果指针的 limb 长度
n	排列数的总数
r	排列数的选择数

警告

n>=r, dst!=NULL, rn>0

返回

返回 dst 的实际 limb 长度

lmmp_nPr_size_()

mp_size_t lmmp_nPr_size_	(	ulong	n,
		ulong	r,
		mp_bitcnt_t *	bits
	)

计算 nPr 排列数的 limb 缓冲区长度

参数

n	排列数的总数
r	排列数的选择数
bits	被修改为 nPr 的2的因子数

警告

r<=n, bits!=NULL

返回

nPr 排列数的 limb 缓冲区长度（比实际长度多）

被这些函数引用 lmmp_odd_nCr_uint_() , 以及 pow_nPr_().

这是这个函数的调用关系图:

lmmp_nthroot_ulong_()

ulong lmmp_nthroot_ulong_	(	ulong	n,
		ulong	root
	)

计算 floor(n^(1/root))

参数

n	被开方数
root	开方次数

返回

floor(n^(1/root))

注解

root=0时，返回0

在文件 nthroot.c 第 55 行定义.

                                               {
    ulong x, currval, base, upper_limit;
 
    if (n == 0 || root == 0)
        return 0;
    if (root == 1)
        return n;
    if (root == 2)
        return lmmp_sqrt_ulong_(n);
    if (root == 3)
        return lmmp_cbrt_ulong_(n);
 
    if (root >= LIMB_BITS || n < (1ULL << root))
        return 1;
 
    /* n <= upper_limit^root */
    upper_limit = max_base[root - 4];
 
    if (upper_limit == 2)
        return upper_limit;
 
    /* upper_limit = 2 for root >= 41 */
    lmmp_debug_assert(root <= 40);
 
    if (root == 4)
        x = sqrt(sqrt(n));
    else
        x = expf(inv_table[root - 5] * logf(n));
 
    base = x;
 
    if (base >= upper_limit)
        base = upper_limit - 1;
 
    currval = pow_n(base, root);
    if (currval == n)
        return base;
 
    while (currval <= n) {
        base++;
        currval = pow_n(base, root);
        if (base == upper_limit)
            break;
    }
 
    while (currval > n) {
        base--;
        currval = pow_n(base, root);
    }
 
    return base;
}

引用了 inv_table, LIMB_BITS, lmmp_cbrt_ulong_(), lmmp_debug_assert, lmmp_sqrt_ulong_(), max_base, n , 以及 pow_n().

函数调用图:

lmmp_odd_factorial_uint_()

mp_size_t lmmp_odd_factorial_uint_	(	mp_ptr	dst,
		mp_size_t	rn,
		uint	n
	)

计算 n! 阶乘的奇数部分

参数

dst	结果指针
rn	结果指针的 limb 长度（factorial_size_ 函数的返回值 - bits/LIMB_BITS）
n	阶乘的阶数

警告

n>0xffff, dst!=NULL, rn>0

返回

返回 dst 的实际 limb 长度

被这些函数引用 lmmp_2factorial_() , 以及 lmmp_odd_factorial_().

这是这个函数的调用关系图:

lmmp_odd_nCr_uint_()

mp_size_t lmmp_odd_nCr_uint_	(	mp_ptr	dst,
		mp_size_t	rn,
		uint	n,
		uint	r
	)

计算 nCr 组合数的奇数部分

参数

dst	结果指针
rn	结果指针的 limb 长度
n	组合数的总数
r	组合数的选择数

返回

返回 dst 的实际 limb 长度

警告

r<=n/2, 0xffff<n, dst!=NULL, rn>0

lmmp_odd_nCr_ushort_()

mp_size_t lmmp_odd_nCr_ushort_	(	mp_ptr	dst,
		mp_size_t	rn,
		uint	n,
		uint	r
	)

计算 nCr 组合数的奇数部分

参数

dst	结果指针
rn	结果指针的 limb 长度
n	组合数的总数
r	组合数的选择数

返回

返回 dst 的实际 limb 长度

警告

r<=n/2, n<=0xffff, dst!=NULL, rn>0

lmmp_odd_nPr_uint_()

mp_size_t lmmp_odd_nPr_uint_	(	mp_ptr	dst,
		mp_size_t	rn,
		ulong	n,
		ulong	r
	)

计算 nPr 排列数的奇数部分

参数

dst	结果指针
rn	结果指针的 limb 长度（nPr_size_ 函数的返回值 - bits/LIMB_BITS）
n	排列数的总数
r	排列数的选择数

警告

0xffffffff>=n>=r, dst!=NULL, rn>0

返回

返回 dst 的实际 limb 长度

被这些函数引用 _odd_nPr_().

这是这个函数的调用关系图:

lmmp_odd_nPr_ulong_()

mp_size_t lmmp_odd_nPr_ulong_	(	mp_ptr	dst,
		mp_size_t	rn,
		ulong	n,
		ulong	r
	)

计算 nPr 排列数的奇数部分

参数

dst	结果指针
rn	结果指针的 limb 长度（nPr_size_ 函数的返回值 - bits/LIMB_BITS）
n	排列数的总数
r	排列数的选择数

警告

n>=r, dst!=NULL, rn>0

返回

返回 dst 的实际 limb 长度

lmmp_odd_nPr_ushort_()

mp_size_t lmmp_odd_nPr_ushort_	(	mp_ptr	dst,
		mp_size_t	rn,
		ulong	n,
		ulong	r
	)

计算 nPr 排列数的奇数部分

参数

dst	结果指针
rn	结果指针的 limb 长度（nPr_size_ 函数的返回值 - bits/LIMB_BITS）
n	排列数的总数
r	排列数的选择数

警告

0xffff>=n>=r, dst!=NULL, rn>0

返回

返回 dst 的实际 limb 长度

被这些函数引用 _odd_nPr_(), lmmp_2factorial_(), lmmp_factorial_(), lmmp_odd_factorial_(), lmmp_odd_multinomial_ushort_() , 以及 lmmp_odd_nCr_ushort_().

这是这个函数的调用关系图:

lmmp_pow_()

mp_size_t lmmp_pow_	(	mp_ptr	dst,
		mp_size_t	rn,
		mp_srcptr	base,
		mp_size_t	n,
		ulong	exp
	)

计算大整数幂 [dst,rn] = [base,n] ^ exp

参数

dst	结果指针
rn	结果指针的 limb 长度
base	底数指针
n	底数指针的 limb 长度
exp	指数

警告

n>0, base[n-1]!=0, sep(dst,base), exp>0

返回

返回 dst 的实际 limb 长度

lmmp_pow_1_()

mp_size_t lmmp_pow_1_	(	mp_ptr	dst,
		mp_size_t	rn,
		mp_limb_t	base,
		ulong	exp
	)

计算幂次方 [dst,rn] = [base,1] ^ exp

参数

dst	结果指针
base	底数
exp	指数

警告

base>=1, exp>0

返回

返回 dst 的实际 limb 长度

被这些函数引用 lmmp_pow_().

这是这个函数的调用关系图:

lmmp_pow_1_size_()

mp_size_t lmmp_pow_1_size_	(	mp_limb_t	base,
		ulong	exp
	)

计算幂次方需要的limb缓冲区长度 base ^ exp

参数

base	底数
exp	指数

警告

exp>0, base>=1

返回

返回值为 base^exp 需要的 limb 缓冲区长度（比实际长度多）

计算幂次方需要的limb缓冲区长度 base ^ exp

This file is part of LAMMP.

LAMMP is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License (LGPL) as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.

This program is distributed WITHOUT ANY WARRANTY.

See https://www.gnu.org/licenses/.

在文件 pow.c 第 25 行定义.

                                                      {
    lmmp_param_assert(base >= 1);
    lmmp_param_assert(exp > 0);
    if (base == 1) {
        return 1;
    } else if (exp <= 2) {
        return 3;
    } else if (exp <= MP_UINT_MAX) {
        /*
        base = b * 2^base_tz
        */
        slong base_tz = lmmp_limb_bits_(base);
        uint32_t b;
        if (base_tz < 32) {
            b = base << (32 - base_tz);
        } else {
            b = base >> (base_tz - 32);
        }
        base_tz = base_tz - 32;
        base_tz *= exp;
        mp_size_t rn = xlog2n_ceil(exp, b);
        rn += base_tz;
        rn = (rn + LIMB_BITS - 1) / LIMB_BITS;
        return rn + 2;
    } else {
        /*
        base = b * 2^base_tz
        */
        slong base_tz = lmmp_limb_bits_(base);
        uint32_t b;
        if (base_tz < 32) {
            b = base << (32 - base_tz);
        } else {
            b = base >> (base_tz - 32);
        }
        base_tz = base_tz - 32;
        base_tz *= exp;
 
        mp_size_t rn;
        /*
        exp = exp' * 2^bits
        exp*log2(base) = exp*log2(b*2^base_tz)
                        = exp*log2(b) + exp*base_tz
                        = exp'*log2(b)*2^bits + exp*base_tz
        */
        mp_bitcnt_t bits = lmmp_limb_bits_(exp);
        bits -= 32;
        exp >>= bits;
        exp++;
        rn = xlog2n_ceil(exp, b) << bits;
        rn += base_tz;
        rn = (rn + LIMB_BITS - 1) / LIMB_BITS;
        return rn + 2;
    }
}

引用了 LIMB_BITS, lmmp_limb_bits_, lmmp_param_assert, MP_UINT_MAX, n , 以及 xlog2n_ceil().

被这些函数引用 lmmp_arith_seqprod_size_(), lmmp_nCr_size_(), lmmp_nPr_size_(), lmmp_pow_size_() , 以及 pow_nPr_().

函数调用图:

这是这个函数的调用关系图:

lmmp_pow_basecase_()

mp_size_t lmmp_pow_basecase_	(	mp_ptr	dst,
		mp_size_t	rn,
		mp_srcptr	base,
		mp_size_t	n,
		ulong	exp
	)

计算奇数次幂算法 [dst,rn] = [base,n] ^ exp

参数

dst	结果指针
rn	结果指针的 limb 缓冲区长度
base	底数指针
n	底数指针的 limb 长度
exp	指数

警告

n>0, base[n-1]!=0, sep(dst,base), [base,n]>1, exp>=3, exp%2==1

返回

返回 dst 的实际 limb 长度

被这些函数引用 lmmp_pow_().

这是这个函数的调用关系图:

lmmp_pow_size_()

mp_size_t lmmp_pow_size_	(	mp_srcptr	base,
		mp_size_t	n,
		ulong	exp
	)

计算幂次方需要的limb缓冲区长度 [base,n] ^ exp

参数

base	底数指针
n	底数 limb 长度
exp	指数

警告

n>0, base[n-1]!=0, [base,n]>1

返回

返回值为 [base,n]^exp 需要的 limb 缓冲区长度（比实际长度多）

在文件 pow.c 第 81 行定义.

                                                                 {
    lmmp_param_assert(n > 0);
    lmmp_param_assert(base[n - 1] != 0);
    if (n == 1) {
        return lmmp_pow_1_size_(base[0], exp);
    }
    if (exp == 1) {
        return n;
    } else if (exp == 2) {
        return n * 2;
    } else {
        /*
        base = b * 2^base_tz
        */
        mp_bitcnt_t base_tz = lmmp_limb_bits_(base[n - 1]);
        uint32_t b;
        if (base_tz < 32) {
            b = base[n - 1] << (32 - base_tz);
            b |= (base[n - 2] >> (LIMB_BITS - 32 + base_tz));
            base_tz = (n - 2) * LIMB_BITS + LIMB_BITS - 32 + base_tz;
        } else if (base_tz == 32) {
            b = base[n - 1];
            base_tz = (n - 1) * LIMB_BITS;
        } else {
            b = base[n - 1] >> (base_tz - 32);
            base_tz = (n - 1) * LIMB_BITS + base_tz - 32;
        }
 
        mp_size_t rn;
        if (exp <= MP_UINT_MAX) {
            rn = exp * base_tz;
            rn += xlog2n_ceil(exp, b);
        } else {
            /*
            exp = exp' * 2^bits
            exp*log2(base) = exp*log2(b*2^base_tz)
                           = exp*log2(b) + exp*base_tz
                           = exp'*log2(b)*2^bits + exp*base_tz
            */
            mp_bitcnt_t bits = lmmp_limb_bits_(exp);
            rn = exp * base_tz;
            bits -= 32;
            exp >>= bits;
            exp++;
            rn += xlog2n_ceil(exp, b) << bits;
        }
        rn = (rn + LIMB_BITS - 1) / LIMB_BITS;
        return rn + 2;
    }
}

引用了 LIMB_BITS, lmmp_limb_bits_, lmmp_param_assert, lmmp_pow_1_size_(), MP_UINT_MAX, n , 以及 xlog2n_ceil().

被这些函数引用 lmmp_pow_().

函数调用图:

这是这个函数的调用关系图:

lmmp_pow_win2_()

mp_size_t lmmp_pow_win2_	(	mp_ptr	dst,
		mp_size_t	rn,
		mp_srcptr	base,
		mp_size_t	n,
		ulong	exp
	)

计算幂次方2比特窗口快速幂算法 [dst,rn] = [base,n] ^ exp

参数

dst	结果指针
rn	结果指针的 limb 长度
base	底数指针
n	底数指针的 limb 长度
exp	指数

警告

n>0, base[n-1]!=0, sep(dst,base), exp>0

返回

返回 dst 的实际 limb 长度

被这些函数引用 lmmp_pow_().

这是这个函数的调用关系图:

lmmp_powmod_uint_()

uint lmmp_powmod_uint_	(	uint	base,
		ulong	exp,
		uint	mod
	)

计算 base^exp 对 mod 取模

参数

base	底数
exp	指数
mod	模数

警告

base < mod, mod > 1

返回

base^exp 对 mod 取模的结果

在文件 is_prime_ulong.c 第 108 行定义.

                                                       {
    lmmp_param_assert(mod > base);
    lmmp_param_assert(mod > 1);
    ulong dst = 1;
    ulong b = base;
    _udiv64_t binv = _udiv64_gen(mod);
    ulong q;
    while (1) {
        if (exp & 1) {
            dst *= b;
            q = _udiv64by64_q_preinv(dst, &binv);
            dst -= q * mod;
        }
        exp >>= 1;
        if (exp == 0)
            break;
        b *= b;
        q = _udiv64by64_q_preinv(b, &binv);
        b -= q * mod;
    }
    return dst;
}

引用了 _udiv64_gen(), _udiv64by64_q_preinv(), lmmp_param_assert , 以及 n.

被这些函数引用 lmmp_powmod_ulong_().

函数调用图:

这是这个函数的调用关系图:

lmmp_powmod_ulong_()

ulong lmmp_powmod_ulong_	(	ulong	base,
		ulong	exp,
		ulong	mod
	)

计算 base^exp 对 mod 取模

参数

base	底数
exp	指数
mod	模数

警告

base < mod, mod > 1

返回

base^exp 对 mod 取模的结果

在文件 is_prime_ulong.c 第 131 行定义.

                                                           {
    if (mod <= MP_UINT_MAX)
        return lmmp_powmod_uint_(base, exp, mod);
    else if (mod <= MONT63_MAX) {
        ulong R2 = mont63_R2(mod);
        ulong m_inv = lmmp_binvert_ulong_(mod);
        m_inv = -m_inv;
        ulong dst = to_mont63(1, R2, mod, m_inv);
        base = to_mont63(base, R2, mod, m_inv);
        while (1) {
            if (exp & 1)
                dst = mont63_mul(dst, base, mod, m_inv);
            exp >>= 1;
            if (exp == 0)
                break;
            base = mont63_mul(base, base, mod, m_inv);
        }
        return from_mont63(dst, mod, m_inv);
    } else {
        ulong R2 = mont64_R2(mod);
        ulong m_inv = lmmp_binvert_ulong_(mod);
        m_inv = -m_inv;
        ulong dst = to_mont64(1, R2, mod, m_inv);
        base = to_mont64(base, R2, mod, m_inv);
        while (1) {
            if (exp & 1)
                dst = mont64_mul(dst, base, mod, m_inv);
            exp >>= 1;
            if (exp == 0)
                break;
            base = mont64_mul(base, base, mod, m_inv);
        }
        return from_mont64(dst, mod, m_inv);
    }
}

引用了 from_mont63(), from_mont64(), lmmp_binvert_ulong_(), lmmp_powmod_uint_(), MONT63_MAX, mont63_mul(), mont63_R2(), mont64_mul(), mont64_R2(), MP_UINT_MAX, n, to_mont63() , 以及 to_mont64().

函数调用图:

lmmp_prev_prime_ulong_()

ulong lmmp_prev_prime_ulong_

(

ulong

n

)

小于等于n的上一个素数

参数

n

起始点（含）

警告

如果 n 小于2，则返回 0

返回

小于等于n的上一个素数，如果n恰好为素数，则返回 n

在文件 is_prime_ulong.c 第 485 行定义.

                                      {
    if (n < ULONG_PRIME_MIN) {
        return 0;
    } else if (n < PRIME_SHORT_TABLE_N) {
        ushort idx = lmmp_prime_cnt16_(n);
        return prime_short_table[idx - 1];
    } else {
        n -= (n % 2 == 0) ? 1 : 0;
        while (1) {
            if (trial_div35711(n)
             || trial_div13(n)
             || trial_div17(n)
             || trial_div19(n)
             || trial_div23(n)
             || trial_div29(n)
             || trial_div31(n)
             || trial_div37(n)
             || trial_div41(n)) {
                n -= 2;
            } else {
                if (lmmp_is_prime_notrial_(n)) {
                    return n;
                } else {
                    n -= 2;
                }
            }
        }
    }
}

引用了 lmmp_is_prime_notrial_(), lmmp_prime_cnt16_(), n, prime_short_table, PRIME_SHORT_TABLE_N, trial_div13(), trial_div17(), trial_div19(), trial_div23(), trial_div29(), trial_div31(), trial_div35711(), trial_div37(), trial_div41() , 以及 ULONG_PRIME_MIN.

函数调用图:

lmmp_primefac_()

mp_size_t lmmp_primefac_	(	mp_ptr	dst,
		mp_size_t	rn,
		uint	n
	)

计算质数阶乘（不超过n的质数累乘）

参数

dst	结果指针
rn	结果指针的 limb 长度
n	质数阶乘的阶数

警告

dst!=NULL, rn>0

返回

dst 的实际 limb 长度

lmmp_primefac_size_()

mp_size_t lmmp_primefac_size_

(

uint

n

)

计算质数阶乘的 limb 缓冲区长度

参数

n

质数阶乘的阶数

返回

质数阶乘的 limb 缓冲区长度

在文件 factorial_extra.c 第 435 行定义.

                                      {
    // Chebyshev's estimate
    // pn# < exp(1.000028*n)
    const double log2 = 0.69314718055994531;  // log(2)
    double l = 1.000028 * (double)n;
    l /= log2;
    mp_size_t rn = (mp_size_t)ceil(l);
    rn = (rn + LIMB_BITS - 1) / LIMB_BITS + 2;  // more two limbs
    return rn;
}

引用了 l, LIMB_BITS , 以及 n.

lmmp_remove_()

mp_size_t lmmp_remove_	(	mp_ptr	np,
		mp_size_t *	nn,
		mp_srcptr	dp,
		mp_size_t	dn
	)

除去[np,nn]中的[dp,dn]的因子

参数

np	被除数指针，将会被修改为除去因子后的数
nn	被除数指针的 limb 长度，将会被修改除去因子后的长度
dp	除数指针
dn	除数指针的 limb 长度

警告

np!=NULL, nn>0, dp!=NULL, dn>0

注解

如果[np,nn]能被[dp,dn]整除，则[np,nn]将被修改为除去因子后的数，nn将被修改为除去因子后的长度。 如果不能被整除，则[np,nn]保持不变，并返回0。

返回

[np,nn]中被[dp,dn]除去的因子的个数，如果不能被整除，则返回0

lmmp_sqrt_ulong_()

ulong lmmp_sqrt_ulong_

(

ulong

a

)

计算算术平方根 floor(sqrt(a))

参数

a

被开方数

返回

floor(sqrt(a))

计算算术平方根 floor(sqrt(a))

This file is part of LAMMP.

LAMMP is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License (LGPL) as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.

This program is distributed WITHOUT ANY WARRANTY.

See https://www.gnu.org/licenses/.

在文件 nthroot.c 第 21 行定义.

                                {
    ulong is;
 
    is = (ulong)sqrt((double)a);
 
    is -= (is * is > a);
    if (is == (1ULL << 32))
        is--;
    return is;
}

引用了 n.

被这些函数引用 lmmp_nthroot_ulong_().

这是这个函数的调用关系图:

lmmp_superfac_()

mp_size_t lmmp_superfac_	(	mp_ptr	dst,
		mp_bitcnt_t	bits,
		mp_size_t	rn,
		ushort	n
	)

计算super阶乘（k!累乘至n）

参数

dst	结果指针
bits	super阶乘的2的因子数
rn	结果指针的 limb 长度
n	super阶乘的阶数

警告

dst!=NULL, rn>0

返回

dst 的实际 limb 长度

lmmp_superfac_size_()

mp_size_t lmmp_superfac_size_	(	ushort	n,
		mp_bitcnt_t *	bits
	)

计算super阶乘的 limb 缓冲区长度

参数

n	super阶乘的阶数
bits	被修改为 super阶乘的2的因子数

警告

bits!=NULL

返回

super阶乘的 limb 缓冲区长度

lmmp_trialdiv_()

ushortp lmmp_trialdiv_	(	mp_srcptr	np,
		mp_size_t	nn,
		ushort	N,
		ushort *	rn
	)

试除法

参数

num	被除数
nn	被除数的 limb 长度
N	试除法尝试的质数最大值
rn	结果指针的 limb 长度

警告

num!=NULL, nn>0, N>2, rn!=NULL

注解

试除法尝试从 2-N 中所有质数进行试除，如果能整除则会插入到返回结果数组中，没有整除的则会返回 NULL。 结果指针请使用 lmmp_free() 函数进行释放。

返回

结果指针，返回不超过N，且能整除[np,nn]的素数（从小到大排列），若没有能够整除的素数，则返回NULL

lmmp_u16_pow_1_()

mp_size_t lmmp_u16_pow_1_	(	mp_ptr	dst,
		mp_size_t	rn,
		ulong	base,
		ulong	exp
	)

计算幂次方 [dst,rn] = [base,1] ^ exp

参数

dst	结果指针
rn	结果指针的 limb 长度
base	底数（16位无符号整数）
exp	指数

警告

0<base<=0xffff, exp>0

返回

返回 dst 的实际 limb 长度

被这些函数引用 _odd_pow_().

这是这个函数的调用关系图:

lmmp_u32_pow_1_()

mp_size_t lmmp_u32_pow_1_	(	mp_ptr	dst,
		mp_size_t	rn,
		ulong	base,
		ulong	exp
	)

计算幂次方 [dst,rn] = [base,1] ^ exp

参数

dst	结果指针
rn	结果指针的 limb 长度
base	底数（32位无符号整数）
exp	指数

警告

0<base<=2^32-1, exp>0

返回

返回 dst 的实际 limb 长度

被这些函数引用 _odd_pow_().

这是这个函数的调用关系图:

lmmp_u4_pow_1_()

mp_size_t lmmp_u4_pow_1_	(	mp_ptr	dst,
		mp_size_t	rn,
		ulong	base,
		ulong	exp
	)

计算幂次方 [dst,rn] = [base,1] ^ exp

参数

dst	结果指针
rn	结果指针的 limb 长度
base	底数（4位无符号整数）
exp	指数

警告

1<=base<=0xf, exp>0

返回

返回 dst 的实际 limb 长度

被这些函数引用 _odd_pow_().

这是这个函数的调用关系图:

lmmp_u64_pow_1_()

mp_size_t lmmp_u64_pow_1_	(	mp_ptr	dst,
		mp_size_t	rn,
		ulong	base,
		ulong	exp
	)

计算幂次方 [dst,rn] = [base,1] ^ exp

参数

dst	结果指针
rn	结果指针的 limb 长度
base	底数（64位无符号整数）
exp	指数

警告

2^32<=base<=2^64-1, exp>0

返回

返回 dst 的实际 limb 长度

lmmp_u8_pow_1_()

mp_size_t lmmp_u8_pow_1_	(	mp_ptr	dst,
		mp_size_t	rn,
		ulong	base,
		ulong	exp
	)

计算幂次方 [dst,rn] = [base,1] ^ exp

参数

dst	结果指针
rn	结果指针的 limb 长度
base	底数（8位无符号整数）
exp	指数

警告

0<base<=0xff, exp>0

返回

返回 dst 的实际 limb 长度

被这些函数引用 _odd_pow_().

这是这个函数的调用关系图: