- std::vector<cln::cl_N> m;
- const int nops = ex_to<numeric>(x1.nops()).to_int();
- for (int i=nops-2; i>=0; i--) {
- m.push_back(ex_to<numeric>(x1.op(i)).to_cl_N());
- }
-
- cln::float_format_t prec = cln::default_float_format;
- cln::cl_N res = cln::complex(cln::cl_float(0, prec), 0);
- cln::cl_N resbuf;
- for (int i=nops; true; i++) {
- // to infinity and beyond ... timewise
- resbuf = res;
- res = res + cln::recip(cln::expt(i,m_1)) * numeric_harmonic(i, m);
- if (cln::zerop(res-resbuf))
- break;
- }
-
- return numeric(res);
+ // buffer for subsums
+ std::vector<cln::cl_N> t(j);
+ cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
+
+ cln::cl_N t0buf;
+ int q = 0;
+ do {
+ t0buf = t[0];
+ q++;
+ t[j-1] = t[j-1] + one / cln::expt(cln::cl_I(q),r[j-1]);
+ for (int k=j-2; k>=0; k--) {
+ t[k] = t[k] + one * t[k+1] / cln::expt(cln::cl_I(q+j-1-k), r[k]);
+ }
+ } while (t[0] != t0buf);