1 // eval_rational_series().
7 #include "cl_LF_tran.h"
12 #include "cln/lfloat.h"
13 #include "cln/integer.h"
15 #include "cln/exception.h"
21 // Evaluates S = sum(N1 <= n < N2, a(n)/b(n) * (p(N1)...p(n))/(q(N1)...q(n)))
22 // and returns P = p(N1)...p(N2-1), Q = q(N1)...q(N2-1), B = B(N1)...B(N2-1)
23 // and T = B*Q*S (all integers). On entry N1 < N2.
24 // P will not be computed if a NULL pointer is passed.
26 static void eval_pqa_series_aux (uintC N1, uintC N2,
27 const cl_pqa_series& args,
28 cl_I* P, cl_I* Q, cl_I* T)
32 throw runtime_exception(); break;
34 if (P) { *P = args.pv[N1]; }
36 *T = args.av[N1] * args.pv[N1];
39 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
41 *Q = args.qv[N1] * args.qv[N1+1];
42 *T = args.qv[N1+1] * args.av[N1] * args.pv[N1]
43 + args.av[N1+1] * p01;
47 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
48 var cl_I p012 = p01 * args.pv[N1+2];
50 var cl_I q12 = args.qv[N1+1] * args.qv[N1+2];
51 *Q = args.qv[N1] * q12;
52 *T = q12 * args.av[N1] * args.pv[N1]
53 + args.qv[N1+2] * args.av[N1+1] * p01
54 + args.av[N1+2] * p012;
58 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
59 var cl_I p012 = p01 * args.pv[N1+2];
60 var cl_I p0123 = p012 * args.pv[N1+3];
61 if (P) { *P = p0123; }
62 var cl_I q23 = args.qv[N1+2] * args.qv[N1+3];
63 var cl_I q123 = args.qv[N1+1] * q23;
64 *Q = args.qv[N1] * q123;
65 *T = q123 * args.av[N1] * args.pv[N1]
66 + q23 * args.av[N1+1] * p01
67 + args.qv[N1+3] * args.av[N1+2] * p012
68 + args.av[N1+3] * p0123;
72 var uintC Nm = (N1+N2)/2; // midpoint
75 eval_pqa_series_aux(N1,Nm,args,&LP,&LQ,<);
76 // Compute right part.
78 eval_pqa_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RT);
79 // Put together partial results.
80 if (P) { *P = LP*RP; }
82 // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT.
89 static void eval_pqsa_series_aux (uintC N1, uintC N2,
90 const cl_pqa_series& args,
91 cl_I* P, cl_I* Q, uintC* QS, cl_I* T)
95 throw runtime_exception(); break;
97 if (P) { *P = args.pv[N1]; }
100 *T = args.av[N1] * args.pv[N1];
103 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
105 *Q = args.qv[N1] * args.qv[N1+1];
106 *QS = args.qsv[N1] + args.qsv[N1+1];
107 *T = ((args.qv[N1+1] * args.av[N1] * args.pv[N1]) << args.qsv[N1+1])
108 + args.av[N1+1] * p01;
112 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
113 var cl_I p012 = p01 * args.pv[N1+2];
114 if (P) { *P = p012; }
115 var cl_I q12 = args.qv[N1+1] * args.qv[N1+2];
116 *Q = args.qv[N1] * q12;
117 *QS = args.qsv[N1] + args.qsv[N1+1] + args.qsv[N1+2];
118 *T = ((q12 * args.av[N1] * args.pv[N1]) << (args.qsv[N1+1] + args.qsv[N1+2]))
119 + ((args.qv[N1+2] * args.av[N1+1] * p01) << args.qsv[N1+2])
120 + args.av[N1+2] * p012;
124 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
125 var cl_I p012 = p01 * args.pv[N1+2];
126 var cl_I p0123 = p012 * args.pv[N1+3];
127 if (P) { *P = p0123; }
128 var cl_I q23 = args.qv[N1+2] * args.qv[N1+3];
129 var cl_I q123 = args.qv[N1+1] * q23;
130 *Q = args.qv[N1] * q123;
131 *QS = args.qsv[N1] + args.qsv[N1+1] + args.qsv[N1+2] + args.qsv[N1+3];
132 *T = ((((((q123 * args.av[N1] * args.pv[N1]) << args.qsv[N1+1])
133 + q23 * args.av[N1+1] * p01) << args.qsv[N1+2])
134 + args.qv[N1+3] * args.av[N1+2] * p012) << args.qsv[N1+3])
135 + args.av[N1+3] * p0123;
139 var uintC Nm = (N1+N2)/2; // midpoint
140 // Compute left part.
143 eval_pqsa_series_aux(N1,Nm,args,&LP,&LQ,&LQS,<);
144 // Compute right part.
147 eval_pqsa_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RQS,&RT);
148 // Put together partial results.
149 if (P) { *P = LP*RP; }
152 // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT.
153 *T = ((RQ*LT) << RQS) + LP*RT;
159 const cl_LF eval_rational_series (uintC N, const cl_pqa_series& args, uintC len)
162 return cl_I_to_LF(0,len);
165 eval_pqa_series_aux(0,N,args,NULL,&Q,&T);
166 return cl_I_to_LF(T,len) / cl_I_to_LF(Q,len);
168 // Precomputation of the shift counts:
169 // Split qv[n] into qv[n]*2^qsv[n].
171 var cl_I* qp = args.qv;
172 var uintC* qsp = args.qsv;
173 for (var uintC n = 0; n < N; n++, qp++, qsp++) {
174 // Pull out maximal power of 2 out of *qp = args.qv[n].
186 eval_pqsa_series_aux(0,N,args,NULL,&Q,&QS,&T);
187 return cl_I_to_LF(T,len) / scale_float(cl_I_to_LF(Q,len),QS);
191 static void eval_pqa_series_aux (uintC N1, uintC N2,
192 cl_pqa_series_stream& args,
193 cl_I* P, cl_I* Q, cl_I* T)
197 throw runtime_exception(); break;
199 var cl_pqa_series_term v0 = args.next(); // [N1]
200 if (P) { *P = v0.p; }
206 var cl_pqa_series_term v0 = args.next(); // [N1]
207 var cl_pqa_series_term v1 = args.next(); // [N1+1]
208 var cl_I p01 = v0.p * v1.p;
211 *T = v1.q * v0.a * v0.p
216 var cl_pqa_series_term v0 = args.next(); // [N1]
217 var cl_pqa_series_term v1 = args.next(); // [N1+1]
218 var cl_pqa_series_term v2 = args.next(); // [N1+2]
219 var cl_I p01 = v0.p * v1.p;
220 var cl_I p012 = p01 * v2.p;
221 if (P) { *P = p012; }
222 var cl_I q12 = v1.q * v2.q;
224 *T = q12 * v0.a * v0.p
230 var cl_pqa_series_term v0 = args.next(); // [N1]
231 var cl_pqa_series_term v1 = args.next(); // [N1+1]
232 var cl_pqa_series_term v2 = args.next(); // [N1+2]
233 var cl_pqa_series_term v3 = args.next(); // [N1+3]
234 var cl_I p01 = v0.p * v1.p;
235 var cl_I p012 = p01 * v2.p;
236 var cl_I p0123 = p012 * v3.p;
237 if (P) { *P = p0123; }
238 var cl_I q23 = v2.q * v3.q;
239 var cl_I q123 = v1.q * q23;
241 *T = q123 * v0.a * v0.p
248 var uintC Nm = (N1+N2)/2; // midpoint
249 // Compute left part.
251 eval_pqa_series_aux(N1,Nm,args,&LP,&LQ,<);
252 // Compute right part.
254 eval_pqa_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RT);
255 // Put together partial results.
256 if (P) { *P = LP*RP; }
258 // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT.
265 const cl_LF eval_rational_series (uintC N, cl_pqa_series_stream& args, uintC len)
268 return cl_I_to_LF(0,len);
270 eval_pqa_series_aux(0,N,args,NULL,&Q,&T);
271 return cl_I_to_LF(T,len) / cl_I_to_LF(Q,len);
274 static void eval_pqa_series_aux (uintC N1, uintC N2,
275 cl_pqa_series_stream& args,
276 cl_R* P, cl_R* Q, cl_R* T,
281 throw runtime_exception(); break;
283 var cl_pqa_series_term v0 = args.next(); // [N1]
284 if (P) { *P = v0.p; }
290 var cl_pqa_series_term v0 = args.next(); // [N1]
291 var cl_pqa_series_term v1 = args.next(); // [N1+1]
292 var cl_I p01 = v0.p * v1.p;
295 *T = v1.q * v0.a * v0.p
300 var cl_pqa_series_term v0 = args.next(); // [N1]
301 var cl_pqa_series_term v1 = args.next(); // [N1+1]
302 var cl_pqa_series_term v2 = args.next(); // [N1+2]
303 var cl_I p01 = v0.p * v1.p;
304 var cl_I p012 = p01 * v2.p;
305 if (P) { *P = p012; }
306 var cl_I q12 = v1.q * v2.q;
308 *T = q12 * v0.a * v0.p
314 var cl_pqa_series_term v0 = args.next(); // [N1]
315 var cl_pqa_series_term v1 = args.next(); // [N1+1]
316 var cl_pqa_series_term v2 = args.next(); // [N1+2]
317 var cl_pqa_series_term v3 = args.next(); // [N1+3]
318 var cl_I p01 = v0.p * v1.p;
319 var cl_I p012 = p01 * v2.p;
320 var cl_I p0123 = p012 * v3.p;
321 if (P) { *P = p0123; }
322 var cl_I q23 = v2.q * v3.q;
323 var cl_I q123 = v1.q * q23;
325 *T = q123 * v0.a * v0.p
332 var uintC Nm = (N1+N2)/2; // midpoint
333 // Compute left part.
335 eval_pqa_series_aux(N1,Nm,args,&LP,&LQ,<,trunclen);
336 // Compute right part.
338 eval_pqa_series_aux(Nm,N2,args,(P?&RP:(cl_R*)0),&RQ,&RT,trunclen);
339 // Put together partial results.
342 truncate_precision(*P,trunclen);
345 truncate_precision(*Q,trunclen);
346 // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT.
348 truncate_precision(*T,trunclen);
354 const cl_LF eval_rational_series (uintC N, cl_pqa_series_stream& args, uintC len, uintC trunclen)
357 return cl_I_to_LF(0,len);
359 eval_pqa_series_aux(0,N,args,NULL,&Q,&T,trunclen);
360 return cl_R_to_LF(T,len) / cl_R_to_LF(Q,len);
362 // Bit complexity (if p(n), q(n), a(n), b(n) have length O(log(n))):