[Closed] Converting Coordinates from OsGrid to Lat Lon

I am trying to create a script that allows me to accurately convert Ordnance Survey Coordinates into Latitude and Longitude for use in a KML exporter. Here is the function I have created so far

function OSGridToLatLong East North =
(

  Local a = 6377563.396, b = 6356256.910              -- Airy 1830 major & minor semi-axes
  Local F0 = 0.9996012717                             -- NatGrid scale factor on central meridian
  Local lat0 = 49*PI/180, lon0 = -2*PI/180  -- NatGrid true origin
  Local N0 = -100000, E0 = 400000                     -- northing & easting of true origin, metres
  Local e2 = 1 - (b*b)/(a*a)                          -- eccentricity squared
  Local n = (a-b)/(a+b), n2 = n*n, n3 = n*n*n

  Local lat=lat0, M=0
  do (
    lat = (North-N0-M)/(a*F0) + lat

    Local Ma = (1 + n + (5/4)*n2 + (5/4)*n3) * (lat-lat0)
    Local Mb = (3*n + 3*n*n + (21/8)*n3) * sin(lat-lat0) * cos(lat+lat0)
    Local Mc = ((15/8)*n2 + (15/8)*n3) * sin(2*(lat-lat0)) * cos(2*(lat+lat0))
    Local Md = (35/24)*n3 * sin(3*(lat-lat0)) * cos(3*(lat+lat0))
    M = b * F0 * (Ma - Mb + Mc - Md)                -- meridional arc

  ) while (North-N0-M >= 0.00001)  -- ie until < 0.01mm

  Local cosLat = cos(lat), sinLat = sin(lat)
  Local nu = a*F0/sqrt(1-e2*sinLat*sinLat)              -- transverse radius of curvature
  Local rho = a*F0*(1-e2)/pow(1-e2*sinLat*sinLat) 1.5  -- meridional radius of curvature
  Local eta2 = nu/rho-1

  Local tanLat = tan(lat)
  Local tan2lat = tanLat*tanLat, tan4lat = tan2lat*tan2lat, tan6lat = tan4lat*tan2lat
  Local secLat = 1/cosLat
  Local nu3 = nu*nu*nu, nu5 = nu3*nu*nu, nu7 = nu5*nu*nu
  Local VII = tanLat/(2*rho*nu)
  Local VIII = tanLat/(24*rho*nu3)*(5+3*tan2lat+eta2-9*tan2lat*eta2)
  Local IX = tanLat/(720*rho*nu5)*(61+90*tan2lat+45*tan4lat)
  Local X = secLat/nu
  Local XI = secLat/(6*nu3)*(nu/rho+2*tan2lat)
  Local XII = secLat/(120*nu5)*(5+28*tan2lat+24*tan4lat)
  Local XIIA = secLat/(5040*nu7)*(61+662*tan2lat+1320*tan4lat+720*tan6lat)

  Local dE = (East-E0), dE2 = dE*dE, dE3 = dE2*dE, dE4 = dE2*dE2, dE5 = dE3*dE2, dE6 = dE4*dE2, dE7 = dE5*dE2
  lat = lat - VII*dE2 + VIII*dE4 - IX*dE6
  Local lon = lon0 + X*dE - XI*dE3 + XII*dE5 - XIIA*dE7


  return  "[" + (formattedPrint (radtodeg lat) format: "12.12g")+ "," + (formattedPrint (radtodeg lon) format: "12.12g") + ",0]"
)

function convert p =
(

	P1 = p
	P1.x = degtorad p.x
	P1.y = degtorad p.y
 
	Local a = 6377563.396d0, b = 6356256.910d0
 
  Local sinPhi = sin(p1.x), cosPhi = cos(p1.x)
  Local sinLambda = sin(p1.y), cosLambda = cos(p1.y)
  Local H = p1.z
 
  Local eSq = (a*a - b*b) / (a*a)
  Local nu = a / sqrt(1 - eSq*sinPhi*sinPhi)
 
  Local x1 = (nu+H) * cosPhi * cosLambda
  Local y1 = (nu+H) * cosPhi * sinLambda
  Local z1 = ((1-eSq)*nu + H) * sinPhi
 
 
   -- apply helmert transform using appropriate params
 
  Local tx = 446.448, ty = -125.157, tz = 542.060
  Local rx = 0.1502/3600 * PI/180  --normalise seconds to radians
  Local ry = 0.2470/3600 * PI/180
  Local rz = 0.8421/3600 * PI/180
  Local s1 = -20.4894/1e6 + 1              -- normalise ppm to (s+1)
 
   --apply transform
  Local x2 = tx + x1*s1 - y1*rz + z1*ry
  Local y2 = ty + x1*rz + y1*s1 - z1*rx
  Local z2 = tz - x1*ry + y1*rx + z1*s1
 
 
   -- convert cartesian to polar coordinates (using ellipse 2)
 
  a = 6378137d0
  b = 6356752.3142d0
  Local precision = 4 / a   --results accurate to around 4 metres
 
  eSq = (a*a - b*b) / (a*a)
  Local p = sqrt(x2*x2 + y2*y2)
  Local phi = atan2 z2 (p*(1-eSq))
	Local phiP = 2*PI
  
	do
	(
		nu = a / sqrt(1 - eSq*sin(phi)*sin(phi))
		phiP = phi
		phi = atan2 (z2 + eSq*nu*sin(phi)) p
	  )	while (abs(phi-phiP) > precision) 
  
  Local lambda = atan2 y2 x2
  H = p/cos(phi) - nu

 
  return #(radtodeg phi, radtodeg lambda, H)
)
convert (execute (OSGridToLatLong 530000 180000))

The maths of which were taken from http://www.movable-type.co.uk/scripts/latlong-gridref.html

However I seem to be experiencing mathematical inaccuracies from the above script? With easting:530000 and Northing:180000 I would expect to return:

easting: 51.503992°
Northing: 0.128329°

(As converted from here http://www.ordnancesurvey.co.uk/gps/transformation )

Instead I seem to return:

easting: 51.7955°
Northing: -0.882482°

Does anybody have any ideas as to what could be causing the inaccuracies? Could it be to do with single-precision floating point?

Any thoughts would be most appreciated.