1 // cl_LF internals, transcendental functions
6 #include "cln/integer.h"
7 #include "cln/integer_ring.h"
8 #include "cln/lfloat.h"
9 #include "float/lfloat/cl_LF.h"
13 // Subroutine for evaluating
14 // sum(0 <= n < N, a(n)/b(n) * (p(0)...p(n))/(q(0)...q(n)))
15 // where all the entries are small integers (ideally polynomials in n).
16 // Some of the factors (a,b,p,q) may be omitted. They are then understood to
17 // be 1. This is fast because it groups factors together before multiplying.
18 // Result will be a cl_LF with len digits.
20 // There are various alternative implementations of the same algorithm that
21 // differ in the way the series is represented and, as a consequence, in memory
24 // 1st implementation (series is precomputed entirely)
26 // Vectors p[0..N-1], q[0..N-1], a[0..N-1], b[0..N-1], N.
27 // Some of the vectors (a,b,p,q) can be a NULL pointer, all of its entries
28 // are then understood to be 1.
29 // If given, a vector qs[0..N-1] which the evaluation routine may use to
30 // split off q[n] into q[n]*2^qs[n]. qs may be NULL, in that case no shift
31 // optimizations will be used. (They are worth it only if a significant
32 // amount of multiplication work can be saved by shifts.)
34 // 2nd implemenation (series is computed on demand, as a stream)
35 // In this alternate implementation the series is not represented as a couple
36 // of arrays, but as a method returning each tuple (p(n),q(n),a(n),b(n))
37 // in turn. This is preferrable if the a(n) are big, in order to avoid too
38 // much memory usage at the same time.
39 // The next() function is called N times and is expected to return
40 // (p(n),q(n),a(n),b(n)) for n=0..N-1 in that order.
42 // 3rd implemenation (series is computed on demand and truncated early)
43 // This is like the second implementation, but it coerces the integer factors
44 // to cl_LF of a given length (trunclen) as soon as the integer factor's size
45 // exceeds the size to store the cl_LF. For this to make sense, trunclen must
46 // not be smaller than len. In practice, this can shave off substantially from
47 // the memory consumption but it also bears a potential for rounding errors.
48 // A minimum trunclen that guarantees correctness must be evaluated on a
49 // case-by-case basis.
51 // As a variation, it is sometimes advantageous to factor out from q[0..N-1]
52 // powers of two and put them back in again at the end of the computation by
53 // left-shifting. These are distinguished by a boolean template parameter.
54 // (Some combinations might not be implemented yet.)
57 // In each of the special cases below, none of (a,b,p,q) can be NULL.
59 struct cl_pqab_series {
65 struct cl_pqab_series_term {
71 struct cl_pqab_series_stream {
72 cl_pqab_series_term (*nextfn)(cl_pqab_series_stream&);
73 cl_pqab_series_term next () { return nextfn(*this); }
75 cl_pqab_series_stream (cl_pqab_series_term (*n)(cl_pqab_series_stream&)) : nextfn (n) {}
77 template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_pqab_series& args, uintC len);
78 template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_pqab_series_stream& args, uintC len);
79 template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_pqab_series_stream& args, uintC len, uintC trunclen);
81 struct cl_pqb_series {
86 struct cl_pqb_series_term {
91 struct cl_pqb_series_stream {
92 cl_pqb_series_term (*nextfn)(cl_pqb_series_stream&);
93 cl_pqb_series_term next () { return nextfn(*this); }
95 cl_pqb_series_stream (cl_pqb_series_term (*n)(cl_pqb_series_stream&)) : nextfn (n) {}
97 template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_pqb_series& args, uintC len);
98 template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_pqb_series_stream& args, uintC len);
99 template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_pqb_series_stream& args, uintC len, uintC trunclen);
101 struct cl_pqa_series {
106 struct cl_pqa_series_term {
111 struct cl_pqa_series_stream {
112 cl_pqa_series_term (*nextfn)(cl_pqa_series_stream&);
113 cl_pqa_series_term next () { return nextfn(*this); }
115 cl_pqa_series_stream (cl_pqa_series_term (*n)(cl_pqa_series_stream&)) : nextfn (n) {}
117 template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_pqa_series& args, uintC len);
118 template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_pqa_series_stream& args, uintC len);
119 template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_pqa_series_stream& args, uintC len, uintC trunclen);
121 struct cl_pq_series {
125 struct cl_pq_series_term {
129 struct cl_pq_series_stream {
130 cl_pq_series_term (*nextfn)(cl_pq_series_stream&);
131 cl_pq_series_term next () { return nextfn(*this); }
133 cl_pq_series_stream (cl_pq_series_term (*n)(cl_pq_series_stream&)) : nextfn (n) {}
135 template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_pq_series& args, uintC len);
136 template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_pq_series_stream& args, uintC len);
137 template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_pq_series_stream& args, uintC len, uintC trunclen);
139 struct cl_pab_series {
144 const cl_LF eval_rational_series (uintC N, const cl_pab_series& args, uintC len);
146 struct cl_pb_series {
150 const cl_LF eval_rational_series (uintC N, const cl_pb_series& args, uintC len);
152 struct cl_pa_series {
156 const cl_LF eval_rational_series (uintC N, const cl_pa_series& args, uintC len);
161 const cl_LF eval_rational_series (uintC N, const cl_p_series& args, uintC len);
163 struct cl_qab_series {
168 template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_qab_series& args, uintC len);
170 struct cl_qb_series {
174 struct cl_qb_series_term {
178 struct cl_qb_series_stream {
179 cl_qb_series_term (*nextfn)(cl_qb_series_stream&);
180 cl_qb_series_term next () { return nextfn(*this); }
182 cl_qb_series_stream (cl_qb_series_term (*n)(cl_qb_series_stream&)) : nextfn (n) {}
184 template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_qb_series& args, uintC len);
185 template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_qb_series_stream& args, uintC len);
187 struct cl_qa_series {
191 template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_qa_series& args, uintC len);
196 struct cl_q_series_term {
199 struct cl_q_series_stream {
200 cl_q_series_term (*nextfn)(cl_q_series_stream&);
201 cl_q_series_term next () { return nextfn(*this); }
203 cl_q_series_stream (cl_q_series_term (*n)(cl_q_series_stream&)) : nextfn (n) {}
205 template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_q_series& args, uintC len);
206 template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_q_series_stream& args, uintC len);
208 struct cl_ab_series {
212 const cl_LF eval_rational_series (uintC N, const cl_ab_series& args, uintC len);
217 const cl_LF eval_rational_series (uintC N, const cl_b_series& args, uintC len);
222 const cl_LF eval_rational_series (uintC N, const cl_a_series& args, uintC len);
226 const cl_LF eval_rational_series (uintC N, const cl__series& args, uintC len);
231 // Evaluates S = sum(N1 <= n < N2, (p(N1)...p(n))/(q(N1)...q(n)))
232 // and U = sum(N1 <= n < N2,
233 // (c(N1)/d(N1)+...+c(n)/d(n))*(p(N1)...p(n))/(q(N1)...q(n)))
235 // P = p(N1)...p(N2-1),
236 // Q = q(N1)...q(N2-1),
238 // C/D = c(N1)/d(N1)+...+c(N2-1)/d(N2-1),
240 // all integers. On entry N1 < N2.
241 struct cl_pqcd_series_term {
248 struct cl_pqcd_series_result {
256 struct cl_pqcd_series_stream {
257 cl_pqcd_series_term (*nextfn)(cl_pqcd_series_stream&);
258 cl_pqcd_series_term next () { return nextfn(*this); }
260 cl_pqcd_series_stream( cl_pqcd_series_term (*n)(cl_pqcd_series_stream&)) : nextfn (n) {}
262 void eval_pqcd_series_aux (uintC N, cl_pqcd_series_term* args, cl_pqcd_series_result<cl_I>& Z, bool rightmost = true);
263 void eval_pqcd_series_aux (uintC N, cl_pqcd_series_stream& args, cl_pqcd_series_result<cl_I>& Z, bool rightmost = true);
264 void eval_pqcd_series_aux (uintC N, cl_pqcd_series_stream& args, cl_pqcd_series_result<cl_R>& Z, uintC trunclen, bool rightmost = true);
265 // Ditto, but returns U/S.
266 const cl_LF eval_pqcd_series (uintC N, cl_pqcd_series_term* args, uintC len);
267 const cl_LF eval_pqcd_series (uintC N, cl_pqcd_series_stream& args, uintC len);
268 const cl_LF eval_pqcd_series (uintC N, cl_pqcd_series_stream& args, uintC len, uintC trunclen);
270 // [Special case c(n)=1.]
272 // Evaluates S = sum(N1 <= n < N2, (p(N1)...p(n))/(q(N1)...q(n)))
273 // and U = sum(N1 <= n < N2, (1/d(N1)+...+1/d(n))*(p(N1)...p(n))/(q(N1)...q(n)))
275 // P = p(N1)...p(N2-1),
276 // Q = q(N1)...q(N2-1),
278 // C/D = 1/d(N1)+...+1/d(N2-1),
280 // all integers. On entry N1 < N2.
281 struct cl_pqd_series_term {
287 struct cl_pqd_series_result {
295 struct cl_pqd_series_stream {
296 cl_pqd_series_term (*nextfn)(cl_pqd_series_stream&);
297 cl_pqd_series_term next () { return nextfn(*this); }
299 cl_pqd_series_stream( cl_pqd_series_term (*n)(cl_pqd_series_stream&)) : nextfn (n) {}
301 void eval_pqd_series_aux (uintC N, cl_pqd_series_term* args, cl_pqd_series_result<cl_I>& Z, bool rightmost = true);
302 void eval_pqd_series_aux (uintC N, cl_pqd_series_stream& args, cl_pqd_series_result<cl_I>& Z, bool rightmost = true);
303 void eval_pqd_series_aux (uintC N, cl_pqd_series_stream& args, cl_pqd_series_result<cl_R>& Z, uintC trunclen, bool rightmost = true);
304 // Ditto, but returns U/S.
305 const cl_LF eval_pqd_series (uintC N, cl_pqd_series_term* args, uintC len);
306 const cl_LF eval_pqd_series (uintC N, cl_pqd_series_stream& args, uintC len);
307 const cl_LF eval_pqd_series (uintC N, cl_pqd_series_stream& args, uintC len, uintC trunclen);
309 // Helper function to divide q by the largest s such that 2^s divides q, returns s.
311 pullout_shiftcount(cl_I& q)
322 // Helper function to convert integer of length > trunclen to long float of
323 // length = trunclen.
325 truncate_precision(cl_R& x, uintC trunclen)
327 if (instanceof(x,cl_I_ring) &&
328 integer_length(the<cl_I>(x))>trunclen*intDsize) {
329 x = cl_I_to_LF(the<cl_I>(x),trunclen);
335 #endif /* _CL_LF_TRAN_H */