1 ```#ifndef NPY_EXTINT128_H_ ``` 2 ```#define NPY_EXTINT128_H_ ``` 3 4 5 ```typedef struct { ``` 6 ``` signed char sign; ``` 7 ``` npy_uint64 lo, hi; ``` 8 ```} npy_extint128_t; ``` 9 10 11 ```/* Integer addition with overflow checking */ ``` 12 ```static NPY_INLINE npy_int64 ``` 13 ```safe_add(npy_int64 a, npy_int64 b, char *overflow_flag) ``` 14 ```{ ``` 15 1 ``` if (a > 0 && b > NPY_MAX_INT64 - a) { ``` 16 1 ``` *overflow_flag = 1; ``` 17 ``` } ``` 18 1 ``` else if (a < 0 && b < NPY_MIN_INT64 - a) { ``` 19 1 ``` *overflow_flag = 1; ``` 20 ``` } ``` 21 1 ``` return a + b; ``` 22 ```} ``` 23 24 25 ```/* Integer subtraction with overflow checking */ ``` 26 ```static NPY_INLINE npy_int64 ``` 27 ```safe_sub(npy_int64 a, npy_int64 b, char *overflow_flag) ``` 28 ```{ ``` 29 1 ``` if (a >= 0 && b < a - NPY_MAX_INT64) { ``` 30 1 ``` *overflow_flag = 1; ``` 31 ``` } ``` 32 1 ``` else if (a < 0 && b > a - NPY_MIN_INT64) { ``` 33 1 ``` *overflow_flag = 1; ``` 34 ``` } ``` 35 1 ``` return a - b; ``` 36 ```} ``` 37 38 39 ```/* Integer multiplication with overflow checking */ ``` 40 ```static NPY_INLINE npy_int64 ``` 41 1 ```safe_mul(npy_int64 a, npy_int64 b, char *overflow_flag) ``` 42 ```{ ``` 43 1 ``` if (a > 0) { ``` 44 1 ``` if (b > NPY_MAX_INT64 / a || b < NPY_MIN_INT64 / a) { ``` 45 1 ``` *overflow_flag = 1; ``` 46 ``` } ``` 47 ``` } ``` 48 1 ``` else if (a < 0) { ``` 49 1 ``` if (b > 0 && a < NPY_MIN_INT64 / b) { ``` 50 1 ``` *overflow_flag = 1; ``` 51 ``` } ``` 52 1 ``` else if (b < 0 && a < NPY_MAX_INT64 / b) { ``` 53 1 ``` *overflow_flag = 1; ``` 54 ``` } ``` 55 ``` } ``` 56 1 ``` return a * b; ``` 57 ```} ``` 58 59 60 ```/* Long integer init */ ``` 61 ```static NPY_INLINE npy_extint128_t ``` 62 ```to_128(npy_int64 x) ``` 63 ```{ ``` 64 ``` npy_extint128_t result; ``` 65 1 ``` result.sign = (x >= 0 ? 1 : -1); ``` 66 1 ``` if (x >= 0) { ``` 67 1 ``` result.lo = x; ``` 68 ``` } ``` 69 ``` else { ``` 70 1 ``` result.lo = (npy_uint64)(-(x + 1)) + 1; ``` 71 ``` } ``` 72 1 ``` result.hi = 0; ``` 73 1 ``` return result; ``` 74 ```} ``` 75 76 77 ```static NPY_INLINE npy_int64 ``` 78 ```to_64(npy_extint128_t x, char *overflow) ``` 79 ```{ ``` 80 1 ``` if (x.hi != 0 || ``` 81 1 ``` (x.sign > 0 && x.lo > NPY_MAX_INT64) || ``` 82 1 ``` (x.sign < 0 && x.lo != 0 && x.lo - 1 > -(NPY_MIN_INT64 + 1))) { ``` 83 1 ``` *overflow = 1; ``` 84 ``` } ``` 85 1 ``` return x.lo * x.sign; ``` 86 ```} ``` 87 88 89 ```/* Long integer multiply */ ``` 90 ```static NPY_INLINE npy_extint128_t ``` 91 1 ```mul_64_64(npy_int64 a, npy_int64 b) ``` 92 ```{ ``` 93 ``` npy_extint128_t x, y, z; ``` 94 ``` npy_uint64 x1, x2, y1, y2, r1, r2, prev; ``` 95 96 1 ``` x = to_128(a); ``` 97 1 ``` y = to_128(b); ``` 98 99 1 ``` x1 = x.lo & 0xffffffff; ``` 100 1 ``` x2 = x.lo >> 32; ``` 101 102 1 ``` y1 = y.lo & 0xffffffff; ``` 103 1 ``` y2 = y.lo >> 32; ``` 104 105 1 ``` r1 = x1*y2; ``` 106 1 ``` r2 = x2*y1; ``` 107 108 1 ``` z.sign = x.sign * y.sign; ``` 109 1 ``` z.hi = x2*y2 + (r1 >> 32) + (r2 >> 32); ``` 110 1 ``` z.lo = x1*y1; ``` 111 112 ``` /* Add with carry */ ``` 113 1 ``` prev = z.lo; ``` 114 1 ``` z.lo += (r1 << 32); ``` 115 1 ``` if (z.lo < prev) { ``` 116 1 ``` ++z.hi; ``` 117 ``` } ``` 118 119 1 ``` prev = z.lo; ``` 120 1 ``` z.lo += (r2 << 32); ``` 121 1 ``` if (z.lo < prev) { ``` 122 1 ``` ++z.hi; ``` 123 ``` } ``` 124 125 1 ``` return z; ``` 126 ```} ``` 127 128 129 ```/* Long integer add */ ``` 130 ```static NPY_INLINE npy_extint128_t ``` 131 1 ```add_128(npy_extint128_t x, npy_extint128_t y, char *overflow) ``` 132 ```{ ``` 133 ``` npy_extint128_t z; ``` 134 135 1 ``` if (x.sign == y.sign) { ``` 136 1 ``` z.sign = x.sign; ``` 137 1 ``` z.hi = x.hi + y.hi; ``` 138 1 ``` if (z.hi < x.hi) { ``` 139 1 ``` *overflow = 1; ``` 140 ``` } ``` 141 1 ``` z.lo = x.lo + y.lo; ``` 142 1 ``` if (z.lo < x.lo) { ``` 143 1 ``` if (z.hi == NPY_MAX_UINT64) { ``` 144 1 ``` *overflow = 1; ``` 145 ``` } ``` 146 1 ``` ++z.hi; ``` 147 ``` } ``` 148 ``` } ``` 149 1 ``` else if (x.hi > y.hi || (x.hi == y.hi && x.lo >= y.lo)) { ``` 150 1 ``` z.sign = x.sign; ``` 151 1 ``` z.hi = x.hi - y.hi; ``` 152 1 ``` z.lo = x.lo; ``` 153 1 ``` z.lo -= y.lo; ``` 154 1 ``` if (z.lo > x.lo) { ``` 155 1 ``` --z.hi; ``` 156 ``` } ``` 157 ``` } ``` 158 ``` else { ``` 159 1 ``` z.sign = y.sign; ``` 160 1 ``` z.hi = y.hi - x.hi; ``` 161 1 ``` z.lo = y.lo; ``` 162 1 ``` z.lo -= x.lo; ``` 163 1 ``` if (z.lo > y.lo) { ``` 164 1 ``` --z.hi; ``` 165 ``` } ``` 166 ``` } ``` 167 168 1 ``` return z; ``` 169 ```} ``` 170 171 172 ```/* Long integer negation */ ``` 173 ```static NPY_INLINE npy_extint128_t ``` 174 ```neg_128(npy_extint128_t x) ``` 175 ```{ ``` 176 1 ``` npy_extint128_t z = x; ``` 177 1 ``` z.sign *= -1; ``` 178 1 ``` return z; ``` 179 ```} ``` 180 181 182 ```static NPY_INLINE npy_extint128_t ``` 183 ```sub_128(npy_extint128_t x, npy_extint128_t y, char *overflow) ``` 184 ```{ ``` 185 1 ``` return add_128(x, neg_128(y), overflow); ``` 186 ```} ``` 187 188 189 ```static NPY_INLINE npy_extint128_t ``` 190 ```shl_128(npy_extint128_t v) ``` 191 ```{ ``` 192 ``` npy_extint128_t z; ``` 193 1 ``` z = v; ``` 194 1 ``` z.hi <<= 1; ``` 195 1 ``` z.hi |= (z.lo & (((npy_uint64)1) << 63)) >> 63; ``` 196 1 ``` z.lo <<= 1; ``` 197 1 ``` return z; ``` 198 ```} ``` 199 200 201 ```static NPY_INLINE npy_extint128_t ``` 202 ```shr_128(npy_extint128_t v) ``` 203 ```{ ``` 204 ``` npy_extint128_t z; ``` 205 1 ``` z = v; ``` 206 1 ``` z.lo >>= 1; ``` 207 1 ``` z.lo |= (z.hi & 0x1) << 63; ``` 208 1 ``` z.hi >>= 1; ``` 209 1 ``` return z; ``` 210 ```} ``` 211 212 ```static NPY_INLINE int ``` 213 1 ```gt_128(npy_extint128_t a, npy_extint128_t b) ``` 214 ```{ ``` 215 1 ``` if (a.sign > 0 && b.sign > 0) { ``` 216 1 ``` return (a.hi > b.hi) || (a.hi == b.hi && a.lo > b.lo); ``` 217 ``` } ``` 218 1 ``` else if (a.sign < 0 && b.sign < 0) { ``` 219 1 ``` return (a.hi < b.hi) || (a.hi == b.hi && a.lo < b.lo); ``` 220 ``` } ``` 221 1 ``` else if (a.sign > 0 && b.sign < 0) { ``` 222 1 ``` return a.hi != 0 || a.lo != 0 || b.hi != 0 || b.lo != 0; ``` 223 ``` } ``` 224 ``` else { ``` 225 ``` return 0; ``` 226 ``` } ``` 227 ```} ``` 228 229 230 ```/* Long integer divide */ ``` 231 ```static NPY_INLINE npy_extint128_t ``` 232 1 ```divmod_128_64(npy_extint128_t x, npy_int64 b, npy_int64 *mod) ``` 233 ```{ ``` 234 ``` npy_extint128_t remainder, pointer, result, divisor; ``` 235 1 ``` char overflow = 0; ``` 236 237 ``` assert(b > 0); ``` 238 239 1 ``` if (b <= 1 || x.hi == 0) { ``` 240 1 ``` result.sign = x.sign; ``` 241 1 ``` result.lo = x.lo / b; ``` 242 1 ``` result.hi = x.hi / b; ``` 243 1 ``` *mod = x.sign * (x.lo % b); ``` 244 1 ``` return result; ``` 245 ``` } ``` 246 247 ``` /* Long division, not the most efficient choice */ ``` 248 1 ``` remainder = x; ``` 249 1 ``` remainder.sign = 1; ``` 250 251 1 ``` divisor.sign = 1; ``` 252 1 ``` divisor.hi = 0; ``` 253 1 ``` divisor.lo = b; ``` 254 255 1 ``` result.sign = 1; ``` 256 1 ``` result.lo = 0; ``` 257 1 ``` result.hi = 0; ``` 258 259 1 ``` pointer.sign = 1; ``` 260 1 ``` pointer.lo = 1; ``` 261 1 ``` pointer.hi = 0; ``` 262 263 1 ``` while ((divisor.hi & (((npy_uint64)1) << 63)) == 0 && ``` 264 1 ``` gt_128(remainder, divisor)) { ``` 265 1 ``` divisor = shl_128(divisor); ``` 266 1 ``` pointer = shl_128(pointer); ``` 267 ``` } ``` 268 269 1 ``` while (pointer.lo || pointer.hi) { ``` 270 1 ``` if (!gt_128(divisor, remainder)) { ``` 271 1 ``` remainder = sub_128(remainder, divisor, &overflow); ``` 272 1 ``` result = add_128(result, pointer, &overflow); ``` 273 ``` } ``` 274 1 ``` divisor = shr_128(divisor); ``` 275 1 ``` pointer = shr_128(pointer); ``` 276 ``` } ``` 277 278 ``` /* Fix signs and return; cannot overflow */ ``` 279 1 ``` result.sign = x.sign; ``` 280 1 ``` *mod = x.sign * remainder.lo; ``` 281 282 1 ``` return result; ``` 283 ```} ``` 284 285 286 ```/* Divide and round down (positive divisor; no overflows) */ ``` 287 ```static NPY_INLINE npy_extint128_t ``` 288 1 ```floordiv_128_64(npy_extint128_t a, npy_int64 b) ``` 289 ```{ ``` 290 ``` npy_extint128_t result; ``` 291 ``` npy_int64 remainder; ``` 292 1 ``` char overflow = 0; ``` 293 ``` assert(b > 0); ``` 294 1 ``` result = divmod_128_64(a, b, &remainder); ``` 295 1 ``` if (a.sign < 0 && remainder != 0) { ``` 296 1 ``` result = sub_128(result, to_128(1), &overflow); ``` 297 ``` } ``` 298 1 ``` return result; ``` 299 ```} ``` 300 301 302 ```/* Divide and round up (positive divisor; no overflows) */ ``` 303 ```static NPY_INLINE npy_extint128_t ``` 304 1 ```ceildiv_128_64(npy_extint128_t a, npy_int64 b) ``` 305 ```{ ``` 306 ``` npy_extint128_t result; ``` 307 ``` npy_int64 remainder; ``` 308 1 ``` char overflow = 0; ``` 309 ``` assert(b > 0); ``` 310 1 ``` result = divmod_128_64(a, b, &remainder); ``` 311 1 ``` if (a.sign > 0 && remainder != 0) { ``` 312 1 ``` result = add_128(result, to_128(1), &overflow); ``` 313 ``` } ``` 314 1 ``` return result; ``` 315 ```} ``` 316 317 ```#endif ```

Read our documentation on viewing source code .