I'm not sure what your problem is. The image I get from OSM looks like this:
And after running your commands I get this:
Which looks like how I would expect a lon/lat dataset to look, which is essentially an equirectangular projection:
The term 'projection' is often used as a synonym for the more correct term, coordinate reference system (CRS), which can include geographic and projected coordinate reference systems.
When a geographic CRS is displayed in two dimensions, its angular units are treated as if they are linear--they're just displayed. They are not displayed using a Mercator projection algorithm. The closest projection algorithm is Plate Carrée. Its forward equations are:
x (easting) = R*lambda (if the central meridian is zero)
y (northing) = R*phi
where
R = radius of the sphere (in your case the semimajor axis of WGS84 would be used)
lambda = longitude in radians
phi = latitude in radians
so you can see that true Plate Carrée coordinates are just scaled compared to displaying latitude-longitude values in degrees.
While the x equation is the same for the spherical Mercator equations, the y equation is different:
y = R ln tan (PI/4 + phi/2)
When you request data from an image-based web service like WMS, usually the layers have been cached (pre-built) in various CRS. The service then publishes which CRS are supported.
Note: Unfortunately, my company (Esri) is guilty of popularizing the usage of 'projection' instead of 'coordinate reference system'. I would just like to state that I started at Esri after that erroneous usage began.
Best Answer
Projection
EPSG:3857
uses the methodEPSG:1024
for Forward and Reverse calculations.The method is fully documented at Geomatics Guidance Note No 7, part 2
There's an mathematical example at page 40 where you can use to cross-reference your code.
17/11/2020 EDIT: The link to the above documentation from times to times, gets broken as things evolve. If that happens please feel free to update it as necessary.