src/float/transcendental/cl_LF_ratseries_stream_pqab.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/abort.h"
  15 #include "cl_LF.h"
  16
  17 namespace cln {
  18
  19 static void eval_pqab_series_aux (uintL N1, uintL N2,
  20                                   cl_pqab_series_stream& args,
  21                                   cl_I* P, cl_I* Q, cl_I* B, cl_I* T)
  22 {
  23         switch (N2 - N1) {
  24         case 0:
  25                 cl_abort(); break;
  26         case 1: {
  27                 var cl_pqab_series_term v0 = args.next(); // [N1]
  28                 if (P) { *P = v0.p; }
  29                 *Q = v0.q;
  30                 *B = v0.b;
  31                 *T = v0.a * v0.p;
  32                 break;
  33                 }
  34         case 2: {
  35                 var cl_pqab_series_term v0 = args.next(); // [N1]
  36                 var cl_pqab_series_term v1 = args.next(); // [N1+1]
  37                 var cl_I p01 = v0.p * v1.p;
  38                 if (P) { *P = p01; }
  39                 *Q = v0.q * v1.q;
  40                 *B = v0.b * v1.b;
  41                 *T = v1.b * v1.q * v0.a * v0.p
  42                    + v0.b * v1.a * p01;
  43                 break;
  44                 }
  45         case 3: {
  46                 var cl_pqab_series_term v0 = args.next(); // [N1]
  47                 var cl_pqab_series_term v1 = args.next(); // [N1+1]
  48                 var cl_pqab_series_term v2 = args.next(); // [N1+2]
  49                 var cl_I p01 = v0.p * v1.p;
  50                 var cl_I p012 = p01 * v2.p;
  51                 if (P) { *P = p012; }
  52                 var cl_I q12 = v1.q * v2.q;
  53                 *Q = v0.q * q12;
  54                 var cl_I b12 = v1.b * v2.b;
  55                 *B = v0.b * b12;
  56                 *T = b12 * q12 * v0.a * v0.p
  57                    + v0.b * (v2.b * v2.q * v1.a * p01
  58                              + v1.b * v2.a * p012);
  59                 break;
  60                 }
  61         case 4: {
  62                 var cl_pqab_series_term v0 = args.next(); // [N1]
  63                 var cl_pqab_series_term v1 = args.next(); // [N1+1]
  64                 var cl_pqab_series_term v2 = args.next(); // [N1+2]
  65                 var cl_pqab_series_term v3 = args.next(); // [N1+3]
  66                 var cl_I p01 = v0.p * v1.p;
  67                 var cl_I p012 = p01 * v2.p;
  68                 var cl_I p0123 = p012 * v3.p;
  69                 if (P) { *P = p0123; }
  70                 var cl_I q23 = v2.q * v3.q;
  71                 var cl_I q123 = v1.q * q23;
  72                 *Q = v0.q * q123;
  73                 var cl_I b01 = v0.b * v1.b;
  74                 var cl_I b23 = v2.b * v3.b;
  75                 *B = b01 * b23;
  76                 *T = b23 * (v1.b * q123 * v0.a * v0.p
  77                             + v0.b * q23 * v1.a * p01)
  78                    + b01 * (v3.b * v3.q * v2.a * p012
  79                             + v2.b * v3.a * p0123);
  80                 break;
  81                 }
  82         default: {
  83                 var uintL Nm = (N1+N2)/2; // midpoint
  84                 // Compute left part.
  85                 var cl_I LP, LQ, LB, LT;
  86                 eval_pqab_series_aux(N1,Nm,args,&LP,&LQ,&LB,&LT);
  87                 // Compute right part.
  88                 var cl_I RP, RQ, RB, RT;
  89                 eval_pqab_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RB,&RT);
  90                 // Put together partial results.
  91                 if (P) { *P = LP*RP; }
  92                 *Q = LQ*RQ;
  93                 *B = LB*RB;
  94                 // S = LS + LP/LQ * RS, so T = RB*RQ*LT + LB*LP*RT.
  95                 *T = RB*RQ*LT + LB*LP*RT;
  96                 break;
  97                 }
  98         }
  99 }
 100
 101 const cl_LF eval_rational_series (uintL N, cl_pqab_series_stream& args, uintC len)
 102 {
 103         if (N==0)
 104                 return cl_I_to_LF(0,len);
 105         var cl_I Q, B, T;
 106         eval_pqab_series_aux(0,N,args,NULL,&Q,&B,&T);
 107         return cl_I_to_LF(T,len) / cl_I_to_LF(B*Q,len);
 108 }
 109
 110 }  // namespace cln