This paper presents a new algorithm for volume-constrained expected compliance minimization of continuum structures with probabilistic loading directions using analytically determined exact objective and gradient functions. The algorithm is based upon the finding that for a particular set of statistical parameters the integration in the expected compliance function can be done symbolically and automatically using symbolic manipulation software. In this study, Mathematica was used to integrate and simplify the highly nonlinear expected compliance function. It will be demonstrated by examples that the result of the symbolic pre-processing step is a simple linear function defined on a particular subset of the inverse stiffness entries which is needed in the compliance computation. The coefficient set of this linear function forms the base of the exact analytical gradient computation used in the optimal solution searching optimality criteria (OC) method to define the steepest descent direction. Naturally, the applied OC method can be replaced by any other appropriate nonlinear solver. Matlab codes of the algorithm for 2D and 3D structures with exact analytical sensitivities have been developed using the topology optimization codes presented by Andreassen et al. [1] for 2D and Liu and Tovar [2] for 3D structures as starting points. Illustrative examples with Mathematica and Matlab codes are presented to demonstrate the essence and viability of the proposed approach and highlight the potential of the automatic symbolic computation in structure optimization.