¡Buenas!

Creo que lo mas intuitivo es trabajar con la funcion inversa y aplicar Bolzano:

x = log_a(b) -> a^x = b -> a^x - b = 0

Por lo tanto se trata de encontrar el cero de la función f(x) = a^x - b

Aplicando un poco de analisis sabemos que a^x es siempre positivo (siempre que a sea positivo), por lo tanto tal cero existirá si y solo si b es positivo.

En estas condiciones tendrás que buscar dos puntos x₀ y₀ en los que la función tenga signo distinto. Así sabrás que el cero de la función se encuentra entre estos dos puntos.

Solo te queda iterar. En cada paso tendrás que buscar el punto medio del intervalo y evaluarlo para saber si la función toma un valor positivo o negativo en dicho punto. Una vez sepas el valor de la función al evaluarlo en el punto medio podrás sustituir la cota superior o inferior que tenias inicialmente por el punto medio.

Este proceso parará cuando:
- La evaluación de la función en el punto medio de exactamente cero.
- La diferencia entre la cota superior e inferior sea menor a un valor que tu consideres un error aceptable para la solución del problema.
- Cuando entre una iteración y la siguiente no cambie el punto medio (esto sucede debido al incremento mínimo que existe en la representación interna de números en coma flotante)

El valor que devolverá la función logaritmo será el punto medio del intervalo en el que han acabado las iteraciones, ya que será el valor de x que ha hecho que la función f(x) = a^x - b se aproxime mas a cero (el valor mas proximo a log_a(b))

¡Saludos!

Implementando de manera rápida lo que comentó @diskontrol que como he podido comprobar solo es válido para logaritmos en base 10, he sacado este código (seguro que es un churro, no me mateis    ). En cuanto lo último comentado por @Puntoinfinito eso solo vale para logaritmos exactos. ¿@do-while al fin y al cabo es ensayo y error no? Si yo tengo log₂4 -> 2^x = 4 vas probando valores enteros hasta que tengas un resultado mayor que 4 o igual y vas reduciendo ... no sé tampoco tengo nivel matemático para hacerlo, implemente la serie de Taylor pero es muy lenta, dejo el código al que @diskontrol se refería:

#include <iostream>
#include <sstream>
#include <string>
 
using namespace std;
 
#define PRECISION 1000
 
#define e 2.7182818284590452353602874713526624977572470936999595749669676277240766303535
 
double ln (double num);
double log10 (double num, std::string &str);
double pow (double base, register int exp);
 
 
int main(int argc , char * argv [])
{
    string str;
    cout <<"Resultado en la precision maxima de de double es: "<< log10(34.4,str) << endl;
    cout << "Resultado truncado a "<< PRECISION << " decimales es: " << endl<< str;
    return 0;
}
 
double pow (double base, register int exp)
{
    double ret = 1;
    while(exp!=0)
    {
        ret = base*ret;
        --exp;
    }
    return ret;
}
double ln (double num)
{
    string unused;
    return (log10(num,unused)/log10(e,unused));
}
double log10 (double num,std::string &str)
{
    double ret = 0;
    double decimal  = 0;
 
    if(num < 10)
        num = pow(num,10);
    else
    {
        while(num >= 10)
        {
            num /= 10;
            ++ret;
        }
        num = pow(num,10);
    }
 
    ostringstream convert;   // stream usado para la conversión.
    convert <<ret;
    str += convert.str();
    str += '.';
    for(register int i = 1 ; i != PRECISION; ++i)
    {
        while(num >= 10)
        {
            num /= 10;
            ++decimal;
        }
        ret += decimal*(1/(pow(10,i)));
        num = pow(num,10);
 
        str +=(char) decimal+48;
        decimal = 0;
    }
    return ret;
}

XD, si y no. Se puede considerar una especie de ensayo y error, pero aplicando la lógica. Después del primer paso conseguimos dos puntos x₀, y₀ en los que la función evaluada en el primero punto, en este caso en concreto, es menor que cero y evaluada en el segundo es mayor que cero:

                                                                      f(y₀)



x₀------------------------------------------------------y₀ (eje OX)


f(x₀)

Como la función es continua, si unes f(x₀) y f(y₀) con una linea, siempre tendrás que atravesar el eje OX. Es decir, en algún punto entre x₀ e y₀ la función sera cero. Ese es el punto que buscamos y sabemos que está entre ambos valores.

Por lo tanto obtenemos el punto z₀ = (x₀ + y₀) / 2 y evaluamos la función en dicho punto. Imagina que sale que la función es positiva en z₀:

                                                                      f(y₀)

                                    f(z₀)

x₀---------------------------z₀--------------------------y₀ (eje OX)


f(x₀)

Sigue sucediendo lo mismo de antes, entre x₀ y z₀ la función se anula, pero ahora tenemos un intervalo mas pequeño que el primero. Como z₀ da un valor positivo, hacemos x₁ := x₀ e y₁ := z₀. Y aplicamos el proceso anterior al intervalo [x₁,y₁] para obtener un punto z₁ que estará todavía mas próximo al cero de la función.

A fin de cuentas se trata de eso. De aproximar en cada paso cada vez mas el cero de la función dividiendo el intervalo de búsqueda en cada iteración. Si consideramos un error e, se tendrá que:

(y₀ - x₀) / 2ⁿ < e si y solo si

y₀ - x₀ < e 2ⁿ si y solo si

(y₀ - x₀) / e < 2ⁿ si y solo si

log₂((y₀ - x₀) / e) < n

Es decir, como mucho tendremos que iterar una vez mas que la parte entera de log₂((y₀ - x₀) / e) para obtener un resultado que se aproxime a la raíz con una precisión mayor que e/2.

¡Saludos!