I am trying to convert coordinates from EPSG:3857 Spherical Mercator to UTM local projection. The problem is that for my use case I cannot use any web interface or library, I should implement the formulas directly in the project.
This is highly nontrivial due to ellipsoidal vs spherical approximations of the earth's surface. If you need this with any degree of accuracy, I would strongly suggest either calling a library such as
On a high level, you can convert from EPSG:3857 back to (lat,lon), more specifically EPSG:4326, pretty easily. This is because both use a spherical approximation for the earth's surface. For standalone mathematical formulae, see for instance Converting latitude, longitude (EPSG:4326) into EPSG:3857? (this goes the other way) or E-telier's answer at https://stackoverflow.com/questions/37523872/converting-coordinates-from-epsg-3857-to-4326
The challenge is then converting (lat,lon) (4326) to UTM. That actually depends on the datum used for UTM in your local area as well as the UTM zone. The UTM zone is determine by longitude and latitude bands (see https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#Latitude_bands for the broad rules as well as for exceptions; you may be able to simplify this if you only work in one specific area that is in a fixed UTM zone).
As to datum, some so-called simplified (but actually quite complicated) formulae if you use the WGS 84 oblate spheroid are listed at https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#Simplified_formulae
However, most implementations of UTM coordinates (on topo maps etc) use a more suitable datum, for instance NAD27 or NAD83, possibly with local or temporal modifications (like NAD83 CSRS in my home of Canada). This is where complexity and the need to use an external library will be most acute.
Editing to add: There is an alternate method of attack if your ROI is very small. You could use ArcGIS or QGIS/GDAL/proj to precisely determine the correct mapping 3857 to (appropriate local) UTM coordinates for the corners of your ROI. Then you could just interpolate for coordinates inside the ROI. The scale factor of EPSG:3857 scales with cos(latitude). If your ROI is small enough that latitude and cos(latitude) is essentially constant, this will be a simple linear transformation. This approach would essentially hard-code the results of a one-time proj- (or similar-) assisted invocation to take care of the difficult spheroid-approximation math. It is, however, of course not easily portable to other ROIs.