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