高精度WGS84与GCJ-02坐标转换 C++版、Java版和Python版

WGS84转GCJ-02坐标系容易y=x+x²+x³+x⁴+...，但是反过来就难了。。原因在于高次代数方程求根。(https://www.bilibili.com/read/cv2121106)(一元n次方程的解法)  牛顿法和一元n次方程

要实现高精度转换，需要进行迭代，而这比较耗时。如果只取前4位，多项式则容易求解，但是精度不够。

有没有较快的方法实现高精度逆变换。

C++：https://www.cnblogs.com/2008nmj/p/14198389.html

以lonGCJ-02,latGCJ-02为观测值，而lonWGS84,latWGS84为未知数，列误差方程并线性化，则线性化后的误差方程为：

v=Ax-l

其中，，，

而，系数矩阵A展开为：

要使得误差v=0，那么需要什么条件：

/**
 *  \brief Covert geodetic coordinate in GCJ-02 coordinate system to geodetic coordinate
 *         in WGS84 coordinate system
 *
 *  \param [in] gcj02lon: longitude in GCJ-02 coordinate system [unit:degree]
 *  \param [in] gcj02lat: latitude in GCJ-02 coordinate system [unit:degree]
 *  \return Returns geodetic coordinate in WGS84 coordinate system
 *  \remark The encryption formula is known,and use the analytical derivation method to 
 *            solve the problem with high precision.
 *    \detail Assuming the encryption formula is 
 *
 *            gcj02lon = Wgs2Gcj(wgs84lon, wgs84lat)
 *            gcj02lat = Wgs2Gcj(wgs84lon, wgs84lat)
 *
 *     In the rectification process, (wgs84lon, wgs84lat) are unknown items. Obviously,
 *   this is a system of nonlinear equations.
 *
 *   The linear formed error functions of forward intersection come from
 *   consideration of a Taylor series expansion.
 *           V = AX - b
 *    here:
 *    V: The residuals of the observed values
 *    A: The jacobian matrix
 *    X: The modification of the unknown items
 *    b: The constant terms of the error functions
 *
 *    Then the error functions written in vector form are:
 *    | V_lon | = | dlongcj_dlonwgs  dlongcj_dlatwgs |  | d_lonwgs | - | l_lon |
 *    | V_lat | = | dlatgcj_dlonwgs  dlatgcj_dlatwgs |  | d_latwgs | - | l_lat |
 *    here:
 *    l_lon = longcj - longcj'                 // the modification of longcj
 *    l_lat = latgcj - latgcj'                 // the modification of latgcj
 *    longcj : the observed longitude in GCJ-02
 *    latgcj : the observed latitude in GCJ-02
 *    longcj' = Wgs2Gcj(wgs84lon',wgs84lat')    // estimated longitude in GCJ-02
 *    latgcj' = Wgs2Gcj(wgs84lon',wgs84lat')    // estimated latitude in GCJ-02
 *    wgs84lon' : estimated longitude in WGS84
 *    wgs84lat' : estimated latitude in WGS84
 *    d_lonwgs : unknown items
 *    d_latwgs : unknown items
 *    wgs84lon = wgs84lon' + d_lonwgs                           // update
 *    wgs84lat = wgs84lat' + d_latwgs
 *
 *      let V = [V_lon V_lat]T = 0, then
 *      d_latwgs = (l_lon * dlatgcj_dlonwgs - l_lat * dlongcj_dlonwgs) /
 *            (dlongcj_dlatwgs * dlatgcj_dlonwgs - dlatgcj_dlatwgs * dlongcj_dlonwgs)
 *      d_lonwgs = (l_lon - dlongcj_dlatwgs * d_latwgs) / dlongcj_dlonwgs
 *
 *    This iterative procedure is repeated until X= [d_lonwgs d_latwgs]T are
 *    sufficiently small.
 *  \time 01:54:46 2020/06/13
 */