Converting latitude, longitude (epsg:4326) into EPSG:3857?

I'm looking to study the math involved to see how a given latitude and longitude in decimal degrees can be converted to EPSG:3857. Can anybody point to a good reference or a possible Ppython/C/C++/Java library source code for checking the math?

from pyproj import Proj, transform

P3857 = Proj(init='epsg:3857')
P4326 = Proj(init='epsg:4326')

x,y = transform(P4326, P3857, lon, lat)

This code snippet is C# and using an array called vertex to hold [x, y]

double smRadius = 6378136.98;
double smRange = smRadius * Math.PI * 2.0;
double smLonToX = smRange / 360.0;
double smRadiansOverDegrees = Math.PI / 180.0;

...

// compute x-map-unit
vertex *= smLonToX;

double y = vertex;

// compute y-map-unit
if (y > 86.0)
{
vertex = smRange;
}
else if (y < -86.0)
{
vertex = -smRange;
}
else
{
y = Math.Log(Math.Tan(y) + (1.0 / Math.Cos(y)), Math.E);
vertex = y * smRadius;
}

Following up on the answer of @Russel at ICS. Today I rewrote his code to convert wgs84 (epsg:4326) coordinates to wgs84/pseudo mercator (epsg:3857) in R:

vertex = list(x_coordinate, y_coordinate)
smRange = smRadius * pi * 2.0
smLonToX = smRange / 360.0
smRadiansOverDegrees = pi / 180.0

vertex[] = vertex[] *smLonToX

y = vertex[]

if (y > 86.0){
vertex[] = smRange
} else if (y < -86.0){
vertex[] = -smRange
} else {
y = y * smRadiansOverDegrees
y = log(tan(y) + (1.0 / cos(y)), base = exp(1))
vertex[] = y * smRadius
}