src/float/transcendental/cl_LF_tran.h

   1 // cl_LF internals, transcendental functions
   2
   3 #ifndef _CL_LF_TRAN_H
   4 #define _CL_LF_TRAN_H
   5
   6 #include "cln/integer.h"
   7 #include "cln/integer_ring.h"
   8 #include "cln/lfloat.h"
   9 #include "float/lfloat/cl_LF.h"
  10
  11 namespace cln {
  12
  13 // Subroutine for evaluating
  14 // sum(0 <= n < N, a(n)/b(n) * (p(0)...p(n))/(q(0)...q(n)))
  15 // where all the entries are small integers (ideally polynomials in n).
  16 // Some of the factors (a,b,p,q) may be omitted. They are then understood to
  17 // be 1. This is fast because it groups factors together before multiplying.
  18 // Result will be a cl_LF with len digits.
  19
  20 // There are various alternative implementations of the same algorithm that
  21 // differ in the way the series is represented and, as a consequence, in memory
  22 // consumption.
  23 //
  24 // 1st implementation (series is precomputed entirely)
  25 // Arguments:
  26 //   Vectors p[0..N-1], q[0..N-1], a[0..N-1], b[0..N-1], N.
  27 //   Some of the vectors (a,b,p,q) can be a NULL pointer, all of its entries
  28 //   are then understood to be 1.
  29 //   If given, a vector qs[0..N-1] which the evaluation routine may use to
  30 //   split off q[n] into q[n]*2^qs[n]. qs may be NULL, in that case no shift
  31 //   optimizations will be used. (They are worth it only if a significant
  32 //   amount of multiplication work can be saved by shifts.)
  33 //
  34 // 2nd implemenation (series is computed on demand, as a stream)
  35 // In this alternate implementation the series is not represented as a couple
  36 // of arrays, but as a method returning each tuple (p(n),q(n),a(n),b(n))
  37 // in turn. This is preferrable if the a(n) are big, in order to avoid too
  38 // much memory usage at the same time.
  39 // The next() function is called N times and is expected to return
  40 // (p(n),q(n),a(n),b(n)) for n=0..N-1 in that order.
  41 //
  42 // 3rd implemenation (series is computed on demand and truncated early)
  43 // This is like the second implementation, but it coerces the integer factors
  44 // to cl_LF of a given length (trunclen) as soon as the integer factor's size
  45 // exceeds the size to store the cl_LF. For this to make sense, trunclen must
  46 // not be smaller than len. In practice, this can shave off substantially from
  47 // the memory consumption but it also bears a potential for rounding errors.
  48 // A minimum trunclen that guarantees correctness must be evaluated on a
  49 // case-by-case basis.
  50 //
  51 // As a variation, it is sometimes advantageous to factor out from q[0..N-1]
  52 // powers of two and put them back in again at the end of the computation by
  53 // left-shifting. These are distinguished by a boolean template parameter.
  54 // (Some combinations might not be implemented yet.)
  55
  56
  57 // In each of the special cases below, none of (a,b,p,q) can be NULL.
  58
  59 struct cl_pqab_series {
  60         const cl_I* pv;
  61               cl_I* qv;
  62         const cl_I* av;
  63         const cl_I* bv;
  64 };
  65 struct cl_pqab_series_term {
  66         cl_I p;
  67         cl_I q;
  68         cl_I a;
  69         cl_I b;
  70 };
  71 struct cl_pqab_series_stream {
  72         cl_pqab_series_term (*nextfn)(cl_pqab_series_stream&);
  73         cl_pqab_series_term next () { return nextfn(*this); }
  74         // Constructor.
  75         cl_pqab_series_stream (cl_pqab_series_term (*n)(cl_pqab_series_stream&)) : nextfn (n) {}
  76 };
  77 template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_pqab_series& args, uintC len);
  78 template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_pqab_series_stream& args, uintC len);
  79 template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_pqab_series_stream& args, uintC len, uintC trunclen);
  80
  81 struct cl_pqb_series {
  82         const cl_I* pv;
  83               cl_I* qv;
  84         const cl_I* bv;
  85 };
  86 struct cl_pqb_series_term {
  87         cl_I p;
  88         cl_I q;
  89         cl_I b;
  90 };
  91 struct cl_pqb_series_stream {
  92         cl_pqb_series_term (*nextfn)(cl_pqb_series_stream&);
  93         cl_pqb_series_term next () { return nextfn(*this); }
  94         // Constructor.
  95         cl_pqb_series_stream (cl_pqb_series_term (*n)(cl_pqb_series_stream&)) : nextfn (n) {}
  96 };
  97 template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_pqb_series& args, uintC len);
  98 template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_pqb_series_stream& args, uintC len);
  99 template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_pqb_series_stream& args, uintC len, uintC trunclen);
 100
 101 struct cl_pqa_series {
 102         const cl_I* pv;
 103               cl_I* qv;
 104         const cl_I* av;
 105 };
 106 struct cl_pqa_series_term {
 107         cl_I p;
 108         cl_I q;
 109         cl_I a;
 110 };
 111 struct cl_pqa_series_stream {
 112         cl_pqa_series_term (*nextfn)(cl_pqa_series_stream&);
 113         cl_pqa_series_term next () { return nextfn(*this); }
 114         // Constructor.
 115         cl_pqa_series_stream (cl_pqa_series_term (*n)(cl_pqa_series_stream&)) : nextfn (n) {}
 116 };
 117 template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_pqa_series& args, uintC len);
 118 template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_pqa_series_stream& args, uintC len);
 119 template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_pqa_series_stream& args, uintC len, uintC trunclen);
 120
 121 struct cl_pq_series {
 122         const cl_I* pv;
 123               cl_I* qv;
 124 };
 125 struct cl_pq_series_term {
 126         cl_I p;
 127         cl_I q;
 128 };
 129 struct cl_pq_series_stream {
 130         cl_pq_series_term (*nextfn)(cl_pq_series_stream&);
 131         cl_pq_series_term next () { return nextfn(*this); }
 132         // Constructor.
 133         cl_pq_series_stream (cl_pq_series_term (*n)(cl_pq_series_stream&)) : nextfn (n) {}
 134 };
 135 template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_pq_series& args, uintC len);
 136 template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_pq_series_stream& args, uintC len);
 137 template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_pq_series_stream& args, uintC len, uintC trunclen);
 138
 139 struct cl_pab_series {
 140         const cl_I* pv;
 141         const cl_I* av;
 142         const cl_I* bv;
 143 };
 144 const cl_LF eval_rational_series (uintC N, const cl_pab_series& args, uintC len);
 145
 146 struct cl_pb_series {
 147         const cl_I* pv;
 148         const cl_I* bv;
 149 };
 150 const cl_LF eval_rational_series (uintC N, const cl_pb_series& args, uintC len);
 151
 152 struct cl_pa_series {
 153         const cl_I* pv;
 154         const cl_I* av;
 155 };
 156 const cl_LF eval_rational_series (uintC N, const cl_pa_series& args, uintC len);
 157
 158 struct cl_p_series {
 159         const cl_I* pv;
 160 };
 161 const cl_LF eval_rational_series (uintC N, const cl_p_series& args, uintC len);
 162
 163 struct cl_qab_series {
 164               cl_I* qv;
 165         const cl_I* av;
 166         const cl_I* bv;
 167 };
 168 template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_qab_series& args, uintC len);
 169
 170 struct cl_qb_series {
 171               cl_I* qv;
 172         const cl_I* bv;
 173 };
 174 struct cl_qb_series_term {
 175         cl_I q;
 176         cl_I b;
 177 };
 178 struct cl_qb_series_stream {
 179         cl_qb_series_term (*nextfn)(cl_qb_series_stream&);
 180         cl_qb_series_term next () { return nextfn(*this); }
 181         // Constructor.
 182         cl_qb_series_stream (cl_qb_series_term (*n)(cl_qb_series_stream&)) : nextfn (n) {}
 183 };
 184 template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_qb_series& args, uintC len);
 185 template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_qb_series_stream& args, uintC len);
 186
 187 struct cl_qa_series {
 188               cl_I* qv;
 189         const cl_I* av;
 190 };
 191 template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_qa_series& args, uintC len);
 192
 193 struct cl_q_series {
 194               cl_I* qv;
 195 };
 196 struct cl_q_series_term {
 197         cl_I q;
 198 };
 199 struct cl_q_series_stream {
 200         cl_q_series_term (*nextfn)(cl_q_series_stream&);
 201         cl_q_series_term next () { return nextfn(*this); }
 202         // Constructor.
 203         cl_q_series_stream (cl_q_series_term (*n)(cl_q_series_stream&)) : nextfn (n) {}
 204 };
 205 template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_q_series& args, uintC len);
 206 template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_q_series_stream& args, uintC len);
 207
 208 struct cl_ab_series {
 209         const cl_I* av;
 210         const cl_I* bv;
 211 };
 212 const cl_LF eval_rational_series (uintC N, const cl_ab_series& args, uintC len);
 213
 214 struct cl_b_series {
 215         const cl_I* bv;
 216 };
 217 const cl_LF eval_rational_series (uintC N, const cl_b_series& args, uintC len);
 218
 219 struct cl_a_series {
 220         const cl_I* av;
 221 };
 222 const cl_LF eval_rational_series (uintC N, const cl_a_series& args, uintC len);
 223
 224 struct cl__series {
 225 };
 226 const cl_LF eval_rational_series (uintC N, const cl__series& args, uintC len);
 227
 228
 229 // [Generalization.]
 230 // Subroutine:
 231 // Evaluates S = sum(N1 <= n < N2, (p(N1)...p(n))/(q(N1)...q(n)))
 232 // and U = sum(N1 <= n < N2,
 233 //             (c(N1)/d(N1)+...+c(n)/d(n))*(p(N1)...p(n))/(q(N1)...q(n)))
 234 // and returns
 235 //     P = p(N1)...p(N2-1),
 236 //     Q = q(N1)...q(N2-1),
 237 //     T = Q*S,
 238 //     C/D = c(N1)/d(N1)+...+c(N2-1)/d(N2-1),
 239 //     V = D*Q*U,
 240 // all integers. On entry N1 < N2.
 241 struct cl_pqcd_series_term {
 242         cl_I p;
 243         cl_I q;
 244         cl_I c;
 245         cl_I d;
 246 };
 247 template<class cl_T>
 248 struct cl_pqcd_series_result {
 249         cl_T P;
 250         cl_T Q;
 251         cl_T T;
 252         cl_T C;
 253         cl_T D;
 254         cl_T V;
 255 };
 256 struct cl_pqcd_series_stream {
 257         cl_pqcd_series_term (*nextfn)(cl_pqcd_series_stream&);
 258         cl_pqcd_series_term next () { return nextfn(*this); }
 259         // Constructor.
 260         cl_pqcd_series_stream( cl_pqcd_series_term (*n)(cl_pqcd_series_stream&)) : nextfn (n) {}
 261 };
 262 void eval_pqcd_series_aux (uintC N, cl_pqcd_series_term* args, cl_pqcd_series_result<cl_I>& Z, bool rightmost = true);
 263 void eval_pqcd_series_aux (uintC N, cl_pqcd_series_stream& args, cl_pqcd_series_result<cl_I>& Z, bool rightmost = true);
 264 void eval_pqcd_series_aux (uintC N, cl_pqcd_series_stream& args, cl_pqcd_series_result<cl_R>& Z, uintC trunclen, bool rightmost = true);
 265 // Ditto, but returns U/S.
 266 const cl_LF eval_pqcd_series (uintC N, cl_pqcd_series_term* args, uintC len);
 267 const cl_LF eval_pqcd_series (uintC N, cl_pqcd_series_stream& args, uintC len);
 268 const cl_LF eval_pqcd_series (uintC N, cl_pqcd_series_stream& args, uintC len, uintC trunclen);
 269
 270 // [Special case c(n)=1.]
 271 // Subroutine:
 272 // Evaluates S = sum(N1 <= n < N2, (p(N1)...p(n))/(q(N1)...q(n)))
 273 // and U = sum(N1 <= n < N2, (1/d(N1)+...+1/d(n))*(p(N1)...p(n))/(q(N1)...q(n)))
 274 // and returns
 275 //     P = p(N1)...p(N2-1),
 276 //     Q = q(N1)...q(N2-1),
 277 //     T = Q*S,
 278 //     C/D = 1/d(N1)+...+1/d(N2-1),
 279 //     V = D*Q*U,
 280 // all integers. On entry N1 < N2.
 281 struct cl_pqd_series_term {
 282         cl_I p;
 283         cl_I q;
 284         cl_I d;
 285 };
 286 template<class cl_T>
 287 struct cl_pqd_series_result {
 288         cl_T P;
 289         cl_T Q;
 290         cl_T T;
 291         cl_T C;
 292         cl_T D;
 293         cl_T V;
 294 };
 295 struct cl_pqd_series_stream {
 296         cl_pqd_series_term (*nextfn)(cl_pqd_series_stream&);
 297         cl_pqd_series_term next () { return nextfn(*this); }
 298         // Constructor.
 299         cl_pqd_series_stream( cl_pqd_series_term (*n)(cl_pqd_series_stream&)) : nextfn (n) {}
 300 };
 301 void eval_pqd_series_aux (uintC N, cl_pqd_series_term* args, cl_pqd_series_result<cl_I>& Z, bool rightmost = true);
 302 void eval_pqd_series_aux (uintC N, cl_pqd_series_stream& args, cl_pqd_series_result<cl_I>& Z, bool rightmost = true);
 303 void eval_pqd_series_aux (uintC N, cl_pqd_series_stream& args, cl_pqd_series_result<cl_R>& Z, uintC trunclen, bool rightmost = true);
 304 // Ditto, but returns U/S.
 305 const cl_LF eval_pqd_series (uintC N, cl_pqd_series_term* args, uintC len);
 306 const cl_LF eval_pqd_series (uintC N, cl_pqd_series_stream& args, uintC len);
 307 const cl_LF eval_pqd_series (uintC N, cl_pqd_series_stream& args, uintC len, uintC trunclen);
 308
 309 // Helper function to divide q by the largest s such that 2^s divides q, returns s.
 310 inline uintC
 311 pullout_shiftcount(cl_I& q)
 312 {
 313         var uintC qs = 0;
 314         if (!zerop(q)) {
 315                 qs = ord2(q);
 316                 if (qs > 0)
 317                         q = q >> qs;
 318         }
 319         return qs;
 320 }
 321
 322 // Helper function to convert integer of length > trunclen to long float of
 323 // length = trunclen.
 324 inline void
 325 truncate_precision(cl_R& x, uintC trunclen)
 326 {
 327         if (instanceof(x,cl_I_ring) &&
 328             integer_length(the<cl_I>(x))>trunclen*intDsize) {
 329                 x = cl_I_to_LF(the<cl_I>(x),trunclen);
 330         }
 331 }
 332
 333 }  // namespace cln
 334
 335 #endif /* _CL_LF_TRAN_H */