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