numth.h 文件参考

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

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

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

浏览源代码.

宏定义
#define	INLINE_   static inline

类型定义
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
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_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_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下的逆元
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)
	计算 [numa,na] 在B^n 下的逆元
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算法）
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 缓冲区长度
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
static 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
static 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_remove_ (mp_ptr np, mp_size_t *nn, mp_srcptr dp, mp_size_t dn)
	除去[np,nn]中的[dp,dn]的因子
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

宏定义说明

INLINE_

#define INLINE_   static inline

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

类型定义说明

schar

typedef int8_t schar

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

scharp

typedef int8_t* scharp

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

sint

typedef int32_t sint

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

sintp

typedef int32_t* sintp

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

slong

typedef int64_t slong

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

slongp

typedef int64_t* slongp

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

sshort

typedef int16_t sshort

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

sshortp

typedef int16_t* sshortp

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

uchar

typedef uint8_t uchar

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

ucharp

typedef uint8_t* ucharp

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

uint

typedef uint32_t uint

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

uintp

typedef uint32_t* uintp

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

ulong

typedef uint64_t ulong

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

ulongp

typedef uint64_t* ulongp

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

ushort

typedef uint16_t ushort

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

ushortp

typedef uint16_t* ushortp

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

函数说明

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

返回

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

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

                                                           {
    double x_m = (double)x / m;
    double log_l = (n + 1) * log(m) + lgamma(x_m + n + 1) - lgamma(x_m);
    double log2_l = log_l / LOG2_;
    mp_size_t rn = ceil(log2_l / LIMB_BITS) + 2; /* more two limbs */
    return rn;
}

引用了 LIMB_BITS , 以及 LOG2_.

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 第 47 行定义.

                                                 {
    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 , 以及 xn.

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

函数调用图:

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

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_binvert_uint_()

uint lmmp_binvert_uint_

(

uint

a

)

计算 a 在2^32下的逆元

参数

a

待求逆元

警告

a%2==1

返回

逆元

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

                                {
    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 */
    return r;
}

引用了 binv_tab , 以及 lmmp_param_assert.

lmmp_binvert_ulong_()

ulong lmmp_binvert_ulong_

(

ulong

a

)

计算 a 在2^64下的逆元

参数

a

待求逆元

警告

a%2==1

返回

逆元

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

                                   {
    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.

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

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

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的因子数

返回

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

lmmp_gcd_11_()

mp_limb_t lmmp_gcd_11_	(	mp_limb_t	u,
		mp_limb_t	v
	)

计算 [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)

返回

dst 的实际 limb 长度

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

参数

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

返回

最大公约数

警告

u!=0, v!=0

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

                                                 {
    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_().

被这些函数引用 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 第 30 行定义.

                                                                   {
    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.

被这些函数引用 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 长度

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

                                                               {
    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_().

被这些函数引用 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 第 146 行定义.

                                                                            {
    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.

函数调用图:

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 长度

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

                                                                                                 {
    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, 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 第 233 行定义.

                                                                                               {
    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, mp_gcd_lehmer_t::m21, lehmer_stack_t::mp, lehmer_stack_t::np, TEMP_B_DECL, TEMP_B_FREE , 以及 lehmer_stack_t::tp.

函数调用图:

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 第 249 行定义.

                                 {
    if (n < PRIME_SHORT_TABLE_N) {
        return lmmp_is_prime_table_(n);
    }
    if (n % 2 == 0)
        return false;
    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(), PRIME_SHORT_TABLE_N , 以及 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 第 365 行定义.

                                   {
    if (n < PRIME_SHORT_TABLE_N) {
        return lmmp_is_prime_table_(n);
    }
    if (n % 2 == 0)
        return false;
    if (trial_div35711(n))
        return false;
    return lmmp_is_prime_notrial_(n);
}

引用了 lmmp_is_prime_notrial_(), lmmp_is_prime_table_(), PRIME_SHORT_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

在文件 nCr.c 第 231 行定义.

                                                                                {
    lmmp_debug_assert(r <= (n / 2));
    mp_size_t shw = bits / LIMB_BITS;
    bits %= LIMB_BITS;
    lmmp_zero(dst, shw);
    if (n <= NCR_SHORT_LIMIT)
        rn = lmmp_odd_nCr_ushort_(dst + shw, rn - shw, n, r);
    else
        rn = lmmp_odd_nCr_uint_(dst + shw, rn - shw, n, r);
    if (bits > 0) {
        dst[shw + rn] = lmmp_shl_(dst + shw, dst + shw, rn, bits);
        rn += shw + 1;
        rn -= dst[rn - 1] == 0;
    } else {
        rn += shw;
    }
    return rn;
}

引用了 LIMB_BITS, lmmp_debug_assert, lmmp_odd_nCr_uint_(), lmmp_odd_nCr_ushort_(), lmmp_shl_(), lmmp_zero , 以及 NCR_SHORT_LIMIT.

函数调用图:

lmmp_nCr_size_()

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

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

参数

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

返回

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 第 435 行定义.

                                      {
    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, 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的因子数

返回

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

被这些函数引用 pow_nPr_().

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

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_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_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_()

static mp_size_t lmmp_pow_1_size_ ( mp_limb_t  base, ulong  exp  )

inlinestatic

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

参数

base	底数
exp	指数

警告

exp>0, base>=1

返回

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

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

                                                              {
    lmmp_param_assert(base >= 1);
    lmmp_param_assert(exp > 0);
    return (ceil((double)(exp)*log2((double)base)) + LIMB_BITS - 1) / LIMB_BITS + 2; /* more two limbs */
}

引用了 LIMB_BITS , 以及 lmmp_param_assert.

被这些函数引用 lmmp_nPr_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	dst 的 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_()

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

inlinestatic

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

参数

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

警告

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

返回

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

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

                                                                         {
    lmmp_param_assert(n > 0);
    lmmp_param_assert(base[n - 1] != 0);
    mp_size_t rn = exp * (n - 1) * LIMB_BITS + ceil((double)exp * log2(base[n - 1]));
    return (rn + LIMB_BITS - 1) / LIMB_BITS + 2; /* more two limbs */
}

引用了 LIMB_BITS , 以及 lmmp_param_assert.

被这些函数引用 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 第 98 行定义.

                                                       {
    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.

被这些函数引用 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 第 121 行定义.

                                                           {
    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, 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 第 465 行定义.

                                      {
    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_(), 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_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

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

                                                                             {
    lmmp_param_assert(np != NULL && *nn > 0);
    lmmp_param_assert(dp != NULL && dn > 0);
 
    mp_srcptr pd_pow[MAX_EXP];
    mp_size_t pn_pow[MAX_EXP];
    // qp as quotient, rp as remainder, divp as divdend
    mp_ptr qp, rp, divp;
    mp_ptr prod;
    mp_size_t divn = *nn, qn;
    mp_size_t ret = 0;
    int i, j;
    pd_pow[0] = dp;
    pn_pow[0] = dn;
 
    TEMP_DECL;
 
    qp = TALLOC_TYPE(*nn, mp_limb_t);
    rp = TALLOC_TYPE(*nn, mp_limb_t);
 
    divp = np;
 
    for (i = 1; i < MAX_EXP; i++) {
        // if qn == 0, means cannot be divided by pd_pow[i - 1]
        if (qn = try_div_(qp, rp, divp, divn, pd_pow[i - 1], pn_pow[i - 1])) {
            divn = qn;
            LMMP_SWAP(qp, divp, mp_ptr);
 
            ret += 1 << (i - 1);
            pn_pow[i] = 2 * pn_pow[i - 1];
            if (divn < pn_pow[i]) {
                ++i;
                break;
            }
            prod = TALLOC_TYPE(pn_pow[i], mp_limb_t);
            lmmp_sqr_(prod, pd_pow[i - 1], pn_pow[i - 1]);
            pn_pow[i] -= (prod[pn_pow[i] - 1] == 0) ? 1 : 0;
            pd_pow[i] = prod;
        } else {
            break;
        }
    }
    for (j = i - 1; j > 0; --j) {
        if (qn = try_div_(qp, rp, divp, divn, pd_pow[j - 1], pn_pow[j - 1])) {
            divn = qn;
            LMMP_SWAP(qp, divp, mp_ptr);
 
            ret += 1 << (j - 1);
        }
    }
    if (qn == 0) {
        // 无法整除
        lmmp_copy(np, divp, divn);
        *nn = divn;
    } else {
        // 整除
        lmmp_copy(np, qp, qn);
        *nn = qn;
    }
 
    TEMP_FREE;
    return ret;
}

引用了 lmmp_copy, lmmp_param_assert, lmmp_sqr_(), LMMP_SWAP, MAX_EXP, TALLOC_TYPE, TEMP_DECL, TEMP_FREE , 以及 try_div_().

函数调用图:

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。

返回

结果指针，返回不超过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	底数
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	底数
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	底数
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	底数
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	底数
exp	指数

警告

0<base<=0xff, exp>0

返回

返回 dst 的实际 limb 长度

被这些函数引用 _odd_pow_().

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