src/float/transcendental/cl_LF_ratsumseries_pqcd_aux.cc

   1 // eval_pqcd_series_aux().
   2
   3 // General includes.
   4 #include "base/cl_sysdep.h"
   5
   6 // Specification.
   7 #include "float/transcendental/cl_LF_tran.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cln/integer.h"
  13 #include "cln/real.h"
  14 #include "cln/exception.h"
  15
  16 namespace cln {
  17
  18 void eval_pqcd_series_aux (uintC N, cl_pqcd_series_term* args, cl_pqcd_series_result<cl_I>& Z, bool rightmost)
  19 {
  20         // N = N2-N1
  21         switch (N) {
  22         case 0:
  23                 throw runtime_exception(); break;
  24         case 1:
  25                 if (!rightmost) { Z.P = args[0].p; }
  26                 Z.Q = args[0].q;
  27                 Z.T = args[0].p;
  28                 if (!rightmost) { Z.C = args[0].c; }
  29                 Z.D = args[0].d;
  30                 Z.V = args[0].c * args[0].p;
  31                 break;
  32         case 2: {
  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;
  37                 Z.T = p0q1;
  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;
  43                 break;
  44                 }
  45         case 3: {
  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;
  58                 break;
  59                 }
  60         default: {
  61                 var uintC Nm = N/2; // midpoint
  62                 // Compute left part.
  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; }
  70                 Z.Q = L.Q * R.Q;
  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; }
  75                 Z.D = L.D * R.D;
  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;
  79                 break;
  80                 }
  81         }
  82 }
  83
  84 void eval_pqcd_series_aux (uintC N, cl_pqcd_series_stream& args, cl_pqcd_series_result<cl_I>& Z, bool rightmost)
  85 {
  86         // N = N2-N1
  87         switch (N) {
  88         case 0:
  89                 throw runtime_exception(); break;
  90         case 1: {
  91                 var cl_pqcd_series_term v0 = args.next(); // [N1]
  92                 if (!rightmost) { Z.P = v0.p; }
  93                 Z.Q = v0.q;
  94                 Z.T = v0.p;
  95                 if (!rightmost) { Z.C = v0.c; }
  96                 Z.D = v0.d;
  97                 Z.V = v0.c * v0.p;
  98                 break;
  99                 }
 100         case 2: {
 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; }
 105                 Z.Q = v0.q * v1.q;
 106                 var cl_I p0q1 = v0.p * v1.q + p01;
 107                 Z.T = p0q1;
 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; }
 111                 Z.D = v0.d * v1.d;
 112                 Z.V = c0d1 * p0q1 + c1d0 * p01;
 113                 break;
 114                 }
 115         case 3: {
 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; }
 129                 Z.D = d01 * v2.d;
 130                 Z.V = v2.d * (v2.q * (c0d1 * p0q1 + c1d0 * p01) + (c0d1 + c1d0) * p012) + v2.c * d01 * p012;
 131                 break;
 132                 }
 133         default: {
 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; }
 143                 Z.Q = L.Q * R.Q;
 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; }
 148                 Z.D = L.D * R.D;
 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;
 152                 break;
 153                 }
 154         }
 155 }
 156
 157 void eval_pqcd_series_aux (uintC N, cl_pqcd_series_stream& args, cl_pqcd_series_result<cl_R>& Z, uintC trunclen, bool rightmost)
 158 {
 159         // N = N2-N1
 160         switch (N) {
 161         case 0:
 162                 throw runtime_exception(); break;
 163         case 1: {
 164                 var cl_pqcd_series_term v0 = args.next(); // [N1]
 165                 if (!rightmost) { Z.P = v0.p; }
 166                 Z.Q = v0.q;
 167                 Z.T = v0.p;
 168                 if (!rightmost) { Z.C = v0.c; }
 169                 Z.D = v0.d;
 170                 Z.V = v0.c * v0.p;
 171                 break;
 172                 }
 173         case 2: {
 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; }
 178                 Z.Q = v0.q * v1.q;
 179                 var cl_I p0q1 = v0.p * v1.q + p01;
 180                 Z.T = p0q1;
 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; }
 184                 Z.D = v0.d * v1.d;
 185                 Z.V = c0d1 * p0q1 + c1d0 * p01;
 186                 break;
 187                 }
 188         case 3: {
 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; }
 202                 Z.D = d01 * v2.d;
 203                 Z.V = v2.d * (v2.q * (c0d1 * p0q1 + c1d0 * p01) + (c0d1 + c1d0) * p012) + v2.c * d01 * p012;
 204                 break;
 205                 }
 206         default: {
 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.
 215                 if (!rightmost) {
 216                         Z.P = L.P * R.P;
 217                         truncate_precision(Z.P,trunclen);
 218                 }
 219                 Z.Q = L.Q * R.Q;
 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);
 225                 if (!rightmost) {
 226                         Z.C = L.C * R.D + L.D * R.C;
 227                         truncate_precision(Z.C,trunclen);
 228                 }
 229                 Z.D = L.D * R.D;
 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);
 235                 break;
 236                 }
 237         }
 238 }
 239
 240 }  // namespace cln