The complementary Monte Carlo and series expansions methods of computer simulations have been described to investigate the critical behaviour of the Ashkin-Teller model in three dimensions. In the first method the invariance of the ratio of the square of the second moment of the order parameter to its fourth moment in the critical region has been exploited and some critical points on the phase boundaries have been calculated in the regions where the continuous transitions are expected. The continuity of the order parameter on the critical lines is verified by a finite size scaling analysis. Large-scale simulations have been performed on SGI Power Challenge XL and L supercomputers using the 64-bit random number generator. The numerically generated series expansions method is described for which the effective algorithm for generation of graphs starting from polygons and based on collapsing the unlinked vertices, is introduced. The new feature of our algorithm is that for each graph we introduce new links between unlinked vertices and we decorate bonds with new vertices, so that more complex graphs in an early stage of the graph generation procedure are obtained. The resulting set of graphs enables the application of the series expansions method and achievement of the precision of allocation of points on the phase diagram comparable to the precision of the Monte Carlo method, i.e. at least 3 decimal digits.