src/float/transcendental/cl_LF_ratseries_pq.cc

   1 // eval_rational_series().
   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/lfloat.h"
  13 #include "cln/integer.h"
  14 #include "cln/real.h"
  15 #include "cln/exception.h"
  16 #include "cl_LF.h"
  17
  18 namespace cln {
  19
  20 // Subroutine.
  21 // Evaluates S = sum(N1 <= n < N2, a(n)/b(n) * (p(N1)...p(n))/(q(N1)...q(n)))
  22 // and returns P = p(N1)...p(N2-1), Q = q(N1)...q(N2-1), B = B(N1)...B(N2-1)
  23 // and T = B*Q*S (all integers). On entry N1 < N2.
  24 // P will not be computed if a NULL pointer is passed.
  25
  26 static void eval_pq_series_aux (uintC N1, uintC N2,
  27                                 const cl_pq_series& args,
  28                                 cl_I* P, cl_I* Q, cl_I* T)
  29 {
  30         switch (N2 - N1) {
  31         case 0:
  32                 throw runtime_exception(); break;
  33         case 1:
  34                 if (P) { *P = args.pv[N1]; }
  35                 *Q = args.qv[N1];
  36                 *T = args.pv[N1];
  37                 break;
  38         case 2: {
  39                 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
  40                 if (P) { *P = p01; }
  41                 *Q = args.qv[N1] * args.qv[N1+1];
  42                 *T = args.qv[N1+1] * args.pv[N1]
  43                    + p01;
  44                 break;
  45                 }
  46         case 3: {
  47                 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
  48                 var cl_I p012 = p01 * args.pv[N1+2];
  49                 if (P) { *P = p012; }
  50                 var cl_I q12 = args.qv[N1+1] * args.qv[N1+2];
  51                 *Q = args.qv[N1] * q12;
  52                 *T = q12 * args.pv[N1]
  53                    + args.qv[N1+2] * p01
  54                    + p012;
  55                 break;
  56                 }
  57         case 4: {
  58                 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
  59                 var cl_I p012 = p01 * args.pv[N1+2];
  60                 var cl_I p0123 = p012 * args.pv[N1+3];
  61                 if (P) { *P = p0123; }
  62                 var cl_I q23 = args.qv[N1+2] * args.qv[N1+3];
  63                 var cl_I q123 = args.qv[N1+1] * q23;
  64                 *Q = args.qv[N1] * q123;
  65                 *T = q123 * args.pv[N1]
  66                    + q23 * p01
  67                    + args.qv[N1+3] * p012
  68                    + p0123;
  69                 break;
  70                 }
  71         default: {
  72                 var uintC Nm = (N1+N2)/2; // midpoint
  73                 // Compute left part.
  74                 var cl_I LP, LQ, LT;
  75                 eval_pq_series_aux(N1,Nm,args,&LP,&LQ,&LT);
  76                 // Compute right part.
  77                 var cl_I RP, RQ, RT;
  78                 eval_pq_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RT);
  79                 // Put together partial results.
  80                 if (P) { *P = LP*RP; }
  81                 *Q = LQ*RQ;
  82                 // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT.
  83                 *T = RQ*LT + LP*RT;
  84                 break;
  85                 }
  86         }
  87 }
  88
  89 static void eval_pqs_series_aux (uintC N1, uintC N2,
  90                                  const cl_pq_series& args,
  91                                  cl_I* P, cl_I* Q, uintC* QS, cl_I* T)
  92 {
  93         switch (N2 - N1) {
  94         case 0:
  95                 throw runtime_exception(); break;
  96         case 1:
  97                 if (P) { *P = args.pv[N1]; }
  98                 *Q = args.qv[N1];
  99                 *QS = args.qsv[N1];
 100                 *T = args.pv[N1];
 101                 break;
 102         case 2: {
 103                 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
 104                 if (P) { *P = p01; }
 105                 *Q = args.qv[N1] * args.qv[N1+1];
 106                 *QS = args.qsv[N1] + args.qsv[N1+1];
 107                 *T = ((args.qv[N1+1] * args.pv[N1]) << args.qsv[N1+1])
 108                    + p01;
 109                 break;
 110                 }
 111         case 3: {
 112                 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
 113                 var cl_I p012 = p01 * args.pv[N1+2];
 114                 if (P) { *P = p012; }
 115                 var cl_I q12 = args.qv[N1+1] * args.qv[N1+2];
 116                 *Q = args.qv[N1] * q12;
 117                 *QS = args.qsv[N1] + args.qsv[N1+1] + args.qsv[N1+2];
 118                 *T = ((q12 * args.pv[N1]) << (args.qsv[N1+1] + args.qsv[N1+2]))
 119                    + ((args.qv[N1+2] * p01) << args.qsv[N1+2])
 120                    + p012;
 121                 break;
 122                 }
 123         case 4: {
 124                 var cl_I p01 = args.pv[N1] * args.pv[N1+1];
 125                 var cl_I p012 = p01 * args.pv[N1+2];
 126                 var cl_I p0123 = p012 * args.pv[N1+3];
 127                 if (P) { *P = p0123; }
 128                 var cl_I q23 = args.qv[N1+2] * args.qv[N1+3];
 129                 var cl_I q123 = args.qv[N1+1] * q23;
 130                 *Q = args.qv[N1] * q123;
 131                 *QS = args.qsv[N1] + args.qsv[N1+1] + args.qsv[N1+2] + args.qsv[N1+3];
 132                 *T = ((((((q123 * args.pv[N1]) << args.qsv[N1+1])
 133                          + q23 * p01) << args.qsv[N1+2])
 134                        + args.qv[N1+3] * p012) << args.qsv[N1+3])
 135                    + p0123;
 136                 break;
 137                 }
 138         default: {
 139                 var uintC Nm = (N1+N2)/2; // midpoint
 140                 // Compute left part.
 141                 var cl_I LP, LQ, LT;
 142                 var uintC LQS;
 143                 eval_pqs_series_aux(N1,Nm,args,&LP,&LQ,&LQS,&LT);
 144                 // Compute right part.
 145                 var cl_I RP, RQ, RT;
 146                 var uintC RQS;
 147                 eval_pqs_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RQS,&RT);
 148                 // Put together partial results.
 149                 if (P) { *P = LP*RP; }
 150                 *Q = LQ*RQ;
 151                 *QS = LQS+RQS;
 152                 // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT.
 153                 *T = ((RQ*LT) << RQS) + LP*RT;
 154                 break;
 155                 }
 156         }
 157 }
 158
 159 const cl_LF eval_rational_series (uintC N, const cl_pq_series& args, uintC len)
 160 {
 161         if (N==0)
 162                 return cl_I_to_LF(0,len);
 163         var cl_I Q, T;
 164         if (!args.qsv) {
 165                 eval_pq_series_aux(0,N,args,NULL,&Q,&T);
 166                 return cl_I_to_LF(T,len) / cl_I_to_LF(Q,len);
 167         } else {
 168                 // Precomputation of the shift counts:
 169                 // Split qv[n] into qv[n]*2^qsv[n].
 170                 {
 171                         var cl_I* qp = args.qv;
 172                         var uintC* qsp = args.qsv;
 173                         for (var uintC n = 0; n < N; n++, qp++, qsp++) {
 174                                 // Pull out maximal power of 2 out of *qp = args.qv[n].
 175                                 var uintC qs = 0;
 176                                 if (!zerop(*qp)) {
 177                                         qs = ord2(*qp);
 178                                         if (qs > 0)
 179                                                 *qp = *qp >> qs;
 180                                 }
 181                                 *qsp = qs;
 182                         }
 183                 }
 184                 // Main computation.
 185                 var uintC QS;
 186                 eval_pqs_series_aux(0,N,args,NULL,&Q,&QS,&T);
 187                 return cl_I_to_LF(T,len) / scale_float(cl_I_to_LF(Q,len),QS);
 188         }
 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 const cl_LF eval_rational_series (uintC N, cl_pq_series_stream& args, uintC len)
 266 {
 267         if (N==0)
 268                 return cl_I_to_LF(0,len);
 269         var cl_I Q, T;
 270         eval_pq_series_aux(0,N,args,NULL,&Q,&T);
 271         return cl_I_to_LF(T,len) / cl_I_to_LF(Q,len);
 272 }
 273
 274 static void eval_pq_series_aux (uintC N1, uintC N2,
 275                                 cl_pq_series_stream& args,
 276                                 cl_R* P, cl_R* Q, cl_R* T,
 277                                 uintC trunclen)
 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                 if (P) { *P = v0.p; }
 285                 *Q = v0.q;
 286                 *T = v0.p;
 287                 break;
 288                 }
 289         case 2: {
 290                 var cl_pq_series_term v0 = args.next(); // [N1]
 291                 var cl_pq_series_term v1 = args.next(); // [N1+1]
 292                 var cl_I p01 = v0.p * v1.p;
 293                 if (P) { *P = p01; }
 294                 *Q = v0.q * v1.q;
 295                 *T = v1.q * v0.p
 296                    + p01;
 297                 break;
 298                 }
 299         case 3: {
 300                 var cl_pq_series_term v0 = args.next(); // [N1]
 301                 var cl_pq_series_term v1 = args.next(); // [N1+1]
 302                 var cl_pq_series_term v2 = args.next(); // [N1+2]
 303                 var cl_I p01 = v0.p * v1.p;
 304                 var cl_I p012 = p01 * v2.p;
 305                 if (P) { *P = p012; }
 306                 var cl_I q12 = v1.q * v2.q;
 307                 *Q = v0.q * q12;
 308                 *T = q12 * v0.p
 309                    + v2.q * p01
 310                    + p012;
 311                 break;
 312                 }
 313         case 4: {
 314                 var cl_pq_series_term v0 = args.next(); // [N1]
 315                 var cl_pq_series_term v1 = args.next(); // [N1+1]
 316                 var cl_pq_series_term v2 = args.next(); // [N1+2]
 317                 var cl_pq_series_term v3 = args.next(); // [N1+3]
 318                 var cl_I p01 = v0.p * v1.p;
 319                 var cl_I p012 = p01 * v2.p;
 320                 var cl_I p0123 = p012 * v3.p;
 321                 if (P) { *P = p0123; }
 322                 var cl_I q23 = v2.q * v3.q;
 323                 var cl_I q123 = v1.q * q23;
 324                 *Q = v0.q * q123;
 325                 *T = q123 * v0.p
 326                    + q23 * p01
 327                    + v3.q * p012
 328                    + p0123;
 329                 break;
 330                 }
 331         default: {
 332                 var uintC Nm = (N1+N2)/2; // midpoint
 333                 // Compute left part.
 334                 var cl_R LP, LQ, LT;
 335                 eval_pq_series_aux(N1,Nm,args,&LP,&LQ,&LT,trunclen);
 336                 // Compute right part.
 337                 var cl_R RP, RQ, RT;
 338                 eval_pq_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RT,trunclen);
 339                 // Put together partial results.
 340                 if (P) {
 341                         *P = LP*RP;
 342                         truncate_precision(*P,trunclen);
 343                 }
 344                 *Q = LQ*RQ;
 345                 truncate_precision(*Q,trunclen);
 346                 // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT.
 347                 *T = RQ*LT + LP*RT;
 348                 truncate_precision(*T,trunclen);
 349                 break;
 350                 }
 351         }
 352 }
 353
 354 const cl_LF eval_rational_series (uintC N, cl_pq_series_stream& args, uintC len, uintC trunclen)
 355 {
 356         if (N==0)
 357                 return cl_I_to_LF(0,len);
 358         var cl_R Q, T;
 359         eval_pq_series_aux(0,N,args,NULL,&Q,&T,trunclen);
 360         return cl_R_to_LF(T,len) / cl_R_to_LF(Q,len);
 361 }
 362 // Bit complexity (if p(n), q(n), a(n), b(n) have length O(log(n))):
 363 // O(log(N)^2*M(N)).
 364
 365 }  // namespace cln