One can generate these formulas systematically using the method described in Arndt's book (as already mentioned by meij in a comment).

I have tabulated what I believe are the best possible n-term formulas for the first n primes up to n = 25, along with the analogous formulas for primitive arctangents up to n = 22: see Table 1 and Table 2 in https://arxiv.org/abs/2207.02501

These are the first 7 sets of denominators:

{3}

{7, 17}

{31, 49, 161}

{251, 449, 4801, 8749}

{351, 1079, 4801, 8749, 19601}

{1574, 4801, 8749, 13311, 21295, 246401}

{8749, 21295, 24751, 28799, 74359, 388961, 672281}

And for the first 25 primes:

{373632043520429, 386624124661501, 473599589105798, 478877529936961, 523367485875499, 543267330048757, 666173153712219, 1433006524150291, 1447605165402271, 1744315135589377, 1796745215731101, 1814660314218751, 2236100361188849, 2767427997467797, 2838712971108351, 3729784979457601, 4573663454608289, 9747977591754401, 11305332448031249, 17431549081705001, 21866103101518721, 34903240221563713, 99913980938200001, 332110803172167361, 19182937474703818751}

The algorithm here is just to start with the largest candidate denominator (the list is finite) and keep adding denominators that preserve linear independence of the exponents vectors. I believe (but have not tried to prove) that this results in the lowest possible Lehmer measure.

Note that for example the 7-term set {28799, 57121, 62425, 74359, 246401, 388961, 672281} mentioned by Douglas Zare does not work because there is a linear dependence: acoth(57121) - acoth(62425) - acoth(672281) = 0.