1 // eval_rational_series<bool>().
4 #include "base/cl_sysdep.h"
7 #include "float/transcendental/cl_LF_tran.h"
12 #include "cln/lfloat.h"
13 #include "cln/integer.h"
15 #include "cln/exception.h"
16 #include "float/lfloat/cl_LF.h"
17 #include "base/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_pq_series_aux (uintC N1, uintC N2,
28 const cl_pq_series& args,
29 cl_I* P, cl_I* Q, cl_I* T)
33 throw runtime_exception(); break;
35 if (P) { *P = 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 *T = args.qv[N1+1] * args.pv[N1]
48 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
49 var cl_I p012 = p01 * args.pv[N1+2];
51 var cl_I q12 = args.qv[N1+1] * args.qv[N1+2];
52 *Q = args.qv[N1] * q12;
53 *T = q12 * args.pv[N1]
59 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
60 var cl_I p012 = p01 * args.pv[N1+2];
61 var cl_I p0123 = p012 * args.pv[N1+3];
62 if (P) { *P = p0123; }
63 var cl_I q23 = args.qv[N1+2] * args.qv[N1+3];
64 var cl_I q123 = args.qv[N1+1] * q23;
65 *Q = args.qv[N1] * q123;
66 *T = q123 * args.pv[N1]
68 + args.qv[N1+3] * p012
73 var uintC Nm = (N1+N2)/2; // midpoint
76 eval_pq_series_aux(N1,Nm,args,&LP,&LQ,<);
77 // Compute right part.
79 eval_pq_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RT);
80 // Put together partial results.
81 if (P) { *P = LP*RP; }
83 // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT.
91 const cl_LF eval_rational_series<false> (uintC N, const cl_pq_series& args, uintC len)
94 return cl_I_to_LF(0,len);
96 eval_pq_series_aux(0,N,args,NULL,&Q,&T);
97 return cl_I_to_LF(T,len) / cl_I_to_LF(Q,len);
100 static void eval_pqs_series_aux (uintC N1, uintC N2,
101 const cl_pq_series& args, const uintC* qsv,
102 cl_I* P, cl_I* Q, uintC* QS, cl_I* T)
106 throw runtime_exception(); break;
108 if (P) { *P = args.pv[N1]; }
114 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
116 *Q = args.qv[N1] * args.qv[N1+1];
117 *QS = qsv[N1] + qsv[N1+1];
118 *T = ((args.qv[N1+1] * args.pv[N1]) << qsv[N1+1])
123 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
124 var cl_I p012 = p01 * args.pv[N1+2];
125 if (P) { *P = p012; }
126 var cl_I q12 = args.qv[N1+1] * args.qv[N1+2];
127 *Q = args.qv[N1] * q12;
128 *QS = qsv[N1] + qsv[N1+1] + qsv[N1+2];
129 *T = ((q12 * args.pv[N1]) << (qsv[N1+1] + qsv[N1+2]))
130 + ((args.qv[N1+2] * p01) << qsv[N1+2])
135 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
136 var cl_I p012 = p01 * args.pv[N1+2];
137 var cl_I p0123 = p012 * args.pv[N1+3];
138 if (P) { *P = p0123; }
139 var cl_I q23 = args.qv[N1+2] * args.qv[N1+3];
140 var cl_I q123 = args.qv[N1+1] * q23;
141 *Q = args.qv[N1] * q123;
142 *QS = qsv[N1] + qsv[N1+1] + qsv[N1+2] + qsv[N1+3];
143 *T = ((((((q123 * args.pv[N1]) << qsv[N1+1])
144 + q23 * p01) << qsv[N1+2])
145 + args.qv[N1+3] * p012) << qsv[N1+3])
150 var uintC Nm = (N1+N2)/2; // midpoint
151 // Compute left part.
154 eval_pqs_series_aux(N1,Nm,args,qsv,&LP,&LQ,&LQS,<);
155 // Compute right part.
158 eval_pqs_series_aux(Nm,N2,args,qsv,(P?&RP:(cl_I*)0),&RQ,&RQS,&RT);
159 // Put together partial results.
160 if (P) { *P = LP*RP; }
163 // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT.
164 *T = ((RQ*LT) << RQS) + LP*RT;
171 const cl_LF eval_rational_series<true> (uintC N, const cl_pq_series& args, uintC len)
174 return cl_I_to_LF(0,len);
176 // Precomputation of the shift counts:
177 // Split qv[n] into qv[n]*2^qsv[n].
179 var uintC* qsv = (uintC*) cl_alloca(N*sizeof(uintC));
180 var cl_I* qp = args.qv;
181 var uintC* qsp = qsv;
182 for (var uintC n = 0; n < N; n++, qp++, qsp++) {
183 *qsp = pullout_shiftcount(*qp);
187 eval_pqs_series_aux(0,N,args,qsv,NULL,&Q,&QS,&T);
188 return cl_I_to_LF(T,len) / scale_float(cl_I_to_LF(Q,len),QS);
191 static void eval_pq_series_aux (uintC N1, uintC N2,
192 cl_pq_series_stream& args,
193 cl_I* P, cl_I* Q, cl_I* T)
197 throw runtime_exception(); break;
199 var cl_pq_series_term v0 = args.next(); // [N1]
200 if (P) { *P = v0.p; }
206 var cl_pq_series_term v0 = args.next(); // [N1]
207 var cl_pq_series_term v1 = args.next(); // [N1+1]
208 var cl_I p01 = v0.p * v1.p;
216 var cl_pq_series_term v0 = args.next(); // [N1]
217 var cl_pq_series_term v1 = args.next(); // [N1+1]
218 var cl_pq_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;
230 var cl_pq_series_term v0 = args.next(); // [N1]
231 var cl_pq_series_term v1 = args.next(); // [N1+1]
232 var cl_pq_series_term v2 = args.next(); // [N1+2]
233 var cl_pq_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;
248 var uintC Nm = (N1+N2)/2; // midpoint
249 // Compute left part.
251 eval_pq_series_aux(N1,Nm,args,&LP,&LQ,<);
252 // Compute right part.
254 eval_pq_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.
266 const cl_LF eval_rational_series<false> (uintC N, cl_pq_series_stream& args, uintC len)
269 return cl_I_to_LF(0,len);
271 eval_pq_series_aux(0,N,args,NULL,&Q,&T);
272 return cl_I_to_LF(T,len) / cl_I_to_LF(Q,len);
275 static void eval_pqs_series_aux (uintC N1, uintC N2,
276 cl_pq_series_stream& args,
277 cl_I* P, cl_I* Q, uintC* QS, cl_I* T)
281 throw runtime_exception(); break;
283 var cl_pq_series_term v0 = args.next(); // [N1]
284 var uintC qs0 = pullout_shiftcount(v0.q);
285 if (P) { *P = v0.p; }
292 var cl_pq_series_term v0 = args.next(); // [N1]
293 var cl_pq_series_term v1 = args.next(); // [N1+1]
294 var uintC qs0 = pullout_shiftcount(v0.q);
295 var uintC qs1 = pullout_shiftcount(v1.q);
296 var cl_I p01 = v0.p * v1.p;
300 *T = ((v1.q * v0.p) << qs1)
305 var cl_pq_series_term v0 = args.next(); // [N1]
306 var cl_pq_series_term v1 = args.next(); // [N1+1]
307 var cl_pq_series_term v2 = args.next(); // [N1+2]
308 var uintC qs0 = pullout_shiftcount(v0.q);
309 var uintC qs1 = pullout_shiftcount(v1.q);
310 var uintC qs2 = pullout_shiftcount(v2.q);
311 var cl_I p01 = v0.p * v1.p;
312 var cl_I p012 = p01 * v2.p;
313 if (P) { *P = p012; }
314 var cl_I q12 = v1.q * v2.q;
316 *QS = qs0 + qs1 + qs2;
317 *T = ((q12 * v0.p) << (qs1 + qs2))
318 + ((v2.q * p01) << qs2)
323 var cl_pq_series_term v0 = args.next(); // [N1]
324 var cl_pq_series_term v1 = args.next(); // [N1+1]
325 var cl_pq_series_term v2 = args.next(); // [N1+2]
326 var cl_pq_series_term v3 = args.next(); // [N1+3]
327 var uintC qs0 = pullout_shiftcount(v0.q);
328 var uintC qs1 = pullout_shiftcount(v1.q);
329 var uintC qs2 = pullout_shiftcount(v2.q);
330 var uintC qs3 = pullout_shiftcount(v3.q);
331 var cl_I p01 = v0.p * v1.p;
332 var cl_I p012 = p01 * v2.p;
333 var cl_I p0123 = p012 * v3.p;
334 if (P) { *P = p0123; }
335 var cl_I q23 = v2.q * v3.q;
336 var cl_I q123 = v1.q * q23;
338 *QS = qs0 + qs1 + qs2 + qs3;
339 *T = ((((((q123 * v0.p) << qs1)
341 + v3.q * p012) << qs3)
346 var uintC Nm = (N1+N2)/2; // midpoint
347 // Compute left part.
350 eval_pqs_series_aux(N1,Nm,args,&LP,&LQ,&LQS,<);
351 // Compute right part.
354 eval_pqs_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RQS,&RT);
355 // Put together partial results.
356 if (P) { *P = LP*RP; }
359 // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT.
360 *T = ((RQ*LT) << RQS) + LP*RT;
367 const cl_LF eval_rational_series<true> (uintC N, cl_pq_series_stream& args, uintC len)
370 return cl_I_to_LF(0,len);
373 eval_pqs_series_aux(0,N,args,NULL,&Q,&QS,&T);
374 return cl_I_to_LF(T,len) / scale_float(cl_I_to_LF(Q,len),QS);
377 static void eval_pq_series_aux (uintC N1, uintC N2,
378 cl_pq_series_stream& args,
379 cl_R* P, cl_R* Q, cl_R* T,
384 throw runtime_exception(); break;
386 var cl_pq_series_term v0 = args.next(); // [N1]
387 if (P) { *P = v0.p; }
393 var cl_pq_series_term v0 = args.next(); // [N1]
394 var cl_pq_series_term v1 = args.next(); // [N1+1]
395 var cl_I p01 = v0.p * v1.p;
403 var cl_pq_series_term v0 = args.next(); // [N1]
404 var cl_pq_series_term v1 = args.next(); // [N1+1]
405 var cl_pq_series_term v2 = args.next(); // [N1+2]
406 var cl_I p01 = v0.p * v1.p;
407 var cl_I p012 = p01 * v2.p;
408 if (P) { *P = p012; }
409 var cl_I q12 = v1.q * v2.q;
417 var cl_pq_series_term v0 = args.next(); // [N1]
418 var cl_pq_series_term v1 = args.next(); // [N1+1]
419 var cl_pq_series_term v2 = args.next(); // [N1+2]
420 var cl_pq_series_term v3 = args.next(); // [N1+3]
421 var cl_I p01 = v0.p * v1.p;
422 var cl_I p012 = p01 * v2.p;
423 var cl_I p0123 = p012 * v3.p;
424 if (P) { *P = p0123; }
425 var cl_I q23 = v2.q * v3.q;
426 var cl_I q123 = v1.q * q23;
435 var uintC Nm = (N1+N2)/2; // midpoint
436 // Compute left part.
438 eval_pq_series_aux(N1,Nm,args,&LP,&LQ,<,trunclen);
439 // Compute right part.
441 eval_pq_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RT,trunclen);
442 // Put together partial results.
445 truncate_precision(*P,trunclen);
448 truncate_precision(*Q,trunclen);
449 // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT.
451 truncate_precision(*T,trunclen);
458 const cl_LF eval_rational_series<false> (uintC N, cl_pq_series_stream& args, uintC len, uintC trunclen)
461 return cl_I_to_LF(0,len);
463 eval_pq_series_aux(0,N,args,NULL,&Q,&T,trunclen);
464 return cl_R_to_LF(T,len) / cl_R_to_LF(Q,len);
466 // Bit complexity (if p(n), q(n), a(n), b(n) have length O(log(n))):