src/float/transcendental/cl_LF_ratseries_pqa.cc

   1 // eval_rational_series<bool>().
   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/lfloat.h"
  13 #include "cln/integer.h"
  14 #include "cln/real.h"
  15 #include "cln/exception.h"
  16 #include "float/lfloat/cl_LF.h"
  17 #include "base/cl_alloca.h"
  18
  19 namespace cln {
  20
  21 // Subroutine.
  22 // Evaluates S = sum(N1 <= n < N2, a(n)/b(n) * (p(N1)...p(n))/(q(N1)...q(n)))
  23 // and returns P = p(N1)...p(N2-1), Q = q(N1)...q(N2-1), B = B(N1)...B(N2-1)
  24 // and T = B*Q*S (all integers). On entry N1 < N2.
  25 // P will not be computed if a NULL pointer is passed.
  26
  27 static void eval_pqa_series_aux (uintC N1, uintC N2,
  28                                  const cl_pqa_series& args,
  29                                  cl_I* P, cl_I* Q, cl_I* T)
  30 {
  31         switch (N2 - N1) {
  32         case 0:
  33                 throw runtime_exception(); break;
  34         case 1:
  35                 if (P) { *P = args.pv[N1]; }
  36                 *Q = args.qv[N1];
  37                 *T = args.av[N1] * args.pv[N1];
  38                 break;
  39         case 2: {
  40                 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
  41                 if (P) { *P = p01; }
  42                 *Q = args.qv[N1] * args.qv[N1+1];
  43                 *T = args.qv[N1+1] * args.av[N1] * args.pv[N1]
  44                    + args.av[N1+1] * p01;
  45                 break;
  46                 }
  47         case 3: {
  48                 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
  49                 var cl_I p012 = p01 * args.pv[N1+2];
  50                 if (P) { *P = p012; }
  51                 var cl_I q12 = args.qv[N1+1] * args.qv[N1+2];
  52                 *Q = args.qv[N1] * q12;
  53                 *T = q12 * args.av[N1] * args.pv[N1]
  54                    + args.qv[N1+2] * args.av[N1+1] * p01
  55                    + args.av[N1+2] * p012;
  56                 break;
  57                 }
  58         case 4: {
  59                 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
  60                 var cl_I p012 = p01 * args.pv[N1+2];
  61                 var cl_I p0123 = p012 * args.pv[N1+3];
  62                 if (P) { *P = p0123; }
  63                 var cl_I q23 = args.qv[N1+2] * args.qv[N1+3];
  64                 var cl_I q123 = args.qv[N1+1] * q23;
  65                 *Q = args.qv[N1] * q123;
  66                 *T = q123 * args.av[N1] * args.pv[N1]
  67                    + q23 * args.av[N1+1] * p01
  68                    + args.qv[N1+3] * args.av[N1+2] * p012
  69                    + args.av[N1+3] * p0123;
  70                 break;
  71                 }
  72         default: {
  73                 var uintC Nm = (N1+N2)/2; // midpoint
  74                 // Compute left part.
  75                 var cl_I LP, LQ, LT;
  76                 eval_pqa_series_aux(N1,Nm,args,&LP,&LQ,&LT);
  77                 // Compute right part.
  78                 var cl_I RP, RQ, RT;
  79                 eval_pqa_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RT);
  80                 // Put together partial results.
  81                 if (P) { *P = LP*RP; }
  82                 *Q = LQ*RQ;
  83                 // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT.
  84                 *T = RQ*LT + LP*RT;
  85                 break;
  86                 }
  87         }
  88 }
  89
  90 template<>
  91 const cl_LF eval_rational_series<false> (uintC N, const cl_pqa_series& args, uintC len)
  92 {
  93         if (N==0)
  94                 return cl_I_to_LF(0,len);
  95         var cl_I Q, T;
  96         eval_pqa_series_aux(0,N,args,NULL,&Q,&T);
  97         return cl_I_to_LF(T,len) / cl_I_to_LF(Q,len);
  98 }
  99
 100 static void eval_pqsa_series_aux (uintC N1, uintC N2,
 101                                   const cl_pqa_series& args, const uintC* qsv,
 102                                   cl_I* P, cl_I* Q, uintC* QS, cl_I* T)
 103 {
 104         switch (N2 - N1) {
 105         case 0:
 106                 throw runtime_exception(); break;
 107         case 1:
 108                 if (P) { *P = args.pv[N1]; }
 109                 *Q = args.qv[N1];
 110                 *QS = qsv[N1];
 111                 *T = args.av[N1] * args.pv[N1];
 112                 break;
 113         case 2: {
 114                 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
 115                 if (P) { *P = p01; }
 116                 *Q = args.qv[N1] * args.qv[N1+1];
 117                 *QS = qsv[N1] + qsv[N1+1];
 118                 *T = ((args.qv[N1+1] * args.av[N1] * args.pv[N1]) << qsv[N1+1])
 119                    + args.av[N1+1] * p01;
 120                 break;
 121                 }
 122         case 3: {
 123                 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
 124                 var cl_I p012 = p01 * args.pv[N1+2];
 125                 if (P) { *P = p012; }
 126                 var cl_I q12 = args.qv[N1+1] * args.qv[N1+2];
 127                 *Q = args.qv[N1] * q12;
 128                 *QS = qsv[N1] + qsv[N1+1] + qsv[N1+2];
 129                 *T = ((q12 * args.av[N1] * args.pv[N1]) << (qsv[N1+1] + qsv[N1+2]))
 130                    + ((args.qv[N1+2] * args.av[N1+1] * p01) << qsv[N1+2])
 131                    + args.av[N1+2] * p012;
 132                 break;
 133                 }
 134         case 4: {
 135                 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
 136                 var cl_I p012 = p01 * args.pv[N1+2];
 137                 var cl_I p0123 = p012 * args.pv[N1+3];
 138                 if (P) { *P = p0123; }
 139                 var cl_I q23 = args.qv[N1+2] * args.qv[N1+3];
 140                 var cl_I q123 = args.qv[N1+1] * q23;
 141                 *Q = args.qv[N1] * q123;
 142                 *QS = qsv[N1] + qsv[N1+1] + qsv[N1+2] + qsv[N1+3];
 143                 *T = ((((((q123 * args.av[N1] * args.pv[N1]) << qsv[N1+1])
 144                          + q23 * args.av[N1+1] * p01) << qsv[N1+2])
 145                        + args.qv[N1+3] * args.av[N1+2] * p012) << qsv[N1+3])
 146                    + args.av[N1+3] * p0123;
 147                 break;
 148                 }
 149         default: {
 150                 var uintC Nm = (N1+N2)/2; // midpoint
 151                 // Compute left part.
 152                 var cl_I LP, LQ, LT;
 153                 var uintC LQS;
 154                 eval_pqsa_series_aux(N1,Nm,args,qsv,&LP,&LQ,&LQS,&LT);
 155                 // Compute right part.
 156                 var cl_I RP, RQ, RT;
 157                 var uintC RQS;
 158                 eval_pqsa_series_aux(Nm,N2,args,qsv,(P?&RP:(cl_I*)0),&RQ,&RQS,&RT);
 159                 // Put together partial results.
 160                 if (P) { *P = LP*RP; }
 161                 *Q = LQ*RQ;
 162                 *QS = LQS+RQS;
 163                 // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT.
 164                 *T = ((RQ*LT) << RQS) + LP*RT;
 165                 break;
 166                 }
 167         }
 168 }
 169
 170 template<>
 171 const cl_LF eval_rational_series<true> (uintC N, const cl_pqa_series& args, uintC len)
 172 {
 173         if (N==0)
 174                 return cl_I_to_LF(0,len);
 175         var cl_I Q, T;
 176         // Precomputation of the shift counts:
 177         // Split qv[n] into qv[n]*2^qsv[n].
 178         CL_ALLOCA_STACK;
 179         var uintC* qsv = (uintC*) cl_alloca(N*sizeof(uintC));
 180         var cl_I* qp = args.qv;
 181         var uintC* qsp = qsv;
 182         for (var uintC n = 0; n < N; n++, qp++, qsp++) {
 183                 *qsp = pullout_shiftcount(*qp);
 184         }
 185         // Main computation.
 186         var uintC QS;
 187         eval_pqsa_series_aux(0,N,args,qsv,NULL,&Q,&QS,&T);
 188         return cl_I_to_LF(T,len) / scale_float(cl_I_to_LF(Q,len),QS);
 189 }
 190
 191 static void eval_pqa_series_aux (uintC N1, uintC N2,
 192                                  cl_pqa_series_stream& args,
 193                                  cl_I* P, cl_I* Q, cl_I* T)
 194 {
 195         switch (N2 - N1) {
 196         case 0:
 197                 throw runtime_exception(); break;
 198         case 1: {
 199                 var cl_pqa_series_term v0 = args.next(); // [N1]
 200                 if (P) { *P = v0.p; }
 201                 *Q = v0.q;
 202                 *T = v0.a * v0.p;
 203                 break;
 204                 }
 205         case 2: {
 206                 var cl_pqa_series_term v0 = args.next(); // [N1]
 207                 var cl_pqa_series_term v1 = args.next(); // [N1+1]
 208                 var cl_I p01 = v0.p * v1.p;
 209                 if (P) { *P = p01; }
 210                 *Q = v0.q * v1.q;
 211                 *T = v1.q * v0.a * v0.p
 212                    + v1.a * p01;
 213                 break;
 214                 }
 215         case 3: {
 216                 var cl_pqa_series_term v0 = args.next(); // [N1]
 217                 var cl_pqa_series_term v1 = args.next(); // [N1+1]
 218                 var cl_pqa_series_term v2 = args.next(); // [N1+2]
 219                 var cl_I p01 = v0.p * v1.p;
 220                 var cl_I p012 = p01 * v2.p;
 221                 if (P) { *P = p012; }
 222                 var cl_I q12 = v1.q * v2.q;
 223                 *Q = v0.q * q12;
 224                 *T = q12 * v0.a * v0.p
 225                    + v2.q * v1.a * p01
 226                    + v2.a * p012;
 227                 break;
 228                 }
 229         case 4: {
 230                 var cl_pqa_series_term v0 = args.next(); // [N1]
 231                 var cl_pqa_series_term v1 = args.next(); // [N1+1]
 232                 var cl_pqa_series_term v2 = args.next(); // [N1+2]
 233                 var cl_pqa_series_term v3 = args.next(); // [N1+3]
 234                 var cl_I p01 = v0.p * v1.p;
 235                 var cl_I p012 = p01 * v2.p;
 236                 var cl_I p0123 = p012 * v3.p;
 237                 if (P) { *P = p0123; }
 238                 var cl_I q23 = v2.q * v3.q;
 239                 var cl_I q123 = v1.q * q23;
 240                 *Q = v0.q * q123;
 241                 *T = q123 * v0.a * v0.p
 242                    + q23 * v1.a * p01
 243                    + v3.q * v2.a * p012
 244                    + v3.a * p0123;
 245                 break;
 246                 }
 247         default: {
 248                 var uintC Nm = (N1+N2)/2; // midpoint
 249                 // Compute left part.
 250                 var cl_I LP, LQ, LT;
 251                 eval_pqa_series_aux(N1,Nm,args,&LP,&LQ,&LT);
 252                 // Compute right part.
 253                 var cl_I RP, RQ, RT;
 254                 eval_pqa_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RT);
 255                 // Put together partial results.
 256                 if (P) { *P = LP*RP; }
 257                 *Q = LQ*RQ;
 258                 // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT.
 259                 *T = RQ*LT + LP*RT;
 260                 break;
 261                 }
 262         }
 263 }
 264
 265 template<>
 266 const cl_LF eval_rational_series<false> (uintC N, cl_pqa_series_stream& args, uintC len)
 267 {
 268         if (N==0)
 269                 return cl_I_to_LF(0,len);
 270         var cl_I Q, T;
 271         eval_pqa_series_aux(0,N,args,NULL,&Q,&T);
 272         return cl_I_to_LF(T,len) / cl_I_to_LF(Q,len);
 273 }
 274
 275 static void eval_pqa_series_aux (uintC N1, uintC N2,
 276                                  cl_pqa_series_stream& args,
 277                                  cl_R* P, cl_R* Q, cl_R* T,
 278                                  uintC trunclen)
 279 {
 280         switch (N2 - N1) {
 281         case 0:
 282                 throw runtime_exception(); break;
 283         case 1: {
 284                 var cl_pqa_series_term v0 = args.next(); // [N1]
 285                 if (P) { *P = v0.p; }
 286                 *Q = v0.q;
 287                 *T = v0.a * v0.p;
 288                 break;
 289                 }
 290         case 2: {
 291                 var cl_pqa_series_term v0 = args.next(); // [N1]
 292                 var cl_pqa_series_term v1 = args.next(); // [N1+1]
 293                 var cl_I p01 = v0.p * v1.p;
 294                 if (P) { *P = p01; }
 295                 *Q = v0.q * v1.q;
 296                 *T = v1.q * v0.a * v0.p
 297                    + v1.a * p01;
 298                 break;
 299                 }
 300         case 3: {
 301                 var cl_pqa_series_term v0 = args.next(); // [N1]
 302                 var cl_pqa_series_term v1 = args.next(); // [N1+1]
 303                 var cl_pqa_series_term v2 = args.next(); // [N1+2]
 304                 var cl_I p01 = v0.p * v1.p;
 305                 var cl_I p012 = p01 * v2.p;
 306                 if (P) { *P = p012; }
 307                 var cl_I q12 = v1.q * v2.q;
 308                 *Q = v0.q * q12;
 309                 *T = q12 * v0.a * v0.p
 310                    + v2.q * v1.a * p01
 311                    + v2.a * p012;
 312                 break;
 313                 }
 314         case 4: {
 315                 var cl_pqa_series_term v0 = args.next(); // [N1]
 316                 var cl_pqa_series_term v1 = args.next(); // [N1+1]
 317                 var cl_pqa_series_term v2 = args.next(); // [N1+2]
 318                 var cl_pqa_series_term v3 = args.next(); // [N1+3]
 319                 var cl_I p01 = v0.p * v1.p;
 320                 var cl_I p012 = p01 * v2.p;
 321                 var cl_I p0123 = p012 * v3.p;
 322                 if (P) { *P = p0123; }
 323                 var cl_I q23 = v2.q * v3.q;
 324                 var cl_I q123 = v1.q * q23;
 325                 *Q = v0.q * q123;
 326                 *T = q123 * v0.a * v0.p
 327                    + q23 * v1.a * p01
 328                    + v3.q * v2.a * p012
 329                    + v3.a * p0123;
 330                 break;
 331                 }
 332         default: {
 333                 var uintC Nm = (N1+N2)/2; // midpoint
 334                 // Compute left part.
 335                 var cl_R LP, LQ, LT;
 336                 eval_pqa_series_aux(N1,Nm,args,&LP,&LQ,&LT,trunclen);
 337                 // Compute right part.
 338                 var cl_R RP, RQ, RT;
 339                 eval_pqa_series_aux(Nm,N2,args,(P?&RP:(cl_R*)0),&RQ,&RT,trunclen);
 340                 // Put together partial results.
 341                 if (P) {
 342                         *P = LP*RP;
 343                         truncate_precision(*P,trunclen);
 344                 }
 345                 *Q = LQ*RQ;
 346                 truncate_precision(*Q,trunclen);
 347                 // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT.
 348                 *T = RQ*LT + LP*RT;
 349                 truncate_precision(*T,trunclen);
 350                 break;
 351                 }
 352         }
 353 }
 354
 355 template<>
 356 const cl_LF eval_rational_series<false> (uintC N, cl_pqa_series_stream& args, uintC len, uintC trunclen)
 357 {
 358         if (N==0)
 359                 return cl_I_to_LF(0,len);
 360         var cl_R Q, T;
 361         eval_pqa_series_aux(0,N,args,NULL,&Q,&T,trunclen);
 362         return cl_R_to_LF(T,len) / cl_R_to_LF(Q,len);
 363 }
 364 // Bit complexity (if p(n), q(n), a(n), b(n) have length O(log(n))):
 365 // O(log(N)^2*M(N)).
 366
 367 }  // namespace cln