Approximating an integral representation of the Number Partition Problem

One can write out an integral whose solution gives the number of solutions to the NP-Complete Number Partition Problem and I'm wondering if anyone has an suggestions or ideas on who to solve or approximate this integral analytically or numerically.

Given $a_0, a_1, \dots, a_{n-1} \in \mathbb{Z}$, the Number Partition Problem asks for a partition of the $a_k$'s such that $\sum_{k=0}^{n-1} \sigma_k a_k = 0$, where $\sigma_k \in ( -1,1 )$.

Consider the following Integral:

$$I(a_0, a_1, \dots, a_{n-1}) = \frac{1}{2 \pi} \int_{-\pi}^{\pi} \prod_{k=0}^{n-1} (e^{ 2 \pi i a_k \theta } + e^{ -2 \pi i a_k \theta } ) d\theta$$

$I(\dots)$ will count the number of solutions for a given instance of $(a_0, a_1, \dots, a_{n-1})$. Solving this integral is worse than NP-Hard (it's #P) so asking for a general solution is out of the question. But can one do any sort of approximation, either analytically or numerically? If you choose the $a_k$'s with some distribution, say uniform on some interval, can you exploit that randomness to help you approximate this integral?

Any ideas would be appreciated.

note: This has been studied by Borgs et all for the NPP Phase Transition and that's where I first saw this integral representation of the Number Partition Problem, but their analysis relies on approximating the family of instances given a uniform distribution on the $a_k$'s rather than trying to solve a particular instance, as I'm trying to do above.

-
Whenever I see a product of cosh's I think of the transfer matrix approach to the 1D Ising model. I wouldn't be surprised if you could map this integral onto something spin-glassy. – Steve Huntsman Sep 10 '10 at 1:41
As you mentioned, the $a_k$ are integers; what you then have is a product of cosines, and if the $a_k$ are large enough, the integrand is awfully oscillatory on the real line. Maybe you can either try finding a contour that bypasses these oscillations, or make use of specialized algorithms (e.g. those by Arieh Iserles) for handling finite oscillatory integrals. – J. M. Sep 10 '10 at 3:34

For any given $\theta$, Koksma's inequality should give you a relation between $(1/n)\sum_{k=0}^{n-1}\log(e^{2\pi ia_k\theta}+e^{-2\pi ia_k\theta})$ and $\int_0^1\log(e^{2\pi ix}+e^{-2\pi ix})dx$. The difference between them should be bounded by something depending on the discrepancy of the set of $a_k\theta$. Another relation is given by the Erdos-Turan inequality. Maybe one of these would be helpful.