src/float/transcendental/cl_LF_cossin.cc

   1 // cl_cossin_ratseries().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_F_tran.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cln/lfloat.h"
  13 #include "cl_LF.h"
  14 #include "cln/integer.h"
  15
  16 namespace cln {
  17
  18 inline const cl_LF_cos_sin_t operator* (const cl_LF_cos_sin_t& a, const cl_LF_cos_sin_t& b)
  19 {
  20         return cl_LF_cos_sin_t(a.cos*b.cos-a.sin*b.sin,a.sin*b.cos+a.cos*b.sin);
  21 }
  22
  23 const cl_LF_cos_sin_t cl_cossin_ratseries (const cl_LF& x)
  24 {
  25         // Similar to expx_ratseries.
  26         var uintC len = TheLfloat(x)->len;
  27         var cl_idecoded_float x_ = integer_decode_float(x);
  28         // x = (-1)^sign * 2^exponent * mantissa
  29         var uintL lq = cl_I_to_UL(- x_.exponent);
  30         var const cl_I& p = x_.mantissa;
  31         // Compute sin(p/2^lq) and cos(p/2^lq) by splitting into pieces.
  32         var cl_boolean first_factor = cl_true;
  33         var cl_LF_cos_sin_t product;
  34         var uintL b1;
  35         var uintL b2;
  36         for (b1 = 0, b2 = 1; b1 < lq; b1 = b2, b2 = 2*b2) {
  37                 // Piece containing bits b1+1..b2 after "decimal point"
  38                 // in the binary representation of (p/2^lq).
  39                 var uintL lqk = (lq >= b2 ? b2 : lq);
  40                 var cl_I pk = ldb(p,cl_byte(lqk-b1,lq-lqk));
  41                 // Compute sin(pk/2^lqk) and cos(pk/2^lqk).
  42                 if (!zerop(pk)) {
  43                         if (minusp(x_.sign)) { pk = -pk; }
  44                         var cl_LF_cos_sin_t factor = cl_cossin_aux(pk,lqk,len);
  45                         if (first_factor) {
  46                                 product = factor;
  47                                 first_factor = cl_false;
  48                         } else
  49                                 product = product * factor;
  50                 }
  51         }
  52         if (first_factor)
  53                 return cl_LF_cos_sin_t(cl_I_to_LF(1,len),cl_I_to_LF(0,len));
  54         else
  55                 return product;
  56 }
  57 // Bit complexity (N = length(x)): O(log(N)^2*M(N)).
  58
  59 }  // namespace cln