Values for $a$ and $b$ as polynomials in $x,y$, when they exist, correspond to factorization of $\frac{x^p+y^2}{x+y}$ over the imaginary field $K_p:=\mathbb Q(\sqrt{-p})$ of the form:
$$\frac{x^p+y^p}{x+y} = (a(x,y)+\sqrt{-p}b(x,y))\cdot (a(x,y)-\sqrt{-p}b(x,y)).$$
For $p\equiv1\pmod4$, we can similarly work in $K_p:=\mathbb Q(\sqrt{p})$ and look for factorization of the form:
$$\frac{x^p+y^p}{x+y} = (a(x,y)+\sqrt{p}b(x,y))\cdot (a(x,y)-\sqrt{p}b(x,y)).$$

[This Sage code](https://sagecell.sagemath.org/?z=eJytkstqwzAQRfeB_MN0UTJKXcdOd6bOMvt2m9phohcqtixkGQwh_14ZN30sSqFUmxHM1bmXGQmpQMtwpBM6ViwXEM9z-jjuoISngYSnYPjeyEYg3uesRnSxbDZbtnbsUM0P9lGsiIfOI-BYu7tJgWMswGYF9b30AahpEHTakjvyTipluJE29NhQexIEoYCQamo608e-fR00BYmM3ZQaOg_idluWGah41YkAY6PxO9_HBM53ItJrFFO87zKj_mBbRlsGV4P1p4P4B_juC1uGwVvwPzJ4AfyQVSz5TZNXcR4zdLmYMropofOmlUdPVkt8SPIsuy46NmxAtTq7S_ESzh__4LJib9N1os8=&lang=sage&interacts=eJyLjgUAARUAuQ==) (updated 2024-06-30) factors $\frac{x^p+1}{x+1}$ over $K_p$ and converts the result into the required form, and reports the (univariate) polynomials $a(x,1)$ and $b(x,1)$. As an example, it reports $a$ and $b$ for all primes below $100$. 

Example for $p=23$:
$$a = (x + 1) \cdot (x^{10} - \frac{3}{2} x^{9} - x^{8} + 5 x^{7} - \frac{17}{2} x^{6} + \frac{21}{2} x^{5} - \frac{17}{2} x^{4} + 5 x^{3} - x^{2} - \frac{3}{2} x + 1),$$
$$b= \frac{1}{2} \cdot (x - 1) \cdot x \cdot (x^{8} + x^{5} - x^{4} + x^{3} + 1).$$