Quantum++  v0.8.8.2
C++11 quantum computing library
number_theory.h
Go to the documentation of this file.
1 /*
2  * Quantum++
3  *
4  * Copyright (c) 2013 - 2016 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  }
62  else // x < 0
63  {
64  result.push_back(std::llround(std::ceil(x)));
65  x = 1 / (x - std::ceil(x));
66  }
67  if (!std::isfinite(x) || std::abs(x) > cut)
68  return result;
69  }
70 
71  return result;
72 }
73 
84 inline double contfrac2x(const std::vector<int>& cf, idx n)
85 {
86  // EXCEPTION CHECKS
87 
88  if (cf.size() == 0)
89  throw Exception("qpp::contfrac2x()", Exception::Type::ZERO_SIZE);
90 
91  if (n == 0)
92  throw Exception("qpp::contfrac2x()", Exception::Type::OUT_OF_RANGE);
93  // END EXCEPTION CHECKS
94 
95  if (n > cf.size())
96  n = cf.size();
97 
98  if (n == 1) // degenerate case, integer
99  return cf[0];
100 
101  double tmp = 1. / cf[n - 1];
102  for (idx i = n - 2; i != 0; --i)
103  {
104  tmp = 1. / (tmp + cf[i]);
105  }
106 
107  return cf[0] + tmp;
108 }
109 
117 inline double contfrac2x(const std::vector<int>& cf)
118 {
119  // EXCEPTION CHECKS
120 
121  if (cf.size() == 0)
122  throw Exception("qpp::contfrac2x()", Exception::Type::ZERO_SIZE);
123  // END EXCEPTION CHECKS
124 
125  if (cf.size() == 1) // degenerate case, integer
126  return cf[0];
127 
128  double tmp = 1. / cf[cf.size() - 1];
129  for (idx i = cf.size() - 2; i != 0; --i)
130  {
131  tmp = 1. / (tmp + cf[i]);
132  }
133 
134  return cf[0] + tmp;
135 }
136 
145 inline bigint gcd(bigint m, bigint n)
146 {
147  // EXCEPTION CHECKS
148 
149  if (m == 0 && n == 0)
150  throw Exception("qpp::gcd()", Exception::Type::OUT_OF_RANGE);
151  // END EXCEPTION CHECKS
152 
153  if (m == 0 || n == 0)
154  return (std::max(std::abs(m), std::abs(n)));
155 
156  bigint result = 1;
157  while (n)
158  {
159  result = n;
160  n = m % result;
161  m = result;
162  }
163 
164  return (result > 0) ? result : -result;
165 }
166 
174 inline bigint gcd(const std::vector<bigint>& ns)
175 {
176  // EXCEPTION CHECKS
177 
178  if (ns.size() == 0)
179  throw Exception("qpp::gcd()", Exception::Type::ZERO_SIZE);
180  // END EXCEPTION CHECKS
181 
182  bigint result = ns[0]; // convention: gcd({n}) = n
183  for (idx i = 1; i < ns.size(); ++i)
184  {
185  result = gcd(result, ns[i]);
186  }
187 
188  return (result > 0) ? result : -result;
189 }
190 
199 inline bigint lcm(bigint m, bigint n)
200 {
201  if (m == 0 && n == 0)
202  throw Exception("qpp::lcm()", Exception::Type::OUT_OF_RANGE);
203 
204  bigint result = m * n / gcd(m, n);
205 
206  return (result > 0) ? result : -result;
207 }
208 
216 inline bigint lcm(const std::vector<bigint>& ns)
217 {
218  // EXCEPTION CHECKS
219 
220  if (ns.size() == 0)
221  throw Exception("qpp::lcm()", Exception::Type::ZERO_SIZE);
222 
223  if (ns.size() == 1) // convention: lcm({n}) = n
224  return ns[0];
225 
226  if (std::find(std::begin(ns), std::end(ns), 0) != std::end(ns))
227  throw Exception("qpp::lcm()", Exception::Type::OUT_OF_RANGE);
228  // END EXCEPTION CHECKS
229 
230  bigint result = ns[0]; // convention: lcm({n}) = n
231 
232  for (idx i = 1; i < ns.size(); ++i)
233  {
234  result = lcm(result, ns[i]);
235  }
236 
237  return (result > 0) ? result : -result;
238 }
239 
246 inline std::vector<idx> invperm(const std::vector<idx>& perm)
247 {
248  // EXCEPTION CHECKS
249 
250  if (!internal::_check_perm(perm))
251  throw Exception("qpp::invperm()", Exception::Type::PERM_INVALID);
252  // END EXCEPTION CHECKS
253 
254  // construct the inverse
255  std::vector<idx> result(perm.size());
256  for (idx i = 0; i < perm.size(); ++i)
257  result[perm[i]] = i;
258 
259  return result;
260 }
261 
270 inline std::vector<idx> compperm(const std::vector<idx>& perm,
271  const std::vector<idx>& sigma)
272 {
273  // EXCEPTION CHECKS
274 
275  if (!internal::_check_perm(perm))
276  throw Exception("qpp::compperm()", Exception::Type::PERM_INVALID);
277  if (!internal::_check_perm(sigma))
278  throw Exception("qpp::compperm()", Exception::Type::PERM_INVALID);
279  if (perm.size() != sigma.size())
280  throw Exception("qpp::compperm()", Exception::Type::PERM_INVALID);
281  // END EXCEPTION CHECKS
282 
283  // construct the composition perm(sigma)
284  std::vector<idx> result(perm.size());
285  for (idx i = 0; i < perm.size(); ++i)
286  result[i] = perm[sigma[i]];
287 
288  return result;
289 }
290 
299 inline std::vector<bigint> factors(bigint n)
300 {
301  if (n < 0)
302  n = -n;
303 
304  // EXCEPTION CHECKS
305 
306  if (n == 0 || n == 1)
307  throw Exception("qpp::factors()", Exception::Type::OUT_OF_RANGE);
308  // END EXCEPTION CHECKS
309 
310  std::vector<bigint> result;
311  bigint d = 2;
312 
313  while (n > 1)
314  {
315  while (n % d == 0)
316  {
317  result.push_back(d);
318  n /= d;
319  }
320  ++d;
321  if (d * d > n) // changes the running time from O(n) to O(sqrt(n))
322  {
323  if (n > 1)
324  {
325  result.push_back(n);
326  }
327  break;
328  }
329  }
330 
331  return result;
332 }
333 
343 inline bool isprime(bigint n)
344 {
345  if (n < 0)
346  n = -n;
347 
348  // EXCEPTION CHECKS
349 
350  if (n == 0 or n == 1)
351  throw Exception("qpp::isprime()", Exception::Type::OUT_OF_RANGE);
352  // END EXCEPTION CHECKS
353 
354  auto facts = factors(n);
355 
356  if (facts.size() == 1)
357  return true;
358 
359  return false;
360 }
361 
374 {
375  // EXCEPTION CHECKS
376 
377  if (a < 0 || n < 0 || p < 1)
378  throw Exception("qpp::modpow()", Exception::Type::OUT_OF_RANGE);
379 
380  if (a == 0 && n == 0)
381  throw Exception("qpp::modpow()", Exception::Type::OUT_OF_RANGE);
382 
383  // END EXCEPTION CHECKS
384 
385  if (a == 0 && n > 0)
386  return 0;
387 
388  if (n == 0 && p == 1)
389  return 0;
390 
391  bigint result = 1;
392 
393  for (; n > 0; n /= 2)
394  {
395  if (n % 2)
396  result = (result * a) % p; // MULTIPLY
397  a = (a * a) % p; // SQUARE
398  }
399 
400  return result;
401 }
402 
413 inline std::tuple<bigint, bigint, bigint> egcd(bigint m, bigint n)
414 {
415  // EXCEPTION CHECKS
416 
417  if (m == 0 && n == 0)
418  throw Exception("qpp::egcd()", Exception::Type::OUT_OF_RANGE);
419 
420  // END EXCEPTION CHECKS
421 
422  bigint a, b, c;
423 
424  bigint q, r, a1, a2, b1, b2;
425 
426  a2 = 1, a1 = 0, b2 = 0, b1 = 1;
427  while (n)
428  {
429  q = m / n, r = m - q * n;
430  a = a2 - q * a1, b = b2 - q * b1;
431  m = n, n = r;
432  a2 = a1, a1 = a, b2 = b1, b1 = b;
433  }
434  c = m, a = a2, b = b2;
435 
436  // correct the signs
437  if( c < 0)
438  {
439  a = -a;
440  b = -b;
441  c = -c;
442  }
443 
444  return std::make_tuple(a, b, c);
445 }
446 
458 {
459  // EXCEPTION CHECKS
460 
461  if (a == 0 || p == 0 || a < 0 || p < 0)
462  throw Exception("qpp::modinv()", Exception::Type::OUT_OF_RANGE);
463 
464  bigint x, y;
465  bigint gcd_ap;
466  std::tie(x, y, gcd_ap) = egcd(p, a);
467 
468  if (gcd_ap != 1)
469  throw Exception("qpp::modinv()", Exception::Type::OUT_OF_RANGE);
470 
471  // END EXCEPTION CHECKS
472 
473  return (y > 0) ? y : y + p;
474 }
475 
476 } /* namespace qpp */
477 
478 #endif /* NUMBER_THEORY_H_ */
std::vector< bigint > factors(bigint n)
Prime factor decomposition.
Definition: number_theory.h:299
bigint lcm(bigint m, bigint n)
Least common multiple of two integers.
Definition: number_theory.h:199
bool isprime(bigint n)
Primality test.
Definition: number_theory.h:343
bigint modpow(bigint a, bigint n, bigint p)
Fast integer power modulo p based on the SQUARE-AND-MULTIPLY algorithm.
Definition: number_theory.h:373
Quantum++ main namespace.
Definition: codes.h:30
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::tuple< bigint, bigint, bigint > egcd(bigint m, bigint n)
Extended greatest common divisor of two integers.
Definition: number_theory.h:413
std::vector< idx > compperm(const std::vector< idx > &perm, const std::vector< idx > &sigma)
Compose permutations.
Definition: number_theory.h:270
double contfrac2x(const std::vector< int > &cf, idx n)
Real representation of a simple continued fraction.
Definition: number_theory.h:84
std::vector< idx > invperm(const std::vector< idx > &perm)
Inverse permutation.
Definition: number_theory.h:246
bool _check_perm(const std::vector< idx > &perm)
Definition: util.h:242
bigint modinv(bigint a, bigint p)
Modular inverse of a mod p.
Definition: number_theory.h:457
Generates custom exceptions, used when validating function parameters.
Definition: exception.h:39
bigint gcd(bigint m, bigint n)
Greatest common divisor of two integers.
Definition: number_theory.h:145
long long int bigint
Big integer.
Definition: types.h:41
std::vector< int > x2contfrac(double x, idx n, idx cut=1e5)
Simple continued fraction expansion.
Definition: number_theory.h:45
std::size_t idx
Non-negative integer index.
Definition: types.h:36