LAMMP 4.2.0
Lamina High-Precision Arithmetic Library
载入中...
搜索中...
未找到
inlines.h
浏览该文件的文档.
1/**
2 * Copyright (C) 2026 HJimmyK(Jericho Knox)
3 *
4 * This file is part of LAMMP.
5 *
6 * LAMMP is free software: you can redistribute it and/or modify it under
7 * the terms of the GNU Lesser General Public License (LGPL) as published
8 * by the Free Software Foundation; either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed WITHOUT ANY WARRANTY.
12 *
13 * See <https://www.gnu.org/licenses/>.
14 */
15
16#ifndef __LAMMP_INLINES_H__
17#define __LAMMP_INLINES_H__
18
19#include "../lmmpn.h"
20#include "mparam.h"
21
22#ifdef _MSC_VER
23#include <intrin.h>
24#endif
25
26static inline int __lmmp_limb_bits_(mp_limb_t x) {
27 int k = 0;
28 while (x) {
29 x >>= 1;
30 k++;
31 }
32 return k;
33}
34
35static inline int __lmmp_limb_popcnt_(mp_limb_t x) {
36#if defined(__GNUC__) || defined(__clang__)
37 mp_limb_t count;
38#if defined(__x86_64__) && defined(USE_ASM)
39 __asm__ volatile("popcnt %1, %0" : "=r"(count) : "r"(x) : "cc");
40#else
41 count = __builtin_popcountll(x);
42#endif
43 return count;
44#elif defined(_MSC_VER) && (defined(_M_X64) || defined(_M_ARM64))
45 return (int)__popcnt64(x);
46#else
47 int k = 0;
48 while (x) {
49 k += x & 1;
50 x >>= 1;
51 }
52 return k;
53#endif
54}
55
56static inline int __lmmp_leading_zeros_(mp_limb_t x) {
57 if (x == 0) return 64;
58#ifdef __GNUC__
59 return __builtin_clzll(x);
60#elif defined(_MSC_VER) && (defined(_M_X64) || defined(_M_ARM64))
61 unsigned long index;
62 _BitScanReverse64(&index, x);
63 return 63 - (int)index;
64#else
65 int n = 0;
66 if (x <= 0x00000000FFFFFFFF) { n += 32; x <<= 32; }
67 if (x <= 0x0000FFFFFFFFFFFF) { n += 16; x <<= 16; }
68 if (x <= 0x00FFFFFFFFFFFFFF) { n += 8; x <<= 8; }
69 if (x <= 0x0FFFFFFFFFFFFFFF) { n += 4; x <<= 4; }
70 if (x <= 0x3FFFFFFFFFFFFFFF) { n += 2; x <<= 2; }
71 if (x <= 0x7FFFFFFFFFFFFFFF) { n += 1; x <<= 1; }
72 return n;
73#endif
74}
75
76static inline int __lmmp_tailing_zeros_(mp_limb_t x) {
77 if (x == 0) return 64;
78#ifdef __GNUC__
79 return __builtin_ctzll(x);
80#elif defined(_MSC_VER) && (defined(_M_X64) || defined(_M_ARM64))
81 unsigned long index;
82 _BitScanForward64(&index, x);
83 return (int)index;
84#else
85 int n = 0;
86 if ((x & 0x00000000FFFFFFFF) == 0) { n += 32; x >>= 32; }
87 if ((x & 0x000000000000FFFF) == 0) { n += 16; x >>= 16; }
88 if ((x & 0x00000000000000FF) == 0) { n += 8; x >>= 8; }
89 if ((x & 0x000000000000000F) == 0) { n += 4; x >>= 4; }
90 if ((x & 0x0000000000000003) == 0) { n += 2; x >>= 2; }
91 if ((x & 0x0000000000000001) == 0) { n += 1; x >>= 1; }
92 return n;
93#endif
94}
95
96static inline void __lmmp_mullh_(mp_limb_t a, mp_limb_t b, mp_ptr restrict dst) {
97#if (defined(__GNUC__) || defined(__clang__)) && defined(__SIZEOF_INT128__)
98 __uint128_t prod = (__uint128_t)a * b;
99 dst[0] = (mp_limb_t)prod;
100 dst[1] = (mp_limb_t)(prod >> 64);
101#elif defined(_MSC_VER) && (defined(_M_X64) || defined(_M_ARM64))
102 dst[0] = _umul128(a, b, dst + 1);
103#else
104 uint64_t ah = a >> 32, bh = b >> 32;
105 a = (uint32_t)a, b = (uint32_t)b;
106 uint64_t r0 = a * b, r1 = a * bh, r2 = ah * b, r3 = ah * bh;
107 r3 += (r1 >> 32) + (r2 >> 32);
108 r1 = (uint32_t)r1, r2 = (uint32_t)r2;
109 r1 += r2;
110 r1 += (r0 >> 32);
111 dst[1] = r3 + (r1 >> 32);
112 dst[0] = (r1 << 32) | (uint32_t)r0;
113#endif
114}
115
117#if (defined(__GNUC__) || defined(__clang__)) && defined(__SIZEOF_INT128__)
118 __uint128_t t = (__uint128_t)a * (__uint128_t)b;
119 return (mp_limb_t)(t >> 64);
120#elif defined(_MSC_VER) && (defined(_M_X64) || defined(_M_ARM64))
121 return __umulh(a, b);
122#else
123 uint64_t ah = a >> 32, bh = b >> 32;
124 a = (uint32_t)a, b = (uint32_t)b;
125 uint64_t r0 = a * b, r1 = a * bh, r2 = ah * b, r3 = ah * bh;
126 r3 += (r1 >> 32) + (r2 >> 32);
127 r1 = (uint32_t)r1, r2 = (uint32_t)r2;
128 r1 += r2;
129 r1 += (r0 >> 32);
130 return r3 + (r1 >> 32);
131#endif
132}
133
134static inline void __lmmp_mul_n_(mp_ptr restrict dst, mp_srcptr restrict numa, mp_srcptr restrict numb, mp_size_t n) {
136 lmmp_mul_basecase_(dst, numa, n, numb, n);
137 else if (n < MUL_TOOM33_THRESHOLD)
138 lmmp_mul_toom22_(dst, numa, n, numb, n);
139 else if (n < MUL_TOOM44_THRESHOLD)
140 lmmp_mul_toom33_(dst, numa, n, numb, n);
141 else if (n < MUL_FFT_THRESHOLD)
142 lmmp_mul_toom44_(dst, numa, n, numb, n);
143 else
144 lmmp_mul_fft_(dst, numa, n, numb, n);
145}
146
147static inline void __lmmp_sqr_(mp_ptr restrict dst, mp_srcptr restrict numa, mp_size_t na) {
148 if (na < MUL_TOOM22_THRESHOLD)
149 lmmp_sqr_basecase_(dst, numa, na);
150 else if (na < MUL_TOOM33_THRESHOLD)
151 lmmp_sqr_toom2_(dst, numa, na);
152 else if (na < MUL_TOOM44_THRESHOLD)
153 lmmp_sqr_toom3_(dst, numa, na);
154 else if (na < MUL_FFT_THRESHOLD)
155 lmmp_sqr_toom4_(dst, numa, na);
156 else
157 lmmp_mul_fft_(dst, numa, na, numa, na);
158}
159
160#define lmmp_leading_zeros_ __lmmp_leading_zeros_
161#define lmmp_tailing_zeros_ __lmmp_tailing_zeros_
162#define lmmp_limb_bits_ __lmmp_limb_bits_
163#define lmmp_limb_popcnt_ __lmmp_limb_popcnt_
164#define lmmp_mullh_ __lmmp_mullh_
165#define lmmp_mulh_ __lmmp_mulh_
166#define lmmp_sqr_ __lmmp_sqr_
167#define lmmp_mul_n_ __lmmp_mul_n_
168
169#endif // __LAMMP_INLINES_H__
#define k
static int __lmmp_limb_bits_(mp_limb_t x)
Copyright (C) 2026 HJimmyK(Jericho Knox)
Definition inlines.h:26
static void __lmmp_mul_n_(mp_ptr restrict dst, mp_srcptr restrict numa, mp_srcptr restrict numb, mp_size_t n)
Definition inlines.h:134
static mp_limb_t __lmmp_mulh_(mp_limb_t a, mp_limb_t b)
Definition inlines.h:116
static int __lmmp_tailing_zeros_(mp_limb_t x)
Definition inlines.h:76
static int __lmmp_limb_popcnt_(mp_limb_t x)
Definition inlines.h:35
static void __lmmp_mullh_(mp_limb_t a, mp_limb_t b, mp_ptr restrict dst)
Definition inlines.h:96
static void __lmmp_sqr_(mp_ptr restrict dst, mp_srcptr restrict numa, mp_size_t na)
Definition inlines.h:147
static int __lmmp_leading_zeros_(mp_limb_t x)
Definition inlines.h:56
mp_limb_t * mp_ptr
Definition lmmp.h:80
uint64_t mp_size_t
Definition lmmp.h:77
const mp_limb_t * mp_srcptr
Definition lmmp.h:81
uint64_t mp_limb_t
Definition lmmp.h:76
void lmmp_mul_toom22_(mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_srcptr numb, mp_size_t nb)
Toom-22乘法运算 [dst,na+nb] = [numa,na] * [numb,nb]
void lmmp_mul_toom44_(mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_srcptr numb, mp_size_t nb)
Toom-44乘法运算 [dst,na+nb] = [numa,na] * [numb,nb]
void lmmp_sqr_basecase_(mp_ptr dst, mp_srcptr numa, mp_size_t na)
基础平方运算 [dst,2*na] = [numa,na]^2
void lmmp_mul_basecase_(mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_srcptr numb, mp_size_t nb)
基础乘法运算 [dst,na+nb] = [numa,na] * [numb,nb]
void lmmp_sqr_toom2_(mp_ptr dst, mp_srcptr numa, mp_size_t na)
Toom-2平方运算 [dst,2*na] = [numa,na]^2
void lmmp_sqr_toom3_(mp_ptr dst, mp_srcptr numa, mp_size_t na)
Toom-3平方运算 [dst,2*na] = [numa,na]^2
void lmmp_sqr_toom4_(mp_ptr pp, mp_srcptr ap, mp_size_t an)
Toom-4平方运算 [dst,2*na] = [numa,na]^2
void lmmp_mul_fft_(mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_srcptr numb, mp_size_t nb)
FFT乘法运算 [dst,na+nb] = [numa,na] * [numb,nb]
Definition mul_fft.c:1095
void lmmp_mul_toom33_(mp_ptr dst, mp_srcptr numa, mp_size_t na, mp_srcptr numb, mp_size_t nb)
Toom-33乘法运算 [dst,na+nb] = [numa,na] * [numb,nb]
#define MUL_FFT_THRESHOLD
Definition mparam.h:54
#define MUL_TOOM22_THRESHOLD
Definition mparam.h:46
#define MUL_TOOM33_THRESHOLD
Definition mparam.h:50
#define MUL_TOOM44_THRESHOLD
Definition mparam.h:52
#define r2
#define r1
#define r3
#define r0
#define t
#define numb
#define n