Quantum++  v0.8
C++11 quantum computing library
number_theory.h
Go to the documentation of this file.
1 /*
2  * Quantum++
3  *
4  * Copyright (c) 2013 - 2015 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  if (n == 0)
48  throw Exception("qpp::x2contfrac()", Exception::Type::OUT_OF_RANGE);
49 
50  std::vector <int> result;
51 
52  for (idx i = 0; i < n; ++i)
53  {
54  result.push_back(std::llround(std::floor(x)));
55  x = 1 / (x - std::floor(x));
56  if (!std::isfinite(x) || x > cut)
57  return result;
58  }
59 
60  return result;
61 }
62 
73 inline double contfrac2x(const std::vector <int>& cf, idx n)
74 {
75  if (cf.size() == 0)
76  throw Exception("qpp::contfrac2x()", Exception::Type::ZERO_SIZE);
77 
78  if (n == 0)
79  throw Exception("qpp::contfrac2x()", Exception::Type::OUT_OF_RANGE);
80 
81  if (n > cf.size())
82  n = cf.size();
83 
84  if (n == 1) // degenerate case, integer
85  return cf[0];
86 
87  double tmp = 1. / cf[n - 1];
88  for (idx i = n - 2; i != 0; --i)
89  {
90  tmp = 1. / (tmp + cf[i]);
91  }
92 
93  return cf[0] + tmp;
94 }
95 
103 inline double contfrac2x(const std::vector <int>& cf)
104 {
105  if (cf.size() == 0)
106  throw Exception("qpp::contfrac2x()", Exception::Type::ZERO_SIZE);
107 
108  if (cf.size() == 1) // degenerate case, integer
109  return cf[0];
110 
111  double tmp = 1. / cf[cf.size() - 1];
112  for (idx i = cf.size() - 2; i != 0; --i)
113  {
114  tmp = 1. / (tmp + cf[i]);
115  }
116 
117  return cf[0] + tmp;
118 }
119 
129 {
130  if (m == 0 && n == 0)
131  throw Exception("qpp::gcd()", Exception::Type::OUT_OF_RANGE);
132 
133  if (m == 0 || n == 0)
134  return (std::max(m, n));
135 
136  ubigint result = 1;
137  while (n)
138  {
139  result = n;
140  n = m % result;
141  m = result;
142  }
143 
144  return result;
145 }
146 
154 inline ubigint gcd(const std::vector <ubigint>& ns)
155 {
156  if (ns.size() == 0)
157  throw Exception("qpp::gcd()", Exception::Type::ZERO_SIZE);
158 
159  ubigint result = ns[0]; // convention: gcd({n}) = n
160  for (idx i = 1; i < ns.size(); ++i)
161  {
162  result = gcd(result, ns[i]);
163  }
164 
165  return result;
166 }
167 
177 {
178  if (m == 0 || n == 0)
179  throw Exception("qpp::lcm()", Exception::Type::OUT_OF_RANGE);
180 
181  return m * n / gcd(m, n);
182 }
183 
191 inline ubigint lcm(const std::vector <ubigint>& ns)
192 {
193  if (ns.size() == 0)
194  throw Exception("qpp::lcm()", Exception::Type::ZERO_SIZE);
195 
196  if (ns.size() == 1) // convention: lcm({n}) = n
197  return ns[0];
198 
199  if (std::find(std::begin(ns), std::end(ns), 0) != std::end(ns))
200  throw Exception("qpp::lcm()", Exception::Type::OUT_OF_RANGE);
201 
202  ubigint prod =
203  std::accumulate(std::begin(ns),
204  std::end(ns),
205  static_cast<ubigint>(1),
206  std::multiplies<ubigint>());
207 
208  return prod / gcd(ns);
209 }
210 
217 inline std::vector <idx> invperm(const std::vector <idx>& perm)
218 {
219  if (!internal::_check_perm(perm))
220  throw Exception("qpp::invperm()", Exception::Type::PERM_INVALID);
221 
222  // construct the inverse
223  std::vector <idx> result(perm.size());
224  for (idx i = 0; i < perm.size(); ++i)
225  result[perm[i]] = i;
226 
227  return result;
228 }
229 
238 inline std::vector <idx> compperm(const std::vector <idx>& perm,
239  const std::vector <idx>& sigma)
240 {
241  if (!internal::_check_perm(perm))
242  throw Exception("qpp::compperm()", Exception::Type::PERM_INVALID);
243  if (!internal::_check_perm(sigma))
244  throw Exception("qpp::compperm()", Exception::Type::PERM_INVALID);
245  if (perm.size() != sigma.size())
246  throw Exception("qpp::compperm()", Exception::Type::PERM_INVALID);
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 <ubigint> factors(ubigint n)
265 {
266  if (n == 0 || n == 1)
267  throw Exception("qpp::factors()", Exception::Type::OUT_OF_RANGE);
268 
269  std::vector <ubigint> result;
270  ubigint d = 2;
271 
272  while (n > 1)
273  {
274  while (n % d == 0)
275  {
276  result.push_back(d);
277  n /= d;
278  }
279  ++d;
280  if (d * d > n) // changes the running time from O(n) to O(sqrt(n))
281  {
282  if (n > 1)
283  {
284  result.push_back(n);
285  }
286  break;
287  }
288  }
289 
290  return result;
291 }
292 
301 inline bool isprime(ubigint n)
302 {
303  if (n == 0 or n == 1)
304  throw Exception("qpp::isprime()", Exception::Type::OUT_OF_RANGE);
305 
306  std::vector <ubigint> facts = factors(n);
307 
308  if (facts.size() == 1)
309  return true;
310 
311  return false;
312 }
313 
325 {
326  if (p == 0 || (a == 0 && n == 0))
327  throw Exception("qpp::modpow()",
329 
330  ubigint result = 1;
331 
332  for (; n > 0; n /= 2)
333  {
334  if (n % 2)
335  result = (result * a) % p;
336  a = (a * a) % p;
337  }
338 
339  return result;
340 }
341 
342 } /* namespace qpp */
343 
344 #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:213
Quantum++ main namespace.
Definition: codes.h:30
ubigint lcm(ubigint m, ubigint n)
Least common multiple of two positive integers.
Definition: number_theory.h:176
ubigint modpow(ubigint a, ubigint n, ubigint p)
Integer power modulo p.
Definition: number_theory.h:324
std::vector< idx > compperm(const std::vector< idx > &perm, const std::vector< idx > &sigma)
Compose permutations.
Definition: number_theory.h:238
double contfrac2x(const std::vector< int > &cf, idx n)
Real representation of a simple continued fraction.
Definition: number_theory.h:73
std::vector< idx > invperm(const std::vector< idx > &perm)
Inverse permutation.
Definition: number_theory.h:217
bool isprime(ubigint n)
Primality test.
Definition: number_theory.h:301
bool _check_perm(const std::vector< idx > &perm)
Definition: util.h:250
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:264
ubigint gcd(ubigint m, ubigint n)
Greatest common divisor of two non-negative integers.
Definition: number_theory.h:128
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