Background:The five-dimensional quadrupole collective model based on energy density functionals (EDFs) has often been employed to treat long-range correlations associated with shape fluctuations in nuclei. Our goal is to derive the collective inertial functions in the collective Hamiltonian by the local quasiparticle random-phase approximation (QRPA) that correctly takes into account time-odd mean-field effects. Currently, a practical framework to perform the QRPA calculation with the modern EDFs on the (β,γ ) deformation space is not available. Purpose: Toward this goal, we develop an efficient numerical method to perform the QRPA calculation on the (β,γ ) deformation space based on the Skyrme EDF. Methods: We use the finite amplitude method (FAM) for the efficient calculation of QRPA strength functions for multipole external fields. We construct a computational code of FAM-QRPA in the three-dimensional Cartesian coordinate space to handle triaxially deformed superfluid nuclei. Introduction. The shape of atomic nuclei is influenced strongly by the quantum nature of nuclear systems. Excitation spectra and their transition probabilities clearly indicate the existence of shape fluctuations and shape coexistence [1], particularly in transitional regions from spherical to deformed shapes in the ground state. The long-lived fission products (LLFP) from uranium fueled reactors, such as Pd and Zr isotopes, are located in transitional regions on the nuclear chart and demonstrate the shape mixing and coexistence. It is important to understand the basic properties of the LLFPs to develop a possible nuclear transmutation method, which is the main target of an ImPACT program "Reduction and Resource Recycling of High-level Radioactive Wastes through Nuclear Transmutation" [2].One of the standard methods of investigating nuclear manybody problems is the nuclear energy density functional (EDF) theory [3]. The nuclear EDF well describes the ground-state properties of atomic nuclei. However, on the mean-field level, it cannot describe shape fluctuations and shape coexistence. We need to go beyond the mean field for a description of such phenomena, including quantum fluctuations associated with the large-amplitude collective motion. If the EDF were constructed as an expectation value of a well-defined Hamiltonian, a possible extension would be the generator coordinate method (GCM) [4][5][6]. However, most of the EDFs are known to have a singular behavior [7,8], which prevents us from the straightforward application of the GCM.