Quantum++  v0.6
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 std::vector<int> x2contfrac(double x, idx n,
46  idx cut = 1e5)
47 {
48  if (n == 0)
49  throw Exception("qpp::x2contfrac()", Exception::Type::OUT_OF_RANGE);
50 
51  std::vector<int> result;
52 
53  for (idx i = 0; i < n; ++i)
54  {
55  result.push_back(std::llround(std::floor(x)));
56  x = 1 / (x - std::floor(x));
57  if (!std::isfinite(x) || x > cut)
58  return result;
59  }
60 
61  return result;
62 }
63 
74 double contfrac2x(const std::vector<int>& cf, idx n)
75 {
76  if (cf.size() == 0)
77  throw Exception("qpp::contfrac2x()", Exception::Type::ZERO_SIZE);
78 
79  if (n == 0)
80  throw Exception("qpp::contfrac2x()", Exception::Type::OUT_OF_RANGE);
81 
82  if (n > cf.size())
83  n = cf.size();
84 
85  if (n == 1) // degenerate case, integer
86  return cf[0];
87 
88  double tmp = 1. / cf[n - 1];
89  for (idx i = n - 2; i != 0; --i)
90  {
91  tmp = 1. / (tmp + cf[i]);
92  }
93 
94  return cf[0] + tmp;
95 }
96 
104 double contfrac2x(const std::vector<int>& cf)
105 {
106  if (cf.size() == 0)
107  throw Exception("qpp::contfrac2x()", Exception::Type::ZERO_SIZE);
108 
109  if (cf.size() == 1) // degenerate case, integer
110  return cf[0];
111 
112  double tmp = 1. / cf[cf.size() - 1];
113  for (idx i = cf.size() - 2; i != 0; --i)
114  {
115  tmp = 1. / (tmp + cf[i]);
116  }
117 
118  return cf[0] + tmp;
119 }
120 
129 unsigned long long int 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 unsigned long long int gcd(const std::vector<unsigned long long int>& ns)
156 {
157  if (ns.size() == 0)
158  throw Exception("qpp::gcd()", Exception::Type::ZERO_SIZE);
159 
160  unsigned long long int result = ns[0]; // convention: gcd({n}) = n
161  for (idx i = 1; i < ns.size(); ++i)
162  {
163  result = gcd(result, ns[i]);
164  }
165 
166  return result;
167 }
168 
177 unsigned long long int lcm(unsigned long long int m, unsigned long long int n)
178 {
179  if (m == 0 || n == 0)
180  throw Exception("qpp::lcm()", Exception::Type::OUT_OF_RANGE);
181 
182  return m * n / gcd(m, n);
183 }
184 
192 unsigned long long int lcm(const std::vector<unsigned long long int>& ns)
193 {
194  if (ns.size() == 0)
195  throw Exception("qpp::lcm()", Exception::Type::ZERO_SIZE);
196 
197  if (ns.size() == 1) // convention: lcm({n}) = n
198  return ns[0];
199 
200  if (std::find(std::begin(ns), std::end(ns), 0) != std::end(ns))
201  throw Exception("qpp::lcm()", Exception::Type::OUT_OF_RANGE);
202 
203  unsigned long long int prod =
204  std::accumulate(std::begin(ns),
205  std::end(ns),
206  static_cast<unsigned long long int>(1),
207  std::multiplies<unsigned long long int>());
208 
209  return prod / gcd(ns);
210 }
211 
218 std::vector<idx> invperm(const std::vector<idx>& perm)
219 {
220  if (!internal::_check_perm(perm))
221  throw Exception("qpp::invperm()", Exception::Type::PERM_INVALID);
222 
223  // construct the inverse
224  std::vector<idx> result(perm.size());
225  for (idx i = 0; i < perm.size(); ++i)
226  result[perm[i]] = i;
227 
228  return result;
229 }
230 
239 std::vector<idx> compperm(const std::vector<idx>& perm,
240  const std::vector<idx>& sigma)
241 {
242  if (!internal::_check_perm(perm))
243  throw Exception("qpp::compperm()", Exception::Type::PERM_INVALID);
244  if (!internal::_check_perm(sigma))
245  throw Exception("qpp::compperm()", Exception::Type::PERM_INVALID);
246  if (perm.size() != sigma.size())
247  throw Exception("qpp::compperm()", Exception::Type::PERM_INVALID);
248 
249  // construct the composition perm(sigma)
250  std::vector<idx> result(perm.size());
251  for (idx i = 0; i < perm.size(); ++i)
252  result[i] = perm[sigma[i]];
253 
254  return result;
255 }
256 
265 std::vector<unsigned long long int> factors(unsigned long long int n)
266 {
267  if (n == 0 || n == 1)
268  throw Exception("qpp::factors()", Exception::Type::OUT_OF_RANGE);
269 
270  std::vector<unsigned long long int> result;
271  unsigned long long int d = 2;
272 
273  while (n > 1)
274  {
275  while (n % d == 0)
276  {
277  result.push_back(d);
278  n /= d;
279  }
280  ++d;
281  if (d * d > n) // changes the running time from O(n) to O(sqrt(n))
282  {
283  if (n > 1)
284  {
285  result.push_back(n);
286  }
287  break;
288  }
289  }
290 
291  return result;
292 }
293 
302 bool isprime(unsigned long long int n)
303 {
304  if (n == 0 or n == 1)
305  throw Exception("qpp::isprime()", Exception::Type::OUT_OF_RANGE);
306 
307  std::vector<unsigned long long int> facts = factors(n);
308 
309  if (facts.size() == 1)
310  return true;
311 
312  return false;
313 }
314 
315 } /* namespace qpp */
316 
317 #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:177
std::vector< unsigned long long int > factors(unsigned long long int n)
Prime factor decomposition.
Definition: number_theory.h:265
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:302
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:239
double contfrac2x(const std::vector< int > &cf, idx n)
Real representation of a simple continued fraction.
Definition: number_theory.h:74
std::vector< idx > invperm(const std::vector< idx > &perm)
Inverse permutation.
Definition: number_theory.h:218
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