Finalize CLN 1.3.7 release.

[cln.git] / src / float / transcendental / cl_LF_ratsumseries_pqd_aux.cc

diff --git a/src/float/transcendental/cl_LF_ratsumseries_pqd_aux.cc b/src/float/transcendental/cl_LF_ratsumseries_pqd_aux.cc
index f6886bd15cbb504ba3325a6699cf026102acdc0e..b5e85857fc228a43f8cd9e1c0e7fc2540cef45af 100644 (file)
--- a/src/float/transcendental/cl_LF_ratsumseries_pqd_aux.cc
+++ b/src/float/transcendental/cl_LF_ratsumseries_pqd_aux.cc
@@ -1,20 +1,21 @@
 // eval_pqd_series_aux().
 
 // General includes.
-#include "cl_sysdep.h"
+#include "base/cl_sysdep.h"
 
 // Specification.
-#include "cl_LF_tran.h"
+#include "float/transcendental/cl_LF_tran.h"
 
 
 // Implementation.
 
 #include "cln/integer.h"
+#include "cln/real.h"
 #include "cln/exception.h"
 
 namespace cln {
 
-void eval_pqd_series_aux (uintC N, cl_pqd_series_term* args, cl_pqd_series_result& Z, cl_boolean rightmost)
+void eval_pqd_series_aux (uintC N, cl_pqd_series_term* args, cl_pqd_series_result<cl_I>& Z, bool rightmost)
 {
        // N = N2-N1
        switch (N) {
@@ -55,10 +56,10 @@ void eval_pqd_series_aux (uintC N, cl_pqd_series_term* args, cl_pqd_series_resul
        default: {
                var uintC Nm = N/2; // midpoint
                // Compute left part.
-               var cl_pqd_series_result L;
-               eval_pqd_series_aux(Nm,args+0,L,cl_false);
+               var cl_pqd_series_result<cl_I> L;
+               eval_pqd_series_aux(Nm,args+0,L,false);
                // Compute right part.
-               var cl_pqd_series_result R;
+               var cl_pqd_series_result<cl_I> R;
                eval_pqd_series_aux(N-Nm,args+Nm,R,rightmost);
                // Put together partial results.
                if (!rightmost) { Z.P = L.P * R.P; }
@@ -76,4 +77,152 @@ void eval_pqd_series_aux (uintC N, cl_pqd_series_term* args, cl_pqd_series_resul
        }
 }
 
+void eval_pqd_series_aux (uintC N, cl_pqd_series_stream& args, cl_pqd_series_result<cl_I>& Z, bool rightmost)
+{
+       // N = N2-N1
+       switch (N) {
+       case 0:
+               throw runtime_exception(); break;
+       case 1: {
+               var cl_pqd_series_term v0 = args.next(); // [N1]
+               if (!rightmost) { Z.P = v0.p; }
+               Z.Q = v0.q;
+               Z.T = v0.p;
+               if (!rightmost) { Z.C = 1; }
+               Z.D = v0.d;
+               Z.V = v0.p;
+               break;
+               }
+       case 2: {
+               var cl_pqd_series_term v0 = args.next(); // [N1]
+               var cl_pqd_series_term v1 = args.next(); // [N1+1]
+               var cl_I p01 = v0.p * v1.p;
+               if (!rightmost) { Z.P = p01; }
+               Z.Q = v0.q * v1.q;
+               var cl_I p0q1 = v0.p * v1.q + p01;
+               Z.T = p0q1;
+               if (!rightmost) { Z.C = v1.d + v0.d; }
+               Z.D = v0.d * v1.d;
+               Z.V = v1.d * p0q1 + v0.d * p01;
+               break;
+               }
+       case 3: {
+               var cl_pqd_series_term v0 = args.next(); // [N1]
+               var cl_pqd_series_term v1 = args.next(); // [N1+1]
+               var cl_pqd_series_term v2 = args.next(); // [N1+2]
+               var cl_I p01 = v0.p * v1.p;
+               var cl_I p012 = p01 * v2.p;
+               if (!rightmost) { Z.P = p012; }
+               Z.Q = v0.q * v1.q * v2.q;
+               var cl_I p0q1 = v0.p * v1.q + p01;
+               Z.T = v2.q * p0q1 + p012;
+               var cl_I d01 = v0.d * v1.d;
+               if (!rightmost) { Z.C = (v1.d + v0.d) * v2.d + d01; }
+               Z.D = d01 * v2.d;
+               Z.V = v2.d * (v2.q * (v1.d * p0q1 + v0.d * p01) + (v1.d + v0.d) * p012) + d01 * p012;
+               break;
+               }
+       default: {
+               var uintC Nm = N/2; // midpoint
+               // Compute left part.
+               var cl_pqd_series_result<cl_I> L;
+               eval_pqd_series_aux(Nm,args,L,false);
+               // Compute right part.
+               var cl_pqd_series_result<cl_I> R;
+               eval_pqd_series_aux(N-Nm,args,R,rightmost);
+               // Put together partial results.
+               if (!rightmost) { Z.P = L.P * R.P; }
+               Z.Q = L.Q * R.Q;
+               // Z.S = L.S + L.P/L.Q*R.S;
+               var cl_I tmp = L.P * R.T;
+               Z.T = R.Q * L.T + tmp;
+               if (!rightmost) { Z.C = L.C * R.D + L.D * R.C; }
+               Z.D = L.D * R.D;
+               // Z.U = L.U + L.C/L.D * L.P/L.Q * R.S + L.P/L.Q * R.U;
+               // Z.V = R.D * R.Q * L.V + R.D * L.C * L.P * R.T + L.D * L.P * R.V;
+               Z.V = R.D * (R.Q * L.V + L.C * tmp) + L.D * L.P * R.V;
+               break;
+               }
+       }
+}
+
+void eval_pqd_series_aux (uintC N, cl_pqd_series_stream& args, cl_pqd_series_result<cl_R>& Z, uintC trunclen, bool rightmost)
+{
+       // N = N2-N1
+       switch (N) {
+       case 0:
+               throw runtime_exception(); break;
+       case 1: {
+               var cl_pqd_series_term v0 = args.next(); // [N1]
+               if (!rightmost) { Z.P = v0.p; }
+               Z.Q = v0.q;
+               Z.T = v0.p;
+               if (!rightmost) { Z.C = 1; }
+               Z.D = v0.d;
+               Z.V = v0.p;
+               break;
+               }
+       case 2: {
+               var cl_pqd_series_term v0 = args.next(); // [N1]
+               var cl_pqd_series_term v1 = args.next(); // [N1+1]
+               var cl_I p01 = v0.p * v1.p;
+               if (!rightmost) { Z.P = p01; }
+               Z.Q = v0.q * v1.q;
+               var cl_I p0q1 = v0.p * v1.q + p01;
+               Z.T = p0q1;
+               if (!rightmost) { Z.C = v1.d + v0.d; }
+               Z.D = v0.d * v1.d;
+               Z.V = v1.d * p0q1 + v0.d * p01;
+               break;
+               }
+       case 3: {
+               var cl_pqd_series_term v0 = args.next(); // [N1]
+               var cl_pqd_series_term v1 = args.next(); // [N1+1]
+               var cl_pqd_series_term v2 = args.next(); // [N1+2]
+               var cl_I p01 = v0.p * v1.p;
+               var cl_I p012 = p01 * v2.p;
+               if (!rightmost) { Z.P = p012; }
+               Z.Q = v0.q * v1.q * v2.q;
+               var cl_I p0q1 = v0.p * v1.q + p01;
+               Z.T = v2.q * p0q1 + p012;
+               var cl_I d01 = v0.d * v1.d;
+               if (!rightmost) { Z.C = (v1.d + v0.d) * v2.d + d01; }
+               Z.D = d01 * v2.d;
+               Z.V = v2.d * (v2.q * (v1.d * p0q1 + v0.d * p01) + (v1.d + v0.d) * p012) + d01 * p012;
+               break;
+               }
+       default: {
+               var uintC Nm = N/2; // midpoint
+               // Compute left part.
+               var cl_pqd_series_result<cl_R> L;
+               eval_pqd_series_aux(Nm,args,L,trunclen,false);
+               // Compute right part.
+               var cl_pqd_series_result<cl_R> R;
+               eval_pqd_series_aux(N-Nm,args,R,trunclen,rightmost);
+               // Put together partial results.
+               if (!rightmost) {
+                       Z.P = L.P * R.P;
+                       truncate_precision(Z.P,trunclen);
+               }
+               Z.Q = L.Q * R.Q;
+               truncate_precision(Z.Q,trunclen);
+               // Z.S = L.S + L.P/L.Q*R.S;
+               var cl_R tmp = L.P * R.T;
+               Z.T = R.Q * L.T + tmp;
+               truncate_precision(Z.T,trunclen);
+               if (!rightmost) {
+                       Z.C = L.C * R.D + L.D * R.C;
+                       truncate_precision(Z.C,trunclen);
+               }
+               Z.D = L.D * R.D;
+               truncate_precision(Z.D,trunclen);
+               // Z.U = L.U + L.C/L.D * L.P/L.Q * R.S + L.P/L.Q * R.U;
+               // Z.V = R.D * R.Q * L.V + R.D * L.C * L.P * R.T + L.D * L.P * R.V;
+               Z.V = R.D * (R.Q * L.V + L.C * tmp) + L.D * L.P * R.V;
+               truncate_precision(Z.V,trunclen);
+               break;
+               }
+       }
+}
+
 }  // namespace cln