1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
|
/* Internal 128 bit int support.
Copyright (C) 2024-2025 Free Software Foundation, Inc.
This file is part of the GNU C Library.
The GNU C Library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
The GNU C Library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, see
<https://www.gnu.org/licenses/>. */
#ifndef _MATH_INT128_H
#define _MATH_INT128_H
/* Limited support for internal 128 bit integer, used on some math
implementations. It uses compiler builtin type if supported, otherwise
it is emulated. Only unsigned and some operations are currently supported:
- u128_t: the 128 bit unsigned type.
- u128_high: return the high part of the number.
- u128_low: return the low part of the number.
- u128_from_u64: create a 128 bit number from a 64 bit one.
- u128_mul: multiply two 128 bit numbers.
- u128_add: add two 128 bit numbers.
- u128_lshift: left shift a number.
- u128_rshift: right shift a number.
*/
#if defined __BITINT_MAXWIDTH__ && __BITINT_MAXWIDTH__ >= 128
typedef unsigned _BitInt(128) u128;
# define __MATH_INT128_BUILTIN_TYPE 1
#elif defined __SIZEOF_INT128__
typedef unsigned __int128 u128;
# define __MATH_INT128_BUILTIN_TYPE 1
#else
# define __MATH_INT128_BUILTIN_TYPE 0
#endif
#if __MATH_INT128_BUILTIN_TYPE
# define u128_high(__x) (uint64_t)((__x) >> 64)
# define u128_low(__x) (uint64_t)(__x)
# define u128_from_u64(__x) (u128)(__x)
# define u128_mul(__x, __y) (__x) * (__y)
# define u128_add(__x, __y) (__x) + (__y)
# define u128_lshift(__x, __y) (__x) << (__y)
# define u128_rshift(__x, __y) (__x) >> (__y)
#else
typedef struct
{
uint64_t low;
uint64_t high;
} u128;
# define u128_high(__x) (__x).high
# define u128_low(__x) (__x).low
# define u128_from_u64(__x) (u128){.low = (__x), .high = 0}
# define MASK32 (UINT64_C(0xffffffff))
static u128 u128_add (u128 x, u128 y)
{
bool carry = x.low + y.low < x.low;
return (u128) { .high = x.high + y.high + carry, .low = x.low + y.low };
}
static u128 u128_lshift (u128 x, unsigned int n)
{
switch (n)
{
case 0: return x;
case 1 ... 63: return (u128) { .high = (x.high << n) | (x.low >> (64 - n)),
.low = x.low << n };
case 64 ...127: return (u128) { .high = x.low << (n - 64), .low = 0};
default: return (u128) { .high = 0, .low = 0 };
}
}
static u128 u128_rshift (u128 x, unsigned int n)
{
switch (n)
{
case 0: return x;
case 1 ... 63: return (u128) { .high = x.high >> n,
.low = (x.high << (64 - n)) | (x.low >> n) };
case 64 ...127: return (u128) { .high = 0, .low = x.high >> (n - 64) };
default: return (u128) { .high = 0, .low = 0 };
}
}
static u128 u128_mul (u128 x, u128 y)
{
if (x.high == 0 && y.high == 0)
{
uint64_t x0 = x.low & MASK32;
uint64_t x1 = x.low >> 32;
uint64_t y0 = y.low & MASK32;
uint64_t y1 = y.low >> 32;
u128 x0y0 = { .high = 0, .low = x0 * y0 };
u128 x0y1 = { .high = 0, .low = x0 * y1 };
u128 x1y0 = { .high = 0, .low = x1 * y0 };
u128 x1y1 = { .high = x1 * y1, .low = 0 };
/* x0y0 + ((x0y1 + x1y0) << 32) + x1y1 */
return u128_add (u128_add (x0y0,
u128_lshift (u128_add (x0y1, x1y0),
32)),
x1y1);
}
else
{
uint64_t x0 = x.low & MASK32;
uint64_t x1 = x.low >> 32;
uint64_t x2 = x.high & MASK32;
uint64_t x3 = x.high >> 32;
uint64_t y0 = y.low & MASK32;
uint64_t y1 = y.low >> 32;
uint64_t y2 = y.high & MASK32;
uint64_t y3 = y.high >> 32;
u128 x0y0 = { .high = 0, .low = x0 * y0 };
u128 x0y1 = { .high = 0, .low = x0 * y1 };
u128 x0y2 = { .high = 0, .low = x0 * y2 };
u128 x0y3 = { .high = 0, .low = x0 * y3 };
u128 x1y0 = { .high = 0, .low = x1 * y0 };
u128 x1y1 = { .high = 0, .low = x1 * y1 };
u128 x1y2 = { .high = 0, .low = x1 * y2 };
u128 x2y0 = { .high = 0, .low = x2 * y0 };
u128 x2y1 = { .high = 0, .low = x2 * y1 };
u128 x3y0 = { .high = 0, .low = x3 * y0 };
/* x0y0 + ((x0y1 + x1y0) << 32) + ((x0y2 + x1y1 + x2y0) << 64) +
((x0y3 + x1y2 + x2y1 + x3y0) << 96) */
u128 r0 = u128_add (x0y0,
u128_lshift (u128_add (x0y1, x1y0),
32));
u128 r1 = u128_add (u128_lshift (u128_add (u128_add (x0y2, x1y1), x2y0),
64),
u128_lshift (u128_add (u128_add (x0y3, x1y2),
u128_add (x2y1, x3y0)),
96));
return u128_add (r0, r1);
}
}
#endif /* __SIZEOF_INT128__ */
#endif
|
