Quantum++  v0.8.6
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  result.push_back(std::llround(std::floor(x)));
58  x = 1 / (x - std::floor(x));
59  if (!std::isfinite(x) || x > cut)
60  return result;
61  }
62 
63  return result;
64 }
65 
76 inline double contfrac2x(const std::vector<int>& cf, idx n)
77 {
78  // EXCEPTION CHECKS
79 
80  if (cf.size() == 0)
81  throw Exception("qpp::contfrac2x()", Exception::Type::ZERO_SIZE);
82 
83  if (n == 0)
84  throw Exception("qpp::contfrac2x()", Exception::Type::OUT_OF_RANGE);
85  // END EXCEPTION CHECKS
86 
87  if (n > cf.size())
88  n = cf.size();
89 
90  if (n == 1) // degenerate case, integer
91  return cf[0];
92 
93  double tmp = 1. / cf[n - 1];
94  for (idx i = n - 2; i != 0; --i)
95  {
96  tmp = 1. / (tmp + cf[i]);
97  }
98 
99  return cf[0] + tmp;
100 }
101 
109 inline double contfrac2x(const std::vector<int>& cf)
110 {
111  // EXCEPTION CHECKS
112 
113  if (cf.size() == 0)
114  throw Exception("qpp::contfrac2x()", Exception::Type::ZERO_SIZE);
115  // END EXCEPTION CHECKS
116 
117  if (cf.size() == 1) // degenerate case, integer
118  return cf[0];
119 
120  double tmp = 1. / cf[cf.size() - 1];
121  for (idx i = cf.size() - 2; i != 0; --i)
122  {
123  tmp = 1. / (tmp + cf[i]);
124  }
125 
126  return cf[0] + tmp;
127 }
128 
138 {
139  // EXCEPTION CHECKS
140 
141  if (m == 0 && n == 0)
142  throw Exception("qpp::gcd()", Exception::Type::OUT_OF_RANGE);
143  // END EXCEPTION CHECKS
144 
145  if (m == 0 || n == 0)
146  return (std::max(m, n));
147 
148  ubigint result = 1;
149  while (n)
150  {
151  result = n;
152  n = m % result;
153  m = result;
154  }
155 
156  return result;
157 }
158 
166 inline ubigint gcd(const std::vector<ubigint>& ns)
167 {
168  // EXCEPTION CHECKS
169 
170  if (ns.size() == 0)
171  throw Exception("qpp::gcd()", Exception::Type::ZERO_SIZE);
172  // END EXCEPTION CHECKS
173 
174  ubigint result = ns[0]; // convention: gcd({n}) = n
175  for (idx i = 1; i < ns.size(); ++i)
176  {
177  result = gcd(result, ns[i]);
178  }
179 
180  return result;
181 }
182 
192 {
193  if (m == 0 || n == 0)
194  throw Exception("qpp::lcm()", Exception::Type::OUT_OF_RANGE);
195 
196  return m * n / gcd(m, n);
197 }
198 
206 inline ubigint lcm(const std::vector<ubigint>& ns)
207 {
208  // EXCEPTION CHECKS
209 
210  if (ns.size() == 0)
211  throw Exception("qpp::lcm()", Exception::Type::ZERO_SIZE);
212 
213  if (ns.size() == 1) // convention: lcm({n}) = n
214  return ns[0];
215 
216  if (std::find(std::begin(ns), std::end(ns), 0) != std::end(ns))
217  throw Exception("qpp::lcm()", Exception::Type::OUT_OF_RANGE);
218  // END EXCEPTION CHECKS
219 
220  ubigint prod =
221  std::accumulate(std::begin(ns),
222  std::end(ns),
223  static_cast<ubigint>(1),
224  std::multiplies<ubigint>());
225 
226  return prod / gcd(ns);
227 }
228 
235 inline std::vector<idx> invperm(const std::vector<idx>& perm)
236 {
237  // EXCEPTION CHECKS
238 
239  if (!internal::_check_perm(perm))
240  throw Exception("qpp::invperm()", Exception::Type::PERM_INVALID);
241  // END EXCEPTION CHECKS
242 
243  // construct the inverse
244  std::vector<idx> result(perm.size());
245  for (idx i = 0; i < perm.size(); ++i)
246  result[perm[i]] = i;
247 
248  return result;
249 }
250 
259 inline std::vector<idx> compperm(const std::vector<idx>& perm,
260  const std::vector<idx>& sigma)
261 {
262  // EXCEPTION CHECKS
263 
264  if (!internal::_check_perm(perm))
265  throw Exception("qpp::compperm()", Exception::Type::PERM_INVALID);
266  if (!internal::_check_perm(sigma))
267  throw Exception("qpp::compperm()", Exception::Type::PERM_INVALID);
268  if (perm.size() != sigma.size())
269  throw Exception("qpp::compperm()", Exception::Type::PERM_INVALID);
270  // END EXCEPTION CHECKS
271 
272  // construct the composition perm(sigma)
273  std::vector<idx> result(perm.size());
274  for (idx i = 0; i < perm.size(); ++i)
275  result[i] = perm[sigma[i]];
276 
277  return result;
278 }
279 
288 inline std::vector<ubigint> factors(ubigint n)
289 {
290  // EXCEPTION CHECKS
291 
292  if (n == 0 || n == 1)
293  throw Exception("qpp::factors()", Exception::Type::OUT_OF_RANGE);
294  // END EXCEPTION CHECKS
295 
296  std::vector<ubigint> result;
297  ubigint d = 2;
298 
299  while (n > 1)
300  {
301  while (n % d == 0)
302  {
303  result.push_back(d);
304  n /= d;
305  }
306  ++d;
307  if (d * d > n) // changes the running time from O(n) to O(sqrt(n))
308  {
309  if (n > 1)
310  {
311  result.push_back(n);
312  }
313  break;
314  }
315  }
316 
317  return result;
318 }
319 
328 inline bool isprime(ubigint n)
329 {
330  // EXCEPTION CHECKS
331 
332  if (n == 0 or n == 1)
333  throw Exception("qpp::isprime()", Exception::Type::OUT_OF_RANGE);
334  // END EXCEPTION CHECKS
335 
336  std::vector<ubigint> facts = factors(n);
337 
338  if (facts.size() == 1)
339  return true;
340 
341  return false;
342 }
343 
355 {
356  // EXCEPTION CHECKS
357 
358  if (p == 0 || (a == 0 && n == 0))
359  throw Exception("qpp::modpow()",
361  // END EXCEPTION CHECKS
362 
363  ubigint result = 1;
364 
365  for (; n > 0; n /= 2)
366  {
367  if (n % 2)
368  result = (result * a) % p;
369  a = (a * a) % p;
370  }
371 
372  return result;
373 }
374 
375 } /* namespace qpp */
376 
377 #endif /* NUMBER_THEORY_H_ */
unsigned long long int ubigint
Non-negative big integer.
Definition: types.h:46
Derived::Scalar prod(const Eigen::MatrixBase< Derived > &A)
Element-wise product of A.
Definition: functions.h:231
Quantum++ main namespace.
Definition: codes.h:30
ubigint lcm(ubigint m, ubigint n)
Least common multiple of two positive integers.
Definition: number_theory.h:191
ubigint modpow(ubigint a, ubigint n, ubigint p)
Integer power modulo p.
Definition: number_theory.h:354
std::vector< idx > compperm(const std::vector< idx > &perm, const std::vector< idx > &sigma)
Compose permutations.
Definition: number_theory.h:259
double contfrac2x(const std::vector< int > &cf, idx n)
Real representation of a simple continued fraction.
Definition: number_theory.h:76
std::vector< idx > invperm(const std::vector< idx > &perm)
Inverse permutation.
Definition: number_theory.h:235
double sigma(const std::vector< double > &prob, const Container &X)
Standard deviation.
Definition: statistics.h:201
bool isprime(ubigint n)
Primality test.
Definition: number_theory.h:328
bool _check_perm(const std::vector< idx > &perm)
Definition: util.h:257
Generates custom exceptions, used when validating function parameters.
Definition: exception.h:39
std::vector< ubigint > factors(ubigint n)
Prime factor decomposition.
Definition: number_theory.h:288
ubigint gcd(ubigint m, ubigint n)
Greatest common divisor of two non-negative integers.
Definition: number_theory.h:137
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