Is it possible to ask Mathematica to give the maximum value of a function over a specific domain?

NMaximize[{f[n], 1 <= n <= 10^6, n ∈ Integers}, n]

(* Out: {63.041, {n -> 2}} *)

You can then get the corresponding symbolic function value using your definition of f[n]:

f[2] // FullSimplify

Exhaustive search is a reasonable method on a small set:

vals = f[N[Range[10^6], 16]]; // AbsoluteTiming
Position[vals, #] -> # &@Max[vals]
(*
  {60.0485, Null}
  {{2}} -> 63.04102181896
*)

Confirm f[2] is the max (with arbitrary-precision, we get evidence that round-off error has not messed up the search):

Nearest[vals, Max[vals], 2]
(*  {63.04102181896, 56.03221153}  *)

MinMax[Precision /@ vals]
(* {3.15971, 14.2151}  <-- results have > 3 digits precision *)

If you're sure round-off is no problem, machine-precision is faster:

vals = f[N[Range[10^6]]]; // AbsoluteTiming
Position[vals, #] -> # &@Max[vals]
(*
  {1.60528, Null}
  {{2}} -> 63.041
*)

Also, using Simplify[f[n], n \[Element] Integers] is 30-35% faster in this case, but that's because f[n] happens to contain terms like Sin[2 n \[Pi]].