series_flint.h
1 #ifndef SYMENGINE_SERIES_FLINT_H
2 #define SYMENGINE_SERIES_FLINT_H
3 
4 #include <symengine/series.h>
5 #include <symengine/expression.h>
6 #include <symengine/symengine_exception.h>
7 
8 #ifdef HAVE_SYMENGINE_FLINT
9 #include <symengine/flint_wrapper.h>
10 
11 namespace SymEngine
12 {
13 
14 using fqp_t = fmpq_poly_wrapper;
15 // Univariate Rational Coefficient Power SeriesBase using Flint
16 class URatPSeriesFlint
17  : public SeriesBase<fqp_t, fmpq_wrapper, URatPSeriesFlint>
18 {
19 public:
20  URatPSeriesFlint(const fqp_t p, const std::string varname,
21  const unsigned degree);
22  IMPLEMENT_TYPEID(SYMENGINE_URATPSERIESFLINT)
23  virtual int compare(const Basic &o) const;
24  virtual hash_t __hash__() const;
25  virtual RCP<const Basic> as_basic() const;
26  virtual umap_int_basic as_dict() const;
27  virtual RCP<const Basic> get_coeff(int) const;
28 
29  static RCP<const URatPSeriesFlint>
30  series(const RCP<const Basic> &t, const std::string &x, unsigned int prec);
31  static fqp_t var(const std::string &s);
32  static fqp_t convert(const integer_class &x);
33  static fqp_t convert(const rational_class &x);
34  static fqp_t convert(const Rational &x);
35  static fqp_t convert(const Integer &x);
36  static fqp_t convert(const Basic &x);
37  static inline fqp_t mul(const fqp_t &s, const fqp_t &r, unsigned prec)
38  {
39  return s.mullow(r, prec);
40  }
41  static fqp_t pow(const fqp_t &s, int n, unsigned prec);
42  static unsigned ldegree(const fqp_t &s);
43  static inline fmpq_wrapper find_cf(const fqp_t &s, const fqp_t &var,
44  unsigned deg)
45  {
46  return s.get_coeff(deg);
47  }
48  static fmpq_wrapper root(fmpq_wrapper &c, unsigned n);
49  static fqp_t diff(const fqp_t &s, const fqp_t &var);
50  static fqp_t integrate(const fqp_t &s, const fqp_t &var);
51  static fqp_t subs(const fqp_t &s, const fqp_t &var, const fqp_t &r,
52  unsigned prec);
53 
54  static inline fqp_t series_invert(const fqp_t &s, const fqp_t &var,
55  unsigned int prec)
56  {
57  SYMENGINE_ASSERT(not s.get_coeff(0).is_zero());
58  return s.inv_series(prec);
59  }
60  static inline fqp_t series_reverse(const fqp_t &s, const fqp_t &var,
61  unsigned int prec)
62  {
63  SYMENGINE_ASSERT(s.get_coeff(0).is_zero()
64  and not s.get_coeff(1).is_zero());
65  return s.revert_series(prec);
66  }
67  static inline fqp_t series_log(const fqp_t &s, const fqp_t &var,
68  unsigned int prec)
69  {
70  SYMENGINE_ASSERT(s.get_coeff(0).is_one());
71  return s.log_series(prec);
72  }
73  static inline fqp_t series_exp(const fqp_t &s, const fqp_t &var,
74  unsigned int prec)
75  {
76  SYMENGINE_ASSERT(s.get_coeff(0).is_zero());
77  return s.exp_series(prec);
78  }
79  static inline fqp_t series_sin(const fqp_t &s, const fqp_t &var,
80  unsigned int prec)
81  {
82  SYMENGINE_ASSERT(s.get_coeff(0).is_zero());
83  return s.sin_series(prec);
84  }
85 
86  static inline fqp_t series_cos(const fqp_t &s, const fqp_t &var,
87  unsigned int prec)
88  {
89  SYMENGINE_ASSERT(s.get_coeff(0).is_zero());
90  return s.cos_series(prec);
91  }
92 
93  static inline fqp_t series_tan(const fqp_t &s, const fqp_t &var,
94  unsigned int prec)
95  {
96  SYMENGINE_ASSERT(s.get_coeff(0).is_zero());
97  return s.tan_series(prec);
98  }
99  static inline fqp_t series_atan(const fqp_t &s, const fqp_t &var,
100  unsigned int prec)
101  {
102  SYMENGINE_ASSERT(s.get_coeff(0).is_zero());
103  return s.atan_series(prec);
104  }
105  static inline fqp_t series_atanh(const fqp_t &s, const fqp_t &var,
106  unsigned int prec)
107  {
108  SYMENGINE_ASSERT(s.get_coeff(0).is_zero());
109  return s.atanh_series(prec);
110  }
111  static inline fqp_t series_asin(const fqp_t &s, const fqp_t &var,
112  unsigned int prec)
113  {
114  SYMENGINE_ASSERT(s.get_coeff(0).is_zero());
115  return s.asin_series(prec);
116  }
117  static inline fqp_t series_asinh(const fqp_t &s, const fqp_t &var,
118  unsigned int prec)
119  {
120  SYMENGINE_ASSERT(s.get_coeff(0).is_zero());
121  return s.asinh_series(prec);
122  }
123  static inline fqp_t series_acos(const fqp_t &s, const fqp_t &var,
124  unsigned int prec)
125  {
126  throw NotImplementedError("acos() not implemented");
127  }
128  static inline fqp_t series_sinh(const fqp_t &s, const fqp_t &var,
129  unsigned int prec)
130  {
131  SYMENGINE_ASSERT(s.get_coeff(0).is_zero());
132  return s.sinh_series(prec);
133  }
134  static inline fqp_t series_cosh(const fqp_t &s, const fqp_t &var,
135  unsigned int prec)
136  {
137  SYMENGINE_ASSERT(s.get_coeff(0).is_zero());
138  return s.cosh_series(prec);
139  }
140  static inline fqp_t series_tanh(const fqp_t &s, const fqp_t &var,
141  unsigned int prec)
142  {
143  SYMENGINE_ASSERT(s.get_coeff(0).is_zero());
144  return s.tanh_series(prec);
145  }
146  static inline fqp_t series_lambertw(const fqp_t &s, const fqp_t &var,
147  unsigned int prec)
148  {
149  SYMENGINE_ASSERT(s.get_coeff(0).is_zero());
150 
151  fqp_t p1;
152  p1.set_zero();
153 
154  auto steps = step_list(prec);
155  for (const auto step : steps) {
156  const fqp_t e(series_exp(p1, var, step));
157  const fqp_t p2(mul(e, p1, step) - s);
158  const fqp_t p3(
159  series_invert(mul(e, fqp_t(p1 + fqp_t(1)), step), var, step));
160  p1 -= mul(p2, p3, step);
161  }
162  return p1;
163  }
164 
165  static inline fqp_t series_nthroot(const fqp_t &s, int n, const fqp_t &var,
166  unsigned int prec)
167  {
168  fqp_t one;
169  one.set_one();
170  if (n == 0)
171  return one;
172  if (n == 1)
173  return s;
174  if (n == -1)
175  return series_invert(s, var, prec);
176 
177  const unsigned ldeg = ldegree(s);
178  if (ldeg % n != 0) {
179  throw NotImplementedError("Puiseux series not implemented.");
180  }
181  fqp_t ss = s;
182  if (ldeg != 0) {
183  ss = s * pow(var, -ldeg, prec);
184  }
185  fmpq_wrapper ct = find_cf(ss, var, 0);
186  bool do_inv = false;
187  if (n < 0) {
188  n = -n;
189  do_inv = true;
190  }
191 
192  fmpq_wrapper ctroot = root(ct, n);
193  fqp_t res_p = one, sn = fqp_t(ss / ct);
194  auto steps = step_list(prec);
195  for (const auto step : steps) {
196  fqp_t t = mul(pow(res_p, n + 1, step), sn, step);
197  res_p += (res_p - t) / n;
198  }
199  if (ldeg != 0) {
200  res_p *= pow(var, ldeg / n, prec);
201  }
202  if (do_inv)
203  return fqp_t(res_p * ctroot);
204  else
205  return fqp_t(series_invert(res_p, var, prec) * ctroot);
206  }
207 };
208 } // namespace SymEngine
209 
210 #endif // HAVE_SYMENGINE_FLINT
211 
212 #endif // SYMENGINE_SERIES_FLINT_H
#define IMPLEMENT_TYPEID(SYMENGINE_ID)
Inline members and functions.
Definition: basic.h:340
Main namespace for SymEngine package.
Definition: add.cpp:19
RCP< const Basic > mul(const RCP< const Basic > &a, const RCP< const Basic > &b)
Multiplication.
Definition: mul.cpp:358
T pow(T... args)