Use the Euler-MacLaurin formula,
$$\sum_{k=1}^\infty F(k)=\int_0^\infty F(k)\,dk+\tfrac{1}{2}[F(\infty)-F(0)]+\int_0^\infty (k-\text{Int}\,[k]-\tfrac{1}{2})F'(k)\,dk.$$
In this case $F(k)=\sqrt{\max\{1 -k^2/x^2,0\}}$, hence$^\ast$
$$f(x)=\sum_{k=1}^\infty F(k)=\frac{\pi x}{4}-\frac{1}{2}+(3-4\delta x)\frac{\sqrt {2\delta x}}{6\sqrt{x}}-\frac{0.075}{\sqrt{x}}-\frac{1}{12x^2}\,\text{Int}\,[x]+{\cal O}(1/x^{3/2}),$$
where $\text{Int}\,[x]$ denotes the integer part of $x$ and $\delta x=x-\text{Int}\,[x]$. (The coefficient $0.075$ is computed numerically -- can it be computed analytically?)
So the equation in the OP,
$$f(x) = - \frac{1}{2} + \frac{\pi}{4}x + \frac{p(x)}{\sqrt{x}},$$
with $p(x+1)=p(x)$ holds, but only up to corrections of order $1/x$.
The plot compares the exact $p(x)$ (blue) with expansion formula above (gold).
$^\ast$ The integral of $k-\text{Int}\,[k]-\tfrac{1}{2}F'(k)$ over an interval $(n,n+1)$ contributes $-\tfrac{1}{12}x^{-2}$ for large $x$, and there are $\text{Int}\,[x]$ of these intervals, hence the term of order $1/x$ that breaks the periodicity of $p(x)$. The interval $(\text{Int}\,[x],x)$ contributes the term $\tfrac{1}{6}(3-4\delta x)\sqrt{2\delta x/x}$ for large $x$. The term $-0.075x^{-1/2}$ comes from the intervals $(n,n+1)$ with $n$ of order $x$. This is the contribution which I was not able to computer analytically.