(There is a k now in the product index.)
Here is a possible solution.
We are using
https://en.wikipedia.org/wiki/Chebyshev_polynomials .
Here is also some code computing the product numerically first, then in an exact manner (by a computation in the cyclotomic field).
sage: prod( [ cos(2*k*pi/2019).n(digits=50) for k in [1..1009] ] ) * 2^1009
-0.99999999999999999999999999999999999999999999999741
sage: K.<zeta> = CyclotomicField(2019)
sage: QQ( prod( [ (zeta^k+zeta^-k)/2 for k in [1..1009] ] ) ).factor()
-1 * 2^-1009