Utilizing four archival X-ray datasets taken with the Hard X-ray Detector onboard Suzaku, timing studies were performed on three magnetars, 1E 1841−045 (observed in 2006), SGR 0501+4516 (2008), and 1RXS J170849.0−400910 (2009 and 2010). Their pulsations were reconfirmed, typically in an energy range of 12–50 keV. The 11.783 s pulses of 1E 1841−045 and those of SGR 0501+4516 at 5.762 s were periodically phase modulated, with a long period of ≈23.4 ks and ≈16.4 ks, respectively. The pulse-phase modulation was also observed, at ≈46.5 ks, from two datasets of 1RXS J170849.0−400910. In all these cases, the modulation amplitude was 6 to 16% of the pulse cycle. Including previously confirmed four objects, this characteristic timing behavior is now detected from seven magnetars in total, and interpreted as a result of free precession of neutron stars that are deformed to an asphericity of ∼10−4. Assuming that the deformation is due to magnetic stress, these magnetars are inferred to harbor toroidal magnetic fields of Bt ∼ 1016 G. By comparing the estimated Bt of these objects with their poloidal dipole field Bd, the Bt/Bd ratio is found to increase with their characteristic age. Therefore, the toroidal fields of magnetars are likely to last longer than their poloidal fields. This explains the presence of some classes of neutron stars that have relatively weak Bd but are suspected to hide strong Bt inside them.