gsl_sum.h 5.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163
  1. /* sum/gsl_sum.h
  2. *
  3. * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
  4. *
  5. * This program is free software; you can redistribute it and/or modify
  6. * it under the terms of the GNU General Public License as published by
  7. * the Free Software Foundation; either version 3 of the License, or (at
  8. * your option) any later version.
  9. *
  10. * This program is distributed in the hope that it will be useful, but
  11. * WITHOUT ANY WARRANTY; without even the implied warranty of
  12. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  13. * General Public License for more details.
  14. *
  15. * You should have received a copy of the GNU General Public License
  16. * along with this program; if not, write to the Free Software
  17. * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18. */
  19. /* Author: G. Jungman */
  20. #ifndef __GSL_SUM_H__
  21. #define __GSL_SUM_H__
  22. #include <stdlib.h>
  23. #undef __BEGIN_DECLS
  24. #undef __END_DECLS
  25. #ifdef __cplusplus
  26. # define __BEGIN_DECLS extern "C" {
  27. # define __END_DECLS }
  28. #else
  29. # define __BEGIN_DECLS /* empty */
  30. # define __END_DECLS /* empty */
  31. #endif
  32. __BEGIN_DECLS
  33. /* Workspace for Levin U Transform with error estimation,
  34. *
  35. * size = number of terms the workspace can handle
  36. * sum_plain = simple sum of series
  37. * q_num = backward diagonal of numerator; length = size
  38. * q_den = backward diagonal of denominator; length = size
  39. * dq_num = table of numerator derivatives; length = size**2
  40. * dq_den = table of denominator derivatives; length = size**2
  41. * dsum = derivative of sum wrt term i; length = size
  42. */
  43. typedef struct
  44. {
  45. size_t size;
  46. size_t i; /* position in array */
  47. size_t terms_used; /* number of calls */
  48. double sum_plain;
  49. double *q_num;
  50. double *q_den;
  51. double *dq_num;
  52. double *dq_den;
  53. double *dsum;
  54. }
  55. gsl_sum_levin_u_workspace;
  56. gsl_sum_levin_u_workspace *gsl_sum_levin_u_alloc (size_t n);
  57. void gsl_sum_levin_u_free (gsl_sum_levin_u_workspace * w);
  58. /* Basic Levin-u acceleration method.
  59. *
  60. * array = array of series elements
  61. * n = size of array
  62. * sum_accel = result of summation acceleration
  63. * err = estimated error
  64. *
  65. * See [Fessler et al., ACM TOMS 9, 346 (1983) and TOMS-602]
  66. */
  67. int gsl_sum_levin_u_accel (const double *array,
  68. const size_t n,
  69. gsl_sum_levin_u_workspace * w,
  70. double *sum_accel, double *abserr);
  71. /* Basic Levin-u acceleration method with constraints on the terms
  72. * used,
  73. *
  74. * array = array of series elements
  75. * n = size of array
  76. * min_terms = minimum number of terms to sum
  77. * max_terms = maximum number of terms to sum
  78. * sum_accel = result of summation acceleration
  79. * err = estimated error
  80. *
  81. * See [Fessler et al., ACM TOMS 9, 346 (1983) and TOMS-602]
  82. */
  83. int gsl_sum_levin_u_minmax (const double *array,
  84. const size_t n,
  85. const size_t min_terms,
  86. const size_t max_terms,
  87. gsl_sum_levin_u_workspace * w,
  88. double *sum_accel, double *abserr);
  89. /* Basic Levin-u step w/o reference to the array of terms.
  90. * We only need to specify the value of the current term
  91. * to execute the step. See TOMS-745.
  92. *
  93. * sum = t0 + ... + t_{n-1} + term; term = t_{n}
  94. *
  95. * term = value of the series term to be added
  96. * n = position of term in series (starting from 0)
  97. * sum_accel = result of summation acceleration
  98. * sum_plain = simple sum of series
  99. */
  100. int
  101. gsl_sum_levin_u_step (const double term,
  102. const size_t n,
  103. const size_t nmax,
  104. gsl_sum_levin_u_workspace * w,
  105. double *sum_accel);
  106. /* The following functions perform the same calculation without
  107. estimating the errors. They require O(N) storage instead of O(N^2).
  108. This may be useful for summing many similar series where the size
  109. of the error has already been estimated reliably and is not
  110. expected to change. */
  111. typedef struct
  112. {
  113. size_t size;
  114. size_t i; /* position in array */
  115. size_t terms_used; /* number of calls */
  116. double sum_plain;
  117. double *q_num;
  118. double *q_den;
  119. double *dsum;
  120. }
  121. gsl_sum_levin_utrunc_workspace;
  122. gsl_sum_levin_utrunc_workspace *gsl_sum_levin_utrunc_alloc (size_t n);
  123. void gsl_sum_levin_utrunc_free (gsl_sum_levin_utrunc_workspace * w);
  124. int gsl_sum_levin_utrunc_accel (const double *array,
  125. const size_t n,
  126. gsl_sum_levin_utrunc_workspace * w,
  127. double *sum_accel, double *abserr_trunc);
  128. int gsl_sum_levin_utrunc_minmax (const double *array,
  129. const size_t n,
  130. const size_t min_terms,
  131. const size_t max_terms,
  132. gsl_sum_levin_utrunc_workspace * w,
  133. double *sum_accel, double *abserr_trunc);
  134. int gsl_sum_levin_utrunc_step (const double term,
  135. const size_t n,
  136. gsl_sum_levin_utrunc_workspace * w,
  137. double *sum_accel);
  138. __END_DECLS
  139. #endif /* __GSL_SUM_H__ */