Логарифмическая производная гамма-функции

[Ekaterinasun]

Подскажите, пожалуйста, как мне реализовать логарифмическую производную гаммы-функции на pascal или delphi в подпрограммах?

[caxap]

Это дигамма-функция. Поищите её разложение в ряд.

[Yuri Gendelman]

Ekaterinasun в сообщении #372515 писал(а):

Подскажите, пожалуйста, как мне реализовать логарифмическую производную гаммы-функции на pascal или delphi в подпрограммах?

Нашел в своем архиве, пользовался последний раз лет ...дцать назад

Логарифмическая производная гамма-функции (пси-функция)
Ссылки:
- Amit D., CACM, 1962, 12, algorithm 147
- М.И.Агеев, В.П.Алик, Ю.И.Марков, Библиотека алгоритмов 101б-150б / М., "Советское Радио", 1978, алгоритм 147б, стр.74-76

Дословный перевод с Алгол-60 на С:

/*
 * Вычисляется
 *                        d ln Г(x+1)     Г'(x+1)
 *   Psi(x) = psi(x+1) = ------------- = ---------
 *                            d x         Г(x+1)
 *
 * x наращивается до величины a. При этом получается погрешность,
 *               1      -8
 * меньшая чем ----- * a
 *              240
 */

#include <math.h>

double Psif(
    double x,    /* аргумент                                        */
    double a,    /* величина, до которой наращивается x, если x < a */
    int *ierr )  /* возвращает: 0, если получен результат            */
                 /*             1, если x - целое и  x < -1          */
{
    static double  pi = 3.141592654;
    double  Psi,y;

    *ierr=0;
    if (x==0.0)
        return (-0.5772156649);
    Psi=0.0;
    if (x<=-1.0)  {
        if (x==floor(x))  {
            *ierr=1;
            return (0.0);
        }
        x=-x-1.0;
        Psi=pi/tan(pi*x);
    }
    while (x<a)  {
        x+=1.0;
        Psi-=1.0/x;
    }
    y=1.0/(x*x);
    Psi+=log(x)+0.5/x-((y/252.0-0.00833333333)*y+0.0833333333)*y;
    return (Psi);
}

[Ekaterinasun]

Спасибо огромное, буду разбираться.Спасибо.

[worm2]

В качестве мелкой придирки к процитированному коду можно сказать, что надо бы побольше значащих цифр в константах взять. Точность типа double — 15-16 цифр, а там в константах записано только 8-9.
Хотя... я всего не знаю, может быть, там просто точность приближения меньше, тогда не стоит...

[Yuri Gendelman]

worm2 в сообщении #375333 писал(а):

надо бы побольше значащих цифр в константах взять.

Код этот был опубликован в 1962 году на языке Algol-60, задолго до появления 64-битного double. При переводе на русский учитывались исправленные в 60-х - начале 70-х ошибки, но текст алгоритмов не переписывался. Мой вариант на С - это тоже перевод. Для реальной работы алгоритмы из справочников оптимизировались для конкретной платформы. В данном случае константы $\pi$, $C$ и $5/6$ хорошо известны с очень большой точностью. Особенно последняя
В данном алгоритме гораздо большая проблема - правильный выбор параметра "a" в зависимости от требуемой точности.

В прошлом посте я забыл привести ссылку на книгу на русском:
Библиотека алгоритмов 101б-150б