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_pqab_series_aux (uintC N1, uintC N2,
27 const cl_pqab_series& args,
28 cl_I* P, cl_I* Q, cl_I* B, cl_I* T)
32 throw runtime_exception(); break;
34 if (P) { *P = args.pv[N1]; }
37 *T = args.av[N1] * args.pv[N1];
40 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
42 *Q = args.qv[N1] * args.qv[N1+1];
43 *B = args.bv[N1] * args.bv[N1+1];
44 *T = args.bv[N1+1] * args.qv[N1+1] * args.av[N1] * args.pv[N1]
45 + args.bv[N1] * args.av[N1+1] * p01;
49 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
50 var cl_I p012 = p01 * args.pv[N1+2];
52 var cl_I q12 = args.qv[N1+1] * args.qv[N1+2];
53 *Q = args.qv[N1] * q12;
54 var cl_I b12 = args.bv[N1+1] * args.bv[N1+2];
55 *B = args.bv[N1] * b12;
56 *T = b12 * q12 * args.av[N1] * args.pv[N1]
57 + args.bv[N1] * (args.bv[N1+2] * args.qv[N1+2] * args.av[N1+1] * p01
58 + args.bv[N1+1] * args.av[N1+2] * p012);
62 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
63 var cl_I p012 = p01 * args.pv[N1+2];
64 var cl_I p0123 = p012 * args.pv[N1+3];
65 if (P) { *P = p0123; }
66 var cl_I q23 = args.qv[N1+2] * args.qv[N1+3];
67 var cl_I q123 = args.qv[N1+1] * q23;
68 *Q = args.qv[N1] * q123;
69 var cl_I b01 = args.bv[N1] * args.bv[N1+1];
70 var cl_I b23 = args.bv[N1+2] * args.bv[N1+3];
72 *T = b23 * (args.bv[N1+1] * q123 * args.av[N1] * args.pv[N1]
73 + args.bv[N1] * q23 * args.av[N1+1] * p01)
74 + b01 * (args.bv[N1+3] * args.qv[N1+3] * args.av[N1+2] * p012
75 + args.bv[N1+2] * args.av[N1+3] * p0123);
79 var uintC Nm = (N1+N2)/2; // midpoint
81 var cl_I LP, LQ, LB, LT;
82 eval_pqab_series_aux(N1,Nm,args,&LP,&LQ,&LB,<);
83 // Compute right part.
84 var cl_I RP, RQ, RB, RT;
85 eval_pqab_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RB,&RT);
86 // Put together partial results.
87 if (P) { *P = LP*RP; }
90 // S = LS + LP/LQ * RS, so T = RB*RQ*LT + LB*LP*RT.
91 *T = RB*RQ*LT + LB*LP*RT;
97 static void eval_pqsab_series_aux (uintC N1, uintC N2,
98 const cl_pqab_series& args,
99 cl_I* P, cl_I* Q, uintC* QS, cl_I* B, cl_I* T)
103 throw runtime_exception(); break;
105 if (P) { *P = args.pv[N1]; }
109 *T = args.av[N1] * args.pv[N1];
112 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
114 *Q = args.qv[N1] * args.qv[N1+1];
115 *QS = args.qsv[N1] + args.qsv[N1+1];
116 *B = args.bv[N1] * args.bv[N1+1];
117 *T = ((args.bv[N1+1] * args.qv[N1+1] * args.av[N1] * args.pv[N1]) << args.qsv[N1+1])
118 + args.bv[N1] * args.av[N1+1] * p01;
122 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
123 var cl_I p012 = p01 * args.pv[N1+2];
124 if (P) { *P = p012; }
125 var cl_I q12 = args.qv[N1+1] * args.qv[N1+2];
126 *Q = args.qv[N1] * q12;
127 *QS = args.qsv[N1] + args.qsv[N1+1] + args.qsv[N1+2];
128 var cl_I b12 = args.bv[N1+1] * args.bv[N1+2];
129 *B = args.bv[N1] * b12;
130 *T = ((b12 * q12 * args.av[N1] * args.pv[N1]) << (args.qsv[N1+1] + args.qsv[N1+2]))
131 + args.bv[N1] * (((args.bv[N1+2] * args.qv[N1+2] * args.av[N1+1] * p01) << args.qsv[N1+2])
132 + args.bv[N1+1] * args.av[N1+2] * p012);
136 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
137 var cl_I p012 = p01 * args.pv[N1+2];
138 var cl_I p0123 = p012 * args.pv[N1+3];
139 if (P) { *P = p0123; }
140 var cl_I q23 = args.qv[N1+2] * args.qv[N1+3];
141 var cl_I q123 = args.qv[N1+1] * q23;
142 *Q = args.qv[N1] * q123;
143 *QS = args.qsv[N1] + args.qsv[N1+1] + args.qsv[N1+2] + args.qsv[N1+3];
144 var cl_I b01 = args.bv[N1] * args.bv[N1+1];
145 var cl_I b23 = args.bv[N1+2] * args.bv[N1+3];
147 *T = ((b23 * (((args.bv[N1+1] * q123 * args.av[N1] * args.pv[N1]) << args.qsv[N1+1])
148 + args.bv[N1] * q23 * args.av[N1+1] * p01)) << (args.qsv[N1+2] + args.qsv[N1+3]))
149 + b01 * (((args.bv[N1+3] * args.qv[N1+3] * args.av[N1+2] * p012) << args.qsv[N1+3])
150 + args.bv[N1+2] * args.av[N1+3] * p0123);
154 var uintC Nm = (N1+N2)/2; // midpoint
155 // Compute left part.
156 var cl_I LP, LQ, LB, LT;
158 eval_pqsab_series_aux(N1,Nm,args,&LP,&LQ,&LQS,&LB,<);
159 // Compute right part.
160 var cl_I RP, RQ, RB, RT;
162 eval_pqsab_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RQS,&RB,&RT);
163 // Put together partial results.
164 if (P) { *P = LP*RP; }
168 // S = LS + LP/LQ * RS, so T = RB*RQ*LT + LB*LP*RT.
169 *T = ((RB*RQ*LT) << RQS) + LB*LP*RT;
175 const cl_LF eval_rational_series (uintC N, const cl_pqab_series& args, uintC len)
178 return cl_I_to_LF(0,len);
181 eval_pqab_series_aux(0,N,args,NULL,&Q,&B,&T);
182 return cl_I_to_LF(T,len) / cl_I_to_LF(B*Q,len);
184 // Precomputation of the shift counts:
185 // Split qv[n] into qv[n]*2^qsv[n].
187 var cl_I* qp = args.qv;
188 var uintC* qsp = args.qsv;
189 for (var uintC n = 0; n < N; n++, qp++, qsp++) {
190 // Pull out maximal power of 2 out of *qp = args.qv[n].
202 eval_pqsab_series_aux(0,N,args,NULL,&Q,&QS,&B,&T);
203 return cl_I_to_LF(T,len) / scale_float(cl_I_to_LF(B*Q,len),QS);
207 static void eval_pqab_series_aux (uintC N1, uintC N2,
208 cl_pqab_series_stream& args,
209 cl_I* P, cl_I* Q, cl_I* B, cl_I* T)
213 throw runtime_exception(); break;
215 var cl_pqab_series_term v0 = args.next(); // [N1]
216 if (P) { *P = v0.p; }
223 var cl_pqab_series_term v0 = args.next(); // [N1]
224 var cl_pqab_series_term v1 = args.next(); // [N1+1]
225 var cl_I p01 = v0.p * v1.p;
229 *T = v1.b * v1.q * v0.a * v0.p
234 var cl_pqab_series_term v0 = args.next(); // [N1]
235 var cl_pqab_series_term v1 = args.next(); // [N1+1]
236 var cl_pqab_series_term v2 = args.next(); // [N1+2]
237 var cl_I p01 = v0.p * v1.p;
238 var cl_I p012 = p01 * v2.p;
239 if (P) { *P = p012; }
240 var cl_I q12 = v1.q * v2.q;
242 var cl_I b12 = v1.b * v2.b;
244 *T = b12 * q12 * v0.a * v0.p
245 + v0.b * (v2.b * v2.q * v1.a * p01
246 + v1.b * v2.a * p012);
250 var cl_pqab_series_term v0 = args.next(); // [N1]
251 var cl_pqab_series_term v1 = args.next(); // [N1+1]
252 var cl_pqab_series_term v2 = args.next(); // [N1+2]
253 var cl_pqab_series_term v3 = args.next(); // [N1+3]
254 var cl_I p01 = v0.p * v1.p;
255 var cl_I p012 = p01 * v2.p;
256 var cl_I p0123 = p012 * v3.p;
257 if (P) { *P = p0123; }
258 var cl_I q23 = v2.q * v3.q;
259 var cl_I q123 = v1.q * q23;
261 var cl_I b01 = v0.b * v1.b;
262 var cl_I b23 = v2.b * v3.b;
264 *T = b23 * (v1.b * q123 * v0.a * v0.p
265 + v0.b * q23 * v1.a * p01)
266 + b01 * (v3.b * v3.q * v2.a * p012
267 + v2.b * v3.a * p0123);
271 var uintC Nm = (N1+N2)/2; // midpoint
272 // Compute left part.
273 var cl_I LP, LQ, LB, LT;
274 eval_pqab_series_aux(N1,Nm,args,&LP,&LQ,&LB,<);
275 // Compute right part.
276 var cl_I RP, RQ, RB, RT;
277 eval_pqab_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RB,&RT);
278 // Put together partial results.
279 if (P) { *P = LP*RP; }
282 // S = LS + LP/LQ * RS, so T = RB*RQ*LT + LB*LP*RT.
283 *T = RB*RQ*LT + LB*LP*RT;
289 const cl_LF eval_rational_series (uintC N, cl_pqab_series_stream& args, uintC len)
292 return cl_I_to_LF(0,len);
294 eval_pqab_series_aux(0,N,args,NULL,&Q,&B,&T);
295 return cl_I_to_LF(T,len) / cl_I_to_LF(B*Q,len);
298 static void eval_pqab_series_aux (uintC N1, uintC N2,
299 cl_pqab_series_stream& args,
300 cl_R* P, cl_R* Q, cl_R* B, cl_R* T,
305 throw runtime_exception(); break;
307 var cl_pqab_series_term v0 = args.next(); // [N1]
308 if (P) { *P = v0.p; }
315 var cl_pqab_series_term v0 = args.next(); // [N1]
316 var cl_pqab_series_term v1 = args.next(); // [N1+1]
317 var cl_I p01 = v0.p * v1.p;
321 *T = v1.b * v1.q * v0.a * v0.p
326 var cl_pqab_series_term v0 = args.next(); // [N1]
327 var cl_pqab_series_term v1 = args.next(); // [N1+1]
328 var cl_pqab_series_term v2 = args.next(); // [N1+2]
329 var cl_I p01 = v0.p * v1.p;
330 var cl_I p012 = p01 * v2.p;
331 if (P) { *P = p012; }
332 var cl_I q12 = v1.q * v2.q;
334 var cl_I b12 = v1.b * v2.b;
336 *T = b12 * q12 * v0.a * v0.p
337 + v0.b * (v2.b * v2.q * v1.a * p01
338 + v1.b * v2.a * p012);
342 var cl_pqab_series_term v0 = args.next(); // [N1]
343 var cl_pqab_series_term v1 = args.next(); // [N1+1]
344 var cl_pqab_series_term v2 = args.next(); // [N1+2]
345 var cl_pqab_series_term v3 = args.next(); // [N1+3]
346 var cl_I p01 = v0.p * v1.p;
347 var cl_I p012 = p01 * v2.p;
348 var cl_I p0123 = p012 * v3.p;
349 if (P) { *P = p0123; }
350 var cl_I q23 = v2.q * v3.q;
351 var cl_I q123 = v1.q * q23;
353 var cl_I b01 = v0.b * v1.b;
354 var cl_I b23 = v2.b * v3.b;
356 *T = b23 * (v1.b * q123 * v0.a * v0.p
357 + v0.b * q23 * v1.a * p01)
358 + b01 * (v3.b * v3.q * v2.a * p012
359 + v2.b * v3.a * p0123);
363 var uintC Nm = (N1+N2)/2; // midpoint
364 // Compute left part.
365 var cl_R LP, LQ, LB, LT;
366 eval_pqab_series_aux(N1,Nm,args,&LP,&LQ,&LB,<,trunclen);
367 // Compute right part.
368 var cl_R RP, RQ, RB, RT;
369 eval_pqab_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RB,&RT,trunclen);
370 // Put together partial results.
373 truncate_precision(*P,trunclen);
376 truncate_precision(*Q,trunclen);
378 truncate_precision(*B,trunclen);
379 // S = LS + LP/LQ * RS, so T = RB*RQ*LT + LB*LP*RT.
380 *T = RB*RQ*LT + LB*LP*RT;
381 truncate_precision(*T,trunclen);
387 const cl_LF eval_rational_series (uintC N, cl_pqab_series_stream& args, uintC len, uintC trunclen)
390 return cl_I_to_LF(0,len);
392 eval_pqab_series_aux(0,N,args,NULL,&Q,&B,&T,trunclen);
393 return cl_R_to_LF(T,len) / cl_R_to_LF(B*Q,len);
395 // Bit complexity (if p(n), q(n), a(n), b(n) have length O(log(n))):