Quantum++  v0.7
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 
128 inline unsigned long long int
129 gcd(unsigned long long int m, unsigned long long int n)
130 {
131  if (m == 0 && n == 0)
132  throw Exception("qpp::gcd()", Exception::Type::OUT_OF_RANGE);
133 
134  if (m == 0 || n == 0)
135  return (std::max(m, n));
136 
137  unsigned long long int result = 1;
138  while (n)
139  {
140  result = n;
141  n = m % result;
142  m = result;
143  }
144 
145  return result;
146 }
147 
155 inline unsigned long long int
156 gcd(const std::vector<unsigned long long int>& ns)
157 {
158  if (ns.size() == 0)
159  throw Exception("qpp::gcd()", Exception::Type::ZERO_SIZE);
160 
161  unsigned long long int result = ns[0]; // convention: gcd({n}) = n
162  for (idx i = 1; i < ns.size(); ++i)
163  {
164  result = gcd(result, ns[i]);
165  }
166 
167  return result;
168 }
169 
178 inline unsigned long long int
179 lcm(unsigned long long int m, unsigned long long int n)
180 {
181  if (m == 0 || n == 0)
182  throw Exception("qpp::lcm()", Exception::Type::OUT_OF_RANGE);
183 
184  return m * n / gcd(m, n);
185 }
186 
194 inline unsigned long long int
195 lcm(const std::vector<unsigned long long int>& ns)
196 {
197  if (ns.size() == 0)
198  throw Exception("qpp::lcm()", Exception::Type::ZERO_SIZE);
199 
200  if (ns.size() == 1) // convention: lcm({n}) = n
201  return ns[0];
202 
203  if (std::find(std::begin(ns), std::end(ns), 0) != std::end(ns))
204  throw Exception("qpp::lcm()", Exception::Type::OUT_OF_RANGE);
205 
206  unsigned long long int prod =
207  std::accumulate(std::begin(ns),
208  std::end(ns),
209  static_cast<unsigned long long int>(1),
210  std::multiplies<unsigned long long int>());
211 
212  return prod / gcd(ns);
213 }
214 
221 inline std::vector<idx> invperm(const std::vector<idx>& perm)
222 {
223  if (!internal::_check_perm(perm))
224  throw Exception("qpp::invperm()", Exception::Type::PERM_INVALID);
225 
226  // construct the inverse
227  std::vector<idx> result(perm.size());
228  for (idx i = 0; i < perm.size(); ++i)
229  result[perm[i]] = i;
230 
231  return result;
232 }
233 
242 inline std::vector<idx> compperm(const std::vector<idx>& perm,
243  const std::vector<idx>& sigma)
244 {
245  if (!internal::_check_perm(perm))
246  throw Exception("qpp::compperm()", Exception::Type::PERM_INVALID);
247  if (!internal::_check_perm(sigma))
248  throw Exception("qpp::compperm()", Exception::Type::PERM_INVALID);
249  if (perm.size() != sigma.size())
250  throw Exception("qpp::compperm()", Exception::Type::PERM_INVALID);
251 
252  // construct the composition perm(sigma)
253  std::vector<idx> result(perm.size());
254  for (idx i = 0; i < perm.size(); ++i)
255  result[i] = perm[sigma[i]];
256 
257  return result;
258 }
259 
268 inline std::vector<unsigned long long int> factors(unsigned long long int n)
269 {
270  if (n == 0 || n == 1)
271  throw Exception("qpp::factors()", Exception::Type::OUT_OF_RANGE);
272 
273  std::vector<unsigned long long int> result;
274  unsigned long long int d = 2;
275 
276  while (n > 1)
277  {
278  while (n % d == 0)
279  {
280  result.push_back(d);
281  n /= d;
282  }
283  ++d;
284  if (d * d > n) // changes the running time from O(n) to O(sqrt(n))
285  {
286  if (n > 1)
287  {
288  result.push_back(n);
289  }
290  break;
291  }
292  }
293 
294  return result;
295 }
296 
305 inline bool isprime(unsigned long long int n)
306 {
307  if (n == 0 or n == 1)
308  throw Exception("qpp::isprime()", Exception::Type::OUT_OF_RANGE);
309 
310  std::vector<unsigned long long int> facts = factors(n);
311 
312  if (facts.size() == 1)
313  return true;
314 
315  return false;
316 }
317 
318 } /* namespace qpp */
319 
320 #endif /* NUMBER_THEORY_H_ */
unsigned long long int lcm(unsigned long long int m, unsigned long long int n)
Least common multiple of two positive integers.
Definition: number_theory.h:179
std::vector< unsigned long long int > factors(unsigned long long int n)
Prime factor decomposition.
Definition: number_theory.h:268
Derived::Scalar prod(const Eigen::MatrixBase< Derived > &A)
Element-wise product of A.
Definition: functions.h:213
bool isprime(unsigned long long int n)
Primality test.
Definition: number_theory.h:305
Quantum++ main namespace.
Definition: codes.h:30
std::vector< idx > compperm(const std::vector< idx > &perm, const std::vector< idx > &sigma)
Compose permutations.
Definition: number_theory.h:242
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:221
bool _check_perm(const std::vector< idx > &perm)
Definition: util.h:256
Generates custom exceptions, used when validating function parameters.
Definition: exception.h:39
unsigned long long int gcd(unsigned long long int m, unsigned long long int n)
Greatest common divisor of two non-negative integers.
Definition: number_theory.h:129
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