src/float/transcendental/cl_LF_ratsumseries_pqd_aux.cc

   1 // eval_pqd_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/real.h"
  14 #include "cln/exception.h"
  15
  16 namespace cln {
  17
  18 void eval_pqd_series_aux (uintC N, cl_pqd_series_term* args, cl_pqd_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 = 1; }
  29                 Z.D = args[0].d;
  30                 Z.V = 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                 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;
  41                 break;
  42                 }
  43         case 3: {
  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;
  54                 break;
  55                 }
  56         default: {
  57                 var uintC Nm = N/2; // midpoint
  58                 // Compute left part.
  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; }
  66                 Z.Q = L.Q * R.Q;
  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; }
  71                 Z.D = L.D * R.D;
  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;
  75                 break;
  76                 }
  77         }
  78 }
  79
  80 void eval_pqd_series_aux (uintC N, cl_pqd_series_stream& args, cl_pqd_series_result<cl_I>& Z, bool rightmost)
  81 {
  82         // N = N2-N1
  83         switch (N) {
  84         case 0:
  85                 throw runtime_exception(); break;
  86         case 1: {
  87                 var cl_pqd_series_term v0 = args.next(); // [N1]
  88                 if (!rightmost) { Z.P = v0.p; }
  89                 Z.Q = v0.q;
  90                 Z.T = v0.p;
  91                 if (!rightmost) { Z.C = 1; }
  92                 Z.D = v0.d;
  93                 Z.V = v0.p;
  94                 break;
  95                 }
  96         case 2: {
  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; }
 101                 Z.Q = v0.q * v1.q;
 102                 var cl_I p0q1 = v0.p * v1.q + p01;
 103                 Z.T = p0q1;
 104                 if (!rightmost) { Z.C = v1.d + v0.d; }
 105                 Z.D = v0.d * v1.d;
 106                 Z.V = v1.d * p0q1 + v0.d * p01;
 107                 break;
 108                 }
 109         case 3: {
 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; }
 121                 Z.D = d01 * v2.d;
 122                 Z.V = v2.d * (v2.q * (v1.d * p0q1 + v0.d * p01) + (v1.d + v0.d) * p012) + d01 * p012;
 123                 break;
 124                 }
 125         default: {
 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; }
 135                 Z.Q = L.Q * R.Q;
 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; }
 140                 Z.D = L.D * R.D;
 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;
 144                 break;
 145                 }
 146         }
 147 }
 148
 149 void eval_pqd_series_aux (uintC N, cl_pqd_series_stream& args, cl_pqd_series_result<cl_R>& Z, uintC trunclen, bool rightmost)
 150 {
 151         // N = N2-N1
 152         switch (N) {
 153         case 0:
 154                 throw runtime_exception(); break;
 155         case 1: {
 156                 var cl_pqd_series_term v0 = args.next(); // [N1]
 157                 if (!rightmost) { Z.P = v0.p; }
 158                 Z.Q = v0.q;
 159                 Z.T = v0.p;
 160                 if (!rightmost) { Z.C = 1; }
 161                 Z.D = v0.d;
 162                 Z.V = v0.p;
 163                 break;
 164                 }
 165         case 2: {
 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; }
 170                 Z.Q = v0.q * v1.q;
 171                 var cl_I p0q1 = v0.p * v1.q + p01;
 172                 Z.T = p0q1;
 173                 if (!rightmost) { Z.C = v1.d + v0.d; }
 174                 Z.D = v0.d * v1.d;
 175                 Z.V = v1.d * p0q1 + v0.d * p01;
 176                 break;
 177                 }
 178         case 3: {
 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; }
 190                 Z.D = d01 * v2.d;
 191                 Z.V = v2.d * (v2.q * (v1.d * p0q1 + v0.d * p01) + (v1.d + v0.d) * p012) + d01 * p012;
 192                 break;
 193                 }
 194         default: {
 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.
 203                 if (!rightmost) {
 204                         Z.P = L.P * R.P;
 205                         truncate_precision(Z.P,trunclen);
 206                 }
 207                 Z.Q = L.Q * R.Q;
 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);
 213                 if (!rightmost) {
 214                         Z.C = L.C * R.D + L.D * R.C;
 215                         truncate_precision(Z.C,trunclen);
 216                 }
 217                 Z.D = L.D * R.D;
 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);
 223                 break;
 224                 }
 225         }
 226 }
 227
 228 }  // namespace cln