Regular code derived from PJ_tmer.c + pj_mlfn.c

testproj.c

// Compile with gcc -O2 testproj.c  -o testproj -lm

#include <math.h>
#include <stdio.h>

#ifndef M_PI
# define M_PI           3.14159265358979323846  /* pi */
#endif
#define HALFPI (M_PI/2)

#define EPS10   1.e-10
#define FC1 1.
#define FC2 .5
#define FC3 .16666666666666666666
#define FC4 .08333333333333333333
#define FC5 .05
#define FC6 .03333333333333333333
#define FC7 .02380952380952380952
#define FC8 .01785714285714285714

#define EN_SIZE 5


#define C00 1.
#define C02 .25
#define C04 .046875
#define C06 .01953125
#define C08 .01068115234375
#define C22 .75
#define C44 .46875
#define C46 .01302083333333333333
#define C48 .00712076822916666666
#define C66 .36458333333333333333
#define C68 .00569661458333333333
#define C88 .3076171875
#define EPS 1e-11
#define MAX_ITER 10


double
pj_mlfn(double phi, double sphi, double cphi, const double *en) {
        cphi *= sphi;
        sphi *= sphi;
        return(en[0] * phi - cphi * (en[1] + sphi*(en[2]
                + sphi*(en[3] + sphi*en[4]))));
}

double
pj_inv_mlfn(double arg, double es, double _1_div_1_m_es, const double *en/*,
            double* sinphinew, double* cosphinew*/) {
        double s, t, phi, k = _1_div_1_m_es;
        int i;

        phi = arg;
        for (i = MAX_ITER; i ; --i) { /* rarely goes over 2 iterations */
                double cosphi = cos(phi);
                s = sin(phi);
                t = 1. - es * s * s;
                t = (pj_mlfn(phi, s, cosphi, en) - arg) * (t * sqrt(t)) * k;
                if (fabs(t) < EPS)
                {
                    //sin(A+B)=sin A cos B + cos A sin B
                    //sin(A-B)=sin A cos B - cos A sin B
                    //cos(A+B)=cos A cos B - sin A sin B
                    //cos(A-B)=cos A cos B + sin A sin B
                    //*sinphinew = s - cosphi * t;
                    //*cosphinew = cosphi + s * t;
                    phi -= t;
                    return phi;
                }
                phi -= t;
        }
        //pj_ctx_set_errno( ctx, -17 );
        return phi;
}

void compute_mlfn_coeffs(double* en, double es)
{
    double t;
    en[0] = C00 - es * (C02 + es * (C04 + es * (C06 + es * C08)));
    en[1] = es * (C22 - es * (C04 + es * (C06 + es * C08)));
    en[2] = (t = es * es) * (C44 - es * (C46 + es * C48));
    en[3] = (t *= es) * (C66 - es * C68);
    en[4] = t * es * C88;
}

typedef struct { double x, y; }     XY;
typedef struct { double lam, phi; } LP;

int main(int argc, char* argv[])
{
    double y0 = 0;
    double x0 = 500000.;
    double lam0 = -3;
    double k0 = 0.9996;
    double phi0 = 0.;
    double a = 6378137;
    double ra = 1./a;
    double invf = 298.257223563;
    double f = 1 / invf;
    double es  = 2 *f - f*f;
    double en[EN_SIZE];
    double ml0, esp;
    int i;
    double acc = 0;
    double inv_k0 = 1. / k0;
    double _1_div_1_m_es = 1. / (1.-es);

    compute_mlfn_coeffs(en, es);
    
    ml0 = pj_mlfn(phi0, sin(phi0), cos(phi0), en);
    esp = es / (1. - es);

#define NITER (100*1000*1000)
    for(i=0;i<NITER;i++)
    {
        LP lp;
        XY xy;

        double n, con, cosphi, d, ds, sinphi, t;
        xy.x = (1000000.*i) / NITER;
        xy.y = (9000000.*i) / NITER;

        xy.x *= ra;
        xy.y *= ra;

        lp.phi = pj_inv_mlfn(ml0 + xy.y * inv_k0, es, _1_div_1_m_es, en/*,
                                &sinphi, &cosphi*/);
        if (fabs(lp.phi) >= HALFPI) {
                lp.phi = xy.y < 0. ? -HALFPI : HALFPI;
                lp.lam = 0.;
        } else {
                double invcosphi;
                sinphi = sin(lp.phi);
                cosphi = cos(lp.phi);
                invcosphi = 1.0 / cosphi;
                t = fabs(cosphi) > 1e-10 ? sinphi * invcosphi : 0.;
                n = esp * cosphi * cosphi;
                d = xy.x * sqrt(con = 1. - es * sinphi * sinphi) * inv_k0;
                con *= t;
                t *= t;
                ds = d * d;
                lp.phi -= (con * ds * _1_div_1_m_es) * FC2 * (1. -
                        ds * FC4 * (5. + t * (3. - 9. *  n) + n * (1. - 4 * n) -
                        ds * FC6 * (61. + t * (90. - 252. * n +
                                45. * t) + 46. * n
                - ds * FC8 * (1385. + t * (3633. + t * (4095. + 1574. * t)) )
                        )));
                lp.lam = d*(FC1 -
                        ds*FC3*( 1. + 2.*t + n -
                        ds*FC5*(5. + t*(28. + 24.*t + 8.*n) + 6.*n
                - ds * FC7 * (61. + t * (662. + t * (1320. + 720. * t)) )
                ))) * invcosphi;
        }
        //if( (i % (NITER/100)) == 0 )
        //    printf("%.16g %.16g --> %.16g %.16g\n", xy.x * a, xy.y * a,
        //           lp.lam / M_PI * 180, lp.phi / M_PI * 180);
        acc += lp.phi + lp.lam;
    }
    
    printf("Average of all longitude and latitudes: %.16g\n", acc / NITER);

    return 0;
}