¿Cómo usar la función polilogaritmo en c ++?

9

¿Hay alguna directiva de preprocesador que pueda usarse para usar la función polylog? ¿O está incluido en cmath? Si es así, ¿lo llamas Li o Polylog?


EDITAR: lo que realmente estoy tratando de hacer es dar un valor analítico para la integral indefinida de la función

x3ex1

que implica funciones de polilogaritmo. Pero si alguien tiene una sugerencia para otra forma de integrar esta función analíticamente, sería bienvenido cualquier idea.


flamencos
fuente
Para fines de búsqueda: la función que se considera en el OP está relacionada con las funciones de Debye . Esta nota puede ser de alguna utilidad.
JM
Otra relación cercana es la integral incompleta de Fermi-Dirac .
hardmath

Respuestas:

7

Hay una biblioteca GPL'd C, ANANT - Algorithms in Analytic Number Theory de Linas Vepstas, que incluye la implementación de multiprecisión del polilogaritmo, basado en GMP .

Desde su archivo README:

Este proyecto contiene implementaciones ad-hoc de diversas funciones analíticas de interés en la teoría de números, incluida la función gamma, la función zeta de Riemann, el polilogaritmo y la función de signo de interrogación de Minkowski. La implementación utiliza la biblioteca Gnu Multi-Precision (GMP) para realizar todas las operaciones de bajo nivel. El código en este documento tiene licencia bajo los términos de la licencia Gnu GPLv3.

La GSL (GNU Scientific Library) aparentemente solo tiene la función de dilogaritmo . Sin embargo, siguiendo una pista de @JM, se encuentra la función Debye que proporciona la integral ulterior (hasta un múltiplo escalar) implementada con doble precisión (consulte GSL 7.10 Órdenes de funciones Debye 1 a 6):

Dn(x)=nxn0xtndtet1


El software de integración simbólica como Mathematica o Maxima proporciona:

0xt3dtet1=6Li4(ex)6xLi3(ex)+3x2Li2(ex)+x3log(1ex)x44π415

El lado izquierdo es obviamente un valor puramente real si , pero los polilogaritmos mostrados tendrán un valor complejo (porque , por lo que la igualdad depende de la cancelación total de partes imaginarias). Podemos evitar la necesidad de aritmética compleja en este caso sustituyendo la expresión:x>0ex>1

0xt3dtet1=6Li4(ex)6xLi3(ex)3x2Li2(ex)x3Li1(ex)+π415

Esto es una mejora porque con los argumentos de pollogaritmo en , los resultados son valores puramente reales. Tenga en cuenta el resultado adecuado cuando es cero, y esto se logra mediante la cancelación entre el término inicial y la constante. Por lo tanto, el error relativo podría ser un problema para pequeños valores positivos de .[0,1]x=0x

Tenga en cuenta que nuestra misteriosa constante es el límite superior limitante en estas integrales (monótonas):π4/15

0t3dtet1=Γ(4)ζ(4)=6π490

Ahora podemos volver a la pregunta del título, ¿Cómo usar la función de polilogaritmo en c ++? Vale la pena señalar que no existe una implementación estándar de funciones de polilogaritmo para C o incluso C ++ . Si el objetivo es evitar cualquier biblioteca adicional para su implementación, lo bastante bien lo lleva a desarrollar sus propias rutinas, tal vez en la línea sugerida por el documento de David C. Wood al que enlaza la Respuesta de GertVdE.

Además de las rutinas de multiprecisión sugeridas en la primera parte de mi Respuesta, hay una biblioteca matemática de doble precisión madura (gratuita) en Cephes por Stephen L. Moshier que implementa versiones reales ( polylog) y complejas ( cpolylog) de las funciones especiales del polilogaritmo. Aunque su precisión depende en parte de las funciones matemáticas estándar subyacentes de C, la documentación fuente de Cephes informa sobre las pruebas y los errores teóricos máximos para las órdenes 1 a 4 en aproximadamente los límites de la doble precisión.

Alternativamente, es posible que desee utilizar otro software para verificar directamente (sin hacer referencia a los policlogaritmos) las rutinas de cuadratura que escribió para su integral. Como esbozo en esta pregunta de Math.SE , la serie de potencia centrada en el origen de la integral tiene una convergencia limitada, pero esto se puede mitigar mediante el uso de una expansión de fracción continua en su lugar.

Para una satisfacción inmediata, recomiendo las rutinas QUADPACK de cuadratura numérica (gratis) incluidas en Maxima , específicamente quad_qag. Por ejemplo, encuentre la integral sobre [0,5] con este comando Maxima:

(%i1) quad_qag(x^3/(%e^x - 1), x, 0, 5, 2);
(%o1) [4.899892158330582,5.4399730923588665*10^-14,21,0]

De los argumentos de entrada, solo el último lleva una explicación. El quinto argumento para quad_qagespecificar qué regla aplicar en la cuadratura adaptativa. Los valores posibles son de 1 a 6 y ofrecen una sofisticación / precisión cada vez mayor. La línea de salida proporciona primero la cuadratura numérica, seguida de una estimación de su error absoluto, el número de subintervalos / pasos utilizados y un código de retorno (aquí cero significa que no se encontraron errores o condiciones especiales).

hardmath
fuente
1
"Incluso C ++", un mejor enlace para funciones especiales sería en.cppreference.com/w/cpp/numeric/special_math , aunque la función aún no está allí. Es sorprendente que ni siquiera esté en la biblioteca Boost.Math boost.org/doc/libs/1_68_0/libs/math/doc/html/special.html .
alfC
1
@alfC: Gracias por el mejor enlace. Usaré eso arriba, para ilustrar la continua falta de soporte estándar para esta función / familia de funciones.
hardmath
4

En primer lugar, debe elegir en función de su aplicación si necesita una aritmética de alta precisión (es decir, ¿estará satisfecho con los resultados de doble precisión de IEEE para las funciones de polylog o necesita una mayor precisión)? Si necesita alta precisión, puede buscar en la familia de herramientas de la biblioteca GMP.

Si no lo hace, puede usar aproximaciones. Algunas investigaciones literarias me señalaron este artículo . Al final del artículo, hay una "tabla de selección": en función de los argumentos de los polylogs que necesita, puede seleccionar una fórmula de aproximación. Pero tenga cuidado de verificar la estabilidad y la precisión.

Si no necesita demasiadas evaluaciones (no en un bucle anidado), simplemente optaría por la cuadratura numérica utilizando el método de doble exponencial.

GertVdE
fuente
La precisión necesaria sería alrededor de 10 ^ -6. ¿Conoces alguna otra forma que no implique el uso de bibliotecas adicionales?
flamingohats
@flamingohats: luego buscaría las expresiones en el documento al que vinculé en mi respuesta. ¿Cuál es su rango para ? x
GertVdE
@flamingohats: ¿Estás diciendo que necesitas aproximar para o algo por el estilo? No entiendo por qué "solo se necesitarán dos evaluaciones", a menos que quiera decir que este es un requisito / ideal adicional para la cuadratura u otra aproximación. Estoy pensando que podemos determinar una aproximación polinómica o racional que proporcione la precisión requerida. x[0,10]0xt3et1dtx[0,10]
hardmath
Lo siento, me confundí. He aproximado la integral de 0.65 a 5.025 usando el método trapezoidal, y necesito una fórmula para encontrar el valor exacto para poder comparar la aproximación con un valor analítico. Sé que esto será aproximado porque es un número de coma flotante, por lo que una precisión de 1e-6 estará bien. Si de alguna manera puedo aprender cómo ingresar la función polylog en el IDE, debería funcionar.
flamingohats
1
@flamingohats: todavía no está claro si necesita un valor de un disparo ( ) o si realmente lo necesita programado en tu código Si es el primero, ¿por qué no usar Sage, Euler Toolbox, ... para verificar una serie de valores cruciales para usted (ya sea utilizando su implementación de Polylog o sus sofisticadas reglas de cuadratura)? Si es el segundo, implemente un método de cuadratura "tradicional" (como lo hizo) y nuevamente use Sage, Euler, ... para probarlo y validarlo. 0.655.205t3et1dt=4.8498308528256668370925
GertVdE
3

j j =Fj(x)jj=12,12,32

Fj(x)=Lij+1(ex)
x
innisfree
fuente