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