Quantum++  v1.0.0-beta2
C++11 quantum computing library
number_theory.h
Go to the documentation of this file.
1 /*
2  * Quantum++
3  *
4  * Copyright (c) 2013 - 2017 Vlad Gheorghiu (vgheorgh@gmail.com)
5  *
6  * This file is part of Quantum++.
7  *
8  * Quantum++ is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * Quantum++ is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with Quantum++. If not, see <http://www.gnu.org/licenses/>.
20  */
21 
27 #ifndef NUMBER_THEORY_H_
28 #define NUMBER_THEORY_H_
29 
30 // number theory functions
31 
32 namespace qpp
33 {
45 inline std::vector<int> x2contfrac(double x, idx N, idx cut = 1e5)
46 {
47  // EXCEPTION CHECKS
48 
49  if (N == 0)
50  throw Exception("qpp::x2contfrac()", Exception::Type::OUT_OF_RANGE);
51  // END EXCEPTION CHECKS
52 
53  std::vector<int> result;
54 
55  for (idx i = 0; i < N; ++i)
56  {
57  if (x > 0)
58  {
59  result.push_back(std::llround(std::floor(x)));
60  x = 1 / (x - std::floor(x));
61  } else // x < 0
62  {
63  result.push_back(std::llround(std::ceil(x)));
64  x = 1 / (x - std::ceil(x));
65  }
66  if (!std::isfinite(x) || std::abs(x) > cut)
67  return result;
68  }
69 
70  return result;
71 }
72 
83 inline double contfrac2x(const std::vector<int>& cf, idx N)
84 {
85  // EXCEPTION CHECKS
86 
87  if (cf.size() == 0)
88  throw Exception("qpp::contfrac2x()", Exception::Type::ZERO_SIZE);
89 
90  if (N == 0)
91  throw Exception("qpp::contfrac2x()", Exception::Type::OUT_OF_RANGE);
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  {
103  tmp = 1. / (tmp + cf[i]);
104  }
105 
106  return cf[0] + tmp;
107 }
108 
116 inline double contfrac2x(const std::vector<int>& cf)
117 {
118  // EXCEPTION CHECKS
119 
120  if (cf.size() == 0)
121  throw Exception("qpp::contfrac2x()", Exception::Type::ZERO_SIZE);
122  // END EXCEPTION CHECKS
123 
124  if (cf.size() == 1) // degenerate case, integer
125  return cf[0];
126 
127  double tmp = 1. / cf[cf.size() - 1];
128  for (idx i = cf.size() - 2; i != 0; --i)
129  {
130  tmp = 1. / (tmp + cf[i]);
131  }
132 
133  return cf[0] + tmp;
134 }
135 
144 inline bigint gcd(bigint a, bigint b)
145 {
146  // EXCEPTION CHECKS
147 
148  if (a == 0 && b == 0)
149  throw Exception("qpp::gcd()", Exception::Type::OUT_OF_RANGE);
150  // END EXCEPTION CHECKS
151 
152  if (a == 0 || b == 0)
153  return (std::max(std::abs(a), std::abs(b)));
154 
155  bigint result = 1;
156  while (b)
157  {
158  result = b;
159  b = a % result;
160  a = result;
161  }
162 
163  return (result > 0) ? result : -result;
164 }
165 
173 inline bigint gcd(const std::vector<bigint>& as)
174 {
175  // EXCEPTION CHECKS
176 
177  if (as.size() == 0)
178  throw Exception("qpp::gcd()", Exception::Type::ZERO_SIZE);
179  // END EXCEPTION CHECKS
180 
181  bigint result = as[0]; // convention: gcd({a}) = a
182  for (idx i = 1; i < as.size(); ++i)
183  {
184  result = gcd(result, as[i]);
185  }
186 
187  return (result > 0) ? result : -result;
188 }
189 
198 inline bigint lcm(bigint a, bigint b)
199 {
200  // EXCEPTION CHECKS
201 
202  if (a == 0 && b == 0)
203  throw Exception("qpp::lcm()", Exception::Type::OUT_OF_RANGE);
204  // END EXCEPTION CHECKS
205 
206  bigint result = a * b / gcd(a, b);
207 
208  return (result > 0) ? result : -result;
209 }
210 
218 inline bigint lcm(const std::vector<bigint>& as)
219 {
220  // EXCEPTION CHECKS
221 
222  if (as.size() == 0)
223  throw Exception("qpp::lcm()", Exception::Type::ZERO_SIZE);
224 
225  if (as.size() == 1) // convention: lcm({a}) = a
226  return as[0];
227 
228  if (std::find(std::begin(as), std::end(as), 0) != std::end(as))
229  throw Exception("qpp::lcm()", Exception::Type::OUT_OF_RANGE);
230  // END EXCEPTION CHECKS
231 
232  bigint result = as[0]; // convention: lcm({n}) = a
233 
234  for (idx i = 1; i < as.size(); ++i)
235  {
236  result = lcm(result, as[i]);
237  }
238 
239  return (result > 0) ? result : -result;
240 }
241 
248 inline std::vector<idx> invperm(const std::vector<idx>& perm)
249 {
250  // EXCEPTION CHECKS
251 
252  if (!internal::check_perm(perm))
253  throw Exception("qpp::invperm()", Exception::Type::PERM_INVALID);
254  // END EXCEPTION CHECKS
255 
256  // construct the inverse
257  std::vector<idx> result(perm.size());
258  for (idx i = 0; i < perm.size(); ++i)
259  result[perm[i]] = i;
260 
261  return result;
262 }
263 
272 inline std::vector<idx> compperm(const std::vector<idx>& perm,
273  const std::vector<idx>& sigma)
274 {
275  // EXCEPTION CHECKS
276 
277  if (!internal::check_perm(perm))
278  throw Exception("qpp::compperm()", Exception::Type::PERM_INVALID);
279  if (!internal::check_perm(sigma))
280  throw Exception("qpp::compperm()", Exception::Type::PERM_INVALID);
281  if (perm.size() != sigma.size())
282  throw Exception("qpp::compperm()", Exception::Type::PERM_INVALID);
283  // END EXCEPTION CHECKS
284 
285  // construct the composition perm(sigma)
286  std::vector<idx> result(perm.size());
287  for (idx i = 0; i < perm.size(); ++i)
288  result[i] = perm[sigma[i]];
289 
290  return result;
291 }
292 
301 inline std::vector<bigint> factors(bigint a)
302 {
303  // flip the sign if necessary
304  a = std::abs(a);
305 
306  // EXCEPTION CHECKS
307 
308  if (a == 0 || a == 1)
309  throw Exception("qpp::factors()", Exception::Type::OUT_OF_RANGE);
310  // END EXCEPTION CHECKS
311 
312  std::vector<bigint> result;
313  bigint d = 2;
314 
315  while (a > 1)
316  {
317  while (a % d == 0)
318  {
319  result.push_back(d);
320  a /= d;
321  }
322  ++d;
323  if (d * d > a) // changes the running time from O(a) to O(sqrt(a))
324  {
325  if (a > 1)
326  {
327  result.push_back(a);
328  }
329  break;
330  }
331  }
332 
333  return result;
334 }
335 
347 {
348  using ubigint = unsigned long long int;
349 
350  // EXCEPTION CHECKS
351 
352  if (p < 1)
353  throw Exception("qpp::modmul()", Exception::Type::OUT_OF_RANGE);
354  // END EXCEPTION CHECKS
355 
356  if (a == 0 || b == 0)
357  return 0;
358 
359  ubigint ua, ub, up;
360 
361  bool is_positive = true;
362  if (a < 0)
363  {
364  ua = -a;
365  is_positive = false;
366  } else
367  ua = a;
368  if (b < 0)
369  {
370  ub = -b;
371  is_positive = false;
372  } else
373  ub = b;
374 
375  if (a < 0 && b < 0)
376  is_positive = true;
377 
378  up = static_cast<ubigint>(p);
379  ua %= up;
380  ub %= up;
381 
382  // the code below is taken from
383  // http://stackoverflow.com/a/18680280/3093378
384  ubigint res = 0;
385  ubigint temp_b;
386 
387  if (ub > ua)
388  std::swap(ua, ub);
389 
390  /* only needed if un may be >= up */
391  if (ub >= up)
392  {
393  if (up > std::numeric_limits<ubigint>::max() / 2u)
394  ub -= up;
395  else
396  ub %= up;
397  }
398 
399  while (ua != 0)
400  {
401  if (ua & 1)
402  {
403  /* add un to res, modulo p, without overflow */
404  /* equiv to if (res + ub >= p), without overflow */
405  if (ub >= up - res)
406  res -= up;
407  res += ub;
408  }
409  ua >>= 1;
410 
411  /* double b, modulo a */
412  temp_b = ub;
413  if (ub >= up - ub) /* equiv to if (2 * ub >= p), without overflow */
414  temp_b -= up;
415  ub += temp_b;
416  }
417 
418  return is_positive ? static_cast<bigint>(res) :
419  static_cast<bigint>(p - res);
420 }
421 
436 {
437  // EXCEPTION CHECKS
438 
439  if (a < 0 || n < 0 || p < 1)
440  throw Exception("qpp::modpow()", Exception::Type::OUT_OF_RANGE);
441 
442  if (a == 0 && n == 0)
443  throw Exception("qpp::modpow()", Exception::Type::OUT_OF_RANGE);
444  // END EXCEPTION CHECKS
445 
446  if (a == 0 && n > 0)
447  return 0;
448 
449  if (n == 0 && p == 1)
450  return 0;
451 
452  bigint result = 1;
453 
454  for (; n > 0; n /= 2)
455  {
456  if (n % 2)
457  result = modmul(result, a, p) % p; // MULTIPLY
458  a = modmul(a, a, p) % p; // SQUARE
459  }
460 
461  return result;
462 }
463 
474 inline std::tuple<bigint, bigint, bigint> egcd(bigint a, bigint b)
475 {
476  // EXCEPTION CHECKS
477 
478  if (a == 0 && b == 0)
479  throw Exception("qpp::egcd()", Exception::Type::OUT_OF_RANGE);
480  // END EXCEPTION CHECKS
481 
482  bigint m, n, c, q, r;
483  bigint m1 = 0, m2 = 1, n1 = 1, n2 = 0;
484 
485  while (b)
486  {
487  q = a / b, r = a - q * b;
488  m = m2 - q * m1, n = n2 - q * n1;
489  a = b, b = r;
490  m2 = m1, m1 = m, n2 = n1, n1 = n;
491  }
492  c = a, m = m2, n = n2;
493 
494  // correct the signs
495  if (c < 0)
496  {
497  m = -m;
498  n = -n;
499  c = -c;
500  }
501 
502  return std::make_tuple(m, n, c);
503 }
504 
516 {
517  // EXCEPTION CHECKS
518 
519  if (a <= 0 || p <= 0)
520  throw Exception("qpp::modinv()", Exception::Type::OUT_OF_RANGE);
521 
522  bigint x, y;
523  bigint gcd_ap;
524  std::tie(x, y, gcd_ap) = egcd(p, a);
525 
526  if (gcd_ap != 1)
527  throw Exception("qpp::modinv()", Exception::Type::OUT_OF_RANGE);
528  // END EXCEPTION CHECKS
529 
530  return (y > 0) ? y : y + p;
531 }
532 
541 inline bool isprime(bigint p, idx k = 80)
542 {
543  p = std::abs(p);
544 
545  // EXCEPTION CHECKS
546 
547  if (p < 2)
548  throw Exception("qpp::isprime()", Exception::Type::OUT_OF_RANGE);
549  // END EXCEPTION CHECKS
550 
551  if (p == 2 || p == 3)
552  return true;
553 
554 // // perform a Fermat primality test
555  bigint x = rand(2, p - 1);
556  if (modpow(x, p - 1, p) != 1)
557  return false;
558 
559  // compute u and r
560  bigint u = 0, r = 1;
561 
562  // write n − 1 as 2^u * r
563  for (bigint i = p - 1; i % 2 == 0; ++u, i /= 2);
564  r = (p - 1) / static_cast<bigint>(std::llround(std::pow(2, u)));
565 
566  // repeat k times
567  for (idx i = 0; i < k; ++i)
568  {
569  // pick a random integer a in the range [2, p − 2]
570  bigint a = rand(2, p - 2);
571 
572  // set z = a^r mod p
573  bigint z = modpow(a, r, p);
574 
575  if (z == 1 || z == p - 1)
576  continue;
577 
578  // repeat u - 1 times
579  bool jump = false;
580  for (idx j = 0; j < static_cast<idx>(u); ++j)
581  {
582  z = (modmul(z, z, p)) % p;
583  if (z == 1)
584  {
585  // composite
586  return false;
587  }
588  if (z == p - 1)
589  {
590  jump = true;
591  break;
592  }
593  }
594  if (jump)
595  continue;
596 
597  return false;
598  }
599 
600  return true;
601 }
602 
613 // A std::optional<bigint> return type would have been awesome here!
614 inline bigint randprime(bigint a, bigint b, idx N = 1000)
615 {
616  // EXCEPTION CHECKS
617 
618  if (a > b)
619  throw qpp::Exception("qpp::randprime()",
621  // END EXCEPTION CHECKS
622 
623  idx i = 0;
624  for (; i < N; ++i)
625  {
626  // select a candidate
627  bigint candidate = rand(a, b);
628  if (std::abs(candidate) < 2)
629  continue;
630  if (std::abs(candidate) == 2)
631  return candidate;
632 
633  // perform a Fermat test
634  bigint x = rand(2, candidate - 1);
635  if (modpow(x, candidate - 1, candidate) != 1)
636  continue; // candidate fails
637 
638  // passed the Fermat test, perform a Miller-Rabin test
639  if (isprime(candidate))
640  return candidate;
641  }
642 
643  if (i == N)
644  throw qpp::Exception("qpp::randprime()", "Prime not found!");
645 
646  return 0; // so we don't get a warning
647 }
648 
649 } /* namespace qpp */
650 
651 #endif /* NUMBER_THEORY_H_ */
bool isprime(bigint p, idx k=80)
Primality test based on the Miller-Rabin's algorithm.
Definition: number_theory.h:541
bigint modpow(bigint a, bigint n, bigint p)
Fast integer power modulo p based on the SQUARE-AND-MULTIPLY algorithm.
Definition: number_theory.h:435
Quantum++ main namespace.
Definition: codes.h:30
bool check_perm(const std::vector< idx > &perm)
Definition: util.h:257
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:207
std::vector< int > x2contfrac(double x, idx N, idx cut=1e5)
Simple continued fraction expansion.
Definition: number_theory.h:45
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:614
std::vector< idx > compperm(const std::vector< idx > &perm, const std::vector< idx > &sigma)
Compose permutations.
Definition: number_theory.h:272
std::vector< idx > invperm(const std::vector< idx > &perm)
Inverse permutation.
Definition: number_theory.h:248
bigint modinv(bigint a, bigint p)
Modular inverse of a mod p.
Definition: number_theory.h:515
Generates custom exceptions, used when validating function parameters.
Definition: exception.h:39
double rand(double a, double b)
Generates a random real number uniformly distributed in the interval [a, b)
Definition: random.h:44
bigint lcm(bigint a, bigint b)
Least common multiple of two integers.
Definition: number_theory.h:198
long long int bigint
Big integer.
Definition: types.h:41
bigint modmul(bigint a, bigint b, bigint p)
Modular multiplication without overflow.
Definition: number_theory.h:346
std::vector< bigint > factors(bigint a)
Prime factor decomposition.
Definition: number_theory.h:301
std::size_t idx
Non-negative integer index.
Definition: types.h:36
std::tuple< bigint, bigint, bigint > egcd(bigint a, bigint b)
Extended greatest common divisor of two integers.
Definition: number_theory.h:474
double contfrac2x(const std::vector< int > &cf, idx N)
Real representation of a simple continued fraction.
Definition: number_theory.h:83
bigint gcd(bigint a, bigint b)
Greatest common divisor of two integers.
Definition: number_theory.h:144