src/float/transcendental/cl_LF_ratsumseries_pqcd_aux.cc

   1 // eval_pqcd_series_aux().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_LF_tran.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cln/integer.h"
  13 #include "cln/abort.h"
  14
  15 namespace cln {
  16
  17 void eval_pqcd_series_aux (uintL N, cl_pqcd_series_term* args, cl_pqcd_series_result& Z, cl_boolean rightmost)
  18 {
  19         // N = N2-N1
  20         switch (N) {
  21         case 0:
  22                 cl_abort(); break;
  23         case 1:
  24                 if (!rightmost) { Z.P = args[0].p; }
  25                 Z.Q = args[0].q;
  26                 Z.T = args[0].p;
  27                 if (!rightmost) { Z.C = args[0].c; }
  28                 Z.D = args[0].d;
  29                 Z.V = args[0].c * args[0].p;
  30                 break;
  31         case 2: {
  32                 var cl_I p01 = args[0].p * args[1].p;
  33                 if (!rightmost) { Z.P = p01; }
  34                 Z.Q = args[0].q * args[1].q;
  35                 var cl_I p0q1 = args[0].p * args[1].q + p01;
  36                 Z.T = p0q1;
  37                 var cl_I c0d1 = args[0].c * args[1].d;
  38                 var cl_I c1d0 = args[1].c * args[0].d;
  39                 if (!rightmost) { Z.C = c0d1 + c1d0; }
  40                 Z.D = args[0].d * args[1].d;
  41                 Z.V = c0d1 * p0q1 + c1d0 * p01;
  42                 break;
  43                 }
  44         case 3: {
  45                 var cl_I p01 = args[0].p * args[1].p;
  46                 var cl_I p012 = p01 * args[2].p;
  47                 if (!rightmost) { Z.P = p012; }
  48                 Z.Q = args[0].q * args[1].q * args[2].q;
  49                 var cl_I p0q1 = args[0].p * args[1].q + p01;
  50                 Z.T = args[2].q * p0q1 + p012;
  51                 var cl_I c0d1 = args[0].c * args[1].d;
  52                 var cl_I c1d0 = args[1].c * args[0].d;
  53                 var cl_I d01 = args[0].d * args[1].d;
  54                 if (!rightmost) { Z.C = (c0d1 + c1d0) * args[2].d + args[2].c * d01; }
  55                 Z.D = d01 * args[2].d;
  56                 Z.V = args[2].d * (args[2].q * (c0d1 * p0q1 + c1d0 * p01) + (c0d1 + c1d0) * p012) + args[2].c * d01 * p012;
  57                 break;
  58                 }
  59         default: {
  60                 var uintL Nm = N/2; // midpoint
  61                 // Compute left part.
  62                 var cl_pqcd_series_result L;
  63                 eval_pqcd_series_aux(Nm,args+0,L,cl_false);
  64                 // Compute right part.
  65                 var cl_pqcd_series_result R;
  66                 eval_pqcd_series_aux(N-Nm,args+Nm,R,rightmost);
  67                 // Put together partial results.
  68                 if (!rightmost) { Z.P = L.P * R.P; }
  69                 Z.Q = L.Q * R.Q;
  70                 // Z.S = L.S + L.P/L.Q*R.S;
  71                 var cl_I tmp = L.P * R.T;
  72                 Z.T = R.Q * L.T + tmp;
  73                 if (!rightmost) { Z.C = L.C * R.D + L.D * R.C; }
  74                 Z.D = L.D * R.D;
  75                 // Z.U = L.U + L.C/L.D * L.P/L.Q * R.S + L.P/L.Q * R.U;
  76                 // Z.V = R.D * R.Q * L.V + R.D * L.C * L.P * R.T + L.D * L.P * R.V;
  77                 Z.V = R.D * (R.Q * L.V + L.C * tmp) + L.D * L.P * R.V;
  78                 break;
  79                 }
  80         }
  81 }
  82
  83 }  // namespace cln