Projection of a point onto a simple convex polytope

You may want to check out Dykstra's projection algorithm (not to be confused with Dijkstra's algorithm). It precisely does what you want: It is an iterative method to compute the projection onto an intersection of convex sets using the projections onto these sets.

Note, that simply projecting onto the sets alone does not work. Consider the hypercube $[-1,1]^2$ and the hyperplane $\{x: x_1+x_2=1\}$. The projection of $(2,1)$ onto the intersection is $(1,0)$. If one first projects onto the cube, then onto the plane yields $(1/2,1/2)$, which is not the wanted projection.

From your description, I believe the optimization problem describing your projection is:

$$ \begin{aligned} \min & \quad \frac{1}{2} ||x - y||_2^2 \\ \text{s.t.} & \quad \ell_i \leq x_i \leq u_i &\quad i = 1, \dots, n \\ & \quad \alpha \leq c^T x \leq \beta \end{aligned} $$

I believe that the best way would be a simple custom interior point method. The objective augmented with log-barriers for the constraints is:

$$ A(x; \mu) = \frac{1}{2} ||x - y||_2^2 - \mu \left[ \sum_{i=1}^n \ln(u_i - x_i) + \sum_{i=1}^n \ln(x_i - \ell_i) + \ln(c^T x - \alpha) + \ln(\beta - c^T x) \right] $$

Fortunately, the Hessian has a nice structure of a diagonal matrix + two rank-one matrices, and can be easily inverted in linear time by applying the well-known Sherman–Morrison formula twice. So you get very quick convergence, and each iteration takes only linear time.