Recently, QED bremsstrahlung in B meson decays into pair of scalars (πs and/or Ks) has become of interest. If experimental acceptance must be taken into account, the PHOTOS Monte Carlo technique is often used in experimental simulations. We will use scalar QED to benchmark PHOTOS, even though this theory is of limited use for complex objects. We present the analytical form of the kernel used in the older versions of PHOTOS, and the new, exact (scalar QED) one. Matrix element and phase-space Jacobians are separated in the final weight, and future extensions based on measurable electromagnetic form-factors are thus possible. The massive phase-space is controlled in the program with no approximations. Thanks to the iterative solution, all leading and next to leading logarithmic terms are properly reproduced by the Monte Carlo simulation. Simultaneously, full differential distributions over the complete multiple-body phase-space are provided. An agreement of better than 0.01% with independent calculations of scalar QED is demonstrated.