Optimization models involving a polynomial objective function and multiple polynomial constraints with discrete variables are often encountered in engineering, management and systems. Treating the non-convex cross-product terms is the key. State-of- the-art methods usually convert such a problem into a 0-1 mixed integer linear programming problem, and, then adopt a branch-and- bound scheme to find an optimal solution. Much effort has been spent on reducing the required numbers of variables and linear constraints as well as on avoiding unbalanced branch-and- bound trees. This talk presents a novel idea of linearizing the discrete cross-product terms in an extremely effective manner. Theoretical analysis shows that the new method significantly reduces the required number of linear constraints from O(h^3 n^3) to O(hn) for a cubic polynomial discrete program with n variables in h possible values. Numerical experiments also confirm a two-order (10^2 times) reduction in computational time for randomly generated problems with millions of variables and constraints.