¿Cuál es el algoritmo más rápido para encontrar el logaritmo natural de un número grande?

Question

Como para $\pi$ tenemos un algoritmo/serie infinita que puede darnos los primeros 50 decimales en unos 3 términos. Así que si no fuera a calcular como $\ln(25551879\cdots)$ (un número entero realmente enorme, muy probablemente un primo), hasta 100 decimales, ¿cuál será el algoritmo que debo utilizar o que se utiliza en todo el mundo y qué eficacia tiene? Sé que la serie de Taylor es bastante lento en su trabajo, por lo que cualquier otro algoritmo en el que se calcula?

Answer 1

Digamos que necesita una tolerancia absoluta de $2^{-m}$ para la respuesta.

Dado un número de la forma $x=a \cdot 2^n$ , $a \in (1/2,1]$ , escriba $\ln(x)=\ln(a)+n\ln(2)$ .

Ahora calcule $\ln(a)$ tomando $m$ términos de la serie Maclaurin de $\ln(1+x)$ con $x=a-1$ y calcular $\ln(2)$ como $-\ln(1/2)$ tomando $m \lceil \log_2(n) \rceil$ términos de la serie Maclaurin de $\ln(1+x)$ con $x=-1/2$ .

Esta forma es un poco quisquillosa en cuanto a trabajar con números decimales frente a números binarios, pero tiene la ventaja de que el $\ln(a)$ converge en el peor de los casos como $2^{-m}$ en lugar de como $(9/10)^m$ como lo hace el enfoque análogo con el decimal. Tiene el inconveniente de que hay que precalcular $\ln(2)$ a una mayor precisión ya que $n$ será mayor, pero eso no importa mucho porque no es un problema "vivo" (siempre que se aplique algún tope al tamaño de la entrada y al tamaño de su recíproco).

Por lo general, no es así como se implementan las funciones de biblioteca en lenguajes de programación como C. Véase, por ejemplo, e_log.c en http://www.netlib.org/fdlibm/ . Esto comienza con una reducción del argumento similar a la que sugerí anteriormente (donde el límite inferior y el límite superior para $a$ difieren en un factor de $2$ ), entonces convierte el problema en $\ln(1+x)=\ln(1+y)-\ln(1-y)$ donde $y=\frac{x}{2+x}$ . Esta conversión conduce a una cierta aceleración de la serie, ya que la serie de la diferencia sólo tiene potencias Impares, y como $|y|<|x|$ . (A partir de aquí se podría proceder con un enfoque de series de Taylor. Si lo hiciera, utilizaría alrededor de $m/4$ términos, debido a las mencionadas cancelaciones y al hecho de que $y$ está en el entorno de $x/2$ teniendo en cuenta que la reducción de argumentos ya se ha realizado).

A continuación, utilizan un polinomio minimax para aproximar $\frac{\ln(1+y)-\ln(1-y)-2y}{y}$ . Este tipo de enfoque es el que suelo ver cuando compruebo el código fuente de las implementaciones rápidas de las funciones de la biblioteca estándar. Los coeficientes de este polinomio minimax fueron probablemente relativamente caros de calcular, pero de nuevo no es un problema "vivo" por lo que su velocidad no importa mucho.

Answer 2

Se trata esencialmente de un debate más profundo sobre la eficacia y la precisión de los distintos métodos.

---

Reducir el argumento para que esté cerca $1$ :

Fundamentalmente, la mayoría de las respuestas persiguen el mismo objetivo: reducir los argumentos a valores pequeños y utilizar la expansión de Taylor para $\ln(x)$ . Hasta ahora hemos visto 3 enfoques:

1) Factorizar una potencia de 2 y utilizar $\ln(a\cdot2^n)=\ln(a)+n\ln(2)$ .

2) Factorizar una potencia de 10 y utilizar $\ln(a\cdot10^n)=\ln(a)+n\ln(10)$ .

3) Reducir por medio del enraizamiento al cuadrado utilizando $\ln(x)=2\ln(\sqrt x)$ .

Podemos notar que el enraizamiento cuadrado del argumento reduce el argumento mucho más rápido que los otros métodos, que dividen el argumento por una constante repetidamente, ya que $\sqrt x\ll x/10<x/2$ para grandes $x$ . Siendo realistas, si tu entrada no tiene más de, digamos, 1000 dígitos, entonces sólo tienes que hacer la raíz cuadrada unas 10 veces en el peor de los casos. Sin embargo, esto tiene el coste de tener que realizar más cálculos para encontrar la raíz cuadrada, que no son fáciles. Por otro lado, realizar las divisiones es extremadamente fácil. Debido a la naturaleza de cómo escribimos/almacenamos los números, todas las divisiones pueden ser calculadas a la vez simplemente moviendo el punto decimal. A continuación, podemos truncar fácilmente hasta la precisión deseada. Con el enraizamiento cuadrado, el error es más difícil de manejar, y como el registro se multiplica por 2 cada vez, el error absoluto se multiplica por 2 también.

Por supuesto, la elección de escribir el argumento como múltiplo de una potencia de 2 o de una potencia de 10 depende de si lo hace un ordenador o un humano. Es probable que prefiera trabajar en base 10.

También está la cuestión adicional de cuál es nuestra gama de $a$ debe ser. Como queremos que sea lo más cercano a 1, podemos hacer algo de álgebra y ver. Para potencias de 2, queremos $a\in(a_0,2a_0]$ tal que $2a_0-1=1-a_0$ . Si se resuelve esto, se obtiene $a\in[\frac23,\frac43]$ . Para las potencias de 10, queremos $a\in[\frac2{11},\frac{20}{11}]$ .

---

Utilizando expansiones en serie:

A partir de aquí podría utilizar la expansión estándar de Taylor para el logaritmo natural:

$$\ln(a)=\sum_{k=1}^\infty\frac{(-1)^{k+1}}k(a-1)^k=(a-1)-\frac{(a-1)^2}2+\frac{(a-1)^3}3-\frac{(a-1)^4}4+\mathcal O(a-1)^5$$

sin embargo esto no converge tan rápido como se podría lograr realizando una expansión de Taylor más cercana a $a$ . Lo anterior viene dado por la definición integral del logaritmo natural y la expansión de Taylor del integrando en $1$ :

$$\ln(a)=\int_1^a\frac{\mathrm dt}t=\sum_{k=0}^\infty(-1)^k\int_1^a(t-1)^k~\mathrm dt$$

Pero podemos mejorar esto ampliando en medio de $1$ y $a$ :

\begin{align}\ln(a)=\int_1^a\frac{\mathrm dt}t&=\sum_{k=0}^\infty(-1)^k\left(\frac2{a+1}\right)^{k+1}\int_1^a\left(t-\frac{a+1}2\right)^k~\mathrm dt\\&=\sum_{k=0}^\infty\frac{(-1)^k}{k+1}\left(\frac{a-1}{a+1}\right)^{k+1}\left(1-(-1)^{k+1}\right)\\&=\sum_{k=0}^\infty\frac2{2k+1}\left(\frac{a-1}{a+1}\right)^{2k+1}\end{align}

Para $a$ cerca de $1$ lo anterior tiene un error aproximado de $\mathcal O((a-1)/2)^{2n+1}$ cuando se utiliza $n$ términos. Una derivación algebraica de lo anterior es proporcionada por Wikipedia pero no muestra realmente lo rápido que este converge. Como estamos dos veces más cerca del límite más lejano de la integral, ganamos un dígito binario adicional de precisión por término. Pero como la mitad de los términos desaparecen, esto significa que básicamente podemos calcular el doble de dígitos por término. Este es el método mencionado por Respuesta de Ian .

Aquí es un programa en Ruby que calcula el logaritmo usando series.

---

Utilizando métodos de búsqueda de raíces:

Mientras que los métodos en serie son realmente agradables y convergen decentemente rápido, Wikipedia proporciona dos métodos más para una evaluación aún más precisa. El primero lo proporciona Eric Towers y supone el cálculo del logaritmo mediante funciones exponenciales. Por supuesto, como el cálculo es barato y la precisión no se ve afectada tanto, recomendaría reducir el argumento para que vuelva a estar cerca de $1$ . Esto significa que $y$ tal y como se define a continuación, se acercará a $0$ permitiendo un cálculo más rápido de la exponencial. Esto también significa que podemos utilizar simplemente $y_0=0$ como nuestra suposición inicial.

$$y=\ln(x)\Rightarrow x=\exp(y)\Rightarrow x-\exp(y)=0$$

En la que podemos aplicar métodos estándar de búsqueda de raíces, como Método de Newton (duplica los dígitos precisos por paso) o Método Halley (triplica los dígitos precisos por paso).

El cálculo de las funciones exponenciales puede realizarse mediante la expansión de Maclaurin:

$$\operatorname{exmp1}(y)=\exp(y)-1=\sum_{n=1}^\infty\frac{y^n}{n!}=y+\frac{y^2}2+\mathcal O(y^3)$$

Desde $y$ está cerca $0$ hay un gran error de punto flotante en el cálculo $\exp(y)$ que tiene un término dominante de $1$ Así que utilizamos $\operatorname{expm1}(y)$ para sortear esto.

También hay que tener en cuenta que como $\Delta y_n\to0$ es más fácil calcular $\exp(\Delta y_n)$ que $\exp(y_{n+1})$ directamente, y utilizar $\exp(y_{n+1})=\exp(\Delta y_n)\exp(y_n)=\exp(y_n)+\exp(y_n)\operatorname{expm1}(\Delta y_n)$ . Esto reduce la exponenciación principal al primer paso, que es trivial ya que $\exp(0)=1$ .

Dejemos que $y_0=0$ y $\operatorname{expy}_0=1$ .

Para el método de Newton, sea $\displaystyle\Delta y_0=x\operatorname{expy}_0-1$ y:

\begin{cases}\Delta y_n=x\operatorname{expy}_n-1,\\\operatorname{expy}_{n+1}=\operatorname{expy}_n+\operatorname{expy}_n\operatorname{expm1}(-\Delta y_n),\\y_{n+1}=y_n+\Delta y_n\end{cases}

Para el método de Halley, dejemos $\displaystyle\Delta y_0=2\cdot\frac{x-\operatorname{expy}_0}{x+\operatorname{expy}_0}$ y:

\begin{cases}\displaystyle\Delta y_n=2\cdot\frac{x-\operatorname{expy}_n}{x+\operatorname{expy}_n},\\\operatorname{expy}_{n+1}=\operatorname{expy}_n+\operatorname{expy}_n\operatorname{expm1}(\Delta y_n),\\y_{n+1}=y_n+\Delta y_n\end{cases}

Aquí es un programa en Ruby que calcula el logaritmo con el método de Newton y aquí es un programa en Ruby que calcula el logaritmo con el método de Halley.

---

Utilizando el agm:

El media aritmética-geométrica es una poderosa herramienta que se puede utilizar aquí para calcular rápidamente el logaritmo, así como $\pi$ y ciertas integrales. Se define como:

$$a_0=a,b_0=b\\a_{n+1}=\frac{a_n+b_n}2,b_{n+1}=\sqrt{a_nb_n}\\M(a,b)=\lim_{n\to\infty}a_n=\lim_{n\to\infty}b_n$$

Según Wikipedia, esto es tan rápido y barato de calcular que se puede utilizar para calcular la función exponencial utilizando logaritmos más rápido que las series que aproximan la función exponencial. Para conseguir $p$ de precisión, tome una $m$ para que $s=x2^m$ es mayor que $2^{p/2}$ . Podemos entonces calcular el logaritmo natural como:

$$\ln(x)=\lim_{m\to\infty}\frac{\pi x2^{m-2}}{2M(x2^{m-2},1)}-m\ln(2)$$

que es esencialmente un replanteamiento de la tasa de convergencia de $M(1,t)$ como $t\to\infty$ . Para este método, ni siquiera es necesario reducir el argumento. De hecho, ¡podemos aplicarlo directamente a los argumentos grandes!

Sin embargo, este método tiene algunos inconvenientes. El cálculo requiere que calculemos algunas raíces cuadradas en flotadores grandes, pero esto se puede manejar con una clase de flotadores especialmente definida y las funciones respectivas.

Alternativamente, por supuesto, uno podría simplemente reducir el argumento para evitar grandes flotadores como antes.

Aquí es un programa en Ruby que calcula el logaritmo utilizando la media aritmética-geométrica.

Answer 3

No sé cuál es la forma más rápida, pero aquí hay un enfoque razonablemente eficiente:

• Puedes dividir números con Newton-Raphson ;
• Una vez que sepas cómo hacerlo, también se pueden sacar raíces cuadradas con Newton-Raphson;
• Puede utilizar $\ln x=2\ln\sqrt{x}$ tantas veces como sea necesario para que el argumento del logaritmo se acerque a $1$ antes de utilizar la serie Maclaurin de $\ln(1+x)$ ;
• Si necesita logaritmos en otra base, utilice $\log_ax=\frac{\ln x}{\ln a}$ .

Answer 4

Método Halley es iterativo y su convergencia es cúbica. Aplicado aquí, invertiríamos para usar la exponencial (que felizmente es su propia derivada): $$ y_{n+1} = y_n+2 \cdot \frac{x- \mathrm{e}^{y_n}}{x + \mathrm{e}^{y_n}} \text{.} $$ Por ejemplo, con $x = 25551879$ y $y_0 = 2$ (es decir, no está cerca), los iterados (todos calculados con intermedios de 15 dígitos) son $2$ , $4.$ , $5.99999$ , $7.99993$ , $9.99946$ , $11.996$ , $13.9708$ , $15.7959$ , $16.9122$ , $17.056$ , $17.0562$ . Mi punto no es que 15 dígitos sean suficientes, sino que el método alcanzó la precisión mostrada en sólo diez pasos.

Te preguntarás, ¿cómo consigo esos exponenciales? Utiliza la serie de potencias. Eso converge rápidamente para cualquier argumento que puedas ver. Como regla general, empieza con el doble de términos que tu argumento, así que para $\mathrm{e}^{17.0562}$ , empezar con $34$ términos en esta serie de Taylor. Esto da $2.5549{\dots}\times 10^{7}$ con error $2355.61{\dots}$ . A continuación, aumenta el número de términos de los exponenciales en (en este caso) $34$ siempre y cuando su estimación para $y_{n+1}$ todavía cambia dentro de la precisión de su objetivo. Cuando eso deje de suceder, tómalo como tu final $y_{n+1}$ y repite el proceso de ampliación de la serie exponencial hasta que tu $y_{n+2}$ se estabiliza. Continúe hasta que obtenga dos $y$ s en una fila que coinciden con su precisión objetivo (más suficientes bits adicionales que no cambian que al menos uno de ellos es un cero para que usted sabe qué manera de redondear el último bit en su respuesta reportada).

Answer 5

Bueno $$ \ln(25551879) = \ln(0.25551879 \times 10^{8}) $$ $$= \ln(0.25551879) + \ln(10^8) $$ $$= 8 \times \ln(10) + \ln(0.25551879) $$

Desde $\ln(10)$ es una constante que puede ser precalculada con un gran número de decimales, sólo necesitamos un método que converja rápidamente para valores inferiores a $1.0$ . No sé si Taylor es lo suficientemente bueno en ese rango restringido o si hay algún otro método mejor.

Esto aborda la cuestión que has planteado acerca de que el argumento es un número grande. En cuanto a la generación de muchos dígitos, hay un montón de buenas respuestas en esta pregunta anterior .