Qrack  9.13
General classical-emulating-quantum development framework
big_integer.hpp
Go to the documentation of this file.
1 //
3 // (C) Daniel Strano and the Qimcifa contributors, 2022, 2023. All rights reserved.
4 //
5 // This header has been adapted for OpenCL and C, from big_integer.c by Andre Azevedo.
6 //
7 // Original file:
8 //
9 // big_integer.c
10 // Description: "Arbitrary"-precision integer
11 // Author: Andre Azevedo <http://github.com/andreazevedo>
12 //
13 // The MIT License (MIT)
14 //
15 // Copyright (c) 2014 Andre Azevedo
16 //
17 // Permission is hereby granted, free of charge, to any person obtaining a copy
18 // of this software and associated documentation files (the "Software"), to deal
19 // in the Software without restriction, including without limitation the rights
20 // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
21 // copies of the Software, and to permit persons to whom the Software is
22 // furnished to do so, subject to the following conditions:
23 //
24 // The above copyright notice and this permission notice shall be included in all
25 // copies or substantial portions of the Software.
26 //
27 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
28 // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
29 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
30 // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
31 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
32 // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
33 // SOFTWARE.
34 
35 #pragma once
36 
37 #include "config.h"
38 
39 #include <cmath>
40 #include <cstdint>
41 
42 #ifdef __SIZEOF_INT128__
43 #define BIG_INTEGER_WORD_BITS 128U
44 #define BIG_INTEGER_WORD_POWER 7U
45 #define BIG_INTEGER_WORD unsigned __int128
46 #define BIG_INTEGER_HALF_WORD uint64_t
47 #define BIG_INTEGER_HALF_WORD_MASK 0xFFFFFFFFFFFFFFFFULL
48 #define BIG_INTEGER_HALF_WORD_MASK_NOT 0xFFFFFFFFFFFFFFFF0000000000000000ULL
49 #else
50 #define BIG_INTEGER_WORD_BITS 64U
51 #define BIG_INTEGER_WORD_POWER 6U
52 #define BIG_INTEGER_WORD uint64_t
53 #define BIG_INTEGER_HALF_WORD uint32_t
54 #define BIG_INTEGER_HALF_WORD_MASK 0xFFFFFFFFULL
55 #define BIG_INTEGER_HALF_WORD_MASK_NOT 0xFFFFFFFF00000000ULL
56 #endif
57 
58 // This can be any power of 2 greater than (or equal to) 64:
59 constexpr size_t BIG_INTEGER_BITS = (1 << QBCAPPOW);
61 
62 // The rest of the constants need to be consistent with the one above:
66 
67 typedef struct BigInteger {
69 
70  inline BigInteger()
71  {
72  // Intentionally left blank.
73  }
74 
75  inline BigInteger(const BigInteger& val)
76  {
77  for (int i = 0; i < BIG_INTEGER_WORD_SIZE; ++i) {
78  bits[i] = val.bits[i];
79  }
80  }
81 
82  inline BigInteger(const BIG_INTEGER_WORD& val)
83  {
84  bits[0] = val;
85  for (int i = 1; i < BIG_INTEGER_WORD_SIZE; ++i) {
86  bits[i] = 0U;
87  }
88  }
89 
90 #ifdef __SIZEOF_INT128__
91  inline explicit operator unsigned __int128() const { return (unsigned __int128)bits[0U]; }
92 #endif
93  inline explicit operator uint64_t() const { return (uint64_t)bits[0U]; }
94  inline explicit operator uint32_t() const { return (uint32_t)bits[0U]; }
96 
97 inline void bi_set_0(BigInteger* p)
98 {
99  for (int i = 0; i < BIG_INTEGER_WORD_SIZE; ++i) {
100  p->bits[i] = 0U;
101  }
102 }
103 
104 inline BigInteger bi_copy(const BigInteger& in)
105 {
106  BigInteger result;
107  for (int i = 0; i < BIG_INTEGER_WORD_SIZE; ++i) {
108  result.bits[i] = in.bits[i];
109  }
110  return result;
111 }
112 
113 inline void bi_copy_ip(const BigInteger& in, BigInteger* out)
114 {
115  for (int i = 0; i < BIG_INTEGER_WORD_SIZE; ++i) {
116  out->bits[i] = in.bits[i];
117  }
118 }
119 
120 inline int bi_compare(const BigInteger& left, const BigInteger& right)
121 {
122  for (int i = BIG_INTEGER_MAX_WORD_INDEX; i >= 0; --i) {
123  if (left.bits[i] > right.bits[i]) {
124  return 1;
125  }
126  if (left.bits[i] < right.bits[i]) {
127  return -1;
128  }
129  }
130 
131  return 0;
132 }
133 
134 inline int bi_compare_0(const BigInteger& left)
135 {
136  for (int i = 0; i < BIG_INTEGER_WORD_SIZE; ++i) {
137  if (left.bits[i]) {
138  return 1;
139  }
140  }
141 
142  return 0;
143 }
144 
145 inline int bi_compare_1(const BigInteger& left)
146 {
147  for (int i = BIG_INTEGER_MAX_WORD_INDEX; i > 0; --i) {
148  if (left.bits[i]) {
149  return 1;
150  }
151  }
152  if (left.bits[0] > 1U) {
153  return 1;
154  }
155  if (left.bits[0] < 1U) {
156  return -1;
157  }
158 
159  return 0;
160 }
161 
162 inline BigInteger operator+(const BigInteger& left, const BigInteger& right)
163 {
164  BigInteger result;
165  result.bits[0U] = 0U;
166  for (int i = 0; i < BIG_INTEGER_MAX_WORD_INDEX; ++i) {
167  result.bits[i] += left.bits[i] + right.bits[i];
168  result.bits[i + 1] = (result.bits[i] < left.bits[i]) ? 1 : 0;
169  }
171 
172  return result;
173 }
174 
175 inline void bi_add_ip(BigInteger* left, const BigInteger& right)
176 {
177  for (int i = 0; i < BIG_INTEGER_MAX_WORD_INDEX; ++i) {
178  BIG_INTEGER_WORD temp = left->bits[i];
179  left->bits[i] += right.bits[i];
180  int j = i;
181  while ((j < BIG_INTEGER_MAX_WORD_INDEX) && (left->bits[j] < temp)) {
182  temp = left->bits[++j]++;
183  }
184  }
186 }
187 
188 inline BigInteger operator-(const BigInteger& left, const BigInteger& right)
189 {
190  BigInteger result;
191  result.bits[0U] = 0U;
192  for (int i = 0; i < BIG_INTEGER_MAX_WORD_INDEX; ++i) {
193  result.bits[i] += left.bits[i] - right.bits[i];
194  result.bits[i + 1] = (result.bits[i] > left.bits[i]) ? -1 : 0;
195  }
197 
198  return result;
199 }
200 
201 inline void bi_sub_ip(BigInteger* left, const BigInteger& right)
202 {
203  for (int i = 0; i < BIG_INTEGER_MAX_WORD_INDEX; ++i) {
204  BIG_INTEGER_WORD temp = left->bits[i];
205  left->bits[i] -= right.bits[i];
206  int j = i;
207  while ((j < BIG_INTEGER_MAX_WORD_INDEX) && (left->bits[j] > temp)) {
208  temp = left->bits[++j]--;
209  }
210  }
212 }
213 
214 inline void bi_increment(BigInteger* pBigInt, const BIG_INTEGER_WORD& value)
215 {
216  BIG_INTEGER_WORD temp = pBigInt->bits[0];
217  pBigInt->bits[0] += value;
218  if (temp <= pBigInt->bits[0]) {
219  return;
220  }
221  for (int i = 1; i < BIG_INTEGER_WORD_SIZE; i++) {
222  temp = pBigInt->bits[i]++;
223  if (temp <= pBigInt->bits[i]) {
224  break;
225  }
226  }
227 }
228 
229 inline void bi_decrement(BigInteger* pBigInt, const BIG_INTEGER_WORD& value)
230 {
231  BIG_INTEGER_WORD temp = pBigInt->bits[0];
232  pBigInt->bits[0] -= value;
233  if (temp >= pBigInt->bits[0]) {
234  return;
235  }
236  for (int i = 0; i < BIG_INTEGER_WORD_SIZE; i++) {
237  temp = pBigInt->bits[i]--;
238  if (temp >= pBigInt->bits[i]) {
239  break;
240  }
241  }
242 }
243 
245 {
246  BigInteger result;
247  for (int i = 0; i < BIG_INTEGER_WORD_SIZE; ++i) {
248  result.bits[i] = a[i];
249  }
250 
251  return result;
252 }
253 
254 inline BigInteger bi_lshift_word(const BigInteger& left, BIG_INTEGER_WORD rightMult)
255 {
256  if (!rightMult) {
257  return left;
258  }
259 
260  BigInteger result = 0;
261  for (int i = rightMult; i < BIG_INTEGER_WORD_SIZE; ++i) {
262  result.bits[i] = left.bits[i - rightMult];
263  }
264 
265  return result;
266 }
267 
268 inline void bi_lshift_word_ip(BigInteger* left, BIG_INTEGER_WORD rightMult)
269 {
270  rightMult &= 63U;
271  if (!rightMult) {
272  return;
273  }
274  for (int i = rightMult; i < BIG_INTEGER_WORD_SIZE; ++i) {
275  left->bits[i] = left->bits[i - rightMult];
276  }
277  for (BIG_INTEGER_WORD i = 0U; i < rightMult; ++i) {
278  left->bits[i] = 0U;
279  }
280 }
281 
282 inline BigInteger bi_rshift_word(const BigInteger& left, const BIG_INTEGER_WORD& rightMult)
283 {
284  if (!rightMult) {
285  return left;
286  }
287 
288  BigInteger result = 0U;
289  for (int i = rightMult; i < BIG_INTEGER_WORD_SIZE; ++i) {
290  result.bits[i - rightMult] = left.bits[i];
291  }
292 
293  return result;
294 }
295 
296 inline void bi_rshift_word_ip(BigInteger* left, const BIG_INTEGER_WORD& rightMult)
297 {
298  if (!rightMult) {
299  return;
300  }
301  for (int i = rightMult; i < BIG_INTEGER_WORD_SIZE; ++i) {
302  left->bits[i - rightMult] = left->bits[i];
303  }
304  for (BIG_INTEGER_WORD i = 0U; i < rightMult; ++i) {
305  left->bits[BIG_INTEGER_MAX_WORD_INDEX - i] = 0U;
306  }
307 }
308 
310 {
311  const int rShift64 = right >> BIG_INTEGER_WORD_POWER;
312  const int rMod = right - (rShift64 << BIG_INTEGER_WORD_POWER);
313 
314  BigInteger result = bi_lshift_word(left, rShift64);
315  if (!rMod) {
316  return result;
317  }
318 
319  const int rModComp = BIG_INTEGER_WORD_BITS - rMod;
320  BIG_INTEGER_WORD carry = 0U;
321  for (int i = 0; i < BIG_INTEGER_WORD_SIZE; ++i) {
322  right = result.bits[i];
323  result.bits[i] = carry | (right << rMod);
324  carry = right >> rModComp;
325  }
326 
327  return result;
328 }
329 
330 inline void bi_lshift_ip(BigInteger* left, BIG_INTEGER_WORD right)
331 {
332  const int rShift64 = right >> BIG_INTEGER_WORD_POWER;
333  const int rMod = right - (rShift64 << BIG_INTEGER_WORD_POWER);
334 
335  bi_lshift_word_ip(left, rShift64);
336  if (!rMod) {
337  return;
338  }
339 
340  const int rModComp = BIG_INTEGER_WORD_BITS - rMod;
341  BIG_INTEGER_WORD carry = 0U;
342  for (int i = 0; i < BIG_INTEGER_WORD_SIZE; ++i) {
343  right = left->bits[i];
344  left->bits[i] = carry | (right << rMod);
345  carry = right >> rModComp;
346  }
347 }
348 
350 {
351  const int rShift64 = right >> BIG_INTEGER_WORD_POWER;
352  const int rMod = right - (rShift64 << BIG_INTEGER_WORD_POWER);
353 
354  BigInteger result = bi_rshift_word(left, rShift64);
355  if (!rMod) {
356  return result;
357  }
358 
359  const int rModComp = BIG_INTEGER_WORD_BITS - rMod;
360  BIG_INTEGER_WORD carry = 0U;
361  for (int i = BIG_INTEGER_MAX_WORD_INDEX; i >= 0; --i) {
362  right = result.bits[i];
363  result.bits[i] = carry | (right >> rMod);
364  carry = right << rModComp;
365  }
366 
367  return result;
368 }
369 
370 inline void bi_rshift_ip(BigInteger* left, BIG_INTEGER_WORD right)
371 {
372  const int rShift64 = right >> BIG_INTEGER_WORD_POWER;
373  const int rMod = right - (rShift64 << BIG_INTEGER_WORD_POWER);
374 
375  bi_rshift_word_ip(left, rShift64);
376  if (!rMod) {
377  return;
378  }
379 
380  const int rModComp = BIG_INTEGER_WORD_BITS - rMod;
381  BIG_INTEGER_WORD carry = 0U;
382  for (int i = BIG_INTEGER_MAX_WORD_INDEX; i >= 0; --i) {
383  right = left->bits[i];
384  left->bits[i] = carry | (right >> rMod);
385  carry = right << rModComp;
386  }
387 }
388 
389 inline int bi_log2(const BigInteger& n)
390 {
391  int pw = 0;
392  BigInteger p = n >> 1U;
393  while (bi_compare_0(p) != 0) {
394  bi_rshift_ip(&p, 1U);
395  ++pw;
396  }
397  return pw;
398 }
399 
400 inline int bi_and_1(const BigInteger& left) { return left.bits[0] & 1; }
401 
402 inline BigInteger operator&(const BigInteger& left, const BigInteger& right)
403 {
404  BigInteger result;
405  for (int i = 0; i < BIG_INTEGER_WORD_SIZE; ++i) {
406  result.bits[i] = left.bits[i] & right.bits[i];
407  }
408 
409  return result;
410 }
411 
412 inline void bi_and_ip(BigInteger* left, const BigInteger& right)
413 {
414  for (int i = 0; i < BIG_INTEGER_WORD_SIZE; ++i) {
415  left->bits[i] &= right.bits[i];
416  }
417 }
418 
419 inline BigInteger operator|(const BigInteger& left, const BigInteger& right)
420 {
421  BigInteger result;
422  for (int i = 0; i < BIG_INTEGER_WORD_SIZE; ++i) {
423  result.bits[i] = left.bits[i] | right.bits[i];
424  }
425 
426  return result;
427 }
428 
429 inline void bi_or_ip(BigInteger* left, const BigInteger& right)
430 {
431  for (int i = 0; i < BIG_INTEGER_WORD_SIZE; ++i) {
432  left->bits[i] |= right.bits[i];
433  }
434 }
435 
436 inline BigInteger operator^(const BigInteger& left, const BigInteger& right)
437 {
438  BigInteger result;
439  for (int i = 0; i < BIG_INTEGER_WORD_SIZE; ++i) {
440  result.bits[i] = left.bits[i] ^ right.bits[i];
441  }
442 
443  return result;
444 }
445 
446 inline void bi_xor_ip(BigInteger* left, const BigInteger& right)
447 {
448  for (int i = 0; i < BIG_INTEGER_WORD_SIZE; ++i) {
449  left->bits[i] ^= right.bits[i];
450  }
451 }
452 
453 inline BigInteger operator~(const BigInteger& left)
454 {
455  BigInteger result;
456  for (int i = 0; i < BIG_INTEGER_WORD_SIZE; ++i) {
457  result.bits[i] = ~(left.bits[i]);
458  }
459 
460  return result;
461 }
462 
463 inline void bi_not_ip(BigInteger* left)
464 {
465  for (int i = 0; i < BIG_INTEGER_WORD_SIZE; ++i) {
466  left->bits[i] = ~(left->bits[i]);
467  }
468 }
469 
470 inline double bi_to_double(const BigInteger& in)
471 {
472  double toRet = 0.0;
473  for (int i = 0; i < BIG_INTEGER_WORD_SIZE; ++i) {
474  if (in.bits[i]) {
475  toRet += in.bits[i] * pow(2.0, BIG_INTEGER_WORD_BITS * i);
476  }
477  }
478 
479  return toRet;
480 }
481 
482 inline bool operator==(const BigInteger& left, const BigInteger& right) { return bi_compare(left, right) == 0; }
483 inline bool operator<(const BigInteger& left, const BigInteger& right) { return bi_compare(left, right) < 0; }
484 inline bool operator<=(const BigInteger& left, const BigInteger& right) { return bi_compare(left, right) <= 0; }
485 inline bool operator>(const BigInteger& left, const BigInteger& right) { return bi_compare(left, right) > 0; }
486 inline bool operator>=(const BigInteger& left, const BigInteger& right) { return bi_compare(left, right) >= 0; }
487 
489 {
490  bi_increment(&a, 1U);
491  return a;
492 }
493 
495 {
496  bi_decrement(&a, 1U);
497  return a;
498 }
499 
505 
506 #if true
511 BigInteger operator*(const BigInteger& left, const BigInteger& right);
512 #else
517 BigInteger operator*(const BigInteger& left, const BigInteger& right);
518 #endif
519 
524 void bi_div_mod_small(
525  const BigInteger& left, BIG_INTEGER_HALF_WORD right, BigInteger* quotient, BIG_INTEGER_HALF_WORD* rmndr);
526 
531 void bi_div_mod(const BigInteger& left, const BigInteger& right, BigInteger* quotient, BigInteger* rmndr);
BigInteger operator~(const BigInteger &left)
Definition: big_integer.hpp:453
void bi_or_ip(BigInteger *left, const BigInteger &right)
Definition: big_integer.hpp:429
constexpr size_t BIG_INTEGER_BITS
Definition: big_integer.hpp:59
constexpr size_t BIG_INTEGER_HALF_WORD_BITS
Definition: big_integer.hpp:63
BigInteger operator--(BigInteger &a)
Definition: big_integer.hpp:494
void bi_decrement(BigInteger *pBigInt, const BIG_INTEGER_WORD &value)
Definition: big_integer.hpp:229
BigInteger operator&(const BigInteger &left, const BigInteger &right)
Definition: big_integer.hpp:402
BigInteger operator>>(const BigInteger &left, BIG_INTEGER_WORD right)
Definition: big_integer.hpp:349
void bi_increment(BigInteger *pBigInt, const BIG_INTEGER_WORD &value)
Definition: big_integer.hpp:214
void bi_div_mod_small(const BigInteger &left, BIG_INTEGER_HALF_WORD right, BigInteger *quotient, BIG_INTEGER_HALF_WORD *rmndr)
"Schoolbook division" (on half words) Complexity - O(x^2)
Definition: big_integer.cpp:179
#define BIG_INTEGER_WORD_POWER
Definition: big_integer.hpp:51
constexpr int BIG_INTEGER_MAX_WORD_INDEX
Definition: big_integer.hpp:65
int bi_and_1(const BigInteger &left)
Definition: big_integer.hpp:400
bool operator>=(const BigInteger &left, const BigInteger &right)
Definition: big_integer.hpp:486
BigInteger operator|(const BigInteger &left, const BigInteger &right)
Definition: big_integer.hpp:419
BigInteger operator+(const BigInteger &left, const BigInteger &right)
Definition: big_integer.hpp:162
void bi_add_ip(BigInteger *left, const BigInteger &right)
Definition: big_integer.hpp:175
bool operator>(const BigInteger &left, const BigInteger &right)
Definition: big_integer.hpp:485
BigInteger bi_rshift_word(const BigInteger &left, const BIG_INTEGER_WORD &rightMult)
Definition: big_integer.hpp:282
BigInteger operator^(const BigInteger &left, const BigInteger &right)
Definition: big_integer.hpp:436
int bi_compare_0(const BigInteger &left)
Definition: big_integer.hpp:134
BigInteger bi_lshift_word(const BigInteger &left, BIG_INTEGER_WORD rightMult)
Definition: big_integer.hpp:254
int bi_compare(const BigInteger &left, const BigInteger &right)
Definition: big_integer.hpp:120
void bi_and_ip(BigInteger *left, const BigInteger &right)
Definition: big_integer.hpp:412
constexpr int BIG_INTEGER_WORD_SIZE
Definition: big_integer.hpp:60
BigInteger operator<<(const BigInteger &left, BIG_INTEGER_WORD right)
Definition: big_integer.hpp:309
#define BIG_INTEGER_WORD_BITS
Definition: big_integer.hpp:50
#define BIG_INTEGER_WORD
Definition: big_integer.hpp:52
void bi_lshift_word_ip(BigInteger *left, BIG_INTEGER_WORD rightMult)
Definition: big_integer.hpp:268
bool operator==(const BigInteger &left, const BigInteger &right)
Definition: big_integer.hpp:482
double bi_to_double(const BigInteger &in)
Definition: big_integer.hpp:470
void bi_copy_ip(const BigInteger &in, BigInteger *out)
Definition: big_integer.hpp:113
BigInteger bi_copy(const BigInteger &in)
Definition: big_integer.hpp:104
void bi_rshift_word_ip(BigInteger *left, const BIG_INTEGER_WORD &rightMult)
Definition: big_integer.hpp:296
void bi_not_ip(BigInteger *left)
Definition: big_integer.hpp:463
int bi_compare_1(const BigInteger &left)
Definition: big_integer.hpp:145
bool operator<(const BigInteger &left, const BigInteger &right)
Definition: big_integer.hpp:483
void bi_set_0(BigInteger *p)
Definition: big_integer.hpp:97
void bi_xor_ip(BigInteger *left, const BigInteger &right)
Definition: big_integer.hpp:446
void bi_sub_ip(BigInteger *left, const BigInteger &right)
Definition: big_integer.hpp:201
BigInteger operator*(const BigInteger &left, BIG_INTEGER_HALF_WORD right)
"Schoolbook multiplication" (on half words) Complexity - O(x^2)
Definition: big_integer.cpp:39
int bi_log2(const BigInteger &n)
Definition: big_integer.hpp:389
constexpr int BIG_INTEGER_HALF_WORD_SIZE
Definition: big_integer.hpp:64
void bi_div_mod(const BigInteger &left, const BigInteger &right, BigInteger *quotient, BigInteger *rmndr)
Adapted from Qrack! (The fundamental algorithm was discovered before.) Complexity - O(log)
Definition: big_integer.cpp:218
BigInteger bi_load(BIG_INTEGER_WORD *a)
Definition: big_integer.hpp:244
void bi_lshift_ip(BigInteger *left, BIG_INTEGER_WORD right)
Definition: big_integer.hpp:330
bool operator<=(const BigInteger &left, const BigInteger &right)
Definition: big_integer.hpp:484
#define BIG_INTEGER_HALF_WORD
Definition: big_integer.hpp:53
struct BigInteger BigInteger
BigInteger operator-(const BigInteger &left, const BigInteger &right)
Definition: big_integer.hpp:188
void bi_rshift_ip(BigInteger *left, BIG_INTEGER_WORD right)
Definition: big_integer.hpp:370
BigInteger operator++(BigInteger &a)
Definition: big_integer.hpp:488
half pow(half x, half y)
Power function.
Definition: half.hpp:3738
MICROSOFT_QUANTUM_DECL void U(_In_ uintq sid, _In_ uintq q, _In_ double theta, _In_ double phi, _In_ double lambda)
(External API) 3-parameter unitary gate
Definition: pinvoke_api.cpp:1562
Definition: big_integer.hpp:67
BigInteger(const BIG_INTEGER_WORD &val)
Definition: big_integer.hpp:82
BIG_INTEGER_WORD bits[BIG_INTEGER_WORD_SIZE]
Definition: big_integer.hpp:68
BigInteger()
Definition: big_integer.hpp:70
BigInteger(const BigInteger &val)
Definition: big_integer.hpp:75