Функция подсчета чисел Бернулли

@iluxa1810 · Регистрация: 03.09.2011

Author24 — интернет-сервис помощи студентам

Стоит задача: Написать функцию, которая считает числа Бернулли.(Формула снизу).
Вроде бы ничего сложного, но затруднение вызывает подсчет факториала(Он же биноминальный коефициент).
Попытался найти примеры реализаций на других языках, но ничего не нашел.

Напишите или функцию на паскале, или напишите словами, как должны происходить итерации.
Спасибо!

$https://www.cyberforum.ru/cgi-bin/latex.cgi?<br /> B_0=1\\\\B_n=\frac{-1}{n+1}\sum_{k=1}^n \binom{n+1}{k+1} B_{n-k},\quad n\in\mathbb{N}.<br />$

@Puporev · 21.09.2014, 17:33

Pascal

uses crt;
 {подсчет числа сочетаний без факториала}
function sochet(n,m:integer):extended;
var i:integer;
    c:extended;
begin
c:=1;
for i:=1 to m do
c:=c*(n-i+1)/i;
sochet:=c
end;
{подсчет чисел Бернулли с рекурсией}
function bernulli(n:integer):extended;
var k:integer;
    s,t:extended;
begin
if n=0 then bernulli:=1
else if n=1 then bernulli:=-0.5
else
 begin
  s:=0;
  for k:=1 to n do
   begin
    t:=bernulli(n-k);{запомним результат и очистим стек памяти}
    s:=s+sochet(n+1,k+1)*t;{вычислим следующее значение суммы}
   end;
  bernulli:=s*(-1)/(n+1);{вычислим число Бернулли}
 end;
if (n>1) and odd(n) then bernulli:=0{если N>1 и нечетное, то =0}
end;
var n:integer;
begin
clrscr;
write('n=');
readln(n);
write(bernulli(n):0:10);
readln
end.

при N>20 считает медленно...
Первые 20 сравнил с этим
https://ru.wikipedia.org/wiki/... 0%BB%D0%B8

@iluxa1810 · 21.09.2014, 18:16 **[ТС]**

Puporev, Спасибо, а можете комментарии к программе подправить? Некоторые комментарии с иероглифами.

@Puporev · 21.09.2014, 18:41

Поправил

ФедосеевПавел · 05.04.2019, 21:35

Подниму тему, рассмотрев её решение с другой стороны.
Возможно, что преподаватель будет ждать решение, сходное с предложенным Puporev, тогда возникнут вопросы о "сверхзнаниях".

О другом алгоритме вычисления чисел Бернулли

Вычисления производятся по реккурентной формуле
$https://www.cyberforum.ru/cgi-bin/latex.cgi?B_0=1$
$https://www.cyberforum.ru/cgi-bin/latex.cgi?B_n=\frac{-1}{n+1}\sum_{k=1}^{n}{C_{n+1}^{k+1}}{B_{n-k}}$
Несложно заметить, что при вычислении отдельного числа B_n все коэффициенты $https://www.cyberforum.ru/cgi-bin/latex.cgi?C_{n+1}^{k+1}$ в формуле являются элементами одной строки треугольника Паскаля. При вычислении B_n+1 будет использоваться уже следующая строка этого треугольника. Казалось бы, что это ничего не даёт. Но, существуют компактные алгоритмы вычисляющие только одну строку треугольника Паскаля на основе предыдущей строки, причём новая строка располагается на месте предыдущей, т.е. бережно относятся к выделяемой памяти
Получить и напечатать первые n строк треугольника Паскаля

Код для Turbo Pascal

Pascal

{
Вычисление чисел Бернулли (Bernoulli)
    числитель   OEIS:A027641
    знаменатель OEIS:A027642
}
program test;
 
type
  real = double;
  integer = longint;
 
const
  BN: array [0..21] of real = (1, -1 / 2, 1 / 6, 0, -1 / 30, 0, 1 / 42, 0, -1 / 30, 0,
    5 / 66, 0, -691 / 2730, 0, 7 / 6, 0, -3617 / 510, 0, 43867 / 798,
    0, -174611 / 330, 0);
 
  procedure Bernoulli(n: integer; var B: array of real);
  var
    C: array [0..33] of integer;
    i, k: integer;
    S: real;
  begin
    C[0] := 1;
    C[1] := 1;
    B[0] := 1;
    for i := 1 to n do
    begin
      C[i + 1] := 0;
      for k := i + 1 downto 1 do
        C[k] := C[k - 1] + C[k];
      S := 0;
      for k := i - 1 downto 0 do
        S := S + C[k] * B[k];
      B[i] := -S / (i + 1);
    end;
  end;
 
const
  n = 21;
var
  B: array [0..n] of real;
  i: integer;
 
begin
  Bernoulli(n, B);
 
  writeln('   Вычисленное          Точное               Погрешность');
  for i := 0 to 21 do
    writeln(B[i]: 20: 15, ' ', BN[i]: 20: 15, ' ', B[i] - BN[i]: 20: 15);
end.

ФедосеевПавел · 07.09.2020, 18:16

Подниму тему ещё раз, пополнив её ещё одним вариантом алгоритмом вычисления чисел Бернулли.

В статье о вычислениях чисел Бернулли на сайте OEIS приводится алгоритм на Python

Python

def SeidelBernoulli(n) :
    D = [0]*(n+2); D[1] = 1; b = True;
    R = [1, -0.5]; h = 1; p = 1; s = -2
    for i in range(2*n-2) :
        if b :
            h += 1; p = 4*p; s = -s; q = s*(p-1)
            for k in range(h,0,-1) : D[k] += D[k+1]
        else :
            for k in range(1,h,1) : D[k] += D[k-1]
            R.append(D[h-1]/q)
            R.append(0);
        b = not b
    return R
 
#дополню несколькими строками для проверки функции
n=10
T=SeidelBernoulli(2*n+1)
for i in range(2*n+1):
    print(T[i], end='\n ')

Он хорошо переносится на Pascal, и имеет преимущество в том, что сами вычисления до момента деления производятся в целых числах и в вычислениях не участвуют ранее полученные числа Бернулли - не накапливается ошибка.

По сравнению с кодом на Python программа на Pascal:

формирует все числа Бернулли, а не только ненулевые
переменная b заменена на функцию odd(i)
цикл по i заполняет числа от 2 до n, как и требует логика задачи вычисления B_n
заполнение R[0] и R[1] выполняется не в начале подпрограммы, а в конце, т.к. при i=1 значение R[1] заполняется нулём

При отладке я экспериментировал с точностью вычислений, поэтому переопределил вещественный тип и определил целый тип для элементов массива (т.к. без определённых ключей компилятор FPC считает размер типа integer равным 16 бит, что недостаточно для данных вычислений). При компиляции в TurboPascal нужно удалить строку с переопределением integer.

Pascal

type
  real = extended;
  int  = int32;
const
  //точные значения чисел Бернулли для сравнения вычислений
  BN: array [0..21] of real = (1, -1 / 2, 1 / 6, 0, -1 / 30, 0, 1 / 42, 0, -1 / 30, 0,
    5 / 66, 0, -691 / 2730, 0, 7 / 6, 0, -3617 / 510, 0, 43867 / 798,
    0, -174611 / 330, 0);
 
  procedure SeidelBernoulli(n: integer; var R: array of real);
  var
    D: array of int;
    p, s, q: int;
    h, i, k: integer;
  begin
    SetLength(D, n + 2);
    D[1] := 1;
    h := 1;
    p := 1;
    s := -2;
    for i := 1 to n do
    begin
      if odd(i) then
      begin
        Inc(h);
        p := 4 * p;
        s := -s;
        q := s * (p - 1);
        for k := h downto 0 + 1 do
          D[k] := D[k] + D[k + 1];
        R[i] := 0;
      end
      else
      begin
        for k := 1 to h - 1 do
          D[k] := D[k] + D[k - 1];
        R[i] := D[h - 1] / q;
      end;
    end;
    //т.к. в цикле при i=1 записывается R[1]:=0
    //то корректируем на точное и инициализируем R[0]
    R[0] := 1;
    if n > 0 then
      R[1] := -0.5;
  end;
 
var
  B: array of real;
  n: integer;
begin
  n := 20;
  SetLength(B, n + 1);
 
  SeidelBernoulli(n, B);
  writeln('Calc': 14, 'Exact': 21, 'Eps': 20);
  for n := 0 to n do
    writeln(B[n]: 20: 15, ' ', BN[n]: 20: 15, ' ', B[n] - BN[n]: 20: 15);
end.

@iluxa1810 13 / 13 / 11 Регистрация: 03.09.2011 Сообщений: 1,026
		1
	Функция подсчета чисел Бернулли 21.09.2014, 15:56. Показов 3778. Ответов 5 Метки нет (Все метки) Стоит задача: Написать функцию, которая считает числа Бернулли.(Формула снизу). Вроде бы ничего сложного, но затруднение вызывает подсчет факториала(Он же биноминальный коефициент). Попытался найти примеры реализаций на других языках, но ничего не нашел. Напишите или функцию на паскале, или напишите словами, как должны происходить итерации. Спасибо! $https://www.cyberforum.ru/cgi-bin/latex.cgi?<br /> B_0=1\\\\B_n=\frac{-1}{n+1}\sum_{k=1}^n \binom{n+1}{k+1} B_{n-k},\quad n\in\mathbb{N}.<br />$ 0

@iluxa1810 13 / 13 / 11 Регистрация: 03.09.2011 Сообщений: 1,026
	21.09.2014, 18:16 [ТС]	3
	Puporev, Спасибо, а можете комментарии к программе подправить? Некоторые комментарии с иероглифами. 0

@Puporev Почетный модератор 64300 / 47595 / 32743 Регистрация: 18.05.2008 Сообщений: 115,181
	21.09.2014, 18:41	4
	Поправил 1

	07.09.2020, 18:16