Suppose $f(x,y) = \prod_{i=1}^n (a_ix+b_iy)$
where $n$ is a constant larger than 500, and $a_i>0$, $b_i>0$ are known coefficient. There is only one global maximum.
What's the most efficient method to find positive $(x,y)$ that maximizes $f(x,y)$ under the constraint that $c_1x+c_2y=1$, where both $c_1$ and $c_2$ are larger than 0.
On
We assume a knowledge of first semester calculus, in particular that polynomial functions are differentiable, and that the maximum or minimum of such functions is attained on a closed bounded interval either at an endpoint or at a critical point. Certain claims are stated without proof, in part to avoid bogging down a discussion about how to do things with the technical analysis of why they work.
The problem can be viewed as a one-dimensional optimization of a polynomial on a closed bounded interval. Eliminate $y = (1-c_1 x)/c_2$ and restrict $0 \le x \le 1/c_1$.
Then $f(x) = \Pi_{i=1}^{n} ((a_i - \frac{b_i c_1}{c_2})x + \frac{b_i}{c_2})$ is a polynomial of degree $n$. Although all $a_i,b_i$ are assumed positive, we cannot tell from the given information which of $a_i - \frac{b_i c_1}{c_2}$ are positive (or even nonzero). Omitting as positive constant multiples any factors $\frac{b_i}{c_2}$ for which the corresponding $a_i - \frac{b_i c_1}{c_2}$ is zero, we can therefore express:
$$ f(x) = C \; \Pi_{i=1}^n (x - r_i) $$
where the $r_i = \frac{-b_i}{(a_i c_2 - b_i c_1)}$ are the roots (possibly with repetition) of $f(x)$ and $C = \Pi_{i=1}^n (a_i - \frac{b_i c_1}{c_2})$.
Depending on whether $C$ is positive or negative, maximizing $f(x)$ on $[0,1/c_1]$ is equivalent either to maximizing or minimizing the monic, fully reducible polynomial $g(x)=\frac{1}{C} f(x)$. The techniques are all but identical for these two cases, differing only in that maximizing requires attention to subintervals where $g(x) \gt 0$ and minimizing attention to subintervals where $g(x) \lt 0$.
For brevity we shall discuss only the case $C \gt 0$, so that maxima of $f(x)$ coincide with maxima of $g(x)$. However we fully discuss potential multiplicity of roots $r_i$, which we broach by sorting the roots and introducing exponents $m_i$ to represent their multiplicities. To reflect this (potential) reduction in counting distinct linear factors, we write:
$$ g(x) = \Pi_{i=1}^{N} (x - r_i)^{m_i} $$
where $n = \sum_{i=1}^{N} m_i \ge N$.
Our first task is to use these sorted roots and their multiplicities to define the subintervals on which $g(x) > 0 $ and may thus contain maxima. In general some roots $r_i$ may lie outside $[0,1/c_1]$. Those $r_i$ falling below $0$ will affect the magnitude of $g(x)$ on the interval, but not the sign of $g(x)$ (all such factors are positive there). Also the collected factors for roots $r_i$ will contribute a sign to $g(x)$ on the interval, but it is the same sign for all points in the interval.
Therefore the interval $[0,1/c_1]$ can be quickly decomposed into $k+1$ subintervals:
$$ 0 \lt r_{i_0} \lt \ldots \lt r_i \lt r_{i+1} \lt \ldots \lt r_{i_k} \lt 1/c_1 $$
in whose interiors the sign of $g(x)$ does not change. If necessary the sign of $g(x)$ in such subintervals could be identified by a single function evaluation in each interior, but after the first subinterval the sign will alternate or not at $r_i$ accordingly as exponent $m_i$ is odd or even.
For the present discussion we are interested only in the subintervals where $g(x)$ is positive (but bear in mind the opposite would be true were $C < 0$). We can say more than this, namely that $g(x)$ is unimodal on such an interval $[r_i,r_{i+1}]$, increasing from zero at $r_i$ to a positive (local) maximum, then decreasing back to zero at $r_{i+1}$.
A standard algorithm for finding the maximum of a unimodal function is golden section search, which provides a robust way to approximate the local maxima in these intervals, similar (in its guaranteed linear rate of convergence) to using bisection to find roots of a continuous function monotone on an interval.
If the simple use of golden section search on each candidate interval is not "efficient" enough (comparing their maxima to determine the global maximum), there are some further considerations to note.
First, convergence of the search can be accelerated after a few golden section iterations by switching to more sophisticated algorithms for finding a root of the derivative $g'(x)$, e.g. Brent's method (which has superlinear convergence) or Newton's method (which has quadratic convergence but thus requires evaluating the second derivative $g''(x)$). Such methods are effective when the derivative $g'(x)$ has a simple root in $(r_i,r_{i+1})$.
Second, some version of alpha-beta pruning can in principle be applied to the detection of which subinterval contains the global maximum. That is, in an ongoing search the largest maximum found so far can be compared to an upper bound computed as a by-product of searching within subsequent intervals, and when it is infeasible that the later intervals can achieve better than what previously was found, those may be discarded.
Finally, the subintervals may be searched in an order which is more likely to discover the global maximum early in the search, and lead to a more rapid disposition of the later subintervals as suggested above.
Let $$ \begin{align*} p_i &= a_i/c_1 > 0,\\ q_i &= b_i/c_2 > 0,\\ u &= c_1x,\\ v &= c_2y. \end{align*} $$ Then your maximization problem is equivalent to minimizing $$ g(u)=-\log f(x,y)=\sum_{i=1}^n -\log [(p_i-q_i)u + q_i];\quad u\in(0,1). $$ One may verify that $g''(u)=\sum_{i=1}^n \frac{(p_i-q_i)^2}{[(p_i-q_i)u + q_i]^2}$, which is positive if $p_i\neq q_i$ for some $i$. Therefore, on $(0,1)$, $g$ is either a strictly convex function when $p_i\neq q_i$ for some $i$, or a constant function otherwise.
Assume that $g$ is strictly convex. If $g'(0)=\sum_i (\frac{p_i}{q_i}-1)<0$ and $g'(1)=\sum_i (1-\frac{q_i}{p_i})>0$, then $g$ has a unique global minimum inside $(0,1)$, otherwise $g$ is strictly increasing/decreasing on $(0,1)$ and hence it has no extremum on the open interval.
When $g$ does have a unique global maximum, as it is univariate, twice differentiable and its derivatives are easy to compute, I think Newton's method is very hard to beat. However, since this is actually constrained optimization, care must be taken so that the iterates do not fall outside $[0,1]$. You may insert a bisection search step or a golden search step (see hardmath's answer) when the Newton iterate overshoots.