Dado que $ x ^ r + y ^ r = z ^ r $ se puede transformar en una ecuación trinomial,
podríamos darle crédito a Lambert, ya que desarrolló funciones para ella en 1758, consulte ejemplo este artículo aquí de Corless, et.al.

Rápidamente tenemos $ x ^ r * (1+ (y / x) ^ r) = x ^ r * ( z / x) ^ r $.
Establezca $ p = ln (y / x) $, $ q = ln (z / x) $ y $ s = e ^ r $. Luego obtenemos $ 1 + s ^ p = s ^ q $.

Configurando más $ w = s ^ q $ y $ t = p / q $, obtenemos la ecuación trinomial
$ - 1 + w = ​​w ^ t $, que tenemos que resolver para $ w $ para un $ t $ dado.

Editar I:
El documento mencionado anteriormente menciona un desarrollo de la serie por Euler,
que se basa en Lambert, que dice lo siguiente:

$ w ^ n = 1 + n * v + \ frac {1} {2} * n * (n + α + β) * v ^ 2 $
$ \ quad \ quad + \ frac {1} {6} * n * (n + α + 2 * β) * (n + 2 * α + β) * v ^ 3 $
$ \ quad \ quad + \ frac {1} {24} * n * (n + α + 3 * β) * (n + 2 * α + 2 * β) * (n + 3 * α + β) * v ^ 4 $
$ \ quad \ quad + \, etc. $

La forma aproximada coincide con la dada por el OP. Pero con respecto a una verificación completa, es decir, hacer coincidir los parámetros α, β, vy colocarlos
correctamente, estaba demasiado ocupado hasta ahora e hice una pregunta secundaria aquí.

Edición II:
La ecuación que resuelve Euler es $ w ^ α-w ^ β = (α-β) * v * w ^ {α + β} PS Podemos
elegir $ α = t-1 $, $ β = t $, $ v = 1 $ y obtenemos $ w ^ {t-1} -w ^ t = (- 1) * w ^ {2 * t-1} $.

Si multiplicamos ambos lados por $ (- 1) * w ^ {- t + 1} $ obtenemos $ -1 + w = ​​w ^ t $ según sea necesario.
Podemos elegir $ n = 1 $ y obtendrá un desarrollo por $ w $ que está a la par
la función $ f $ del OP.