Quantum++  v1.0-rc3
A modern C++11 quantum computing library
number_theory.h
Go to the documentation of this file.
1 /*
2  * This file is part of Quantum++.
3  *
4  * MIT License
5  *
6  * Copyright (c) 2013 - 2018 Vlad Gheorghiu (vgheorgh@gmail.com)
7  *
8  * Permission is hereby granted, free of charge, to any person obtaining a copy
9  * of this software and associated documentation files (the "Software"), to deal
10  * in the Software without restriction, including without limitation the rights
11  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12  * copies of the Software, and to permit persons to whom the Software is
13  * furnished to do so, subject to the following conditions:
14  *
15  * The above copyright notice and this permission notice shall be included in
16  * all copies or substantial portions of the Software.
17  *
18  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
24  * SOFTWARE.
25  */
26 
32 #ifndef NUMBER_THEORY_H_
33 #define NUMBER_THEORY_H_
34 
35 namespace qpp {
47 inline std::vector<int> x2contfrac(double x, idx N, idx cut = 1e5) {
48  // EXCEPTION CHECKS
49 
50  if (N == 0)
51  throw exception::OutOfRange("qpp::x2contfrac()");
52  // END EXCEPTION CHECKS
53 
54  std::vector<int> result;
55 
56  for (idx i = 0; i < N; ++i) {
57  if (x > 0) {
58  result.push_back(static_cast<int>(std::llround(std::floor(x))));
59  x = 1 / (x - std::floor(x));
60  } else // x < 0
61  {
62  result.push_back(static_cast<int>(std::llround(std::ceil(x))));
63  x = 1 / (x - std::ceil(x));
64  }
65  if (!std::isfinite(x) || std::abs(x) > cut)
66  return result;
67  }
68 
69  return result;
70 }
71 
84 inline double contfrac2x(const std::vector<int>& cf, idx N = idx(-1)) {
85  // EXCEPTION CHECKS
86 
87  if (cf.size() == 0)
88  throw exception::ZeroSize("qpp::contfrac2x()");
89 
90  if (N == 0)
91  throw exception::OutOfRange("qpp::contfrac2x()");
92  // END EXCEPTION CHECKS
93 
94  if (N > cf.size())
95  N = cf.size();
96 
97  if (N == 1) // degenerate case, integer
98  return cf[0];
99 
100  double tmp = 1. / cf[N - 1];
101  for (idx i = N - 2; i != 0; --i) {
102  tmp = 1. / (tmp + cf[i]);
103  }
104 
105  return cf[0] + tmp;
106 }
107 
116 inline bigint gcd(bigint a, bigint b) {
117  // EXCEPTION CHECKS
118 
119  if (a == 0 && b == 0)
120  throw exception::OutOfRange("qpp::gcd()");
121  // END EXCEPTION CHECKS
122 
123  if (a == 0 || b == 0)
124  return (std::max(std::abs(a), std::abs(b)));
125 
126  bigint result = 1;
127  while (b) {
128  result = b;
129  b = a % result;
130  a = result;
131  }
132 
133  return (result > 0) ? result : -result;
134 }
135 
143 inline bigint gcd(const std::vector<bigint>& as) {
144  // EXCEPTION CHECKS
145 
146  if (as.size() == 0)
147  throw exception::ZeroSize("qpp::gcd()");
148  // END EXCEPTION CHECKS
149 
150  bigint result = as[0]; // convention: gcd({a}) = a
151  for (idx i = 1; i < as.size(); ++i) {
152  result = gcd(result, as[i]);
153  }
154 
155  return (result > 0) ? result : -result;
156 }
157 
166 inline bigint lcm(bigint a, bigint b) {
167  // EXCEPTION CHECKS
168 
169  if (a == 0 && b == 0)
170  throw exception::OutOfRange("qpp::lcm()");
171  // END EXCEPTION CHECKS
172 
173  bigint result = a * b / gcd(a, b);
174 
175  return (result > 0) ? result : -result;
176 }
177 
185 inline bigint lcm(const std::vector<bigint>& as) {
186  // EXCEPTION CHECKS
187 
188  if (as.size() == 0)
189  throw exception::ZeroSize("qpp::lcm()");
190 
191  if (as.size() == 1) // convention: lcm({a}) = a
192  return as[0];
193 
194  if (std::find(std::begin(as), std::end(as), 0) != std::end(as))
195  throw exception::OutOfRange("qpp::lcm()");
196  // END EXCEPTION CHECKS
197 
198  bigint result = as[0]; // convention: lcm({n}) = a
199 
200  for (idx i = 1; i < as.size(); ++i) {
201  result = lcm(result, as[i]);
202  }
203 
204  return (result > 0) ? result : -result;
205 }
206 
213 inline std::vector<idx> invperm(const std::vector<idx>& perm) {
214  // EXCEPTION CHECKS
215 
216  if (!internal::check_perm(perm))
217  throw exception::PermInvalid("qpp::invperm()");
218  // END EXCEPTION CHECKS
219 
220  // construct the inverse
221  std::vector<idx> result(perm.size());
222  for (idx i = 0; i < perm.size(); ++i)
223  result[perm[i]] = i;
224 
225  return result;
226 }
227 
236 inline std::vector<idx> compperm(const std::vector<idx>& perm,
237  const std::vector<idx>& sigma) {
238  // EXCEPTION CHECKS
239 
240  if (!internal::check_perm(perm))
241  throw exception::PermInvalid("qpp::compperm()");
242  if (!internal::check_perm(sigma))
243  throw exception::PermInvalid("qpp::compperm()");
244  if (perm.size() != sigma.size())
245  throw exception::PermInvalid("qpp::compperm()");
246  // END EXCEPTION CHECKS
247 
248  // construct the composition perm(sigma)
249  std::vector<idx> result(perm.size());
250  for (idx i = 0; i < perm.size(); ++i)
251  result[i] = perm[sigma[i]];
252 
253  return result;
254 }
255 
264 inline std::vector<bigint> factors(bigint a) {
265  // flip the sign if necessary
266  a = std::abs(a);
267 
268  // EXCEPTION CHECKS
269 
270  if (a == 0 || a == 1)
271  throw exception::OutOfRange("qpp::factors()");
272  // END EXCEPTION CHECKS
273 
274  std::vector<bigint> result;
275  bigint d = 2;
276 
277  while (a > 1) {
278  while (a % d == 0) {
279  result.push_back(d);
280  a /= d;
281  }
282  ++d;
283  if (d * d > a) // changes the running time from O(a) to O(sqrt(a))
284  {
285  if (a > 1) {
286  result.push_back(a);
287  }
288  break;
289  }
290  }
291 
292  return result;
293 }
294 
305 inline bigint modmul(bigint a, bigint b, bigint p) {
306  using ubigint = unsigned long long int;
307 
308  // EXCEPTION CHECKS
309 
310  if (p < 1)
311  throw exception::OutOfRange("qpp::modmul()");
312  // END EXCEPTION CHECKS
313 
314  if (a == 0 || b == 0)
315  return 0;
316 
317  ubigint ua, ub, up;
318 
319  bool is_positive = true;
320  if (a < 0) {
321  ua = -a;
322  is_positive = false;
323  } else
324  ua = a;
325  if (b < 0) {
326  ub = -b;
327  is_positive = false;
328  } else
329  ub = b;
330 
331  if (a < 0 && b < 0)
332  is_positive = true;
333 
334  up = static_cast<ubigint>(p);
335  ua %= up;
336  ub %= up;
337 
338  // the code below is taken from
339  // http://stackoverflow.com/a/18680280/3093378
340  ubigint res = 0;
341  ubigint temp_b;
342 
343  if (ub > ua)
344  std::swap(ua, ub);
345 
346  /* only needed if un may be >= up */
347  if (ub >= up) {
348  if (up > std::numeric_limits<ubigint>::max() / 2u)
349  ub -= up;
350  else
351  ub %= up;
352  }
353 
354  while (ua != 0) {
355  if (ua & 1) {
356  /* add un to res, modulo p, without overflow */
357  /* equiv to if (res + ub >= p), without overflow */
358  if (ub >= up - res)
359  res -= up;
360  res += ub;
361  }
362  ua >>= 1;
363 
364  /* double b, modulo a */
365  temp_b = ub;
366  if (ub >= up - ub) /* equiv to if (2 * ub >= p), without overflow */
367  temp_b -= up;
368  ub += temp_b;
369  }
370 
371  return is_positive ? static_cast<bigint>(res)
372  : static_cast<bigint>(p - res);
373 }
374 
388 inline bigint modpow(bigint a, bigint n, bigint p) {
389  // EXCEPTION CHECKS
390 
391  if (a < 0 || n < 0 || p < 1)
392  throw exception::OutOfRange("qpp::modpow()");
393 
394  if (a == 0 && n == 0)
395  throw exception::OutOfRange("qpp::modpow()");
396  // END EXCEPTION CHECKS
397 
398  if (a == 0 && n > 0)
399  return 0;
400 
401  if (n == 0 && p == 1)
402  return 0;
403 
404  bigint result = 1;
405 
406  for (; n > 0; n /= 2) {
407  if (n % 2)
408  result = modmul(result, a, p) % p; // MULTIPLY
409  a = modmul(a, a, p) % p; // SQUARE
410  }
411 
412  return result;
413 }
414 
425 inline std::tuple<bigint, bigint, bigint> egcd(bigint a, bigint b) {
426  // EXCEPTION CHECKS
427 
428  if (a == 0 && b == 0)
429  throw exception::OutOfRange("qpp::egcd()");
430  // END EXCEPTION CHECKS
431 
432  bigint m, n, c, q, r;
433  bigint m1 = 0, m2 = 1, n1 = 1, n2 = 0;
434 
435  while (b) {
436  q = a / b, r = a - q * b;
437  m = m2 - q * m1, n = n2 - q * n1;
438  a = b, b = r;
439  m2 = m1, m1 = m, n2 = n1, n1 = n;
440  }
441  c = a, m = m2, n = n2;
442 
443  // correct the signs
444  if (c < 0) {
445  m = -m;
446  n = -n;
447  c = -c;
448  }
449 
450  return std::make_tuple(m, n, c);
451 }
452 
463 inline bigint modinv(bigint a, bigint p) {
464  // EXCEPTION CHECKS
465 
466  if (a <= 0 || p <= 0)
467  throw exception::OutOfRange("qpp::modinv()");
468 
469  bigint x, y;
470  bigint gcd_ap;
471  std::tie(x, y, gcd_ap) = egcd(p, a);
472 
473  if (gcd_ap != 1)
474  throw exception::OutOfRange("qpp::modinv()");
475  // END EXCEPTION CHECKS
476 
477  return (y > 0) ? y : y + p;
478 }
479 
488 inline bool isprime(bigint p, idx k = 80) {
489  p = std::abs(p);
490 
491  // EXCEPTION CHECKS
492 
493  if (p < 2)
494  throw exception::OutOfRange("qpp::isprime()");
495  // END EXCEPTION CHECKS
496 
497  if (p == 2 || p == 3)
498  return true;
499 
500  // // perform a Fermat primality test
501  bigint x = rand(2, p - 1);
502  if (modpow(x, p - 1, p) != 1)
503  return false;
504 
505  // compute u and r
506  bigint u = 0, r = 1;
507 
508  // write n − 1 as 2^u * r
509  for (bigint i = p - 1; i % 2 == 0; ++u, i /= 2)
510  ;
511  r = (p - 1) / static_cast<bigint>(std::llround(std::pow(2, u)));
512 
513  // repeat k times
514  for (idx i = 0; i < k; ++i) {
515  // pick a random integer a in the range [2, p − 2]
516  bigint a = rand(2, p - 2);
517 
518  // set z = a^r mod p
519  bigint z = modpow(a, r, p);
520 
521  if (z == 1 || z == p - 1)
522  continue;
523 
524  // repeat u - 1 times
525  bool jump = false;
526  for (idx j = 0; j < static_cast<idx>(u); ++j) {
527  z = (modmul(z, z, p)) % p;
528  if (z == 1) {
529  // composite
530  return false;
531  }
532  if (z == p - 1) {
533  jump = true;
534  break;
535  }
536  }
537  if (jump)
538  continue;
539 
540  return false;
541  }
542 
543  return true;
544 }
545 
556 // A std::optional<bigint> return type would have been awesome here!
557 inline bigint randprime(bigint a, bigint b, idx N = 1000) {
558  // EXCEPTION CHECKS
559 
560  if (a > b)
561  throw exception::OutOfRange("qpp::randprime()");
562  // END EXCEPTION CHECKS
563 
564  idx i = 0;
565  for (; i < N; ++i) {
566  // select a candidate
567  bigint candidate = rand(a, b);
568  if (std::abs(candidate) < 2)
569  continue;
570  if (std::abs(candidate) == 2)
571  return candidate;
572 
573  // perform a Fermat test
574  bigint x = rand(2, candidate - 1);
575  if (modpow(x, candidate - 1, candidate) != 1)
576  continue; // candidate fails
577 
578  // passed the Fermat test, perform a Miller-Rabin test
579  if (isprime(candidate))
580  return candidate;
581  }
582 
583  if (i == N)
584  throw exception::CustomException("qpp::randprime()",
585  "Prime not found!");
586 
587  return 0; // so we don't get a warning
588 }
589 
590 } /* namespace qpp */
591 
592 #endif /* NUMBER_THEORY_H_ */
Custom exception.
Definition: exception.h:569
bool isprime(bigint p, idx k=80)
Primality test based on the Miller-Rabin&#39;s algorithm.
Definition: number_theory.h:488
bigint modpow(bigint a, bigint n, bigint p)
Fast integer power modulo p based on the SQUARE-AND-MULTIPLY algorithm.
Definition: number_theory.h:388
double contfrac2x(const std::vector< int > &cf, idx N=idx(-1))
Real representation of a simple continued fraction.
Definition: number_theory.h:84
Quantum++ main namespace.
Definition: codes.h:35
bool check_perm(const std::vector< idx > &perm)
Definition: util.h:259
double sigma(const std::vector< double > &prob, const Container &X, typename std::enable_if< is_iterable< Container >::value >::type *=nullptr)
Standard deviation.
Definition: statistics.h:197
std::vector< int > x2contfrac(double x, idx N, idx cut=1e5)
Simple continued fraction expansion.
Definition: number_theory.h:47
bigint randprime(bigint a, bigint b, idx N=1000)
Generates a random big prime uniformly distributed in the interval [a, b].
Definition: number_theory.h:557
Invalid permutation exception.
Definition: exception.h:378
std::vector< idx > compperm(const std::vector< idx > &perm, const std::vector< idx > &sigma)
Compose permutations.
Definition: number_theory.h:236
std::vector< idx > invperm(const std::vector< idx > &perm)
Inverse permutation.
Definition: number_theory.h:213
bigint modinv(bigint a, bigint p)
Modular inverse of a mod p.
Definition: number_theory.h:463
double rand(double a, double b)
Generates a random real number uniformly distributed in the interval [a, b)
Definition: random.h:45
Parameter out of range exception.
Definition: exception.h:513
bigint lcm(bigint a, bigint b)
Least common multiple of two integers.
Definition: number_theory.h:166
long long int bigint
Big integer.
Definition: types.h:44
bigint modmul(bigint a, bigint b, bigint p)
Modular multiplication without overflow.
Definition: number_theory.h:305
std::vector< bigint > factors(bigint a)
Prime factor decomposition.
Definition: number_theory.h:264
std::size_t idx
Non-negative integer index.
Definition: types.h:39
std::tuple< bigint, bigint, bigint > egcd(bigint a, bigint b)
Extended greatest common divisor of two integers.
Definition: number_theory.h:425
bigint gcd(bigint a, bigint b)
Greatest common divisor of two integers.
Definition: number_theory.h:116
Object has zero size exception.
Definition: exception.h:132