- cl_RA tmp; // used to store the sum
- for (int i=highest_result+1; i<=n; ++i) {
- // the first two elements:
- tmp = cl_I(-2*i-1)/cl_I(2);
- // accumulate the remaining elements:
- for (int j=0; j<i; ++j)
- tmp = tmp + ::binomial(2*i+3,j*2+2)*results[j];
- // divide by -(nn+1) and store result:
- results.push_back(tmp/cl_I(-2*i-3));
+ // let CLN set up the integer ring we are working in
+ ::cl_univpoly_integer_ring myring = ::cl_find_univpoly_ring(::cl_I_ring);
+
+ // store nonvanishing Bernoulli numbers and the last tangent poly here
+ static vector< ::cl_RA > results;
+ static int highest_result = 0;
+ static ::cl_UP_I T_n = myring->create(1);
+
+ // index of results vector
+ int m = (nn.to_int()-4)/2;
+ // look if the result is precomputed
+ if (m < highest_result)
+ return numeric(results[m]);
+ // reserve the whole chunk we'll need
+ if (results.capacity() < (unsigned)(m+1))
+ results.reserve(m+1);
+
+ // create the generating polynomial T_gen = 1+x^2
+ ::cl_UP_I T_gen = myring->create(2);
+ ::set_coeff(T_gen, 0, cl_I(1));
+ ::set_coeff(T_gen, 2, cl_I(1));
+ ::finalize(T_gen);
+ // T_n will be iterated, start with T_1(x) == 1+x^2 == T_gen
+ if (highest_result==0)
+ T_n = T_gen;
+ // iterate to the n'th tangent polynomial
+ for (int n=highest_result; n<=m; ++n) {
+ // recur tangent polynomial twice
+ for (int i=0; i<2; ++i)
+ T_n = ::deriv(T_n)*T_gen;
+ // fetch T_{n-1}
+ ::cl_I T = ::coeff(T_n,0);
+ // compute B_n from the T_{n-1}:
+ ::cl_RA B = T * (n%2 ? ::cl_I(1) : ::cl_I(-1));
+ B = B * 2*(n+2) / ((::cl_I(1)<<4*(n+2))-(::cl_I(1)<<2*(n+2)));
+ results.push_back(B);
+ ++highest_result;