src/float/transcendental/cl_LF_ratseries_pq.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_pq_series_aux (uintC N1, uintC N2,
  28                                 const cl_pq_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.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.pv[N1]
  44                    + 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.pv[N1]
  54                    + args.qv[N1+2] * p01
  55                    + 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.pv[N1]
  67                    + q23 * p01
  68                    + args.qv[N1+3] * p012
  69                    + 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_pq_series_aux(N1,Nm,args,&LP,&LQ,&LT);
  77                 // Compute right part.
  78                 var cl_I RP, RQ, RT;
  79                 eval_pq_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_pq_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_pq_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_pqs_series_aux (uintC N1, uintC N2,
 101                                  const cl_pq_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.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.pv[N1]) << qsv[N1+1])
 119                    + 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.pv[N1]) << (qsv[N1+1] + qsv[N1+2]))
 130                    + ((args.qv[N1+2] * p01) << qsv[N1+2])
 131                    + 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.pv[N1]) << qsv[N1+1])
 144                          + q23 * p01) << qsv[N1+2])
 145                        + args.qv[N1+3] * p012) << qsv[N1+3])
 146                    + 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_pqs_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_pqs_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_pq_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_pqs_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_pq_series_aux (uintC N1, uintC N2,
 192                                 cl_pq_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_pq_series_term v0 = args.next(); // [N1]
 200                 if (P) { *P = v0.p; }
 201                 *Q = v0.q;
 202                 *T = v0.p;
 203                 break;
 204                 }
 205         case 2: {
 206                 var cl_pq_series_term v0 = args.next(); // [N1]
 207                 var cl_pq_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.p
 212                    + p01;
 213                 break;
 214                 }
 215         case 3: {
 216                 var cl_pq_series_term v0 = args.next(); // [N1]
 217                 var cl_pq_series_term v1 = args.next(); // [N1+1]
 218                 var cl_pq_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.p
 225                    + v2.q * p01
 226                    + p012;
 227                 break;
 228                 }
 229         case 4: {
 230                 var cl_pq_series_term v0 = args.next(); // [N1]
 231                 var cl_pq_series_term v1 = args.next(); // [N1+1]
 232                 var cl_pq_series_term v2 = args.next(); // [N1+2]
 233                 var cl_pq_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.p
 242                    + q23 * p01
 243                    + v3.q * p012
 244                    + 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_pq_series_aux(N1,Nm,args,&LP,&LQ,&LT);
 252                 // Compute right part.
 253                 var cl_I RP, RQ, RT;
 254                 eval_pq_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_pq_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_pq_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_pqs_series_aux (uintC N1, uintC N2,
 276                                  cl_pq_series_stream& args,
 277                                  cl_I* P, cl_I* Q, uintC* QS, cl_I* T)
 278 {
 279         switch (N2 - N1) {
 280         case 0:
 281                 throw runtime_exception(); break;
 282         case 1: {
 283                 var cl_pq_series_term v0 = args.next(); // [N1]
 284                 var uintC qs0 = pullout_shiftcount(v0.q);
 285                 if (P) { *P = v0.p; }
 286                 *Q = v0.q;
 287                 *QS = qs0;
 288                 *T = v0.p;
 289                 break;
 290                 }
 291         case 2: {
 292                 var cl_pq_series_term v0 = args.next(); // [N1]
 293                 var cl_pq_series_term v1 = args.next(); // [N1+1]
 294                 var uintC qs0 = pullout_shiftcount(v0.q);
 295                 var uintC qs1 = pullout_shiftcount(v1.q);
 296                 var cl_I p01 = v0.p * v1.p;
 297                 if (P) { *P = p01; }
 298                 *Q = v0.q * v1.q;
 299                 *QS = qs0 + qs1;
 300                 *T = ((v1.q * v0.p) << qs1)
 301                    + p01;
 302                 break;
 303                 }
 304         case 3: {
 305                 var cl_pq_series_term v0 = args.next(); // [N1]
 306                 var cl_pq_series_term v1 = args.next(); // [N1+1]
 307                 var cl_pq_series_term v2 = args.next(); // [N1+2]
 308                 var uintC qs0 = pullout_shiftcount(v0.q);
 309                 var uintC qs1 = pullout_shiftcount(v1.q);
 310                 var uintC qs2 = pullout_shiftcount(v2.q);
 311                 var cl_I p01 = v0.p * v1.p;
 312                 var cl_I p012 = p01 * v2.p;
 313                 if (P) { *P = p012; }
 314                 var cl_I q12 = v1.q * v2.q;
 315                 *Q = v0.q * q12;
 316                 *QS = qs0 + qs1 + qs2;
 317                 *T = ((q12 * v0.p) << (qs1 + qs2))
 318                    + ((v2.q * p01) << qs2)
 319                    + p012;
 320                 break;
 321                 }
 322         case 4: {
 323                 var cl_pq_series_term v0 = args.next(); // [N1]
 324                 var cl_pq_series_term v1 = args.next(); // [N1+1]
 325                 var cl_pq_series_term v2 = args.next(); // [N1+2]
 326                 var cl_pq_series_term v3 = args.next(); // [N1+3]
 327                 var uintC qs0 = pullout_shiftcount(v0.q);
 328                 var uintC qs1 = pullout_shiftcount(v1.q);
 329                 var uintC qs2 = pullout_shiftcount(v2.q);
 330                 var uintC qs3 = pullout_shiftcount(v3.q);
 331                 var cl_I p01 = v0.p * v1.p;
 332                 var cl_I p012 = p01 * v2.p;
 333                 var cl_I p0123 = p012 * v3.p;
 334                 if (P) { *P = p0123; }
 335                 var cl_I q23 = v2.q * v3.q;
 336                 var cl_I q123 = v1.q * q23;
 337                 *Q = v0.q * q123;
 338                 *QS = qs0 + qs1 + qs2 + qs3;
 339                 *T = ((((((q123 * v0.p) << qs1)
 340                          + q23 * p01) << qs2)
 341                        + v3.q * p012) << qs3)
 342                    + p0123;
 343                 break;
 344                 }
 345         default: {
 346                 var uintC Nm = (N1+N2)/2; // midpoint
 347                 // Compute left part.
 348                 var cl_I LP, LQ, LT;
 349                 var uintC LQS;
 350                 eval_pqs_series_aux(N1,Nm,args,&LP,&LQ,&LQS,&LT);
 351                 // Compute right part.
 352                 var cl_I RP, RQ, RT;
 353                 var uintC RQS;
 354                 eval_pqs_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RQS,&RT);
 355                 // Put together partial results.
 356                 if (P) { *P = LP*RP; }
 357                 *Q = LQ*RQ;
 358                 *QS = LQS+RQS;
 359                 // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT.
 360                 *T = ((RQ*LT) << RQS) + LP*RT;
 361                 break;
 362                 }
 363         }
 364 }
 365
 366 template<>
 367 const cl_LF eval_rational_series<true> (uintC N, cl_pq_series_stream& args, uintC len)
 368 {
 369         if (N==0)
 370                 return cl_I_to_LF(0,len);
 371         var cl_I Q, T;
 372         var uintC QS;
 373         eval_pqs_series_aux(0,N,args,NULL,&Q,&QS,&T);
 374         return cl_I_to_LF(T,len) / scale_float(cl_I_to_LF(Q,len),QS);
 375 }
 376
 377 static void eval_pq_series_aux (uintC N1, uintC N2,
 378                                 cl_pq_series_stream& args,
 379                                 cl_R* P, cl_R* Q, cl_R* T,
 380                                 uintC trunclen)
 381 {
 382         switch (N2 - N1) {
 383         case 0:
 384                 throw runtime_exception(); break;
 385         case 1: {
 386                 var cl_pq_series_term v0 = args.next(); // [N1]
 387                 if (P) { *P = v0.p; }
 388                 *Q = v0.q;
 389                 *T = v0.p;
 390                 break;
 391                 }
 392         case 2: {
 393                 var cl_pq_series_term v0 = args.next(); // [N1]
 394                 var cl_pq_series_term v1 = args.next(); // [N1+1]
 395                 var cl_I p01 = v0.p * v1.p;
 396                 if (P) { *P = p01; }
 397                 *Q = v0.q * v1.q;
 398                 *T = v1.q * v0.p
 399                    + p01;
 400                 break;
 401                 }
 402         case 3: {
 403                 var cl_pq_series_term v0 = args.next(); // [N1]
 404                 var cl_pq_series_term v1 = args.next(); // [N1+1]
 405                 var cl_pq_series_term v2 = args.next(); // [N1+2]
 406                 var cl_I p01 = v0.p * v1.p;
 407                 var cl_I p012 = p01 * v2.p;
 408                 if (P) { *P = p012; }
 409                 var cl_I q12 = v1.q * v2.q;
 410                 *Q = v0.q * q12;
 411                 *T = q12 * v0.p
 412                    + v2.q * p01
 413                    + p012;
 414                 break;
 415                 }
 416         case 4: {
 417                 var cl_pq_series_term v0 = args.next(); // [N1]
 418                 var cl_pq_series_term v1 = args.next(); // [N1+1]
 419                 var cl_pq_series_term v2 = args.next(); // [N1+2]
 420                 var cl_pq_series_term v3 = args.next(); // [N1+3]
 421                 var cl_I p01 = v0.p * v1.p;
 422                 var cl_I p012 = p01 * v2.p;
 423                 var cl_I p0123 = p012 * v3.p;
 424                 if (P) { *P = p0123; }
 425                 var cl_I q23 = v2.q * v3.q;
 426                 var cl_I q123 = v1.q * q23;
 427                 *Q = v0.q * q123;
 428                 *T = q123 * v0.p
 429                    + q23 * p01
 430                    + v3.q * p012
 431                    + p0123;
 432                 break;
 433                 }
 434         default: {
 435                 var uintC Nm = (N1+N2)/2; // midpoint
 436                 // Compute left part.
 437                 var cl_R LP, LQ, LT;
 438                 eval_pq_series_aux(N1,Nm,args,&LP,&LQ,&LT,trunclen);
 439                 // Compute right part.
 440                 var cl_R RP, RQ, RT;
 441                 eval_pq_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RT,trunclen);
 442                 // Put together partial results.
 443                 if (P) {
 444                         *P = LP*RP;
 445                         truncate_precision(*P,trunclen);
 446                 }
 447                 *Q = LQ*RQ;
 448                 truncate_precision(*Q,trunclen);
 449                 // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT.
 450                 *T = RQ*LT + LP*RT;
 451                 truncate_precision(*T,trunclen);
 452                 break;
 453                 }
 454         }
 455 }
 456
 457 template<>
 458 const cl_LF eval_rational_series<false> (uintC N, cl_pq_series_stream& args, uintC len, uintC trunclen)
 459 {
 460         if (N==0)
 461                 return cl_I_to_LF(0,len);
 462         var cl_R Q, T;
 463         eval_pq_series_aux(0,N,args,NULL,&Q,&T,trunclen);
 464         return cl_R_to_LF(T,len) / cl_R_to_LF(Q,len);
 465 }
 466 // Bit complexity (if p(n), q(n), a(n), b(n) have length O(log(n))):
 467 // O(log(N)^2*M(N)).
 468
 469 }  // namespace cln