This article provides a library written in C++ for evaluating multicenter integrals of polarization operators between Cartesian Gaussians. CUDA kernels that take advantage of the acceleration provided by graphical processing units (GPU) are also included. For convenience the integral routines can be exported as a python module. This library thus complements existing open source quantum chemistry integral packages for other one-and two-electron operators, such as GEN1INT and libint. Outline: The reader is assumed to be familiar with the notation in the original articles by Schwerdtfeger and Silberbach. Primed equations refer to the first article (Ref. ). After defining the quantities of interest, some space is devoted to the subcase 2.2, which contained a mistake that was fixed only recently. Then the efficient and numerically stable calculation of certain auxiliary functions is explained in detail. In section 4.1 an upper bound on the polarization integrals based on the Cauchy-Schwarz inequality is given that allows negligible integrals to be efficiently prescreened. Finally, section 5 explains how to use the library. Expressions for nuclear gradients of the polarization integrals are provided in the appendix.