Loading...
Searching...
No Matches
ntheory_funcs.cpp
1#include <symengine/ntheory.h>
2#include <symengine/ntheory_funcs.h>
3#include <symengine/prime_sieve.h>
4
5namespace SymEngine
6{
7
8PrimePi::PrimePi(const RCP<const Basic> &arg) : OneArgFunction(arg)
9{
10 SYMENGINE_ASSIGN_TYPEID()
11 SYMENGINE_ASSERT(is_canonical(arg))
12}
13
14bool PrimePi::is_canonical(const RCP<const Basic> &arg) const
15{
16 if (is_a_Number(*arg) or is_a<Constant>(*arg)) {
17 return false;
18 } else {
19 return true;
20 }
21}
22
23RCP<const Basic> PrimePi::create(const RCP<const Basic> &arg) const
24{
25 return primepi(arg);
26}
27
28RCP<const Basic> primepi(const RCP<const Basic> &arg)
29{
30 if (is_a_Number(*arg)) {
31 if (is_a<NaN>(*arg)) {
32 return arg;
33 } else if (is_a<Infty>(*arg)) {
34 if (down_cast<const Infty &>(*arg).is_negative_infinity()) {
35 return make_rcp<const Integer>(integer_class(0));
36 } else {
37 return arg;
38 }
39 } else if (down_cast<const Number &>(*arg).is_complex()) {
40 throw SymEngineException("Complex can't be passed to primepi!");
41 } else if (down_cast<const Number &>(*arg).is_negative()) {
42 return make_rcp<const Integer>(integer_class(0));
43 }
44 }
45 if (is_a_Number(*arg) or is_a<Constant>(*arg)) {
46 unsigned int num
47 = (unsigned int)down_cast<const Integer &>(*SymEngine::floor(arg))
48 .as_uint();
49 Sieve::iterator pi(num);
50 unsigned long int p = 0;
51 while ((pi.next_prime()) <= num) {
52 p++;
53 }
54 return make_rcp<const Integer>(integer_class(p));
55 }
56 return make_rcp<const PrimePi>(arg);
57}
58
59Primorial::Primorial(const RCP<const Basic> &arg) : OneArgFunction(arg)
60{
61 SYMENGINE_ASSIGN_TYPEID()
62 SYMENGINE_ASSERT(is_canonical(arg))
63}
64
65bool Primorial::is_canonical(const RCP<const Basic> &arg) const
66{
67 if (is_a_Number(*arg) or is_a<Constant>(*arg)) {
68 return false;
69 } else {
70 return true;
71 }
72}
73
74RCP<const Basic> Primorial::create(const RCP<const Basic> &arg) const
75{
76 return primorial(arg);
77}
78
79RCP<const Basic> primorial(const RCP<const Basic> &arg)
80{
81 if (is_a_Number(*arg)) {
82 if (is_a<NaN>(*arg)) {
83 return arg;
84 }
85 if (down_cast<const Number &>(*arg).is_positive()) {
86 if (is_a<Infty>(*arg)) {
87 return arg;
88 }
89 } else {
90 throw SymEngineException(
91 "Only positive numbers are allowed for primorial!");
92 }
93 }
94 if (is_a_Number(*arg) or is_a<Constant>(*arg)) {
95 unsigned long n
96 = down_cast<const Integer &>(*SymEngine::floor(arg)).as_uint();
97 return make_rcp<const Integer>(mp_primorial(n));
98 }
99 return make_rcp<const Primorial>(arg);
100}
101
110RCP<const Basic> polygonal_number(const RCP<const Basic> &s,
111 const RCP<const Basic> &n)
112{
113 if (is_a_Number(*s)) {
114 if (not is_a<Integer>(*s)
115 or not down_cast<const Integer &>(*sub(s, integer(2)))
116 .is_positive()) {
117 throw DomainError("The number of sides of the polygon must be an "
118 "integer greater than 2");
119 }
120 }
121 if (is_a_Number(*n)) {
122 if (not is_a<Integer>(*n)
123 or not down_cast<const Integer &>(*n).is_positive()) {
124 throw DomainError("n must be an integer greater than 0");
125 }
126 }
127
128 if (is_a_Number(*s) and is_a_Number(*n)) {
129 // Optimized numeric calculation
130 auto s_int = down_cast<const Integer &>(*s).as_integer_class();
131 auto n_int = down_cast<const Integer &>(*n).as_integer_class();
132 auto res = mp_polygonal_number(s_int, n_int);
133 return make_rcp<const Integer>(res);
134 }
135
136 RCP<const Integer> m1 = integer(-1);
137 RCP<const Integer> m2 = integer(-2);
138 RCP<const Integer> p2 = integer(2);
139 RCP<const Integer> p4 = integer(4);
140 RCP<const Basic> x = div(
141 add(mul(add(s, m2), pow(n, p2)), mul(add(p4, mul(m1, s)), n)), p2);
142 return x;
143}
144
153RCP<const Basic> principal_polygonal_root(const RCP<const Basic> &s,
154 const RCP<const Basic> &x)
155{
156 if (is_a_Number(*s)) {
157 if (not is_a<Integer>(*s)
158 or not down_cast<const Integer &>(*sub(s, integer(2)))
159 .is_positive()) {
160 throw DomainError("The number of sides of the polygon must be an "
161 "integer greater than 2");
162 }
163 }
164 if (is_a_Number(*x)) {
165 if (not is_a<Integer>(*x)
166 or not down_cast<const Integer &>(*x).is_positive()) {
167 throw DomainError("x must be an integer greater than 0");
168 }
169 }
170
171 if (is_a_Number(*s) and is_a_Number(*x)) {
172 // Optimized numeric calculation
173 auto s_int = down_cast<const Integer &>(*s).as_integer_class();
174 auto x_int = down_cast<const Integer &>(*x).as_integer_class();
175 auto res = mp_principal_polygonal_root(s_int, x_int);
176 return make_rcp<const Integer>(res);
177 }
178
179 RCP<const Integer> m2 = integer(-2);
180 RCP<const Integer> m4 = integer(-4);
181 RCP<const Integer> p2 = integer(2);
182 RCP<const Integer> p8 = integer(8);
183 RCP<const Basic> root
184 = sqrt(add(mul(mul(p8, add(s, m2)), x), pow(add(s, m4), p2)));
185 RCP<const Basic> n = div(add(root, add(s, m4)), mul(p2, add(s, m2)));
186 return n;
187}
188
189} // namespace SymEngine
RCP< const Basic > create(const RCP< const Basic > &arg) const override
Method to construct classes with canonicalization.
Main namespace for SymEngine package.
Definition: add.cpp:19
bool is_a_Number(const Basic &b)
Definition: number.h:130
integer_class mp_polygonal_number(const integer_class &s, const integer_class &n)
Numeric calculation of the n:th s-gonal number.
Definition: ntheory.cpp:1648
RCP< const Basic > div(const RCP< const Basic > &a, const RCP< const Basic > &b)
Division.
Definition: mul.cpp:431
RCP< const Basic > principal_polygonal_root(const RCP< const Basic > &s, const RCP< const Basic > &x)
The principal s-gonal root of x.
RCP< const Basic > sub(const RCP< const Basic > &a, const RCP< const Basic > &b)
Substracts b from a.
Definition: add.cpp:495
RCP< const Basic > mul(const RCP< const Basic > &a, const RCP< const Basic > &b)
Multiplication.
Definition: mul.cpp:352
RCP< const Basic > floor(const RCP< const Basic > &arg)
Canonicalize Floor:
Definition: functions.cpp:611
RCP< const Basic > add(const RCP< const Basic > &a, const RCP< const Basic > &b)
Adds two objects (safely).
Definition: add.cpp:425
integer_class mp_principal_polygonal_root(const integer_class &s, const integer_class &x)
Numeric calculation of the principal s-gonal root of x.
Definition: ntheory.cpp:1665
std::enable_if< std::is_integral< T >::value, RCP< constInteger > >::type integer(T i)
Definition: integer.h:197
RCP< const Basic > polygonal_number(const RCP< const Basic > &s, const RCP< const Basic > &n)
The n:th s-gonal number.