<jats:p>Using GPGPU techniques and multi-precision calculation we developed the code to study QCD phase transition line in the canonical approach. The canonical approach is a powerful tool to investigate sign problem in Lattice QCD. The central part of the canonical approach is the fugacity expansion of the grand canonical partition functions. Canonical partition functions <jats:italic>Z<jats:sub>n</jats:sub></jats:italic>(<jats:italic>T</jats:italic>) are coefficients of this expansion. Using various methods we study properties of <jats:italic>Z<jats:sub>n</jats:sub></jats:italic>(<jats:italic>T</jats:italic>). At the last step we perform cubic spline for temperature dependence of <jats:italic>Z<jats:sub>n</jats:sub></jats:italic>(<jats:italic>T</jats:italic>) at fixed <jats:italic>n</jats:italic> and compute baryon number susceptibility <jats:italic>χB</jats:italic>/<jats:italic>T</jats:italic><jats:sup>2</jats:sup> as function of temperature. After that we compute numerically ∂χ/∂T and restore crossover line in QCD phase diagram. We use improved Wilson fermions and Iwasaki gauge action on the 16<jats:sup>3</jats:sup> × 4 lattice with <jats:italic>m<jats:sub>π</jats:sub></jats:italic>/<jats:italic>m<jats:sub>ρ</jats:sub></jats:italic> = 0.8 as a sandbox to check the canonical approach. In this framework we obtain coefficient in parametrization of crossover line <jats:italic>T<jats:sub>c</jats:sub></jats:italic>(<jats:italic>µ</jats:italic><jats:sup>2</jats:sup><jats:italic><jats:sub>B</jats:sub></jats:italic>) = <jats:italic>T<jats:sub>c</jats:sub></jats:italic>(<jats:italic>C</jats:italic>−<jats:italic>ĸµ</jats:italic><jats:sup>2</jats:sup><jats:italic><jats:sub>B</jats:sub></jats:italic>/<jats:italic>T</jats:italic><jats:sup>2</jats:sup><jats:italic><jats:sub>c</jats:sub></jats:italic>) with <jats:italic>ĸ</jats:italic> = −0.0453 ± 0.0099.</jats:p>