How to fix significant inprecision in converting RDS to ETRS89 with RDNAPTRANS2018 implementation in C#?

Question

I am attempting to convert RDS coordinates to ETRS89. I will have to do this thousands of times a minute for which reason I am creating the implementation myself.

To do this I have gotten the RDNAPTRANS2018 and implemented the formula in C#. This code is at the bottom of the question. However, when executing the code I find each coordinate has about ~100-1000 meters difference from the actual RDS coordinates. An example:

//Given RDS
RD-X: 74402
RD-Y: 448670

//Calculated ETRS89
POINT (52.018499953545344 4.209861310427128)

//RDS of calculated ETRS89
RD-X: 74186
RD-Y: 448449

I have considered there may be an issue with rounding or imprecision but I was not able to find it and have confirmed all constant values used in the equation.

For your consideration the pictures below are the two parts of the formula I am attempting to addapt. The full paper can be requested at this link: (https://formulieren.kadaster.nl/aanvragen_rdnaptrans) and is free of charge.

And due to its relevance chapter 2.4.1:

Finally my current implementation:

// RD parameters
    private const double k = 0.9999079;
    private const double x0 = 155000;
    private const double y0 = 463000;

    private readonly double phi0 = 52.156161.ToRadians(); // Convert degrees to radians
    private readonly double lambda0 = 5.387639.ToRadians(); // Convert degrees to radians
    private const double R = 6369347.914861309;  // Radius of the Earth in meters
    private readonly double epsilon = 0.000000001.ToRadians();

    // Bessel ellipsoid parameters
    private const double a = 6377397.155;
    private const double f = 1 / 299.1528128;
    private readonly double e = Math.Sqrt(f * (2 - f));

    private Point ParseLocationFromMessage(int xRD, int yRD)
    {
        // Inverse oblique stereographic projection from RD to sphere
        double r = Math.Sqrt(Math.Pow(xRD - x0, 2) + Math.Pow(yRD - y0, 2));
        double sinAlpha = (xRD - x0) / r;
        double cosAlpha = (yRD - y0) / r;
        double psi = 2 * Math.Atan(r / (2 * k * R));


        //Try cos/sin alpha cos / sin
        double xBar, yBar, zBar;
        if (xRD != x0 || yRD != y0)
        {
            xBar = Math.Cos(phi0) * Math.Cos(psi) - cosAlpha * Math.Sin(phi0) * Math.Sin(psi);
            yBar = sinAlpha * Math.Sin(psi);
            zBar = cosAlpha * Math.Cos(phi0) * Math.Sin(psi) + Math.Sin(phi0) * Math.Cos(psi);
        }
        else
        {
            xBar = Math.Cos(phi0);
            yBar = 0;
            zBar = Math.Sin(phi0);
        }

        double latitude = Math.Asin(zBar);
        double longitude;
        if (xBar > 0)
        {
            longitude = lambda0 + Math.Atan(yBar / xBar);
        }
        else if (xBar < 0 && xRD >= x0)
        {
            longitude = lambda0 + Math.Atan(yBar / xBar) + 180.0.ToRadians();
        }
        else if (xBar < 0 && xRD < x0)
        {
            longitude = lambda0 + Math.Atan(yBar / xBar) - 180.0.ToRadians();
        }
        else if (xBar == 0 && xRD > x0)
        {
            longitude = lambda0 + 90.0.ToRadians();
        }
        else if (xBar == 0 && xRD < x0)
        {
            longitude = lambda0 - 90.0.ToRadians();
        }
        else
        {
            longitude = lambda0;
        }

        double n = Math.Sqrt(1 + (Math.Pow(e, 2) * Math.Pow(Math.Cos(phi0), 4) / (1 - Math.Pow(e, 2))));
        double Rn = a / Math.Sqrt(1 - Math.Pow(e, 2) * Math.Pow(Math.Sin(phi0), 2));
        double Rm = Rn * (1 - Math.Pow(e, 2)) / 1 - Math.Pow(e, 2) * Math.Pow(Math.Sin(phi0), 2);
        double largePhi0 = Math.Atan(Math.Sqrt(Rm) / Math.Sqrt(Rn) * Math.Tan(phi0));
        double w0 = Math.Log(Math.Tan((largePhi0 + 90.0.ToRadians()) / 2));
        double q0 = Math.Log(Math.Tan((phi0 + 90.0.ToRadians()) / 2)) - (e / 2 * Math.Log(1 + e * Math.Sin(phi0) / (1 - e * Math.Sin(phi0))));
        double m = w0 - n * q0;
        double w = Math.Log(Math.Tan((latitude + 90.0.ToRadians()) / 2));
        double q = (w - m) / n;

        double phi_i = 0;
        double phi_i_next = 0;
        while (Math.Abs(phi_i_next - phi_i) < epsilon)
        {
            phi_i = phi_i_next;

            if (latitude == -90.0.ToRadians() || latitude == 90.0.ToRadians())
            {
                phi_i_next = latitude;
            }
            else if (-90.0.ToRadians() < latitude && latitude < 90.0.ToRadians())
            {
                phi_i_next = 2 * Math.Atan(Math.Exp(q + e / 2 * Math.Log((1 + e * Math.Sin(phi_i)) / (1 - e * Math.Sin(phi_i))))) - 90.0.ToRadians();
            }
        }

        double lambda_n = ((longitude - lambda0) / n) + lambda0;
        double lambda_wrapped = lambda_n.ToDegrees() % 360;
        if (lambda_wrapped > 180) lambda_wrapped -= 360;
        else if (lambda_wrapped <= -180) lambda_wrapped += 360;

        double latitudeFinal = phi_i_next.ToDegrees();
        double longitudeFinal = lambda_wrapped;
        Point point = new Point(new Coordinate(latitudeFinal, longitudeFinal));
        Console.WriteLine(point.ToString());
        Console.ReadLine();
        return point;
    }

Any suggestion as to where I may be wrong or imprecise or inaccurate or all three at the same time are very welcome. Thank you!