1 // eval_rational_series().
7 #include "cl_LF_tran.h"
12 #include "cl_lfloat.h"
13 #include "cl_integer.h"
18 // Evaluates S = sum(N1 <= n < N2, a(n)/b(n) * (p(N1)...p(n))/(q(N1)...q(n)))
19 // and returns P = p(N1)...p(N2-1), Q = q(N1)...q(N2-1), B = B(N1)...B(N2-1)
20 // and T = B*Q*S (all integers). On entry N1 < N2.
21 // P will not be computed if a NULL pointer is passed.
23 static void eval_pq_series_aux (uintL N1, uintL N2,
24 const cl_pq_series& args,
25 cl_I* P, cl_I* Q, cl_I* T)
31 if (P) { *P = args.pv[N1]; }
36 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
38 *Q = args.qv[N1] * args.qv[N1+1];
39 *T = args.qv[N1+1] * args.pv[N1]
44 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
45 var cl_I p012 = p01 * args.pv[N1+2];
47 var cl_I q12 = args.qv[N1+1] * args.qv[N1+2];
48 *Q = args.qv[N1] * q12;
49 *T = q12 * args.pv[N1]
55 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
56 var cl_I p012 = p01 * args.pv[N1+2];
57 var cl_I p0123 = p012 * args.pv[N1+3];
58 if (P) { *P = p0123; }
59 var cl_I q23 = args.qv[N1+2] * args.qv[N1+3];
60 var cl_I q123 = args.qv[N1+1] * q23;
61 *Q = args.qv[N1] * q123;
62 *T = q123 * args.pv[N1]
64 + args.qv[N1+3] * p012
69 var uintL Nm = (N1+N2)/2; // midpoint
72 eval_pq_series_aux(N1,Nm,args,&LP,&LQ,<);
73 // Compute right part.
75 eval_pq_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RT);
76 // Put together partial results.
77 if (P) { *P = LP*RP; }
79 // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT.
86 static void eval_pqs_series_aux (uintL N1, uintL N2,
87 const cl_pq_series& args,
88 cl_I* P, cl_I* Q, uintL* QS, cl_I* T)
94 if (P) { *P = args.pv[N1]; }
100 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
102 *Q = args.qv[N1] * args.qv[N1+1];
103 *QS = args.qsv[N1] + args.qsv[N1+1];
104 *T = ((args.qv[N1+1] * args.pv[N1]) << args.qsv[N1+1])
109 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
110 var cl_I p012 = p01 * args.pv[N1+2];
111 if (P) { *P = p012; }
112 var cl_I q12 = args.qv[N1+1] * args.qv[N1+2];
113 *Q = args.qv[N1] * q12;
114 *QS = args.qsv[N1] + args.qsv[N1+1] + args.qsv[N1+2];
115 *T = ((q12 * args.pv[N1]) << (args.qsv[N1+1] + args.qsv[N1+2]))
116 + ((args.qv[N1+2] * p01) << args.qsv[N1+2])
121 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
122 var cl_I p012 = p01 * args.pv[N1+2];
123 var cl_I p0123 = p012 * args.pv[N1+3];
124 if (P) { *P = p0123; }
125 var cl_I q23 = args.qv[N1+2] * args.qv[N1+3];
126 var cl_I q123 = args.qv[N1+1] * q23;
127 *Q = args.qv[N1] * q123;
128 *QS = args.qsv[N1] + args.qsv[N1+1] + args.qsv[N1+2] + args.qsv[N1+3];
129 *T = ((((((q123 * args.pv[N1]) << args.qsv[N1+1])
130 + q23 * p01) << args.qsv[N1+2])
131 + args.qv[N1+3] * p012) << args.qsv[N1+3])
136 var uintL Nm = (N1+N2)/2; // midpoint
137 // Compute left part.
140 eval_pqs_series_aux(N1,Nm,args,&LP,&LQ,&LQS,<);
141 // Compute right part.
144 eval_pqs_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RQS,&RT);
145 // Put together partial results.
146 if (P) { *P = LP*RP; }
149 // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT.
150 *T = ((RQ*LT) << RQS) + LP*RT;
156 const cl_LF eval_rational_series (uintL N, const cl_pq_series& args, uintC len)
159 return cl_I_to_LF(0,len);
162 eval_pq_series_aux(0,N,args,NULL,&Q,&T);
163 return cl_I_to_LF(T,len) / cl_I_to_LF(Q,len);
165 // Precomputation of the shift counts:
166 // Split qv[n] into qv[n]*2^qsv[n].
168 var cl_I* qp = args.qv;
169 var uintL* qsp = args.qsv;
170 for (var uintL n = 0; n < N; n++, qp++, qsp++) {
171 // Pull out maximal power of 2 out of *qp = args.qv[n].
183 eval_pqs_series_aux(0,N,args,NULL,&Q,&QS,&T);
184 return cl_I_to_LF(T,len) / scale_float(cl_I_to_LF(Q,len),QS);
187 // Bit complexity (if p(n), q(n), a(n), b(n) have length O(log(n))):