1 // eval_rational_series<bool>().
2
3 // General includes.
4 #include "base/cl_sysdep.h"
5
6 // Specification.
7 #include "float/transcendental/cl_LF_tran.h"
8
9
10 // Implementation.
11
12 #include "cln/lfloat.h"
13 #include "cln/integer.h"
14 #include "cln/real.h"
15 #include "cln/exception.h"
16 #include "float/lfloat/cl_LF.h"
17 #include "base/cl_alloca.h"
18
19 namespace cln {
20
21 // Subroutine.
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.
26
eval_pqb_series_aux(uintC N1,uintC N2,const cl_pqb_series & args,cl_I * P,cl_I * Q,cl_I * B,cl_I * T)27 static void eval_pqb_series_aux (uintC N1, uintC N2,
28 const cl_pqb_series& args,
29 cl_I* P, cl_I* Q, cl_I* B, cl_I* T)
30 {
31 switch (N2 - N1) {
32 case 0:
33 throw runtime_exception(); break;
34 case 1:
35 if (P) { *P = args.pv[N1]; }
36 *Q = args.qv[N1];
37 *B = args.bv[N1];
38 *T = args.pv[N1];
39 break;
40 case 2: {
41 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
42 if (P) { *P = p01; }
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.pv[N1]
46 + args.bv[N1] * p01;
47 break;
48 }
49 case 3: {
50 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
51 var cl_I p012 = p01 * args.pv[N1+2];
52 if (P) { *P = p012; }
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.pv[N1]
58 + args.bv[N1] * (args.bv[N1+2] * args.qv[N1+2] * p01
59 + args.bv[N1+1] * p012);
60 break;
61 }
62 case 4: {
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];
72 *B = b01 * b23;
73 *T = b23 * (args.bv[N1+1] * q123 * args.pv[N1]
74 + args.bv[N1] * q23 * p01)
75 + b01 * (args.bv[N1+3] * args.qv[N1+3] * p012
76 + args.bv[N1+2] * p0123);
77 break;
78 }
79 default: {
80 var uintC Nm = (N1+N2)/2; // midpoint
81 // Compute left part.
82 var cl_I LP, LQ, LB, LT;
83 eval_pqb_series_aux(N1,Nm,args,&LP,&LQ,&LB,<);
84 // Compute right part.
85 var cl_I RP, RQ, RB, RT;
86 eval_pqb_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RB,&RT);
87 // Put together partial results.
88 if (P) { *P = LP*RP; }
89 *Q = LQ*RQ;
90 *B = LB*RB;
91 // S = LS + LP/LQ * RS, so T = RB*RQ*LT + LB*LP*RT.
92 *T = RB*RQ*LT + LB*LP*RT;
93 break;
94 }
95 }
96 }
97
98 template<>
eval_rational_series(uintC N,const cl_pqb_series & args,uintC len)99 const cl_LF eval_rational_series<false> (uintC N, const cl_pqb_series& args, uintC len)
100 {
101 if (N==0)
102 return cl_I_to_LF(0,len);
103 var cl_I Q, B, T;
104 eval_pqb_series_aux(0,N,args,NULL,&Q,&B,&T);
105 return cl_I_to_LF(T,len) / cl_I_to_LF(B*Q,len);
106 }
107
eval_pqsb_series_aux(uintC N1,uintC N2,const cl_pqb_series & args,const uintC * qsv,cl_I * P,cl_I * Q,uintC * QS,cl_I * B,cl_I * T)108 static void eval_pqsb_series_aux (uintC N1, uintC N2,
109 const cl_pqb_series& args, const uintC* qsv,
110 cl_I* P, cl_I* Q, uintC* QS, cl_I* B, cl_I* T)
111 {
112 switch (N2 - N1) {
113 case 0:
114 throw runtime_exception(); break;
115 case 1:
116 if (P) { *P = args.pv[N1]; }
117 *Q = args.qv[N1];
118 *QS = qsv[N1];
119 *B = args.bv[N1];
120 *T = args.pv[N1];
121 break;
122 case 2: {
123 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
124 if (P) { *P = p01; }
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.pv[N1]) << qsv[N1+1])
129 + args.bv[N1] * p01;
130 break;
131 }
132 case 3: {
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.pv[N1]) << (qsv[N1+1] + qsv[N1+2]))
142 + args.bv[N1] * (((args.bv[N1+2] * args.qv[N1+2] * p01) << qsv[N1+2])
143 + args.bv[N1+1] * p012);
144 break;
145 }
146 case 4: {
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];
157 *B = b01 * b23;
158 *T = ((b23 * (((args.bv[N1+1] * q123 * args.pv[N1]) << qsv[N1+1])
159 + args.bv[N1] * q23 * p01)) << (qsv[N1+2] + qsv[N1+3]))
160 + b01 * (((args.bv[N1+3] * args.qv[N1+3] * p012) << qsv[N1+3])
161 + args.bv[N1+2] * p0123);
162 break;
163 }
164 default: {
165 var uintC Nm = (N1+N2)/2; // midpoint
166 // Compute left part.
167 var cl_I LP, LQ, LB, LT;
168 var uintC LQS;
169 eval_pqsb_series_aux(N1,Nm,args,qsv,&LP,&LQ,&LQS,&LB,<);
170 // Compute right part.
171 var cl_I RP, RQ, RB, RT;
172 var uintC RQS;
173 eval_pqsb_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; }
176 *Q = LQ*RQ;
177 *QS = LQS+RQS;
178 *B = LB*RB;
179 // S = LS + LP/LQ * RS, so T = RB*RQ*LT + LB*LP*RT.
180 *T = ((RB*RQ*LT) << RQS) + LB*LP*RT;
181 break;
182 }
183 }
184 }
185
186 template<>
eval_rational_series(uintC N,const cl_pqb_series & args,uintC len)187 const cl_LF eval_rational_series<true> (uintC N, const cl_pqb_series& args, uintC len)
188 {
189 if (N==0)
190 return cl_I_to_LF(0,len);
191 var cl_I Q, B, T;
192 // Precomputation of the shift counts:
193 // Split qv[n] into qv[n]*2^qsv[n].
194 CL_ALLOCA_STACK;
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);
200 }
201 // Main computation.
202 var uintC QS;
203 eval_pqsb_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);
205 }
206
eval_pqb_series_aux(uintC N1,uintC N2,cl_pqb_series_stream & args,cl_I * P,cl_I * Q,cl_I * B,cl_I * T)207 static void eval_pqb_series_aux (uintC N1, uintC N2,
208 cl_pqb_series_stream& args,
209 cl_I* P, cl_I* Q, cl_I* B, cl_I* T)
210 {
211 switch (N2 - N1) {
212 case 0:
213 throw runtime_exception(); break;
214 case 1: {
215 var cl_pqb_series_term v0 = args.next(); // [N1]
216 if (P) { *P = v0.p; }
217 *Q = v0.q;
218 *B = v0.b;
219 *T = v0.p;
220 break;
221 }
222 case 2: {
223 var cl_pqb_series_term v0 = args.next(); // [N1]
224 var cl_pqb_series_term v1 = args.next(); // [N1+1]
225 var cl_I p01 = v0.p * v1.p;
226 if (P) { *P = p01; }
227 *Q = v0.q * v1.q;
228 *B = v0.b * v1.b;
229 *T = v1.b * v1.q * v0.p
230 + v0.b * p01;
231 break;
232 }
233 case 3: {
234 var cl_pqb_series_term v0 = args.next(); // [N1]
235 var cl_pqb_series_term v1 = args.next(); // [N1+1]
236 var cl_pqb_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;
241 *Q = v0.q * q12;
242 var cl_I b12 = v1.b * v2.b;
243 *B = v0.b * b12;
244 *T = b12 * q12 * v0.p
245 + v0.b * (v2.b * v2.q * p01
246 + v1.b * p012);
247 break;
248 }
249 case 4: {
250 var cl_pqb_series_term v0 = args.next(); // [N1]
251 var cl_pqb_series_term v1 = args.next(); // [N1+1]
252 var cl_pqb_series_term v2 = args.next(); // [N1+2]
253 var cl_pqb_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;
260 *Q = v0.q * q123;
261 var cl_I b01 = v0.b * v1.b;
262 var cl_I b23 = v2.b * v3.b;
263 *B = b01 * b23;
264 *T = b23 * (v1.b * q123 * v0.p
265 + v0.b * q23 * p01)
266 + b01 * (v3.b * v3.q * p012
267 + v2.b * p0123);
268 break;
269 }
270 default: {
271 var uintC Nm = (N1+N2)/2; // midpoint
272 // Compute left part.
273 var cl_I LP, LQ, LB, LT;
274 eval_pqb_series_aux(N1,Nm,args,&LP,&LQ,&LB,<);
275 // Compute right part.
276 var cl_I RP, RQ, RB, RT;
277 eval_pqb_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RB,&RT);
278 // Put together partial results.
279 if (P) { *P = LP*RP; }
280 *Q = LQ*RQ;
281 *B = LB*RB;
282 // S = LS + LP/LQ * RS, so T = RB*RQ*LT + LB*LP*RT.
283 *T = RB*RQ*LT + LB*LP*RT;
284 break;
285 }
286 }
287 }
288
289 template<>
eval_rational_series(uintC N,cl_pqb_series_stream & args,uintC len)290 const cl_LF eval_rational_series<false> (uintC N, cl_pqb_series_stream& args, uintC len)
291 {
292 if (N==0)
293 return cl_I_to_LF(0,len);
294 var cl_I Q, B, T;
295 eval_pqb_series_aux(0,N,args,NULL,&Q,&B,&T);
296 return cl_I_to_LF(T,len) / cl_I_to_LF(B*Q,len);
297 }
298
eval_pqb_series_aux(uintC N1,uintC N2,cl_pqb_series_stream & args,cl_R * P,cl_R * Q,cl_R * B,cl_R * T,uintC trunclen)299 static void eval_pqb_series_aux (uintC N1, uintC N2,
300 cl_pqb_series_stream& args,
301 cl_R* P, cl_R* Q, cl_R* B, cl_R* T,
302 uintC trunclen)
303 {
304 switch (N2 - N1) {
305 case 0:
306 throw runtime_exception(); break;
307 case 1: {
308 var cl_pqb_series_term v0 = args.next(); // [N1]
309 if (P) { *P = v0.p; }
310 *Q = v0.q;
311 *B = v0.b;
312 *T = v0.p;
313 break;
314 }
315 case 2: {
316 var cl_pqb_series_term v0 = args.next(); // [N1]
317 var cl_pqb_series_term v1 = args.next(); // [N1+1]
318 var cl_I p01 = v0.p * v1.p;
319 if (P) { *P = p01; }
320 *Q = v0.q * v1.q;
321 *B = v0.b * v1.b;
322 *T = v1.b * v1.q * v0.p
323 + v0.b * p01;
324 break;
325 }
326 case 3: {
327 var cl_pqb_series_term v0 = args.next(); // [N1]
328 var cl_pqb_series_term v1 = args.next(); // [N1+1]
329 var cl_pqb_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;
334 *Q = v0.q * q12;
335 var cl_I b12 = v1.b * v2.b;
336 *B = v0.b * b12;
337 *T = b12 * q12 * v0.p
338 + v0.b * (v2.b * v2.q * p01
339 + v1.b * p012);
340 break;
341 }
342 case 4: {
343 var cl_pqb_series_term v0 = args.next(); // [N1]
344 var cl_pqb_series_term v1 = args.next(); // [N1+1]
345 var cl_pqb_series_term v2 = args.next(); // [N1+2]
346 var cl_pqb_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;
353 *Q = v0.q * q123;
354 var cl_I b01 = v0.b * v1.b;
355 var cl_I b23 = v2.b * v3.b;
356 *B = b01 * b23;
357 *T = b23 * (v1.b * q123 * v0.p
358 + v0.b * q23 * p01)
359 + b01 * (v3.b * v3.q * p012
360 + v2.b * p0123);
361 break;
362 }
363 default: {
364 var uintC Nm = (N1+N2)/2; // midpoint
365 // Compute left part.
366 var cl_R LP, LQ, LB, LT;
367 eval_pqb_series_aux(N1,Nm,args,&LP,&LQ,&LB,<,trunclen);
368 // Compute right part.
369 var cl_R RP, RQ, RB, RT;
370 eval_pqb_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RB,&RT,trunclen);
371 // Put together partial results.
372 if (P) {
373 *P = LP*RP;
374 truncate_precision(*P,trunclen);
375 }
376 *Q = LQ*RQ;
377 truncate_precision(*Q,trunclen);
378 *B = LB*RB;
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);
383 break;
384 }
385 }
386 }
387
388 template<>
eval_rational_series(uintC N,cl_pqb_series_stream & args,uintC len,uintC trunclen)389 const cl_LF eval_rational_series<false> (uintC N, cl_pqb_series_stream& args, uintC len, uintC trunclen)
390 {
391 if (N==0)
392 return cl_I_to_LF(0,len);
393 var cl_R Q, B, T;
394 eval_pqb_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);
396 }
397
398 // Bit complexity (if p(n), q(n), a(n), b(n) have length O(log(n))):
399 // O(log(N)^2*M(N)).
400
401 } // namespace cln
402