Численно рассчитать комбинации факториалов и полиномов - PullRequest
Новая
       72

Численно рассчитать комбинации факториалов и полиномов

4 голосов
/ 16 марта 2020

Я пытаюсь написать короткую подпрограмму C ++ для вычисления следующей функции F (i, j, z) для заданных целых чисел j> i (обычно они l ie между 0 и 100) и комплексного числа z (ограничено | z | <100), где L - <a href="https://mathworld.wolfram.com/AssociatedLaguerrePolynomial.html" rel="nofollow noreferrer"> связанные полиномы Лагерра :

enter image description here

Проблема в том, что я хочу эту функцию быть вызываемым изнутри ядра CUDA (то есть с атрибутом __device__). Поэтому стандартные библиотечные / Boost / et c функции исключены, если только они не достаточно просты, чтобы их можно было повторно реализовать самостоятельно - это особенно относится к полиномам Лагерра, существующим в Boost и C ++ 17. Независимо от того, удастся ли мне обернуть какую-либо стандартную функцию для полиномов Лагерра, у меня все еще будет аналогичный предварительный фактор для вычисления формы (z ^ j / j!).

Вопрос: Как Могу ли я сделать относительно простую реализацию такой функции, не внося значительную числовую нестабильность?

Моя идея до сих пор состоит в том, чтобы вычислять L и его предварительный множитель независимо. Предварительный коэффициент я вычислю сначала с 0 до ji и вычислю (z ^ 1 * z ^ 2/2 * ... * z ^ (j-1) / (ji)!). Затем я вычислю оставшийся фактор exp (- | z | ^ 2/2) * (ji)! * sqrt (i! / j!) (либо аналогичным образом, либо через Gamma-функцию, которая реализована в математике CUDA). Идея состоит в том, чтобы найти минимальный алгоритм для вычисления связанного полинома Лагерра, если мне не удастся обернуть реализацию, например, в Boost или GNU C ++.

Правка / примечание: Выражение для F фактически взрывается численно для некоторых значений i / j. Он был неверно получен в источнике, где я его получил, и вместо этого индексы связанных полиномов Лагерра должны быть L_i ^ (ji). Это не лишает законной силы подходы, предложенные в ответах / комментариях.

Ответы [ 2 ]

2 голосов
/ 16 марта 2020

Я рекомендую найти рекуррентное соотношение для коэффициентов полинома Лагерра: 

C(k+1) = g(k)C(k)
g(k) = C(k+1) / C(k)
g(k) = -z * (j - k) / ((j - i + k + 1) * (k + 1)) //Verify this yourself :)

Это позволяет вам избежать большинства факториалов при вычислении полинома.

После этого я буду следовать Идея Северина выполнять вычисления в логарифмах, чтобы не перегружать двойной диапазон с плавающей запятой: 

log(F) = log(sqrt(i!/j!)) - |z|^2 + (j-i) * log(-z) + log(L(|z|^2))
log(L) = log((2*j - i)!) + log(sum) // where the summation is computed using the recurrence relation above

и использовать тот факт, что: 

log(a!) = sum(k=1..a, log(k))

, а также: 

log(z) = log(|z|) + I * arg(z) for complex z
log(-z) = log(|z|) + I * arg(-z)
log(-z) = log(|z|) - I * arg(z)

для части log(sqrt(i!/j!)), которую я бы сделал (предполагая, что j> = i): 

  log(sqrt(i!/j!))
= 0.5 * (log(i!) - log(j!))
= -0.5 * sum(k==i+1..j, log(k))

Я не пробовал это, так что определенно могут быть небольшие ошибки здесь и там. Этот ответ больше о технике, а не о готовности к копированию

1 голос
/ 16 марта 2020

Ну, что вам нужно сделать, это логарифмировать его

Предполагая натуральный логарифм,

q = log (z ^ j / j!) = Log (z ^ j) - log ( j!) = j * log (z) - log (Gamma (j + 1))

Первый термин прост, второй - стандартная функция C ++ lgamma (x) (или вы можете использовать GSL).

вычислить значение q и вернуть cexp(q)

В этом методе вы также можете сложить показатель степени

...