Implementación rápida y precisa de doble precisión de la función gamma incompleta

10

¿Cuál es la forma más moderna de implementar funciones especiales de doble precisión? Necesito la siguiente integral:

Fmetro(t)=0 01tu2metromi-ttu2retu=γ(metro+12,t)2tmetro+12
parametro=0 0,1,2,...yt>0 0, que puede escribirse en términos de la función gamma incompleta inferior. Aquí está mi implementación de Fortran y C:

https://gist.github.com/3764427

que usa la expansión en serie, resume los términos hasta la precisión dada y luego usa las relaciones de recursión para obtener de manera eficiente valores para menor metro. Lo probé bien y obtengo una precisión de 1e-15 para todos los valores de parámetros que necesito, vea los comentarios de la versión Fortran para más detalles.

¿Hay una mejor manera de implementarlo? Aquí hay una implementación de la función gamma en gfortran:

https://github.com/mirrors/gcc/blob/master/libgfortran/intrinsics/c99_functions.c#L1781

está usando la aproximación de funciones racionales en lugar de resumir algunas series infinitas que estoy haciendo. Creo que es un mejor enfoque, porque uno debería obtener una precisión uniforme. ¿Hay alguna forma canónica de abordar estas cosas, o hay que encontrar un algoritmo especial para cada función especial?

Actualización 1 :

Basado en los comentarios, aquí está la implementación usando SLATEC:

https://gist.github.com/3767621

reproduce valores de mi propia función, aproximadamente en el nivel de precisión 1e-15. Sin embargo, noté un problema que para t = 1e-6 ym = 50, el término 2 es igual a 1e-303 y para una "m" más alta simplemente comienza a dar respuestas incorrectas. Mi función no tiene este problema, porque utilizo una serie de relaciones de expansión / recurrencia directamente paraFm. Aquí hay un ejemplo de un valor correcto:tmetro+12Fmetro

,F100(1e-6)=4.97511945200351715E-003

pero no puedo obtener esto usando SLATEC porque el denominador explota. Como puede ver, el valor real de es agradable y pequeño.Fmetro

Actualización 2 :

Para evitar el problema anterior, se puede usar la función dgamit( función Gamma incompleta de Tricomi), entonces F(m, t) = dgamit(m+0.5_dp, t) * gamma(m+0.5_dp) / 2, ya no hay problema con , pero desafortunadamente explota para m 172 . Sin embargo, esto podría ser lo suficientemente alta m para mis propósitos.tgamma(m+0.5_dp)metro172metro

Ondřej Čertík
fuente
2
¿Por qué codificar su propia función? GSL, cephes y SLATEC lo implementan.
Geoff Oxberry
He actualizado la pregunta de por qué no uso SLATEC.
Ondřej Čertík
@ OndřejČertík ¡Has descubierto un error aparentemente! ¡Votó su pregunta!
Ali
Ali --- no es un error en SLATEC, pero en realidad necesito dividir entre t m + 1γ(z,X) para obtener un valor paraFm(t). Por lo tanto, el método numérico que funciona paraγ(z,x)podría no funcionar tan bien paraFm(t). tmetro+12Fmetro(t)γ(z,X)Fmetro(t)
Ondřej Čertík
@ OndřejČertík OK, lo siento, mi error, no revisé su código antes de hacer mi comentario.
Ali

Respuestas:

9

La integral en cuestión también se conoce como la función Boys, después del químico británico Samuel Francis Boys, quien introdujo su uso a principios de la década de 1950. Hace unos años, necesitaba calcular esta función con doble precisión, lo más rápido posible pero con precisión. I logrado alcanzar un error relativo del orden de a través de todo el dominio de entrada.10-15

En general, es ventajoso utilizar diferentes aproximaciones para argumentos pequeños y grandes, donde el cambio óptimo entre "grande" y "pequeño" se determina mejor experimentalmente, y en general es una función de metro . Para mi código, definí argumentos "pequeños" como aquellos que satisfacen la condición .unametro+112

Para argumentos grandes, calculo

Fmetro(una)=12γ(metro+12,una)×pag×pag,  pag=una-12(metro+12)

Este orden de operaciones evita el desbordamiento prematuro. Como aquí solo necesitamos la función gamma incompleta inferior de las órdenes de medio entero en lugar de una función gamma incompleta inferior completamente general, es ventajoso calcular desde una perspectiva de rendimiento

γ(metro+12,una)=Γ(metro+12)-Γ(metro+12,una)

utilizando valores tabulados de y computaciónΓ(m+1Γ(metro+12)acuerdo con esta respuesta, evitando cuidadosamente el problema de la cancelación sustractiva mediante el uso de una operación de suma múltiple fusionada. Una posible optimización adicional es observar que para a lo suficientemente grandea,γ(m+1Γ(metro+12,una)unadentro de una precisión de punto flotante dada.γ(metro+12,una)=Γ(metro+12)

Para pequeños argumentos, comencé con una expansión en serie para la función gamma incompleta inferior de

A. Erdelyi, W. Magnus, F. Oberhettinger y FG Tricomi, "Funciones trascendentales superiores, vol. 2". Nueva York, NY: McGraw-Hill 1953

y lo modificó para calcular la función Boys siguiente manera (truncando la serie cuando el término es lo suficientemente pequeño para una precisión dada):Fmetro(una)

Fmetro(una)=121metro+12Exp(-una)(1+norte=1unanorte(1+metro+12)× ... ×(norte+metro+12))

metro=0 0,1,2,3F0 0(una)=π4 4unamirF(una)mirFERFerferff

metro=1,2,3una<212Fmetro(una)=12una((2metro-1)Fmetro-1(una)-Exp(-una))

unametrometroFmetro-1=12metro-1(2una Fmetro(una)+Exp(-una))

njuffa
fuente
Gracias @njuffa por la gran respuesta. Si crea su código para este código abierto, creo que sería muy útil para muchas personas.
Ondřej Čertík
1
Actualmente, una implementación de CUDA del algoritmo descrito está disponible para su descarga gratuita desde el sitio web del desarrollador de NVIDIA (requiere registro gratuito como desarrollador de CUDA, aprobación generalmente dentro de un día hábil). El código está bajo una licencia BSD, que debería ser compatible con casi cualquier tipo de proyecto.
njuffa
4

Echaría un vistazo al libro de Abramowicz & Stegun, o la nueva revisión que NIST ha publicado hace un par de años y creo que está disponible en línea. También discuten formas de implementar cosas de manera estable.

Wolfgang Bangerth
fuente
Estaba usando esto: dlmf.nist.gov/8 , cuando lo implementé, pero ese es probablemente otro recurso. El capítulo 5 de Recetas numéricas también tiene información interesante, pero solo aplicable a funciones de una variable.
Ondřej Čertík
No creo que encuentres nada mucho más reciente que su referencia de 2001; SLATEC será más viejo que eso.
Geoff Oxberry
1

No parece ser el estado de la técnica , pero SLATEC a ofertas NETLIB "1400 de propósito general rutinas matemáticas y estadísticas." El Gamma incompleto está disponible en las funciones especiales aquí. .

Implementar tales funciones lleva mucho tiempo y es propenso a errores, por lo que no lo haría yo a menos que sea absolutamente necesario. SLATEC ha existido desde hace bastante tiempo y se usa ampliamente, al menos en función de los recuentos de descargas , por lo que esperaría que la implementación esté madura.

Ali
fuente