LAMMP 4.2.0
Lamina High-Precision Arithmetic Library
载入中...
搜索中...
未找到
宏定义 | 函数
mul_toom44.c 文件参考
#include "../../../include/lammp/impl/mparam.h"
#include "../../../include/lammp/impl/toom_interp.h"
+ mul_toom44.c 的引用(Include)关系图:

浏览源代码.

宏定义

#define a0   numa
 
#define a1   (numa + n)
 
#define a2   (numa + 2 * n)
 
#define a3   (numa + 3 * n)
 
#define amx   (dst + n + 1) /* n+1 */
 
#define apx   dst /* n+1 */
 
#define b0   numb
 
#define b1   (numb + n)
 
#define b2   (numb + 2 * n)
 
#define b3   (numb + 3 * n)
 
#define bmx   (dst + 2 * n + 2) /* n+1 */
 
#define bpx   (dst + 4 * n + 2) /* n+1 */
 
#define lmmp_mul_n_(dst, numa, numb, n)
 Copyright (C) 2026 HJimmyK(Jericho Knox)
 
#define tp   (scratch + 8 * n + 5)
 
#define v0   dst /* 2n */
 
#define v1   (dst + 2 * n) /* 2n+1 */
 
#define v2   scratch /* 2n+1 */
 
#define vh   (scratch + 4 * n + 2) /* 2n+1 */
 
#define vinf   (dst + 6 * n) /* s+t */
 
#define vm1   (scratch + 6 * n + 3) /* 2n+1 */
 
#define vm2   (scratch + 2 * n + 1) /* 2n+1 */
 

函数

void lmmp_mul_toom44_ (mp_ptr restrict dst, mp_srcptr restrict numa, mp_size_t na, mp_srcptr restrict numb, mp_size_t nb)
 

宏定义说明

◆ a0

#define a0   numa

◆ a1

#define a1   (numa + n)

◆ a2

#define a2   (numa + 2 * n)

◆ a3

#define a3   (numa + 3 * n)

◆ amx

#define amx   (dst + n + 1) /* n+1 */

◆ apx

#define apx   dst /* n+1 */

◆ b0

#define b0   numb

◆ b1

#define b1   (numb + n)

◆ b2

#define b2   (numb + 2 * n)

◆ b3

#define b3   (numb + 3 * n)

◆ bmx

#define bmx   (dst + 2 * n + 2) /* n+1 */

◆ bpx

#define bpx   (dst + 4 * n + 2) /* n+1 */

◆ lmmp_mul_n_

#define lmmp_mul_n_ (   dst,
  numa,
  numb,
  n 
)
值:
if ((n) < MUL_TOOM22_THRESHOLD) \
lmmp_mul_basecase_((dst), (numa), (n), (numb), (n)); \
else if ((n) < MUL_TOOM33_THRESHOLD) \
lmmp_mul_toom22_((dst), (numa), (n), (numb), (n)); \
else if ((n) < MUL_TOOM44_THRESHOLD) \
lmmp_mul_toom33_((dst), (numa), (n), (numb), (n)); \
else \
lmmp_mul_toom44_((dst), (numa), (n), (numb), (n))
MUL_TOOM22_THRESHOLD
#define MUL_TOOM22_THRESHOLD
Definition mparam.h:46
MUL_TOOM33_THRESHOLD
#define MUL_TOOM33_THRESHOLD
Definition mparam.h:50
MUL_TOOM44_THRESHOLD
#define MUL_TOOM44_THRESHOLD
Definition mparam.h:52
numb
#define numb
n
#define n

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

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

49 {
50 lmmp_param_assert(na >= nb);
51 lmmp_param_assert(4 * na <= 5 * nb);
52 mp_size_t n, s, t;
53 mp_limb_t cy;
54 enum toom7_flags flags;
55
56#define a0 numa
57#define a1 (numa + n)
58#define a2 (numa + 2 * n)
59#define a3 (numa + 3 * n)
60#define b0 numb
61#define b1 (numb + n)
62#define b2 (numb + 2 * n)
63#define b3 (numb + 3 * n)
64
65 lmmp_debug_assert(na >= nb);
66
67 n = (na + 3) >> 2;
68 TEMP_S_DECL;
69 mp_ptr restrict scratch = SALLOC_TYPE(8 * n + 8, mp_limb_t);
70
71 s = na - 3 * n;
72 t = nb - 3 * n;
73
74 lmmp_debug_assert(0 < s && s <= n);
75 lmmp_debug_assert(0 < t && t <= n);
76 lmmp_debug_assert(s >= t);
77
78 /* NOTE: The multiplications to v2, vm2, vh and vm1 overwrites the
79 * following limb, so these must be computed in order, and we need a
80 * one limb gap to tp. */
81#define v0 dst /* 2n */
82#define v1 (dst + 2 * n) /* 2n+1 */
83#define vinf (dst + 6 * n) /* s+t */
84#define v2 scratch /* 2n+1 */
85#define vm2 (scratch + 2 * n + 1) /* 2n+1 */
86#define vh (scratch + 4 * n + 2) /* 2n+1 */
87#define vm1 (scratch + 6 * n + 3) /* 2n+1 */
88#define tp (scratch + 8 * n + 5)
89
90 /* apx and bpx must not overlap with v1 */
91#define apx dst /* n+1 */
92#define amx (dst + n + 1) /* n+1 */
93#define bmx (dst + 2 * n + 2) /* n+1 */
94#define bpx (dst + 4 * n + 2) /* n+1 */
95
96 /* Total scratch need: 8*n + 5 + scratch for recursive calls. This
97 gives roughly 32 n/3 + log term. */
98
99 /* Compute apx = a0 + 2 a1 + 4 a2 + 8 a3 and amx = a0 - 2 a1 + 4 a2 - 8 a3. */
100 flags = (enum toom7_flags)(toom7_w1_neg & lmmp_toom_eval_dgr3_pm2_(apx, amx, numa, n, s, tp));
101
102 /* Compute bpx = b0 + 2 b1 + 4 b2 + 8 b3 and bmx = b0 - 2 b1 + 4 b2 - 8 b3. */
103 flags = (enum toom7_flags)(flags ^ (toom7_w1_neg & lmmp_toom_eval_dgr3_pm2_(bpx, bmx, numb, n, t, tp)));
104
105 lmmp_mul_n_(v2, apx, bpx, n + 1); /* v2, 2n+1 limbs */
106 lmmp_mul_n_(vm2, amx, bmx, n + 1); /* vm2, 2n+1 limbs */
107
108 /* Compute apx = 8 a0 + 4 a1 + 2 a2 + a3 = (((2*a0 + a1) * 2 + a2) * 2 + a3 */
109
110 cy = lmmp_addshl1_n_(apx, a1, a0, n);
111 cy = 2 * cy + lmmp_addshl1_n_(apx, a2, apx, n);
112 if (s < n) {
113 mp_limb_t cy2;
114 cy2 = lmmp_addshl1_n_(apx, a3, apx, s);
115 apx[n] = 2 * cy + lmmp_shl_(apx + s, apx + s, n - s, 1);
116 lmmp_inc_1(apx + s, cy2);
117 } else
118 apx[n] = 2 * cy + lmmp_addshl1_n_(apx, a3, apx, n);
119
120
121 /* Compute bpx = 8 b0 + 4 b1 + 2 b2 + b3 = (((2*b0 + b1) * 2 + b2) * 2 + b3 */
122
123 cy = lmmp_addshl1_n_(bpx, b1, b0, n);
124 cy = 2 * cy + lmmp_addshl1_n_(bpx, b2, bpx, n);
125 if (t < n) {
126 mp_limb_t cy2;
127 cy2 = lmmp_addshl1_n_(bpx, b3, bpx, t);
128 bpx[n] = 2 * cy + lmmp_shl_(bpx + t, bpx + t, n - t, 1);
129 lmmp_inc_1(bpx + t, cy2);
130 } else
131 bpx[n] = 2 * cy + lmmp_addshl1_n_(bpx, b3, bpx, n);
132
133 lmmp_debug_assert(apx[n] < 15);
134 lmmp_debug_assert(bpx[n] < 15);
135
136 lmmp_mul_n_(vh, apx, bpx, n + 1); /* vh, 2n+1 limbs */
137
138 /* Compute apx = a0 + a1 + a2 + a3 and amx = a0 - a1 + a2 - a3. */
139 flags = (enum toom7_flags)(flags | (toom7_w3_neg & lmmp_toom_eval_dgr3_pm1_(apx, amx, numa, n, s, tp)));
140
141 /* Compute bpx = b0 + b1 + b2 + b3 and bmx = b0 - b1 + b2 - b3. */
142 flags = (enum toom7_flags)(flags ^ (toom7_w3_neg & lmmp_toom_eval_dgr3_pm1_(bpx, bmx, numb, n, t, tp)));
143
144 lmmp_mul_n_(vm1, amx, bmx, n + 1); /* vm1, 2n+1 limbs */
145 /* Clobbers amx, bmx. */
146 lmmp_mul_n_(v1, apx, bpx, n + 1); /* v1, 2n+1 limbs */
147
148 lmmp_mul_n_(v0, a0, b0, n);
149 if (s > t)
150 lmmp_mul_(vinf, a3, s, b3, t);
151 else
152 lmmp_mul_n_(vinf, a3, b3, s);
153
154 lmmp_toom_interp7_(dst, n, flags, vm2, vm1, v2, vh, s + t, tp);
155
156 TEMP_S_FREE;
157}
mp_ptr
mp_limb_t * mp_ptr
Definition lmmp.h:80
mp_size_t
uint64_t mp_size_t
Definition lmmp.h:77
lmmp_debug_assert
#define lmmp_debug_assert(x)
Definition lmmp.h:390
mp_limb_t
uint64_t mp_limb_t
Definition lmmp.h:76
lmmp_param_assert
#define lmmp_param_assert(x)
Definition lmmp.h:401
lmmp_mul_
void lmmp_mul_(mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_srcptr numb, mp_size_t nb)
不等长大数乘法操作 [dst,na+nb] = [numa,na] * [numb,nb]
lmmp_addshl1_n_
mp_limb_t lmmp_addshl1_n_(mp_ptr dst, mp_srcptr numa, mp_srcptr numb, mp_size_t n)
加法结合左移1位操作 [dst,n] = [numa,n] + ([numb,n] << 1)
Definition shl.c:66
lmmp_shl_
mp_limb_t lmmp_shl_(mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_size_t shl)
大数左移操作 [dst,na] = [numa,na]<<shl，dst的低shl位填充0
Definition shl.c:19
lmmp_inc_1
#define lmmp_inc_1(p, inc)
大数加指定值宏（预期无进位）
Definition lmmpn.h:950
t
#define t
s
#define s
b0
#define b0
v0
#define v0
a3
#define a3
lmmp_mul_n_
#define lmmp_mul_n_(dst, numa, numb, n)
Copyright (C) 2026 HJimmyK(Jericho Knox)
Definition mul_toom44.c:21
b1
#define b1
v2
#define v2
vm1
#define vm1
apx
#define apx
vh
#define vh
a2
#define a2
a0
#define a0
tp
#define tp
a1
#define a1
bmx
#define bmx
b3
#define b3
b2
#define b2
vinf
#define vinf
bpx
#define bpx
amx
#define amx
v1
#define v1
vm2
#define vm2
scratch
#define scratch
SALLOC_TYPE
#define SALLOC_TYPE(n, type)
Definition tmp_alloc.h:144
TEMP_S_DECL
#define TEMP_S_DECL
Definition tmp_alloc.h:133
TEMP_S_FREE
#define TEMP_S_FREE
Definition tmp_alloc.h:166
toom7_flags
toom7_flags
Definition toom_interp.h:24
toom7_w1_neg
@ toom7_w1_neg
Definition toom_interp.h:24
toom7_w3_neg
@ toom7_w3_neg
Definition toom_interp.h:24
lmmp_toom_interp7_
void lmmp_toom_interp7_(mp_ptr dst, mp_size_t n, enum toom7_flags flags, mp_ptr w1, mp_ptr w3, mp_ptr w4, mp_ptr w5, mp_size_t w6n, mp_ptr tp)
Toom插值计算（7点插值）：用于Toom-44、Toom-53、Toom-62 乘法算法
Definition mul_toom_interp7.c:55
lmmp_toom_eval_dgr3_pm2_
int lmmp_toom_eval_dgr3_pm2_(mp_ptr xp2, mp_ptr xm2, mp_srcptr xp, mp_size_t n, mp_size_t x3n, mp_ptr tp)
Toom-3 专用：3次多项式在 x = +2 和 x = -2 处求值 计算 P(+2) 和 P(-2)，其中 P(x) 是一个3次多项式（4段系数）。
Definition mul_toom_eval.c:39
lmmp_toom_eval_dgr3_pm1_
int lmmp_toom_eval_dgr3_pm1_(mp_ptr xp1, mp_ptr xm1, mp_srcptr xp, mp_size_t n, mp_size_t x3n, mp_ptr tp)
Toom-3 专用：3次多项式在 x = +1 和 x = -1 处求值 计算 P(+1) 和 P(-1)，其中 P(x) 是一个3次多项式（4段系数）。
Definition mul_toom_eval.c:19

◆ tp

#define tp   (scratch + 8 * n + 5)

◆ v0

#define v0   dst /* 2n */

◆ v1

#define v1   (dst + 2 * n) /* 2n+1 */

◆ v2

#define v2   scratch /* 2n+1 */

◆ vh

#define vh   (scratch + 4 * n + 2) /* 2n+1 */

◆ vinf

#define vinf   (dst + 6 * n) /* s+t */

◆ vm1

#define vm1   (scratch + 6 * n + 3) /* 2n+1 */

◆ vm2

#define vm2   (scratch + 2 * n + 1) /* 2n+1 */

函数说明

◆ lmmp_mul_toom44_()

void lmmp_mul_toom44_ ( mp_ptr restrict  dst,
mp_srcptr restrict  numa,
mp_size_t  na,
mp_srcptr restrict  numb,
mp_size_t  nb 
)

在文件 mul_toom44.c 第 49 行定义.

49 {
50 lmmp_param_assert(na >= nb);
51 lmmp_param_assert(4 * na <= 5 * nb);
52 mp_size_t n, s, t;
53 mp_limb_t cy;
54 enum toom7_flags flags;
55
56#define a0 numa
57#define a1 (numa + n)
58#define a2 (numa + 2 * n)
59#define a3 (numa + 3 * n)
60#define b0 numb
61#define b1 (numb + n)
62#define b2 (numb + 2 * n)
63#define b3 (numb + 3 * n)
64
65 lmmp_debug_assert(na >= nb);
66
67 n = (na + 3) >> 2;
68 TEMP_S_DECL;
69 mp_ptr restrict scratch = SALLOC_TYPE(8 * n + 8, mp_limb_t);
70
71 s = na - 3 * n;
72 t = nb - 3 * n;
73
74 lmmp_debug_assert(0 < s && s <= n);
75 lmmp_debug_assert(0 < t && t <= n);
76 lmmp_debug_assert(s >= t);
77
78 /* NOTE: The multiplications to v2, vm2, vh and vm1 overwrites the
79 * following limb, so these must be computed in order, and we need a
80 * one limb gap to tp. */
81#define v0 dst /* 2n */
82#define v1 (dst + 2 * n) /* 2n+1 */
83#define vinf (dst + 6 * n) /* s+t */
84#define v2 scratch /* 2n+1 */
85#define vm2 (scratch + 2 * n + 1) /* 2n+1 */
86#define vh (scratch + 4 * n + 2) /* 2n+1 */
87#define vm1 (scratch + 6 * n + 3) /* 2n+1 */
88#define tp (scratch + 8 * n + 5)
89
90 /* apx and bpx must not overlap with v1 */
91#define apx dst /* n+1 */
92#define amx (dst + n + 1) /* n+1 */
93#define bmx (dst + 2 * n + 2) /* n+1 */
94#define bpx (dst + 4 * n + 2) /* n+1 */
95
96 /* Total scratch need: 8*n + 5 + scratch for recursive calls. This
97 gives roughly 32 n/3 + log term. */
98
99 /* Compute apx = a0 + 2 a1 + 4 a2 + 8 a3 and amx = a0 - 2 a1 + 4 a2 - 8 a3. */
100 flags = (enum toom7_flags)(toom7_w1_neg & lmmp_toom_eval_dgr3_pm2_(apx, amx, numa, n, s, tp));
101
102 /* Compute bpx = b0 + 2 b1 + 4 b2 + 8 b3 and bmx = b0 - 2 b1 + 4 b2 - 8 b3. */
103 flags = (enum toom7_flags)(flags ^ (toom7_w1_neg & lmmp_toom_eval_dgr3_pm2_(bpx, bmx, numb, n, t, tp)));
104
105 lmmp_mul_n_(v2, apx, bpx, n + 1); /* v2, 2n+1 limbs */
106 lmmp_mul_n_(vm2, amx, bmx, n + 1); /* vm2, 2n+1 limbs */
107
108 /* Compute apx = 8 a0 + 4 a1 + 2 a2 + a3 = (((2*a0 + a1) * 2 + a2) * 2 + a3 */
109
110 cy = lmmp_addshl1_n_(apx, a1, a0, n);
111 cy = 2 * cy + lmmp_addshl1_n_(apx, a2, apx, n);
112 if (s < n) {
113 mp_limb_t cy2;
114 cy2 = lmmp_addshl1_n_(apx, a3, apx, s);
115 apx[n] = 2 * cy + lmmp_shl_(apx + s, apx + s, n - s, 1);
116 lmmp_inc_1(apx + s, cy2);
117 } else
118 apx[n] = 2 * cy + lmmp_addshl1_n_(apx, a3, apx, n);
119
120
121 /* Compute bpx = 8 b0 + 4 b1 + 2 b2 + b3 = (((2*b0 + b1) * 2 + b2) * 2 + b3 */
122
123 cy = lmmp_addshl1_n_(bpx, b1, b0, n);
124 cy = 2 * cy + lmmp_addshl1_n_(bpx, b2, bpx, n);
125 if (t < n) {
126 mp_limb_t cy2;
127 cy2 = lmmp_addshl1_n_(bpx, b3, bpx, t);
128 bpx[n] = 2 * cy + lmmp_shl_(bpx + t, bpx + t, n - t, 1);
129 lmmp_inc_1(bpx + t, cy2);
130 } else
131 bpx[n] = 2 * cy + lmmp_addshl1_n_(bpx, b3, bpx, n);
132
133 lmmp_debug_assert(apx[n] < 15);
134 lmmp_debug_assert(bpx[n] < 15);
135
136 lmmp_mul_n_(vh, apx, bpx, n + 1); /* vh, 2n+1 limbs */
137
138 /* Compute apx = a0 + a1 + a2 + a3 and amx = a0 - a1 + a2 - a3. */
139 flags = (enum toom7_flags)(flags | (toom7_w3_neg & lmmp_toom_eval_dgr3_pm1_(apx, amx, numa, n, s, tp)));
140
141 /* Compute bpx = b0 + b1 + b2 + b3 and bmx = b0 - b1 + b2 - b3. */
142 flags = (enum toom7_flags)(flags ^ (toom7_w3_neg & lmmp_toom_eval_dgr3_pm1_(bpx, bmx, numb, n, t, tp)));
143
144 lmmp_mul_n_(vm1, amx, bmx, n + 1); /* vm1, 2n+1 limbs */
145 /* Clobbers amx, bmx. */
146 lmmp_mul_n_(v1, apx, bpx, n + 1); /* v1, 2n+1 limbs */
147
148 lmmp_mul_n_(v0, a0, b0, n);
149 if (s > t)
150 lmmp_mul_(vinf, a3, s, b3, t);
151 else
152 lmmp_mul_n_(vinf, a3, b3, s);
153
154 lmmp_toom_interp7_(dst, n, flags, vm2, vm1, v2, vh, s + t, tp);
155
156 TEMP_S_FREE;
157}

引用了 a0, a1, a2, a3, amx, apx, b0, b1, b2, b3, bmx, bpx, lmmp_addshl1_n_(), lmmp_debug_assert, lmmp_inc_1, lmmp_mul_(), lmmp_mul_n_, lmmp_param_assert, lmmp_shl_(), lmmp_toom_eval_dgr3_pm1_(), lmmp_toom_eval_dgr3_pm2_(), lmmp_toom_interp7_(), n, numb, s, SALLOC_TYPE, scratch, t, TEMP_S_DECL, TEMP_S_FREE, toom7_w1_neg, toom7_w3_neg, tp, v0, v1, v2, vh, vinf, vm1 , 以及 vm2.

+ 函数调用图: