1 // eval_pqd_series_aux().
4 #include "base/cl_sysdep.h"
7 #include "float/transcendental/cl_LF_tran.h"
12 #include "cln/integer.h"
14 #include "cln/exception.h"
18 void eval_pqd_series_aux (uintC N, cl_pqd_series_term* args, cl_pqd_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 = 1; }
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 if (!rightmost) { Z.C = args[1].d + args[0].d; }
39 Z.D = args[0].d * args[1].d;
40 Z.V = args[1].d * p0q1 + args[0].d * p01;
44 var cl_I p01 = args[0].p * args[1].p;
45 var cl_I p012 = p01 * args[2].p;
46 if (!rightmost) { Z.P = p012; }
47 Z.Q = args[0].q * args[1].q * args[2].q;
48 var cl_I p0q1 = args[0].p * args[1].q + p01;
49 Z.T = args[2].q * p0q1 + p012;
50 var cl_I d01 = args[0].d * args[1].d;
51 if (!rightmost) { Z.C = (args[1].d + args[0].d) * args[2].d + d01; }
52 Z.D = d01 * args[2].d;
53 Z.V = args[2].d * (args[2].q * (args[1].d * p0q1 + args[0].d * p01) + (args[1].d + args[0].d) * p012) + d01 * p012;
57 var uintC Nm = N/2; // midpoint
59 var cl_pqd_series_result<cl_I> L;
60 eval_pqd_series_aux(Nm,args+0,L,false);
61 // Compute right part.
62 var cl_pqd_series_result<cl_I> R;
63 eval_pqd_series_aux(N-Nm,args+Nm,R,rightmost);
64 // Put together partial results.
65 if (!rightmost) { Z.P = L.P * R.P; }
67 // Z.S = L.S + L.P/L.Q*R.S;
68 var cl_I tmp = L.P * R.T;
69 Z.T = R.Q * L.T + tmp;
70 if (!rightmost) { Z.C = L.C * R.D + L.D * R.C; }
72 // Z.U = L.U + L.C/L.D * L.P/L.Q * R.S + L.P/L.Q * R.U;
73 // Z.V = R.D * R.Q * L.V + R.D * L.C * L.P * R.T + L.D * L.P * R.V;
74 Z.V = R.D * (R.Q * L.V + L.C * tmp) + L.D * L.P * R.V;
80 void eval_pqd_series_aux (uintC N, cl_pqd_series_stream& args, cl_pqd_series_result<cl_I>& Z, bool rightmost)
85 throw runtime_exception(); break;
87 var cl_pqd_series_term v0 = args.next(); // [N1]
88 if (!rightmost) { Z.P = v0.p; }
91 if (!rightmost) { Z.C = 1; }
97 var cl_pqd_series_term v0 = args.next(); // [N1]
98 var cl_pqd_series_term v1 = args.next(); // [N1+1]
99 var cl_I p01 = v0.p * v1.p;
100 if (!rightmost) { Z.P = p01; }
102 var cl_I p0q1 = v0.p * v1.q + p01;
104 if (!rightmost) { Z.C = v1.d + v0.d; }
106 Z.V = v1.d * p0q1 + v0.d * p01;
110 var cl_pqd_series_term v0 = args.next(); // [N1]
111 var cl_pqd_series_term v1 = args.next(); // [N1+1]
112 var cl_pqd_series_term v2 = args.next(); // [N1+2]
113 var cl_I p01 = v0.p * v1.p;
114 var cl_I p012 = p01 * v2.p;
115 if (!rightmost) { Z.P = p012; }
116 Z.Q = v0.q * v1.q * v2.q;
117 var cl_I p0q1 = v0.p * v1.q + p01;
118 Z.T = v2.q * p0q1 + p012;
119 var cl_I d01 = v0.d * v1.d;
120 if (!rightmost) { Z.C = (v1.d + v0.d) * v2.d + d01; }
122 Z.V = v2.d * (v2.q * (v1.d * p0q1 + v0.d * p01) + (v1.d + v0.d) * p012) + d01 * p012;
126 var uintC Nm = N/2; // midpoint
127 // Compute left part.
128 var cl_pqd_series_result<cl_I> L;
129 eval_pqd_series_aux(Nm,args,L,false);
130 // Compute right part.
131 var cl_pqd_series_result<cl_I> R;
132 eval_pqd_series_aux(N-Nm,args,R,rightmost);
133 // Put together partial results.
134 if (!rightmost) { Z.P = L.P * R.P; }
136 // Z.S = L.S + L.P/L.Q*R.S;
137 var cl_I tmp = L.P * R.T;
138 Z.T = R.Q * L.T + tmp;
139 if (!rightmost) { Z.C = L.C * R.D + L.D * R.C; }
141 // Z.U = L.U + L.C/L.D * L.P/L.Q * R.S + L.P/L.Q * R.U;
142 // Z.V = R.D * R.Q * L.V + R.D * L.C * L.P * R.T + L.D * L.P * R.V;
143 Z.V = R.D * (R.Q * L.V + L.C * tmp) + L.D * L.P * R.V;
149 void eval_pqd_series_aux (uintC N, cl_pqd_series_stream& args, cl_pqd_series_result<cl_R>& Z, uintC trunclen, bool rightmost)
154 throw runtime_exception(); break;
156 var cl_pqd_series_term v0 = args.next(); // [N1]
157 if (!rightmost) { Z.P = v0.p; }
160 if (!rightmost) { Z.C = 1; }
166 var cl_pqd_series_term v0 = args.next(); // [N1]
167 var cl_pqd_series_term v1 = args.next(); // [N1+1]
168 var cl_I p01 = v0.p * v1.p;
169 if (!rightmost) { Z.P = p01; }
171 var cl_I p0q1 = v0.p * v1.q + p01;
173 if (!rightmost) { Z.C = v1.d + v0.d; }
175 Z.V = v1.d * p0q1 + v0.d * p01;
179 var cl_pqd_series_term v0 = args.next(); // [N1]
180 var cl_pqd_series_term v1 = args.next(); // [N1+1]
181 var cl_pqd_series_term v2 = args.next(); // [N1+2]
182 var cl_I p01 = v0.p * v1.p;
183 var cl_I p012 = p01 * v2.p;
184 if (!rightmost) { Z.P = p012; }
185 Z.Q = v0.q * v1.q * v2.q;
186 var cl_I p0q1 = v0.p * v1.q + p01;
187 Z.T = v2.q * p0q1 + p012;
188 var cl_I d01 = v0.d * v1.d;
189 if (!rightmost) { Z.C = (v1.d + v0.d) * v2.d + d01; }
191 Z.V = v2.d * (v2.q * (v1.d * p0q1 + v0.d * p01) + (v1.d + v0.d) * p012) + d01 * p012;
195 var uintC Nm = N/2; // midpoint
196 // Compute left part.
197 var cl_pqd_series_result<cl_R> L;
198 eval_pqd_series_aux(Nm,args,L,trunclen,false);
199 // Compute right part.
200 var cl_pqd_series_result<cl_R> R;
201 eval_pqd_series_aux(N-Nm,args,R,trunclen,rightmost);
202 // Put together partial results.
205 truncate_precision(Z.P,trunclen);
208 truncate_precision(Z.Q,trunclen);
209 // Z.S = L.S + L.P/L.Q*R.S;
210 var cl_R tmp = L.P * R.T;
211 Z.T = R.Q * L.T + tmp;
212 truncate_precision(Z.T,trunclen);
214 Z.C = L.C * R.D + L.D * R.C;
215 truncate_precision(Z.C,trunclen);
218 truncate_precision(Z.D,trunclen);
219 // Z.U = L.U + L.C/L.D * L.P/L.Q * R.S + L.P/L.Q * R.U;
220 // Z.V = R.D * R.Q * L.V + R.D * L.C * L.P * R.T + L.D * L.P * R.V;
221 Z.V = R.D * (R.Q * L.V + L.C * tmp) + L.D * L.P * R.V;
222 truncate_precision(Z.V,trunclen);