1 // eval_rational_series<bool>().
7 #include "cl_LF_tran.h"
12 #include "cln/lfloat.h"
13 #include "cln/integer.h"
15 #include "cln/exception.h"
17 #include "cl_alloca.h"
22 // Evaluates S = sum(N1 <= n < N2, a(n)/b(n) * (p(N1)...p(n))/(q(N1)...q(n)))
23 // and returns P = p(N1)...p(N2-1), Q = q(N1)...q(N2-1), B = B(N1)...B(N2-1)
24 // and T = B*Q*S (all integers). On entry N1 < N2.
25 // P will not be computed if a NULL pointer is passed.
27 static void eval_pqab_series_aux (uintC N1, uintC N2,
28 const cl_pqab_series& args,
29 cl_I* P, cl_I* Q, cl_I* B, cl_I* T)
33 throw runtime_exception(); break;
35 if (P) { *P = args.pv[N1]; }
38 *T = args.av[N1] * args.pv[N1];
41 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
43 *Q = args.qv[N1] * args.qv[N1+1];
44 *B = args.bv[N1] * args.bv[N1+1];
45 *T = args.bv[N1+1] * args.qv[N1+1] * args.av[N1] * args.pv[N1]
46 + args.bv[N1] * args.av[N1+1] * p01;
50 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
51 var cl_I p012 = p01 * args.pv[N1+2];
53 var cl_I q12 = args.qv[N1+1] * args.qv[N1+2];
54 *Q = args.qv[N1] * q12;
55 var cl_I b12 = args.bv[N1+1] * args.bv[N1+2];
56 *B = args.bv[N1] * b12;
57 *T = b12 * q12 * args.av[N1] * args.pv[N1]
58 + args.bv[N1] * (args.bv[N1+2] * args.qv[N1+2] * args.av[N1+1] * p01
59 + args.bv[N1+1] * args.av[N1+2] * p012);
63 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
64 var cl_I p012 = p01 * args.pv[N1+2];
65 var cl_I p0123 = p012 * args.pv[N1+3];
66 if (P) { *P = p0123; }
67 var cl_I q23 = args.qv[N1+2] * args.qv[N1+3];
68 var cl_I q123 = args.qv[N1+1] * q23;
69 *Q = args.qv[N1] * q123;
70 var cl_I b01 = args.bv[N1] * args.bv[N1+1];
71 var cl_I b23 = args.bv[N1+2] * args.bv[N1+3];
73 *T = b23 * (args.bv[N1+1] * q123 * args.av[N1] * args.pv[N1]
74 + args.bv[N1] * q23 * args.av[N1+1] * p01)
75 + b01 * (args.bv[N1+3] * args.qv[N1+3] * args.av[N1+2] * p012
76 + args.bv[N1+2] * args.av[N1+3] * p0123);
80 var uintC Nm = (N1+N2)/2; // midpoint
82 var cl_I LP, LQ, LB, LT;
83 eval_pqab_series_aux(N1,Nm,args,&LP,&LQ,&LB,<);
84 // Compute right part.
85 var cl_I RP, RQ, RB, RT;
86 eval_pqab_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RB,&RT);
87 // Put together partial results.
88 if (P) { *P = LP*RP; }
91 // S = LS + LP/LQ * RS, so T = RB*RQ*LT + LB*LP*RT.
92 *T = RB*RQ*LT + LB*LP*RT;
99 const cl_LF eval_rational_series<false> (uintC N, const cl_pqab_series& args, uintC len)
102 return cl_I_to_LF(0,len);
104 eval_pqab_series_aux(0,N,args,NULL,&Q,&B,&T);
105 return cl_I_to_LF(T,len) / cl_I_to_LF(B*Q,len);
108 static void eval_pqsab_series_aux (uintC N1, uintC N2,
109 const cl_pqab_series& args, const uintC* qsv,
110 cl_I* P, cl_I* Q, uintC* QS, cl_I* B, cl_I* T)
114 throw runtime_exception(); break;
116 if (P) { *P = args.pv[N1]; }
120 *T = args.av[N1] * args.pv[N1];
123 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
125 *Q = args.qv[N1] * args.qv[N1+1];
126 *QS = qsv[N1] + qsv[N1+1];
127 *B = args.bv[N1] * args.bv[N1+1];
128 *T = ((args.bv[N1+1] * args.qv[N1+1] * args.av[N1] * args.pv[N1]) << qsv[N1+1])
129 + args.bv[N1] * args.av[N1+1] * p01;
133 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
134 var cl_I p012 = p01 * args.pv[N1+2];
135 if (P) { *P = p012; }
136 var cl_I q12 = args.qv[N1+1] * args.qv[N1+2];
137 *Q = args.qv[N1] * q12;
138 *QS = qsv[N1] + qsv[N1+1] + qsv[N1+2];
139 var cl_I b12 = args.bv[N1+1] * args.bv[N1+2];
140 *B = args.bv[N1] * b12;
141 *T = ((b12 * q12 * args.av[N1] * args.pv[N1]) << (qsv[N1+1] + qsv[N1+2]))
142 + args.bv[N1] * (((args.bv[N1+2] * args.qv[N1+2] * args.av[N1+1] * p01) << qsv[N1+2])
143 + args.bv[N1+1] * args.av[N1+2] * p012);
147 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
148 var cl_I p012 = p01 * args.pv[N1+2];
149 var cl_I p0123 = p012 * args.pv[N1+3];
150 if (P) { *P = p0123; }
151 var cl_I q23 = args.qv[N1+2] * args.qv[N1+3];
152 var cl_I q123 = args.qv[N1+1] * q23;
153 *Q = args.qv[N1] * q123;
154 *QS = qsv[N1] + qsv[N1+1] + qsv[N1+2] + qsv[N1+3];
155 var cl_I b01 = args.bv[N1] * args.bv[N1+1];
156 var cl_I b23 = args.bv[N1+2] * args.bv[N1+3];
158 *T = ((b23 * (((args.bv[N1+1] * q123 * args.av[N1] * args.pv[N1]) << qsv[N1+1])
159 + args.bv[N1] * q23 * args.av[N1+1] * p01)) << (qsv[N1+2] + qsv[N1+3]))
160 + b01 * (((args.bv[N1+3] * args.qv[N1+3] * args.av[N1+2] * p012) << qsv[N1+3])
161 + args.bv[N1+2] * args.av[N1+3] * p0123);
165 var uintC Nm = (N1+N2)/2; // midpoint
166 // Compute left part.
167 var cl_I LP, LQ, LB, LT;
169 eval_pqsab_series_aux(N1,Nm,args,qsv,&LP,&LQ,&LQS,&LB,<);
170 // Compute right part.
171 var cl_I RP, RQ, RB, RT;
173 eval_pqsab_series_aux(Nm,N2,args,qsv,(P?&RP:(cl_I*)0),&RQ,&RQS,&RB,&RT);
174 // Put together partial results.
175 if (P) { *P = LP*RP; }
179 // S = LS + LP/LQ * RS, so T = RB*RQ*LT + LB*LP*RT.
180 *T = ((RB*RQ*LT) << RQS) + LB*LP*RT;
187 const cl_LF eval_rational_series<true> (uintC N, const cl_pqab_series& args, uintC len)
190 return cl_I_to_LF(0,len);
192 // Precomputation of the shift counts:
193 // Split qv[n] into qv[n]*2^qsv[n].
195 var uintC* qsv = (uintC*) cl_alloca(N*sizeof(uintC));
196 var cl_I* qp = args.qv;
197 var uintC* qsp = qsv;
198 for (var uintC n = 0; n < N; n++, qp++, qsp++) {
199 *qsp = pullout_shiftcount(*qp);
203 eval_pqsab_series_aux(0,N,args,qsv,NULL,&Q,&QS,&B,&T);
204 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;
290 const cl_LF eval_rational_series<false> (uintC N, cl_pqab_series_stream& args, uintC len)
293 return cl_I_to_LF(0,len);
295 eval_pqab_series_aux(0,N,args,NULL,&Q,&B,&T);
296 return cl_I_to_LF(T,len) / cl_I_to_LF(B*Q,len);
299 static void eval_pqab_series_aux (uintC N1, uintC N2,
300 cl_pqab_series_stream& args,
301 cl_R* P, cl_R* Q, cl_R* B, cl_R* T,
306 throw runtime_exception(); break;
308 var cl_pqab_series_term v0 = args.next(); // [N1]
309 if (P) { *P = v0.p; }
316 var cl_pqab_series_term v0 = args.next(); // [N1]
317 var cl_pqab_series_term v1 = args.next(); // [N1+1]
318 var cl_I p01 = v0.p * v1.p;
322 *T = v1.b * v1.q * v0.a * v0.p
327 var cl_pqab_series_term v0 = args.next(); // [N1]
328 var cl_pqab_series_term v1 = args.next(); // [N1+1]
329 var cl_pqab_series_term v2 = args.next(); // [N1+2]
330 var cl_I p01 = v0.p * v1.p;
331 var cl_I p012 = p01 * v2.p;
332 if (P) { *P = p012; }
333 var cl_I q12 = v1.q * v2.q;
335 var cl_I b12 = v1.b * v2.b;
337 *T = b12 * q12 * v0.a * v0.p
338 + v0.b * (v2.b * v2.q * v1.a * p01
339 + v1.b * v2.a * p012);
343 var cl_pqab_series_term v0 = args.next(); // [N1]
344 var cl_pqab_series_term v1 = args.next(); // [N1+1]
345 var cl_pqab_series_term v2 = args.next(); // [N1+2]
346 var cl_pqab_series_term v3 = args.next(); // [N1+3]
347 var cl_I p01 = v0.p * v1.p;
348 var cl_I p012 = p01 * v2.p;
349 var cl_I p0123 = p012 * v3.p;
350 if (P) { *P = p0123; }
351 var cl_I q23 = v2.q * v3.q;
352 var cl_I q123 = v1.q * q23;
354 var cl_I b01 = v0.b * v1.b;
355 var cl_I b23 = v2.b * v3.b;
357 *T = b23 * (v1.b * q123 * v0.a * v0.p
358 + v0.b * q23 * v1.a * p01)
359 + b01 * (v3.b * v3.q * v2.a * p012
360 + v2.b * v3.a * p0123);
364 var uintC Nm = (N1+N2)/2; // midpoint
365 // Compute left part.
366 var cl_R LP, LQ, LB, LT;
367 eval_pqab_series_aux(N1,Nm,args,&LP,&LQ,&LB,<,trunclen);
368 // Compute right part.
369 var cl_R RP, RQ, RB, RT;
370 eval_pqab_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RB,&RT,trunclen);
371 // Put together partial results.
374 truncate_precision(*P,trunclen);
377 truncate_precision(*Q,trunclen);
379 truncate_precision(*B,trunclen);
380 // S = LS + LP/LQ * RS, so T = RB*RQ*LT + LB*LP*RT.
381 *T = RB*RQ*LT + LB*LP*RT;
382 truncate_precision(*T,trunclen);
389 const cl_LF eval_rational_series<false> (uintC N, cl_pqab_series_stream& args, uintC len, uintC trunclen)
392 return cl_I_to_LF(0,len);
394 eval_pqab_series_aux(0,N,args,NULL,&Q,&B,&T,trunclen);
395 return cl_R_to_LF(T,len) / cl_R_to_LF(B*Q,len);
397 // Bit complexity (if p(n), q(n), a(n), b(n) have length O(log(n))):