← Back to branch summary

~ubuntu-branches/ubuntu/saucy/deal.ii/saucy

« back to all changes in this revision

Viewing changes to contrib/boost/include/boost/math/special_functions/zeta.hpp

  • Committer: Package Import Robot
  • Author(s): "Adam C. Powell, IV", Adam C. Powell, IV, Christophe Trophime
  • Date: 2012-02-21 06:57:30 UTC
  • mfrom: (3.1.7 sid)
  • Revision ID: package-import@ubuntu.com-20120221065730-r2iz70lg557wcd2e
Tags: 7.1.0-1
http://bugs.debian.org/652057
http://bugs.debian.org/652058
http://bugs.debian.org/637655
http://bugs.debian.org/659245
[ Adam C. Powell, IV ]
* New upstream (closes: #652057).
* Updated to use PETSc and SLEPc 3.2, and forward-ported all patches.
* Removed NetCDF Build-Depends because it uses serial HDF5.
* Made Sacado cmath patch work with new configure.
* Moved -dev package symlink lines in rules to arch all section.

[ Christophe Trophime ]
* debian/rules:
   - add dh_strip --dbg-package to generate dbg package (closes: #652058)
   - add .install files to simplify rules
* Add support for mumps, arpack (closes: #637655)
* Add patch for slepc 3.2 (closes: #659245)

Show diffs side-by-side

added added

removed removed

Lines of Context:
contrib/boost/include/boost/math/special_functions/zeta.hpp
1
 
//  Copyright John Maddock 2007.
2
 
//  Use, modification and distribution are subject to the
3
 
//  Boost Software License, Version 1.0. (See accompanying file
4
 
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5
 
 
6
 
#ifndef BOOST_MATH_ZETA_HPP
7
 
#define BOOST_MATH_ZETA_HPP
8
 
 
9
 
#ifdef _MSC_VER
10
 
#pragma once
11
 
#endif
12
 
 
13
 
#include <boost/math/tools/precision.hpp>
14
 
#include <boost/math/tools/series.hpp>
15
 
#include <boost/math/policies/error_handling.hpp>
16
 
#include <boost/math/special_functions/gamma.hpp>
17
 
#include <boost/math/special_functions/sin_pi.hpp>
18
 
 
19
 
namespace boost{ namespace math{ namespace detail{
20
 
 
21
 
template <class T, class Policy>
22
 
struct zeta_series_cache_size
23
 
{
24
 
   //
25
 
   // Work how large to make our cache size when evaluating the series 
26
 
   // evaluation:  normally this is just large enough for the series
27
 
   // to have converged, but for arbitrary precision types we need a 
28
 
   // really large cache to achieve reasonable precision in a reasonable 
29
 
   // time.  This is important when constructing rational approximations
30
 
   // to zeta for example.
31
 
   //
32
 
   typedef typename boost::math::policies::precision<T,Policy>::type precision_type;
33
 
   typedef typename mpl::if_<
34
 
      mpl::less_equal<precision_type, mpl::int_<0> >,
35
 
      mpl::int_<5000>,
36
 
      typename mpl::if_<
37
 
         mpl::less_equal<precision_type, mpl::int_<64> >,
38
 
         mpl::int_<70>,
39
 
         typename mpl::if_<
40
 
            mpl::less_equal<precision_type, mpl::int_<113> >,
41
 
            mpl::int_<100>,
42
 
            mpl::int_<5000>
43
 
         >::type
44
 
      >::type
45
 
   >::type type;
46
 
};
47
 
 
48
 
template <class T, class Policy>
49
 
T zeta_series_imp(T s, T sc, const Policy&)
50
 
{
51
 
   //
52
 
   // Series evaluation from:
53
 
   // Havil, J. Gamma: Exploring Euler's Constant. 
54
 
   // Princeton, NJ: Princeton University Press, 2003.
55
 
   //
56
 
   // See also http://mathworld.wolfram.com/RiemannZetaFunction.html
57
 
   //
58
 
   BOOST_MATH_STD_USING
59
 
   T sum = 0;
60
 
   T mult = 0.5;
61
 
   T change;
62
 
   typedef typename zeta_series_cache_size<T,Policy>::type cache_size;
63
 
   T powers[cache_size::value] = { 0, };
64
 
   unsigned n = 0;
65
 
   do{
66
 
      T binom = -static_cast<T>(n);
67
 
      T nested_sum = 1;
68
 
      if(n < sizeof(powers) / sizeof(powers[0]))
69
 
         powers[n] = pow(static_cast<T>(n + 1), -s);
70
 
      for(unsigned k = 1; k <= n; ++k)
71
 
      {
72
 
         T p;
73
 
         if(k < sizeof(powers) / sizeof(powers[0]))
74
 
         {
75
 
            p = powers[k];
76
 
            //p = pow(k + 1, -s);
77
 
         }
78
 
         else
79
 
            p = pow(static_cast<T>(k + 1), -s);
80
 
         nested_sum += binom * p;
81
 
        binom *= (k - static_cast<T>(n)) / (k + 1);
82
 
      }
83
 
      change = mult * nested_sum;
84
 
      sum += change;
85
 
      mult /= 2;
86
 
      ++n;
87
 
   }while(fabs(change / sum) > tools::epsilon<T>());
88
 
 
89
 
   return sum * 1 / -boost::math::powm1(T(2), sc);
90
 
}
91
 
//
92
 
// Classical p-series:
93
 
//
94
 
template <class T>
95
 
struct zeta_series2
96
 
{
97
 
   typedef T result_type;
98
 
   zeta_series2(T _s) : s(-_s), k(1){}
99
 
   T operator()()
100
 
   {
101
 
      BOOST_MATH_STD_USING
102
 
      return pow(static_cast<T>(k++), s);
103
 
   }
104
 
private:
105
 
   T s;
106
 
   unsigned k;
107
 
};
108
 
 
109
 
template <class T, class Policy>
110
 
inline T zeta_series2_imp(T s, const Policy& pol)
111
 
{
112
 
   boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();;
113
 
   zeta_series2<T> f(s);
114
 
   T result = tools::sum_series(
115
 
      f, 
116
 
      policies::get_epsilon<T, Policy>(),
117
 
      max_iter);
118
 
   policies::check_series_iterations("boost::math::zeta_series2<%1%>(%1%)", max_iter, pol);
119
 
   return result;
120
 
}
121
 
 
122
 
template <class T, class Policy>
123
 
T zeta_imp_prec(T s, T sc, const Policy& pol, const mpl::int_<0>&)
124
 
{
125
 
   BOOST_MATH_STD_USING
126
 
   T result; 
127
 
   //
128
 
   // Only use power series if it will converge in 100 
129
 
   // iterations or less: the more iterations it consumes
130
 
   // the slower convergence becomes so we have to be very 
131
 
   // careful in it's usage.
132
 
   //
133
 
   if (s > -log(tools::epsilon<T>()) / 4.5)
134
 
      result = detail::zeta_series2_imp(s, pol);
135
 
   else
136
 
      result = detail::zeta_series_imp(s, sc, pol);
137
 
   return result;
138
 
}
139
 
 
140
 
template <class T, class Policy>
141
 
inline T zeta_imp_prec(T s, T sc, const Policy&, const mpl::int_<53>&)
142
 
{
143
 
   BOOST_MATH_STD_USING
144
 
   T result;
145
 
   if(s < 1)
146
 
   {
147
 
      // Rational Approximation
148
 
      // Maximum Deviation Found:                     2.020e-18
149
 
      // Expected Error Term:                         -2.020e-18
150
 
      // Max error found at double precision:         3.994987e-17
151
 
      static const T P[6] = {    
152
 
         0.24339294433593750202L,
153
 
         -0.49092470516353571651L,
154
 
         0.0557616214776046784287L,
155
 
         -0.00320912498879085894856L,
156
 
         0.000451534528645796438704L,
157
 
         -0.933241270357061460782e-5L,
158
 
        };
159
 
      static const T Q[6] = {    
160
 
         1L,
161
 
         -0.279960334310344432495L,
162
 
         0.0419676223309986037706L,
163
 
         -0.00413421406552171059003L,
164
 
         0.00024978985622317935355L,
165
 
         -0.101855788418564031874e-4L,
166
 
      };
167
 
      result = tools::evaluate_polynomial(P, sc) / tools::evaluate_polynomial(Q, sc);
168
 
      result -= 1.2433929443359375F;
169
 
      result += (sc);
170
 
      result /= (sc);
171
 
   }
172
 
   else if(s <= 2)
173
 
   {
174
 
      // Maximum Deviation Found:        9.007e-20
175
 
      // Expected Error Term:            9.007e-20
176
 
      static const T P[6] = {    
177
 
         0.577215664901532860516,
178
 
         0.243210646940107164097,
179
 
         0.0417364673988216497593,
180
 
         0.00390252087072843288378,
181
 
         0.000249606367151877175456,
182
 
         0.110108440976732897969e-4,
183
 
      };
184
 
      static const T Q[6] = {    
185
 
         1,
186
 
         0.295201277126631761737,
187
 
         0.043460910607305495864,
188
 
         0.00434930582085826330659,
189
 
         0.000255784226140488490982,
190
 
         0.10991819782396112081e-4,
191
 
      };
192
 
      result = tools::evaluate_polynomial(P, -sc) / tools::evaluate_polynomial(Q, -sc);
193
 
      result += 1 / (-sc);
194
 
   }
195
 
   else if(s <= 4)
196
 
   {
197
 
      // Maximum Deviation Found:          5.946e-22
198
 
      // Expected Error Term:              -5.946e-22
199
 
      static const float Y = 0.6986598968505859375;
200
 
      static const T P[6] = {    
201
 
         -0.0537258300023595030676,
202
 
         0.0445163473292365591906,
203
 
         0.0128677673534519952905,
204
 
         0.00097541770457391752726,
205
 
         0.769875101573654070925e-4,
206
 
         0.328032510000383084155e-5,
207
 
      };
208
 
      static const T Q[7] = {    
209
 
         1,
210
 
         0.33383194553034051422,
211
 
         0.0487798431291407621462,
212
 
         0.00479039708573558490716,
213
 
         0.000270776703956336357707,
214
 
         0.106951867532057341359e-4,
215
 
         0.236276623974978646399e-7,
216
 
      };
217
 
      result = tools::evaluate_polynomial(P, s - 2) / tools::evaluate_polynomial(Q, s - 2);
218
 
      result += Y + 1 / (-sc);
219
 
   }
220
 
   else if(s <= 7)
221
 
   {
222
 
      // Maximum Deviation Found:                     2.955e-17
223
 
      // Expected Error Term:                         2.955e-17
224
 
      // Max error found at double precision:         2.009135e-16
225
 
 
226
 
      static const T P[6] = {    
227
 
         -2.49710190602259410021,
228
 
         -2.60013301809475665334,
229
 
         -0.939260435377109939261,
230
 
         -0.138448617995741530935,
231
 
         -0.00701721240549802377623,
232
 
         -0.229257310594893932383e-4,
233
 
      };
234
 
      static const T Q[9] = {    
235
 
         1,
236
 
         0.706039025937745133628,
237
 
         0.15739599649558626358,
238
 
         0.0106117950976845084417,
239
 
         -0.36910273311764618902e-4,
240
 
         0.493409563927590008943e-5,
241
 
         -0.234055487025287216506e-6,
242
 
         0.718833729365459760664e-8,
243
 
         -0.1129200113474947419e-9,
244
 
      };
245
 
      result = tools::evaluate_polynomial(P, s - 4) / tools::evaluate_polynomial(Q, s - 4);
246
 
      result = 1 + exp(result);
247
 
   }
248
 
   else if(s < 15)
249
 
   {
250
 
      // Maximum Deviation Found:                     7.117e-16
251
 
      // Expected Error Term:                         7.117e-16
252
 
      // Max error found at double precision:         9.387771e-16
253
 
      static const T P[7] = {    
254
 
         -4.78558028495135619286,
255
 
         -1.89197364881972536382,
256
 
         -0.211407134874412820099,
257
 
         -0.000189204758260076688518,
258
 
         0.00115140923889178742086,
259
 
         0.639949204213164496988e-4,
260
 
         0.139348932445324888343e-5,
261
 
        };
262
 
      static const T Q[9] = {    
263
 
         1,
264
 
         0.244345337378188557777,
265
 
         0.00873370754492288653669,
266
 
         -0.00117592765334434471562,
267
 
         -0.743743682899933180415e-4,
268
 
         -0.21750464515767984778e-5,
269
 
         0.471001264003076486547e-8,
270
 
         -0.833378440625385520576e-10,
271
 
         0.699841545204845636531e-12,
272
 
        };
273
 
      result = tools::evaluate_polynomial(P, s - 7) / tools::evaluate_polynomial(Q, s - 7);
274
 
      result = 1 + exp(result);
275
 
   }
276
 
   else if(s < 36)
277
 
   {
278
 
      // Max error in interpolated form:             1.668e-17
279
 
      // Max error found at long double precision:   1.669714e-17
280
 
      static const T P[8] = {    
281
 
         -10.3948950573308896825,
282
 
         -2.85827219671106697179,
283
 
         -0.347728266539245787271,
284
 
         -0.0251156064655346341766,
285
 
         -0.00119459173416968685689,
286
 
         -0.382529323507967522614e-4,
287
 
         -0.785523633796723466968e-6,
288
 
         -0.821465709095465524192e-8,
289
 
      };
290
 
      static const T Q[10] = {    
291
 
         1,
292
 
         0.208196333572671890965,
293
 
         0.0195687657317205033485,
294
 
         0.00111079638102485921877,
295
 
         0.408507746266039256231e-4,
296
 
         0.955561123065693483991e-6,
297
 
         0.118507153474022900583e-7,
298
 
         0.222609483627352615142e-14,
299
 
      };
300
 
      result = tools::evaluate_polynomial(P, s - 15) / tools::evaluate_polynomial(Q, s - 15);
301
 
      result = 1 + exp(result);
302
 
   }
303
 
   else if(s < 56)
304
 
   {
305
 
      result = 1 + pow(T(2), -s);
306
 
   }
307
 
   else
308
 
   {
309
 
      result = 1;
310
 
   }
311
 
   return result;
312
 
}
313
 
 
314
 
template <class T, class Policy>
315
 
T zeta_imp_prec(T s, T sc, const Policy&, const mpl::int_<64>&)
316
 
{
317
 
   BOOST_MATH_STD_USING
318
 
   T result;
319
 
   if(s < 1)
320
 
   {
321
 
      // Rational Approximation
322
 
      // Maximum Deviation Found:                     3.099e-20
323
 
      // Expected Error Term:                         3.099e-20
324
 
      // Max error found at long double precision:    5.890498e-20
325
 
      static const T P[6] = {    
326
 
         0.243392944335937499969L,
327
 
         -0.496837806864865688082L,
328
 
         0.0680008039723709987107L,
329
 
         -0.00511620413006619942112L,
330
 
         0.000455369899250053003335L,
331
 
         -0.279496685273033761927e-4L,
332
 
        };
333
 
      static const T Q[7] = {    
334
 
         1L,
335
 
         -0.30425480068225790522L,
336
 
         0.050052748580371598736L,
337
 
         -0.00519355671064700627862L,
338
 
         0.000360623385771198350257L,
339
 
         -0.159600883054550987633e-4L,
340
 
         0.339770279812410586032e-6L,
341
 
      };
342
 
      result = tools::evaluate_polynomial(P, sc) / tools::evaluate_polynomial(Q, sc);
343
 
      result -= 1.2433929443359375F;
344
 
      result += (sc);
345
 
      result /= (sc);
346
 
   }
347
 
   else if(s <= 2)
348
 
   {
349
 
      // Maximum Deviation Found:                     1.059e-21
350
 
      // Expected Error Term:                         1.059e-21
351
 
      // Max error found at long double precision:    1.626303e-19
352
 
 
353
 
      static const T P[6] = {    
354
 
         0.577215664901532860605L,
355
 
         0.222537368917162139445L,
356
 
         0.0356286324033215682729L,
357
 
         0.00304465292366350081446L,
358
 
         0.000178102511649069421904L,
359
 
         0.700867470265983665042e-5L,
360
 
      };
361
 
      static const T Q[7] = {    
362
 
         1L,
363
 
         0.259385759149531030085L,
364
 
         0.0373974962106091316854L,
365
 
         0.00332735159183332820617L,
366
 
         0.000188690420706998606469L,
367
 
         0.635994377921861930071e-5L,
368
 
         0.226583954978371199405e-7L,
369
 
      };
370
 
      result = tools::evaluate_polynomial(P, -sc) / tools::evaluate_polynomial(Q, -sc);
371
 
      result += 1 / (-sc);
372
 
   }
373
 
   else if(s <= 4)
374
 
   {
375
 
      // Maximum Deviation Found:          5.946e-22
376
 
      // Expected Error Term:              -5.946e-22
377
 
      static const float Y = 0.6986598968505859375;
378
 
      static const T P[7] = {    
379
 
         -0.053725830002359501027L,
380
 
         0.0470551187571475844778L,
381
 
         0.0101339410415759517471L,
382
 
         0.00100240326666092854528L,
383
 
         0.685027119098122814867e-4L,
384
 
         0.390972820219765942117e-5L,
385
 
         0.540319769113543934483e-7L,
386
 
      };
387
 
      static const T Q[8] = {    
388
 
         1,
389
 
         0.286577739726542730421L,
390
 
         0.0447355811517733225843L,
391
 
         0.00430125107610252363302L,
392
 
         0.000284956969089786662045L,
393
 
         0.116188101609848411329e-4L,
394
 
         0.278090318191657278204e-6L,
395
 
         -0.19683620233222028478e-8L,
396
 
      };
397
 
      result = tools::evaluate_polynomial(P, s - 2) / tools::evaluate_polynomial(Q, s - 2);
398
 
      result += Y + 1 / (-sc);
399
 
   }
400
 
   else if(s <= 7)
401
 
   {
402
 
      // Max error found at long double precision: 8.132216e-19
403
 
      static const T P[8] = {    
404
 
         -2.49710190602259407065L,
405
 
         -3.36664913245960625334L,
406
 
         -1.77180020623777595452L,
407
 
         -0.464717885249654313933L,
408
 
         -0.0643694921293579472583L,
409
 
         -0.00464265386202805715487L,
410
 
         -0.000165556579779704340166L,
411
 
         -0.252884970740994069582e-5L,
412
 
      };
413
 
      static const T Q[9] = {    
414
 
         1,
415
 
         1.01300131390690459085L,
416
 
         0.387898115758643503827L,
417
 
         0.0695071490045701135188L,
418
 
         0.00586908595251442839291L,
419
 
         0.000217752974064612188616L,
420
 
         0.397626583349419011731e-5L,
421
 
         -0.927884739284359700764e-8L,
422
 
         0.119810501805618894381e-9L,
423
 
      };
424
 
      result = tools::evaluate_polynomial(P, s - 4) / tools::evaluate_polynomial(Q, s - 4);
425
 
      result = 1 + exp(result);
426
 
   }
427
 
   else if(s < 15)
428
 
   {
429
 
      // Max error in interpolated form:              1.133e-18
430
 
      // Max error found at long double precision:    2.183198e-18
431
 
      static const T P[9] = {    
432
 
         -4.78558028495135548083L,
433
 
         -3.23873322238609358947L,
434
 
         -0.892338582881021799922L,
435
 
         -0.131326296217965913809L,
436
 
         -0.0115651591773783712996L,
437
 
         -0.000657728968362695775205L,
438
 
         -0.252051328129449973047e-4L,
439
 
         -0.626503445372641798925e-6L,
440
 
         -0.815696314790853893484e-8L,
441
 
        };
442
 
      static const T Q[9] = {    
443
 
         1,
444
 
         0.525765665400123515036L,
445
 
         0.10852641753657122787L,
446
 
         0.0115669945375362045249L,
447
 
         0.000732896513858274091966L,
448
 
         0.30683952282420248448e-4L,
449
 
         0.819649214609633126119e-6L,
450
 
         0.117957556472335968146e-7L,
451
 
         -0.193432300973017671137e-12L,
452
 
        };
453
 
      result = tools::evaluate_polynomial(P, s - 7) / tools::evaluate_polynomial(Q, s - 7);
454
 
      result = 1 + exp(result);
455
 
   }
456
 
   else if(s < 42)
457
 
   {
458
 
      // Max error in interpolated form:             1.668e-17
459
 
      // Max error found at long double precision:   1.669714e-17
460
 
      static const T P[9] = {    
461
 
         -10.3948950573308861781L,
462
 
         -2.82646012777913950108L,
463
 
         -0.342144362739570333665L,
464
 
         -0.0249285145498722647472L,
465
 
         -0.00122493108848097114118L,
466
 
         -0.423055371192592850196e-4L,
467
 
         -0.1025215577185967488e-5L,
468
 
         -0.165096762663509467061e-7L,
469
 
         -0.145392555873022044329e-9L,
470
 
      };
471
 
      static const T Q[10] = {    
472
 
         1,
473
 
         0.205135978585281988052L,
474
 
         0.0192359357875879453602L,
475
 
         0.00111496452029715514119L,
476
 
         0.434928449016693986857e-4L,
477
 
         0.116911068726610725891e-5L,
478
 
         0.206704342290235237475e-7L,
479
 
         0.209772836100827647474e-9L,
480
 
         -0.939798249922234703384e-16L,
481
 
         0.264584017421245080294e-18L,
482
 
      };
483
 
      result = tools::evaluate_polynomial(P, s - 15) / tools::evaluate_polynomial(Q, s - 15);
484
 
      result = 1 + exp(result);
485
 
   }
486
 
   else if(s < 63)
487
 
   {
488
 
      result = 1 + pow(T(2), -s);
489
 
   }
490
 
   else
491
 
   {
492
 
      result = 1;
493
 
   }
494
 
   return result;
495
 
}
496
 
 
497
 
template <class T, class Policy>
498
 
T zeta_imp_prec(T s, T sc, const Policy& pol, const mpl::int_<113>&)
499
 
{
500
 
   BOOST_MATH_STD_USING
501
 
   T result;
502
 
   if(s < 1)
503
 
   {
504
 
      // Rational Approximation
505
 
      // Maximum Deviation Found:                     9.493e-37
506
 
      // Expected Error Term:                         9.492e-37
507
 
      // Max error found at long double precision:    7.281332e-31
508
 
 
509
 
      static const T P[10] = {    
510
 
         -1L,
511
 
         -0.0353008629988648122808504280990313668L,
512
 
         0.0107795651204927743049369868548706909L,
513
 
         0.000523961870530500751114866884685172975L,
514
 
         -0.661805838304910731947595897966487515e-4L,
515
 
         -0.658932670403818558510656304189164638e-5L,
516
 
         -0.103437265642266106533814021041010453e-6L,
517
 
         0.116818787212666457105375746642927737e-7L,
518
 
         0.660690993901506912123512551294239036e-9L,
519
 
         0.113103113698388531428914333768142527e-10L,
520
 
        };
521
 
      static const T Q[11] = {    
522
 
         1L,
523
 
         -0.387483472099602327112637481818565459L,
524
 
         0.0802265315091063135271497708694776875L,
525
 
         -0.0110727276164171919280036408995078164L,
526
 
         0.00112552716946286252000434849173787243L,
527
 
         -0.874554160748626916455655180296834352e-4L,
528
 
         0.530097847491828379568636739662278322e-5L,
529
 
         -0.248461553590496154705565904497247452e-6L,
530
 
         0.881834921354014787309644951507523899e-8L,
531
 
         -0.217062446168217797598596496310953025e-9L,
532
 
         0.315823200002384492377987848307151168e-11L,
533
 
      };
534
 
      result = tools::evaluate_polynomial(P, sc) / tools::evaluate_polynomial(Q, sc);
535
 
      result += (sc);
536
 
      result /= (sc);
537
 
   }
538
 
   else if(s <= 2)
539
 
   {
540
 
      // Maximum Deviation Found:                     1.616e-37
541
 
      // Expected Error Term:                         -1.615e-37
542
 
 
543
 
      static const T P[10] = {    
544
 
         0.577215664901532860606512090082402431L,
545
 
         0.255597968739771510415479842335906308L,
546
 
         0.0494056503552807274142218876983542205L,
547
 
         0.00551372778611700965268920983472292325L,
548
 
         0.00043667616723970574871427830895192731L,
549
 
         0.268562259154821957743669387915239528e-4L,
550
 
         0.109249633923016310141743084480436612e-5L,
551
 
         0.273895554345300227466534378753023924e-7L,
552
 
         0.583103205551702720149237384027795038e-9L,
553
 
         -0.835774625259919268768735944711219256e-11L,
554
 
      };
555
 
      static const T Q[11] = {    
556
 
         1L,
557
 
         0.316661751179735502065583176348292881L,
558
 
         0.0540401806533507064453851182728635272L,
559
 
         0.00598621274107420237785899476374043797L,
560
 
         0.000474907812321704156213038740142079615L,
561
 
         0.272125421722314389581695715835862418e-4L,
562
 
         0.112649552156479800925522445229212933e-5L,
563
 
         0.301838975502992622733000078063330461e-7L,
564
 
         0.422960728687211282539769943184270106e-9L,
565
 
         -0.377105263588822468076813329270698909e-11L,
566
 
         -0.581926559304525152432462127383600681e-13L,
567
 
      };
568
 
      result = tools::evaluate_polynomial(P, -sc) / tools::evaluate_polynomial(Q, -sc);
569
 
      result += 1 / (-sc);
570
 
   }
571
 
   else if(s <= 4)
572
 
   {
573
 
      // Maximum Deviation Found:                     1.891e-36
574
 
      // Expected Error Term:                         -1.891e-36
575
 
      // Max error found: 2.171527e-35
576
 
 
577
 
      static const float Y = 0.6986598968505859375;
578
 
      static const T P[11] = {    
579
 
         -0.0537258300023595010275848333539748089L,
580
 
         0.0429086930802630159457448174466342553L,
581
 
         0.0136148228754303412510213395034056857L,
582
 
         0.00190231601036042925183751238033763915L,
583
 
         0.000186880390916311438818302549192456581L,
584
 
         0.145347370745893262394287982691323657e-4L,
585
 
         0.805843276446813106414036600485884885e-6L,
586
 
         0.340818159286739137503297172091882574e-7L,
587
 
         0.115762357488748996526167305116837246e-8L,
588
 
         0.231904754577648077579913403645767214e-10L,
589
 
         0.340169592866058506675897646629036044e-12L,
590
 
      };
591
 
      static const T Q[12] = {    
592
 
         1L,
593
 
         0.363755247765087100018556983050520554L,
594
 
         0.0696581979014242539385695131258321598L,
595
 
         0.00882208914484611029571547753782014817L,
596
 
         0.000815405623261946661762236085660996718L,
597
 
         0.571366167062457197282642344940445452e-4L,
598
 
         0.309278269271853502353954062051797838e-5L,
599
 
         0.12822982083479010834070516053794262e-6L,
600
 
         0.397876357325018976733953479182110033e-8L,
601
 
         0.8484432107648683277598472295289279e-10L,
602
 
         0.105677416606909614301995218444080615e-11L,
603
 
         0.547223964564003701979951154093005354e-15L,
604
 
      };
605
 
      result = tools::evaluate_polynomial(P, s - 2) / tools::evaluate_polynomial(Q, s - 2);
606
 
      result += Y + 1 / (-sc);
607
 
   }
608
 
   else if(s <= 6)
609
 
   {
610
 
      // Max error in interpolated form:             1.510e-37
611
 
      // Max error found at long double precision:   2.769266e-34
612
 
 
613
 
      static const T Y = 3.28348541259765625F;
614
 
 
615
 
      static const T P[13] = {    
616
 
         0.786383506575062179339611614117697622L,
617
 
         0.495766593395271370974685959652073976L,
618
 
         -0.409116737851754766422360889037532228L,
619
 
         -0.57340744006238263817895456842655987L,
620
 
         -0.280479899797421910694892949057963111L,
621
 
         -0.0753148409447590257157585696212649869L,
622
 
         -0.0122934003684672788499099362823748632L,
623
 
         -0.00126148398446193639247961370266962927L,
624
 
         -0.828465038179772939844657040917364896e-4L,
625
 
         -0.361008916706050977143208468690645684e-5L,
626
 
         -0.109879825497910544424797771195928112e-6L,
627
 
         -0.214539416789686920918063075528797059e-8L,
628
 
         -0.15090220092460596872172844424267351e-10L,
629
 
      };
630
 
      static const T Q[14] = {    
631
 
         1L,
632
 
         1.69490865837142338462982225731926485L,
633
 
         1.22697696630994080733321401255942464L,
634
 
         0.495409420862526540074366618006341533L,
635
 
         0.122368084916843823462872905024259633L,
636
 
         0.0191412993625268971656513890888208623L,
637
 
         0.00191401538628980617753082598351559642L,
638
 
         0.000123318142456272424148930280876444459L,
639
 
         0.531945488232526067889835342277595709e-5L,
640
 
         0.161843184071894368337068779669116236e-6L,
641
 
         0.305796079600152506743828859577462778e-8L,
642
 
         0.233582592298450202680170811044408894e-10L,
643
 
         -0.275363878344548055574209713637734269e-13L,
644
 
         0.221564186807357535475441900517843892e-15L,
645
 
      };
646
 
      result = tools::evaluate_polynomial(P, s - 4) / tools::evaluate_polynomial(Q, s - 4);
647
 
      result -= Y;
648
 
      result = 1 + exp(result);
649
 
   }
650
 
   else if(s < 10)
651
 
   {
652
 
      // Max error in interpolated form:             1.999e-34
653
 
      // Max error found at long double precision:   2.156186e-33
654
 
 
655
 
      static const T P[13] = {    
656
 
         -4.0545627381873738086704293881227365L,
657
 
         -4.70088348734699134347906176097717782L,
658
 
         -2.36921550900925512951976617607678789L,
659
 
         -0.684322583796369508367726293719322866L,
660
 
         -0.126026534540165129870721937592996324L,
661
 
         -0.015636903921778316147260572008619549L,
662
 
         -0.00135442294754728549644376325814460807L,
663
 
         -0.842793965853572134365031384646117061e-4L,
664
 
         -0.385602133791111663372015460784978351e-5L,
665
 
         -0.130458500394692067189883214401478539e-6L,
666
 
         -0.315861074947230418778143153383660035e-8L,
667
 
         -0.500334720512030826996373077844707164e-10L,
668
 
         -0.420204769185233365849253969097184005e-12L,
669
 
        };
670
 
      static const T Q[14] = {    
671
 
         1L,
672
 
         0.97663511666410096104783358493318814L,
673
 
         0.40878780231201806504987368939673249L,
674
 
         0.0963890666609396058945084107597727252L,
675
 
         0.0142207619090854604824116070866614505L,
676
 
         0.00139010220902667918476773423995750877L,
677
 
         0.940669540194694997889636696089994734e-4L,
678
 
         0.458220848507517004399292480807026602e-5L,
679
 
         0.16345521617741789012782420625435495e-6L,
680
 
         0.414007452533083304371566316901024114e-8L,
681
 
         0.68701473543366328016953742622661377e-10L,
682
 
         0.603461891080716585087883971886075863e-12L,
683
 
         0.294670713571839023181857795866134957e-16L,
684
 
         -0.147003914536437243143096875069813451e-18L,
685
 
        };
686
 
      result = tools::evaluate_polynomial(P, s - 6) / tools::evaluate_polynomial(Q, s - 6);
687
 
      result = 1 + exp(result);
688
 
   }
689
 
   else if(s < 17)
690
 
   {
691
 
      // Max error in interpolated form:             1.641e-32
692
 
      // Max error found at long double precision:   1.696121e-32
693
 
      static const T P[13] = {    
694
 
         -6.91319491921722925920883787894829678L,
695
 
         -3.65491257639481960248690596951049048L,
696
 
         -0.813557553449954526442644544105257881L,
697
 
         -0.0994317301685870959473658713841138083L,
698
 
         -0.00726896610245676520248617014211734906L,
699
 
         -0.000317253318715075854811266230916762929L,
700
 
         -0.66851422826636750855184211580127133e-5L,
701
 
         0.879464154730985406003332577806849971e-7L,
702
 
         0.113838903158254250631678791998294628e-7L,
703
 
         0.379184410304927316385211327537817583e-9L,
704
 
         0.612992858643904887150527613446403867e-11L,
705
 
         0.347873737198164757035457841688594788e-13L,
706
 
         -0.289187187441625868404494665572279364e-15L,
707
 
        };
708
 
      static const T Q[14] = {    
709
 
         1L,
710
 
         0.427310044448071818775721584949868806L,
711
 
         0.074602514873055756201435421385243062L,
712
 
         0.00688651562174480772901425121653945942L,
713
 
         0.000360174847635115036351323894321880445L,
714
 
         0.973556847713307543918865405758248777e-5L,
715
 
         -0.853455848314516117964634714780874197e-8L,
716
 
         -0.118203513654855112421673192194622826e-7L,
717
 
         -0.462521662511754117095006543363328159e-9L,
718
 
         -0.834212591919475633107355719369463143e-11L,
719
 
         -0.5354594751002702935740220218582929e-13L,
720
 
         0.406451690742991192964889603000756203e-15L,
721
 
         0.887948682401000153828241615760146728e-19L,
722
 
         -0.34980761098820347103967203948619072e-21L,
723
 
        };
724
 
      result = tools::evaluate_polynomial(P, s - 10) / tools::evaluate_polynomial(Q, s - 10);
725
 
      result = 1 + exp(result);
726
 
   }
727
 
   else if(s < 30)
728
 
   {
729
 
      // Max error in interpolated form:             1.563e-31
730
 
      // Max error found at long double precision:   1.562725e-31
731
 
 
732
 
      static const T P[13] = {    
733
 
         -11.7824798233959252791987402769438322L,
734
 
         -4.36131215284987731928174218354118102L,
735
 
         -0.732260980060982349410898496846972204L,
736
 
         -0.0744985185694913074484248803015717388L,
737
 
         -0.00517228281320594683022294996292250527L,
738
 
         -0.000260897206152101522569969046299309939L,
739
 
         -0.989553462123121764865178453128769948e-5L,
740
 
         -0.286916799741891410827712096608826167e-6L,
741
 
         -0.637262477796046963617949532211619729e-8L,
742
 
         -0.106796831465628373325491288787760494e-9L,
743
 
         -0.129343095511091870860498356205376823e-11L,
744
 
         -0.102397936697965977221267881716672084e-13L,
745
 
         -0.402663128248642002351627980255756363e-16L,
746
 
      };
747
 
      static const T Q[14] = {    
748
 
         1L,
749
 
         0.311288325355705609096155335186466508L,
750
 
         0.0438318468940415543546769437752132748L,
751
 
         0.00374396349183199548610264222242269536L,
752
 
         0.000218707451200585197339671707189281302L,
753
 
         0.927578767487930747532953583797351219e-5L,
754
 
         0.294145760625753561951137473484889639e-6L,
755
 
         0.704618586690874460082739479535985395e-8L,
756
 
         0.126333332872897336219649130062221257e-9L,
757
 
         0.16317315713773503718315435769352765e-11L,
758
 
         0.137846712823719515148344938160275695e-13L,
759
 
         0.580975420554224366450994232723910583e-16L,
760
 
         -0.291354445847552426900293580511392459e-22L,
761
 
         0.73614324724785855925025452085443636e-25L,
762
 
      };
763
 
      result = tools::evaluate_polynomial(P, s - 17) / tools::evaluate_polynomial(Q, s - 17);
764
 
      result = 1 + exp(result);
765
 
   }
766
 
   else if(s < 74)
767
 
   {
768
 
      // Max error in interpolated form:             2.311e-27
769
 
      // Max error found at long double precision:   2.297544e-27
770
 
      static const T P[14] = {    
771
 
         -20.7944102007844314586649688802236072L,
772
 
         -4.95759941987499442499908748130192187L,
773
 
         -0.563290752832461751889194629200298688L,
774
 
         -0.0406197001137935911912457120706122877L,
775
 
         -0.0020846534789473022216888863613422293L,
776
 
         -0.808095978462109173749395599401375667e-4L,
777
 
         -0.244706022206249301640890603610060959e-5L,
778
 
         -0.589477682919645930544382616501666572e-7L,
779
 
         -0.113699573675553496343617442433027672e-8L,
780
 
         -0.174767860183598149649901223128011828e-10L,
781
 
         -0.210051620306761367764549971980026474e-12L,
782
 
         -0.189187969537370950337212675466400599e-14L,
783
 
         -0.116313253429564048145641663778121898e-16L,
784
 
         -0.376708747782400769427057630528578187e-19L,
785
 
      };
786
 
      static const T Q[16] = {    
787
 
         1L,
788
 
         0.205076752981410805177554569784219717L,
789
 
         0.0202526722696670378999575738524540269L,
790
 
         0.001278305290005994980069466658219057L,
791
 
         0.576404779858501791742255670403304787e-4L,
792
 
         0.196477049872253010859712483984252067e-5L,
793
 
         0.521863830500876189501054079974475762e-7L,
794
 
         0.109524209196868135198775445228552059e-8L,
795
 
         0.181698713448644481083966260949267825e-10L,
796
 
         0.234793316975091282090312036524695562e-12L,
797
 
         0.227490441461460571047545264251399048e-14L,
798
 
         0.151500292036937400913870642638520668e-16L,
799
 
         0.543475775154780935815530649335936121e-19L,
800
 
         0.241647013434111434636554455083309352e-28L,
801
 
         -0.557103423021951053707162364713587374e-31L,
802
 
         0.618708773442584843384712258199645166e-34L,
803
 
      };
804
 
      result = tools::evaluate_polynomial(P, s - 30) / tools::evaluate_polynomial(Q, s - 30);
805
 
      result = 1 + exp(result);
806
 
   }
807
 
   else if(s < 117)
808
 
   {
809
 
      result = 1 + pow(T(2), -s);
810
 
   }
811
 
   else
812
 
   {
813
 
      result = 1;
814
 
   }
815
 
   return result;
816
 
}
817
 
 
818
 
template <class T, class Policy, class Tag>
819
 
T zeta_imp(T s, T sc, const Policy& pol, const Tag& tag)
820
 
{
821
 
   BOOST_MATH_STD_USING
822
 
   if(s == 1)
823
 
      return policies::raise_pole_error<T>(
824
 
         "boost::math::zeta<%1%>", 
825
 
         "Evaluation of zeta function at pole %1%", 
826
 
         s, pol);
827
 
   T result;
828
 
   if(s == 0)
829
 
   {
830
 
      result = -0.5;
831
 
   }
832
 
   else if(s < 0)
833
 
   {
834
 
      std::swap(s, sc);
835
 
      if(floor(sc/2) == sc/2)
836
 
         result = 0;
837
 
      else
838
 
      {
839
 
         result = boost::math::sin_pi(0.5f * sc, pol)
840
 
            * 2 * pow(2 * constants::pi<T>(), -s) 
841
 
            * boost::math::tgamma(s, pol) 
842
 
            * zeta_imp(s, sc, pol, tag);
843
 
      }
844
 
   }
845
 
   else
846
 
   {
847
 
      result = zeta_imp_prec(s, sc, pol, tag);
848
 
   }
849
 
   return result;
850
 
}
851
 
 
852
 
} // detail
853
 
 
854
 
template <class T, class Policy>
855
 
inline typename tools::promote_args<T>::type zeta(T s, const Policy&)
856
 
{
857
 
   typedef typename tools::promote_args<T>::type result_type;
858
 
   typedef typename policies::evaluation<result_type, Policy>::type value_type;
859
 
   typedef typename policies::precision<result_type, Policy>::type precision_type;
860
 
   typedef typename policies::normalise<
861
 
      Policy, 
862
 
      policies::promote_float<false>, 
863
 
      policies::promote_double<false>, 
864
 
      policies::discrete_quantile<>,
865
 
      policies::assert_undefined<> >::type forwarding_policy;
866
 
   typedef typename mpl::if_<
867
 
      mpl::less_equal<precision_type, mpl::int_<0> >,
868
 
      mpl::int_<0>,
869
 
      typename mpl::if_<
870
 
         mpl::less_equal<precision_type, mpl::int_<53> >,
871
 
         mpl::int_<53>,  // double
872
 
         typename mpl::if_<
873
 
            mpl::less_equal<precision_type, mpl::int_<64> >,
874
 
            mpl::int_<64>, // 80-bit long double
875
 
            typename mpl::if_<
876
 
               mpl::less_equal<precision_type, mpl::int_<113> >,
877
 
               mpl::int_<113>, // 128-bit long double
878
 
               mpl::int_<0> // too many bits, use generic version.
879
 
            >::type
880
 
         >::type
881
 
      >::type
882
 
   >::type tag_type;
883
 
   //typedef mpl::int_<0> tag_type;
884
 
 
885
 
   return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::zeta_imp(
886
 
      static_cast<value_type>(s),
887
 
      static_cast<value_type>(1 - static_cast<value_type>(s)),
888
 
      forwarding_policy(),
889
 
      tag_type()), "boost::math::zeta<%1%>(%1%)");
890
 
}
891
 
 
892
 
template <class T>
893
 
inline typename tools::promote_args<T>::type zeta(T s)
894
 
{
895
 
   return zeta(s, policies::policy<>());
896
 
}
897
 
 
898
 
}} // namespaces
899
 
 
900
 
#endif // BOOST_MATH_ZETA_HPP
901
 
 
902
 
 
903