1 // eval_pqcd_series_aux().
7 #include "cl_LF_tran.h"
12 #include "cln/integer.h"
14 #include "cln/exception.h"
18 void eval_pqcd_series_aux (uintC N, cl_pqcd_series_term* args, cl_pqcd_series_result<cl_I>& Z, bool rightmost)
23 throw runtime_exception(); break;
25 if (!rightmost) { Z.P = args[0].p; }
28 if (!rightmost) { Z.C = args[0].c; }
30 Z.V = args[0].c * args[0].p;
33 var cl_I p01 = args[0].p * args[1].p;
34 if (!rightmost) { Z.P = p01; }
35 Z.Q = args[0].q * args[1].q;
36 var cl_I p0q1 = args[0].p * args[1].q + p01;
38 var cl_I c0d1 = args[0].c * args[1].d;
39 var cl_I c1d0 = args[1].c * args[0].d;
40 if (!rightmost) { Z.C = c0d1 + c1d0; }
41 Z.D = args[0].d * args[1].d;
42 Z.V = c0d1 * p0q1 + c1d0 * p01;
46 var cl_I p01 = args[0].p * args[1].p;
47 var cl_I p012 = p01 * args[2].p;
48 if (!rightmost) { Z.P = p012; }
49 Z.Q = args[0].q * args[1].q * args[2].q;
50 var cl_I p0q1 = args[0].p * args[1].q + p01;
51 Z.T = args[2].q * p0q1 + p012;
52 var cl_I c0d1 = args[0].c * args[1].d;
53 var cl_I c1d0 = args[1].c * args[0].d;
54 var cl_I d01 = args[0].d * args[1].d;
55 if (!rightmost) { Z.C = (c0d1 + c1d0) * args[2].d + args[2].c * d01; }
56 Z.D = d01 * args[2].d;
57 Z.V = args[2].d * (args[2].q * (c0d1 * p0q1 + c1d0 * p01) + (c0d1 + c1d0) * p012) + args[2].c * d01 * p012;
61 var uintC Nm = N/2; // midpoint
63 var cl_pqcd_series_result<cl_I> L;
64 eval_pqcd_series_aux(Nm,args+0,L,false);
65 // Compute right part.
66 var cl_pqcd_series_result<cl_I> R;
67 eval_pqcd_series_aux(N-Nm,args+Nm,R,rightmost);
68 // Put together partial results.
69 if (!rightmost) { Z.P = L.P * R.P; }
71 // Z.S = L.S + L.P/L.Q*R.S;
72 var cl_I tmp = L.P * R.T;
73 Z.T = R.Q * L.T + tmp;
74 if (!rightmost) { Z.C = L.C * R.D + L.D * R.C; }
76 // Z.U = L.U + L.C/L.D * L.P/L.Q * R.S + L.P/L.Q * R.U;
77 // Z.V = R.D * R.Q * L.V + R.D * L.C * L.P * R.T + L.D * L.P * R.V;
78 Z.V = R.D * (R.Q * L.V + L.C * tmp) + L.D * L.P * R.V;
84 void eval_pqcd_series_aux (uintC N, cl_pqcd_series_stream& args, cl_pqcd_series_result<cl_I>& Z, bool rightmost)
89 throw runtime_exception(); break;
91 var cl_pqcd_series_term v0 = args.next(); // [N1]
92 if (!rightmost) { Z.P = v0.p; }
95 if (!rightmost) { Z.C = v0.c; }
101 var cl_pqcd_series_term v0 = args.next(); // [N1]
102 var cl_pqcd_series_term v1 = args.next(); // [N1+1]
103 var cl_I p01 = v0.p * v1.p;
104 if (!rightmost) { Z.P = p01; }
106 var cl_I p0q1 = v0.p * v1.q + p01;
108 var cl_I c0d1 = v0.c * v1.d;
109 var cl_I c1d0 = v1.c * v0.d;
110 if (!rightmost) { Z.C = c0d1 + c1d0; }
112 Z.V = c0d1 * p0q1 + c1d0 * p01;
116 var cl_pqcd_series_term v0 = args.next(); // [N1]
117 var cl_pqcd_series_term v1 = args.next(); // [N1+1]
118 var cl_pqcd_series_term v2 = args.next(); // [N1+2]
119 var cl_I p01 = v0.p * v1.p;
120 var cl_I p012 = p01 * v2.p;
121 if (!rightmost) { Z.P = p012; }
122 Z.Q = v0.q * v1.q * v2.q;
123 var cl_I p0q1 = v0.p * v1.q + p01;
124 Z.T = v2.q * p0q1 + p012;
125 var cl_I c0d1 = v0.c * v1.d;
126 var cl_I c1d0 = v1.c * v0.d;
127 var cl_I d01 = v0.d * v1.d;
128 if (!rightmost) { Z.C = (c0d1 + c1d0) * v2.d + v2.c * d01; }
130 Z.V = v2.d * (v2.q * (c0d1 * p0q1 + c1d0 * p01) + (c0d1 + c1d0) * p012) + v2.c * d01 * p012;
134 var uintC Nm = N/2; // midpoint
135 // Compute left part.
136 var cl_pqcd_series_result<cl_I> L;
137 eval_pqcd_series_aux(Nm,args,L,false);
138 // Compute right part.
139 var cl_pqcd_series_result<cl_I> R;
140 eval_pqcd_series_aux(N-Nm,args,R,rightmost);
141 // Put together partial results.
142 if (!rightmost) { Z.P = L.P * R.P; }
144 // Z.S = L.S + L.P/L.Q*R.S;
145 var cl_I tmp = L.P * R.T;
146 Z.T = R.Q * L.T + tmp;
147 if (!rightmost) { Z.C = L.C * R.D + L.D * R.C; }
149 // Z.U = L.U + L.C/L.D * L.P/L.Q * R.S + L.P/L.Q * R.U;
150 // Z.V = R.D * R.Q * L.V + R.D * L.C * L.P * R.T + L.D * L.P * R.V;
151 Z.V = R.D * (R.Q * L.V + L.C * tmp) + L.D * L.P * R.V;
157 void eval_pqcd_series_aux (uintC N, cl_pqcd_series_stream& args, cl_pqcd_series_result<cl_R>& Z, uintC trunclen, bool rightmost)
162 throw runtime_exception(); break;
164 var cl_pqcd_series_term v0 = args.next(); // [N1]
165 if (!rightmost) { Z.P = v0.p; }
168 if (!rightmost) { Z.C = v0.c; }
174 var cl_pqcd_series_term v0 = args.next(); // [N1]
175 var cl_pqcd_series_term v1 = args.next(); // [N1+1]
176 var cl_I p01 = v0.p * v1.p;
177 if (!rightmost) { Z.P = p01; }
179 var cl_I p0q1 = v0.p * v1.q + p01;
181 var cl_I c0d1 = v0.c * v1.d;
182 var cl_I c1d0 = v1.c * v0.d;
183 if (!rightmost) { Z.C = c0d1 + c1d0; }
185 Z.V = c0d1 * p0q1 + c1d0 * p01;
189 var cl_pqcd_series_term v0 = args.next(); // [N1]
190 var cl_pqcd_series_term v1 = args.next(); // [N1+1]
191 var cl_pqcd_series_term v2 = args.next(); // [N1+2]
192 var cl_I p01 = v0.p * v1.p;
193 var cl_I p012 = p01 * v2.p;
194 if (!rightmost) { Z.P = p012; }
195 Z.Q = v0.q * v1.q * v2.q;
196 var cl_I p0q1 = v0.p * v1.q + p01;
197 Z.T = v2.q * p0q1 + p012;
198 var cl_I c0d1 = v0.c * v1.d;
199 var cl_I c1d0 = v1.c * v0.d;
200 var cl_I d01 = v0.d * v1.d;
201 if (!rightmost) { Z.C = (c0d1 + c1d0) * v2.d + v2.c * d01; }
203 Z.V = v2.d * (v2.q * (c0d1 * p0q1 + c1d0 * p01) + (c0d1 + c1d0) * p012) + v2.c * d01 * p012;
207 var uintC Nm = N/2; // midpoint
208 // Compute left part.
209 var cl_pqcd_series_result<cl_R> L;
210 eval_pqcd_series_aux(Nm,args,L,trunclen,false);
211 // Compute right part.
212 var cl_pqcd_series_result<cl_R> R;
213 eval_pqcd_series_aux(N-Nm,args,R,trunclen,rightmost);
214 // Put together partial results.
217 truncate_precision(Z.P,trunclen);
220 truncate_precision(Z.Q,trunclen);
221 // Z.S = L.S + L.P/L.Q*R.S;
222 var cl_R tmp = L.P * R.T;
223 Z.T = R.Q * L.T + tmp;
224 truncate_precision(Z.T,trunclen);
226 Z.C = L.C * R.D + L.D * R.C;
227 truncate_precision(Z.C,trunclen);
230 truncate_precision(Z.D,trunclen);
231 // Z.U = L.U + L.C/L.D * L.P/L.Q * R.S + L.P/L.Q * R.U;
232 // Z.V = R.D * R.Q * L.V + R.D * L.C * L.P * R.T + L.D * L.P * R.V;
233 Z.V = R.D * (R.Q * L.V + L.C * tmp) + L.D * L.P * R.V;
234 truncate_precision(Z.V,trunclen);