Efficient modeling of eddy current testing (ECT) signals is needed in many areas of industry. Design of probes may be improved and interpretation of experimental signals better understood by using dedicated numerical simulation tools, if they are computationally effective and accurate, yet remain simple enough to be applied at an end-user level. A boundary element method (BEM), dedicated to the numerical simulation of ECT signals due to complex narrow cracks within a planar multilayered structure (PMS) and presenting arbitrary orientations, is investigated. The theoretical formulation relies on the calculation of the dyadic Green operator, associated to the PMS, via appropriate vector wave function expansions. Then, the use of the discrete complex image method followed by the application of the generalized pencil of function method is proposed for efficient computation of this operator. Results of the validation of the complete model by comparison with the experimental data acquired in the laboratory-controlled conditions and with data computed by a finite element code are discussed.