Quantum++  v1.0-rc2
A modern 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 namespace qpp
31 {
43 inline std::vector<int> x2contfrac(double x, idx N, idx cut = 1e5)
44 {
45  // EXCEPTION CHECKS
46 
47  if (N == 0)
48  throw exception::OutOfRange("qpp::x2contfrac()");
49  // END EXCEPTION CHECKS
50 
51  std::vector<int> result;
52 
53  for (idx i = 0; i < N; ++i)
54  {
55  if (x > 0)
56  {
57  result.push_back(static_cast<int>(std::llround(std::floor(x))));
58  x = 1 / (x - std::floor(x));
59  } else // x < 0
60  {
61  result.push_back(static_cast<int>(std::llround(std::ceil(x))));
62  x = 1 / (x - std::ceil(x));
63  }
64  if (!std::isfinite(x) || std::abs(x) > cut)
65  return result;
66  }
67 
68  return result;
69 }
70 
83 inline double contfrac2x(const std::vector<int>& cf, idx N = idx(-1))
84 {
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  {
103  tmp = 1. / (tmp + cf[i]);
104  }
105 
106  return cf[0] + tmp;
107 }
108 
117 inline bigint gcd(bigint a, bigint b)
118 {
119  // EXCEPTION CHECKS
120 
121  if (a == 0 && b == 0)
122  throw exception::OutOfRange("qpp::gcd()");
123  // END EXCEPTION CHECKS
124 
125  if (a == 0 || b == 0)
126  return (std::max(std::abs(a), std::abs(b)));
127 
128  bigint result = 1;
129  while (b)
130  {
131  result = b;
132  b = a % result;
133  a = result;
134  }
135 
136  return (result > 0) ? result : -result;
137 }
138 
146 inline bigint gcd(const std::vector<bigint>& as)
147 {
148  // EXCEPTION CHECKS
149 
150  if (as.size() == 0)
151  throw exception::ZeroSize("qpp::gcd()");
152  // END EXCEPTION CHECKS
153 
154  bigint result = as[0]; // convention: gcd({a}) = a
155  for (idx i = 1; i < as.size(); ++i)
156  {
157  result = gcd(result, as[i]);
158  }
159 
160  return (result > 0) ? result : -result;
161 }
162 
171 inline bigint lcm(bigint a, bigint b)
172 {
173  // EXCEPTION CHECKS
174 
175  if (a == 0 && b == 0)
176  throw exception::OutOfRange("qpp::lcm()");
177  // END EXCEPTION CHECKS
178 
179  bigint result = a * b / gcd(a, b);
180 
181  return (result > 0) ? result : -result;
182 }
183 
191 inline bigint lcm(const std::vector<bigint>& as)
192 {
193  // EXCEPTION CHECKS
194 
195  if (as.size() == 0)
196  throw exception::ZeroSize("qpp::lcm()");
197 
198  if (as.size() == 1) // convention: lcm({a}) = a
199  return as[0];
200 
201  if (std::find(std::begin(as), std::end(as), 0) != std::end(as))
202  throw exception::OutOfRange("qpp::lcm()");
203  // END EXCEPTION CHECKS
204 
205  bigint result = as[0]; // convention: lcm({n}) = a
206 
207  for (idx i = 1; i < as.size(); ++i)
208  {
209  result = lcm(result, as[i]);
210  }
211 
212  return (result > 0) ? result : -result;
213 }
214 
221 inline std::vector<idx> invperm(const std::vector<idx>& perm)
222 {
223  // EXCEPTION CHECKS
224 
225  if (!internal::check_perm(perm))
226  throw exception::PermInvalid("qpp::invperm()");
227  // END EXCEPTION CHECKS
228 
229  // construct the inverse
230  std::vector<idx> result(perm.size());
231  for (idx i = 0; i < perm.size(); ++i)
232  result[perm[i]] = i;
233 
234  return result;
235 }
236 
245 inline std::vector<idx> compperm(const std::vector<idx>& perm,
246  const std::vector<idx>& sigma)
247 {
248  // EXCEPTION CHECKS
249 
250  if (!internal::check_perm(perm))
251  throw exception::PermInvalid("qpp::compperm()");
252  if (!internal::check_perm(sigma))
253  throw exception::PermInvalid("qpp::compperm()");
254  if (perm.size() != sigma.size())
255  throw exception::PermInvalid("qpp::compperm()");
256  // END EXCEPTION CHECKS
257 
258  // construct the composition perm(sigma)
259  std::vector<idx> result(perm.size());
260  for (idx i = 0; i < perm.size(); ++i)
261  result[i] = perm[sigma[i]];
262 
263  return result;
264 }
265 
274 inline std::vector<bigint> factors(bigint a)
275 {
276  // flip the sign if necessary
277  a = std::abs(a);
278 
279  // EXCEPTION CHECKS
280 
281  if (a == 0 || a == 1)
282  throw exception::OutOfRange("qpp::factors()");
283  // END EXCEPTION CHECKS
284 
285  std::vector<bigint> result;
286  bigint d = 2;
287 
288  while (a > 1)
289  {
290  while (a % d == 0)
291  {
292  result.push_back(d);
293  a /= d;
294  }
295  ++d;
296  if (d * d > a) // changes the running time from O(a) to O(sqrt(a))
297  {
298  if (a > 1)
299  {
300  result.push_back(a);
301  }
302  break;
303  }
304  }
305 
306  return result;
307 }
308 
320 {
321  using ubigint = unsigned long long int;
322 
323  // EXCEPTION CHECKS
324 
325  if (p < 1)
326  throw exception::OutOfRange("qpp::modmul()");
327  // END EXCEPTION CHECKS
328 
329  if (a == 0 || b == 0)
330  return 0;
331 
332  ubigint ua, ub, up;
333 
334  bool is_positive = true;
335  if (a < 0)
336  {
337  ua = -a;
338  is_positive = false;
339  } else
340  ua = a;
341  if (b < 0)
342  {
343  ub = -b;
344  is_positive = false;
345  } else
346  ub = b;
347 
348  if (a < 0 && b < 0)
349  is_positive = true;
350 
351  up = static_cast<ubigint>(p);
352  ua %= up;
353  ub %= up;
354 
355  // the code below is taken from
356  // http://stackoverflow.com/a/18680280/3093378
357  ubigint res = 0;
358  ubigint temp_b;
359 
360  if (ub > ua)
361  std::swap(ua, ub);
362 
363  /* only needed if un may be >= up */
364  if (ub >= up)
365  {
366  if (up > std::numeric_limits<ubigint>::max() / 2u)
367  ub -= up;
368  else
369  ub %= up;
370  }
371 
372  while (ua != 0)
373  {
374  if (ua & 1)
375  {
376  /* add un to res, modulo p, without overflow */
377  /* equiv to if (res + ub >= p), without overflow */
378  if (ub >= up - res)
379  res -= up;
380  res += ub;
381  }
382  ua >>= 1;
383 
384  /* double b, modulo a */
385  temp_b = ub;
386  if (ub >= up - ub) /* equiv to if (2 * ub >= p), without overflow */
387  temp_b -= up;
388  ub += temp_b;
389  }
390 
391  return is_positive ? static_cast<bigint>(res) :
392  static_cast<bigint>(p - res);
393 }
394 
409 {
410  // EXCEPTION CHECKS
411 
412  if (a < 0 || n < 0 || p < 1)
413  throw exception::OutOfRange("qpp::modpow()");
414 
415  if (a == 0 && n == 0)
416  throw exception::OutOfRange("qpp::modpow()");
417  // END EXCEPTION CHECKS
418 
419  if (a == 0 && n > 0)
420  return 0;
421 
422  if (n == 0 && p == 1)
423  return 0;
424 
425  bigint result = 1;
426 
427  for (; n > 0; n /= 2)
428  {
429  if (n % 2)
430  result = modmul(result, a, p) % p; // MULTIPLY
431  a = modmul(a, a, p) % p; // SQUARE
432  }
433 
434  return result;
435 }
436 
447 inline std::tuple<bigint, bigint, bigint> egcd(bigint a, bigint b)
448 {
449  // EXCEPTION CHECKS
450 
451  if (a == 0 && b == 0)
452  throw exception::OutOfRange("qpp::egcd()");
453  // END EXCEPTION CHECKS
454 
455  bigint m, n, c, q, r;
456  bigint m1 = 0, m2 = 1, n1 = 1, n2 = 0;
457 
458  while (b)
459  {
460  q = a / b, r = a - q * b;
461  m = m2 - q * m1, n = n2 - q * n1;
462  a = b, b = r;
463  m2 = m1, m1 = m, n2 = n1, n1 = n;
464  }
465  c = a, m = m2, n = n2;
466 
467  // correct the signs
468  if (c < 0)
469  {
470  m = -m;
471  n = -n;
472  c = -c;
473  }
474 
475  return std::make_tuple(m, n, c);
476 }
477 
489 {
490  // EXCEPTION CHECKS
491 
492  if (a <= 0 || p <= 0)
493  throw exception::OutOfRange("qpp::modinv()");
494 
495  bigint x, y;
496  bigint gcd_ap;
497  std::tie(x, y, gcd_ap) = egcd(p, a);
498 
499  if (gcd_ap != 1)
500  throw exception::OutOfRange("qpp::modinv()");
501  // END EXCEPTION CHECKS
502 
503  return (y > 0) ? y : y + p;
504 }
505 
514 inline bool isprime(bigint p, idx k = 80)
515 {
516  p = std::abs(p);
517 
518  // EXCEPTION CHECKS
519 
520  if (p < 2)
521  throw exception::OutOfRange("qpp::isprime()");
522  // END EXCEPTION CHECKS
523 
524  if (p == 2 || p == 3)
525  return true;
526 
527 // // perform a Fermat primality test
528  bigint x = rand(2, p - 1);
529  if (modpow(x, p - 1, p) != 1)
530  return false;
531 
532  // compute u and r
533  bigint u = 0, r = 1;
534 
535  // write n − 1 as 2^u * r
536  for (bigint i = p - 1; i % 2 == 0; ++u, i /= 2);
537  r = (p - 1) / static_cast<bigint>(std::llround(std::pow(2, u)));
538 
539  // repeat k times
540  for (idx i = 0; i < k; ++i)
541  {
542  // pick a random integer a in the range [2, p − 2]
543  bigint a = rand(2, p - 2);
544 
545  // set z = a^r mod p
546  bigint z = modpow(a, r, p);
547 
548  if (z == 1 || z == p - 1)
549  continue;
550 
551  // repeat u - 1 times
552  bool jump = false;
553  for (idx j = 0; j < static_cast<idx>(u); ++j)
554  {
555  z = (modmul(z, z, p)) % p;
556  if (z == 1)
557  {
558  // composite
559  return false;
560  }
561  if (z == p - 1)
562  {
563  jump = true;
564  break;
565  }
566  }
567  if (jump)
568  continue;
569 
570  return false;
571  }
572 
573  return true;
574 }
575 
586 // A std::optional<bigint> return type would have been awesome here!
587 inline bigint randprime(bigint a, bigint b, idx N = 1000)
588 {
589  // EXCEPTION CHECKS
590 
591  if (a > b)
592  throw exception::OutOfRange("qpp::randprime()");
593  // END EXCEPTION CHECKS
594 
595  idx i = 0;
596  for (; i < N; ++i)
597  {
598  // select a candidate
599  bigint candidate = rand(a, b);
600  if (std::abs(candidate) < 2)
601  continue;
602  if (std::abs(candidate) == 2)
603  return candidate;
604 
605  // perform a Fermat test
606  bigint x = rand(2, candidate - 1);
607  if (modpow(x, candidate - 1, candidate) != 1)
608  continue; // candidate fails
609 
610  // passed the Fermat test, perform a Miller-Rabin test
611  if (isprime(candidate))
612  return candidate;
613  }
614 
615  if (i == N)
616  throw exception::CustomException("qpp::randprime()", "Prime not found!");
617 
618  return 0; // so we don't get a warning
619 }
620 
621 } /* namespace qpp */
622 
623 #endif /* NUMBER_THEORY_H_ */
Custom exception.
Definition: exception.h:636
bool isprime(bigint p, idx k=80)
Primality test based on the Miller-Rabin&#39;s algorithm.
Definition: number_theory.h:514
bigint modpow(bigint a, bigint n, bigint p)
Fast integer power modulo p based on the SQUARE-AND-MULTIPLY algorithm.
Definition: number_theory.h:408
double contfrac2x(const std::vector< int > &cf, idx N=idx(-1))
Real representation of a simple continued fraction.
Definition: number_theory.h:83
Quantum++ main namespace.
Definition: codes.h:30
bool check_perm(const std::vector< idx > &perm)
Definition: util.h:278
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:204
std::vector< int > x2contfrac(double x, idx N, idx cut=1e5)
Simple continued fraction expansion.
Definition: number_theory.h:43
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:587
Invalid permutation exception.
Definition: exception.h:412
std::vector< idx > compperm(const std::vector< idx > &perm, const std::vector< idx > &sigma)
Compose permutations.
Definition: number_theory.h:245
std::vector< idx > invperm(const std::vector< idx > &perm)
Inverse permutation.
Definition: number_theory.h:221
bigint modinv(bigint a, bigint p)
Modular inverse of a mod p.
Definition: number_theory.h:488
double rand(double a, double b)
Generates a random real number uniformly distributed in the interval [a, b)
Definition: random.h:41
Parameter out of range exception.
Definition: exception.h:567
bigint lcm(bigint a, bigint b)
Least common multiple of two integers.
Definition: number_theory.h:171
long long int bigint
Big integer.
Definition: types.h:40
bigint modmul(bigint a, bigint b, bigint p)
Modular multiplication without overflow.
Definition: number_theory.h:319
std::vector< bigint > factors(bigint a)
Prime factor decomposition.
Definition: number_theory.h:274
std::size_t idx
Non-negative integer index.
Definition: types.h:35
std::tuple< bigint, bigint, bigint > egcd(bigint a, bigint b)
Extended greatest common divisor of two integers.
Definition: number_theory.h:447
bigint gcd(bigint a, bigint b)
Greatest common divisor of two integers.
Definition: number_theory.h:117
Object has zero size exception.
Definition: exception.h:134