I now believe that my question (and suggestion that proof $(1)$ should have become standard before Lacroix) relied on the misconception that tangent was easier to differentiate than arctangent. In fact the calculus of inverse trigonometric functions took off earlier, as has been explained by G. Eneström (1905), C. Boyer (1947), or V. J. Katz (1987):
it was quite common [at first] to deal with what we call the arcsine function rather than the sine.
C. Wilson (2001, 2007) concurs and stresses that our trig functions with periodic graphs weren’t much seen or differentiated until Euler “found” them to solve 2nd order linear differential equations (1741); so much so that he still wrote in (1749, p.15):
as this way of operating is not yet commonplace, it will be apropos to warn that the differentials of the formulas $\sin.𝜑$ : $\cos.𝜑$ : $\mathrm{tang}.𝜑$ : $\cot.𝜑$ are $d𝜑\,\cos.𝜑$ : $-d𝜑\,\sin.𝜑$ : $\smash[b]{\frac{d𝜑}{\cos.𝜑\,^2}}$ & $\smash[b]{-\frac{d𝜑}{\sin.𝜑\,^2}}$
— and e.g. in (1796, p.163) L’Huilier still computed $\tan'$ from $\arctan'$ rather than vice versa.
In this vein, $d𝜃 = \frac{dt}{1+t^2}$ was not proved like $(1)$, but by a differential triangle argument similar to $(*)$ but simpler and attributed to Cotes (Aestimatio errorum, 1722): in modern notation, parametrize the unit circle with $(x,y)=\frac{(1,\,t)}{\sqrt{1+t^2}}$ and obtain
$$
d\theta
=\sqrt{\smash{dx^2+dy^2}\vphantom{a^2}}
=\frac{dt}{1+t^2}
\tag6
$$
($=\mathrm{CE}$ in Cotes’ figure, which became standard in many books even before his own — the list could almost be described as “everyone but Euler”):
1708 Charles-René Reyneau §590 fig. 41
1718 John Craig pp.52–54
1722 Roger Cotes (posthumous) Lemma II
1730 Edmund Stone p.63 fig. 13
1736 James Hodgson p.230
1736 John Muller §247 fig. 153
1737 Thomas Simpson §143
1742 Colin MacLaurin §195 fig. 52
1743 William Emerson pp.171–172 fig. 76
1748 Maria Agnesi p.639 fig. 4
1749 Charles Walmesley (credits Cotes) pp.3,53 fig. 10
1749 William Emerson p.29 fig. 6
1750 Thomas Simpson §142
1754 Louis-Antoine de Bougainville p.24 fig. 9
1761 Abraham Kästner §299 fig. 18
1765 Jean Le Rond D’Alembert p.640 fig. 25
1767 Étienne Bézout p.146 fig. 46
1768 Thomas Le Seur & François Jacquier p.63 fig. 7
1774 Jean Saury pp.25,63 fig. 3
1779 Samuel Horsley pp.298–299
1786 Simon L’Huilier pp.103–104 fig. 20
1795 Simon L’Huilier §76 fig. 17
As to our usual proof $(1)$, it appears before Lacroix in 18-year-old Legendre’s Theses mathematicæ (1770), then in a book by their common teacher J.-F. Marie and several others:
1770 Adrien-Marie Legendre pp.10,16
1772 Joseph-François Marie §904
1777 Jacques-Antoine Joseph Cousin p.81
1781 Claude Bertrand p.140
1781 Louis Lefèvre-Gineau p.31
1795 Simon L’Huilier §76
1797 Sylvestre-François Lacroix pp.113–114
1801 Joseph-Louis Lagrange p.81
1810 Sylvestre-François Lacroix pp.lii,203–204
As noted by the OP we can replace $f$ by $af(bx)$ for suitable $a,b\in\mathbb{R}$ so that wlog we can take $c=1$ and ensure that $\sup f=-\inf f=1$.
Firstly we note that $f(z)$ is infinitely differentiable on $\mathbb{R}$ so we can form the taylor series at 0, $f(z)=\sum_{i=0}^{\infty}\frac{f^{(i)}(0)}{i!}z^i$. Since $|f^{(i)}(0)|\leq 1$ for all $i\geq 0$ we have $\sum_{i=0}^{\infty}|\frac{f^{(i)}(0)}{i!}z^i|=\sum_{i=0}^{\infty}\frac{|f^{(i)}(0)|}{i!}|z|^i\leq \sum_{i=0}^{\infty}\frac{|z|^i}{i!}$ which converges for all $z$ to $e^{|z|}$.
Hence $F(z)=\sum_{i=0}^{\infty}\frac{a_i}{i!}z^i$ is an absolutely convergent series defining an entire function on $\mathbb{C}$ agreeing with $f$ on $\mathbb{R}$ s.t. $\sup f^{(k)}=-\inf f^{(k)}=1$ for all $k\in \mathbb{Z}_{\geq0}$ and $|F(z)|\leq e^{|z|}$ for all $z\in \mathbb{C}$.
We now determine the form of $F$ given the preceding conditions.
First note that Bernstein proved the following (see Rahman and Tariq$^1$) as an extension of his related inequality for polynomials:
Theorem Let $g$ be an entire function of exponential type $\tau>0$ such that $|g(x)|\leq M$ on the real axis. Then $$\sup_{-\infty<x<+\infty}|g^{'}(x)|\leq M\tau.$$
Now $|F(z)|\leq e^{|z|}$ for all $z\in \mathbb{C}$ and therefore we know that $F$ is of exponential type 1. In addition $|F(x)|=|f(x)|\leq 1$ on the real axis.
We are therefore interested in the conditions of equality in the above.
Fortunately in their book "Analytic Theory of Polynomials"$^2$ Rahman and Schmeisser prove (Theorem 14.1.7) that equality holds in the above if and only if $g(z)=ae^{i\tau z}+be^{-i\tau z}$ where $|a|+|b|=M$. (You can access the relevant pages 513-514 on google books)
$|F(x)|\leq 1$ for $x\in\mathbb{R}$ and hence in the above theorem we may take $M=1$, $\tau=1$ and $g=F$ since $|F|$ is of exponential type 1. Also $\sup_{x\in\mathbb{R}}|F^{'}(x)|=1$ which means we have the case of equality.
Thus $F(z)=ae^{iz}+be^{-iz}$ where $|a|+|b|=M=1$.
On the real axis $F$ is real valued and agrees with $f$ so we must have $f(x)=ae^{ix}+be^{-ix}$, $x\in \mathbb{R}$ with $a=\bar{b}$. Setting $a=c+id$ we obtain $f(x)=2c\cos x-2d\sin x$ with $|a|=|b|=2\sqrt{c^2+d^2}=1$. Rewriting with $C=2c$, $D=-2d$, we obtain
$$f(x)=C\cos x+D\sin x$$
for some constants $C,D\in \mathbb{R}$, $C^2+D^2=1$.
This proves the OP's conjecture. Note that the condition $C^2+D^2=1$ is due to the bound we imposed on $f$ due to our normalisation which is not part of the definition of $S(c)$.
1 Rahman, Q. I.; Tariq, Q. M., On Bernstein’s inequality for entire functions of exponential type, J. Math. Anal. Appl. 359, No. 1, 168-180 (2009). ZBL1168.30002.
2 Rahman, Q. I.; Schmeisser, G., Analytic theory of polynomials, London Mathematical Society Monographs. New Series 26. Oxford: Oxford University Press (ISBN 0-19-853493-0/hbk). xiv, 742 p. (2002). ZBL1072.30006.
Best Answer
First we put it in the notation of Mathematica $K(k)$ is $K(k^2)$. So our function will be $$f(k)= k K(k^2)\sinh\Bigl(\frac{\pi}{2}\frac{K(1-k^2)}{K(k^2)}\Bigr).$$ Now we change variables (W486, Whittaker, Watson p.~486) . $$k=\frac{\vartheta_2^2(q)}{\vartheta_3^2(q)} \quad (*)$$ Where$\newcommand\Z{\mathbb{Z}}$ $$\vartheta_2(q)=2q^{\frac14}(1+q^2+q^6+\cdots)= \sum_{n\in\Z}q^{(n-\frac12)^2}= 2q^{\frac14}\prod_{n=1}^\infty\{(1-q^{2n})(1+q^{2n})^2\}$$ $$\vartheta_3(q)=1+2q+2q^4+2q^9+\cdots=\sum_{n\in\Z}q^{n^2}= \prod_{n=1}^\infty\{(1-q^{2n})(1+q^{2n-1})^2\}$$ $$\vartheta_4(q)=1-2q+2q^4-2q^9+\cdots=\sum_{n\in\Z}(-1)^nq^{n^2}= \prod_{n=1}^\infty\{(1-q^{2n})(1-q^{2n-1})^2\}$$ The function of $q$ in (*) is differentiable and increasing on $(0,1)$ it is $0$ in $0$ and $1$ in $1$. (we shall write $\vartheta_j$ to denote $\vartheta_j(q)$). Since (W467) $$\vartheta_2^4+\vartheta_4^4=\vartheta_3^4$$ we have $$1-k^2=1-\frac{\vartheta_2^4}{\vartheta_3^4}=\frac{\vartheta_4^4}{\vartheta_3^4}$$ The interesting thing about this change of variables is that $$K(k^2)=K\Bigl(\frac{\vartheta_2^4}{\vartheta_3^4}\Bigr)=\frac{\pi}{2}\vartheta_3^2,\qquad K(1-k^2)=K\Bigl(\frac{\vartheta_4^4}{\vartheta_3^4}\Bigr)=\frac{\log(1/q)}{2}\vartheta_3^2.$$
Now our function is $$k K(k^2)\sinh\Bigl(\frac{\pi}{2}\frac{K(1-k^2)}{K(k^2)}\Bigr) =\frac{\pi}{4}\Bigl(\frac{1}{\sqrt{q}}-\sqrt{q}\Bigr) \vartheta_2^2$$ that must be decreasing in $q$.
We add to the above some comments:
We need to show that $f(q):=(1-q)\vartheta_2^2/4\sqrt{q}$ is decreasing for $0 < q < 1 $. But
$$f(q)=(1-q)\prod_{n=1}^\infty (1-q^{2n})^2(1+q^{2n})^4= (1-q)\prod_{n=1}^\infty (1-q^{4n})^2(1+q^{2n})^2$$ This is the same to prove that the logarithmic derivative is negative $$-\frac{1}{1-q}-\sum_{n=1}^\infty\frac{8nq^{4n-1}}{1-q^{4n}}+\sum_{n=1}^\infty \frac{4nq^{2n-1}}{1+q^{2n}}$$ multiply this by $q>0$ and expand in series $$-\sum_{m=1}^\infty q^m-\sum_{n=1}^\infty \sum_{k=1}^\infty 8nq^{4nk}+ \sum_{n=1}^\infty\sum_{k=1}^\infty (-1)^{k+1}4nq^{2nk}$$
To show that this is negative observe that the only positive terms are those in the third sum with $k=2j+1$ odd, we will pair this term with that in the first sum corresponding to the same $n$ and $k=j$, these two terms are $$-8nq^{4nj}+4nq^{2n(2j+1)}=-4nq^{4nj}(2-q^{2n})<0.$$ This leaves only the terms with $j=0$ without pair. The remaining positive terms adds to $$\sum_{n=1}^\infty 4nq^{2n}=\frac{4q^2}{(1-q^2)^2}.$$ These terms we compensate with the first sum $$-\frac{q}{1-q}+\frac{4q^2}{(1-q^2)^2}=-\frac{q(1+q)(1-q^2)-4q^2}{(1-q^2)^2}$$ This is negative for $ 0 < q < 0.295598 $.
There is an intrinsic difficulty to treat larger values of $q$.
I propose to use the modularity of the theta function:
We have the equality considering $\vartheta_j$ as functions of $q$ $$\vartheta_2(e^{-\frac{\pi}{x}})=\sqrt{x}\vartheta_4(e^{-\pi x})$$ It follows that putting $q=e^{-\frac{\pi}{x}}$ $$\frac{1}{4}\Bigl(\frac{1}{\sqrt{q}}-\sqrt{q}\Bigr)\vartheta_2^2(e^{-\pi/x})= \frac{x}{2}\sinh\frac{\pi}{2x}\, \vartheta_4^2(e^{-\pi x}).$$ We must show this function is decreasing Since $$\vartheta_4(e^{-\pi x}) = \prod_{n=1}^\infty (1-e^{-2\pi n x})(1-e^{-\pi(2n-1)x})^2$$ is almost $1$ for $x$ near infinite, we only have to show that for $x$ large $\frac{x}{2}\sinh\frac{\pi}{2x}$ is decreasing
In this way we are proving that our function decrease when the initial $q$ is near 1. This was the difficult part before.
This strategy must be sufficient.