Approximate solutions to non polynomial equations

Just an idea how to solve the problem without specifying d1[p],d2[p] (Perturbation theory, MMA version <12 )

First assume the solution p can be approximated by p -> p0 + p1 edr + p2 edr^2 + p3 edr^3 (* edr=e/r *)

Substituting into the equation and series expansion

Clear[d1,d2]
ser = Series[-F + Sqrt[p^2 + m^2] + (edr d1[p]) + (edr^2 d2[p]) /. 
p ->  p0 + p1 edr + p2 edr^2 + p3 edr^3, {edr, 0, 2}] // Normal
(*-F + Sqrt[m^2 + p0^2] + edr ((p0 p1)/Sqrt[m^2 + p0^2] + d1[p0]) +edr^2 ((m^2 p1^2)/(2 (m^2 + p0^2)^(3/2)) + (m^2 p0 p2)/(m^2 + p0^2)^(3/2) + (p0^3 p2)/(m^2 + p0^2)^(3/2) + d1[p0] + p1 Derivative[1][d1][p0])*)

The series coefficients can be solved for p0,p1,p2

CoefficientList[ser, edr] // FullSimplify[#, F > m > 0] &
Simplify[Solve[% == 0, {p0, p1, p2}], F >m> 0]
(*{{
p0 -> -Sqrt[F^2 - m^2], 
p1 -> (F d1[-Sqrt[F^2 - m^2]])/Sqrt[F^2 - m^2], 
p2 -> (1/(2 (F^2 - m^2)^2))(m^2 Sqrt[F^2 - m^2] d1[-Sqrt[F^2 - m^2]]^2 +2 F (F^2 - m^2)^(3/2) d2[-Sqrt[F^2 - m^2]] +2 F^2 (F^2 - m^2) d1[-Sqrt[F^2 - m^2]]Derivative[1][d1][-Sqrt[F^2 - m^2]])
}, {
p0 -> Sqrt[F^2 - m^2], 
p1 -> -((F d1[Sqrt[F^2 - m^2]])/Sqrt[F^2 - m^2]), 
p2 -> -(1/(2 (F^2 - m^2)^2))(m^2 Sqrt[F^2 - m^2] d1[Sqrt[F^2 - m^2]]^2 +2 F (F^2 - m^2)^(3/2) d2[Sqrt[F^2 - m^2]] +2 F^2 (-F^2 + m^2) d1[Sqrt[F^2 - m^2]] Derivative[1][d1][Sqrt[F^2 - m^2]])
}}*)

I'll make an example of the functions $d_i(p)$:

d1[p_] = Sin[p];
d2[p_] = p^2 - 7;

Series-expansion of the solution for small $e$: I'll use $E=$F because the symbol E is already in use. The order of the series-expansion is set to 2 here but can be increased at will,

Assuming[F > 0, 
  AsymptoticSolve[F == Sqrt[p^2 + m^2] + (e d1[p])/r + (e^2 d2[p])/r^2,
    {p, Sqrt[F^2 - m^2]}, {e, 0, 2}] // FullSimplify]

{{p -> (1/( 4 ((F - m) (F + m))^(3/2) r^2))(-e^2 (m^2 + 4 F (F - m) (F + m) (-7 + F^2 - m^2)) + 4 (F^2 - m^2)^2 r^2 + e^2 m^2 Cos[2 Sqrt[F^2 - m^2]] + 4 e F (-F^2 r + m^2 r + e F Sqrt[F^2 - m^2] Cos[Sqrt[F^2 - m^2]]) Sin[Sqrt[ F^2 - m^2]])}}

This is the solution that, for $e=0$, passes through $p_0=\sqrt{F^2-m^2}$. There is also the other branch that passes through $p_0=-\sqrt{F^2-m^2}$; choose the branch in the AsymptoticSolve input.

For the more complex problem in the comments, it's sufficient to extend the assumptions a bit:

Assuming[ma > 0 && mb > 0 && F > ma + mb, 
  AsymptoticSolve[
    F == Sqrt[p^2 + ma^2] + Sqrt[p^2 + mb^2] + (e d1[p])/r + (e^2 d2[p])/r^2,
    {p, Sqrt[(F - ma - mb) (F + ma - mb) (F - ma + mb) (F + ma + mb)]/(2 F)},
    {e, 0, 2}] // Refine]

(lengthy output)