How to minimize floating point inaccuracy in rotat

I have a program which visualizes translation and rotation on a 3D space. I have no problem with the main code but after hours of testing i found out that the function which causes inaccuracies is the rotation on arbitrary axis.

My object contains array of points, and three vectors which is its local axes. I use rotation function with these axes. The function works well when my position is at the origin, but once the point is translated to other than the origin it shows slightly wrong rotation until the point is out of place.

I have used this website as my guide, please jump to the bottom of the page at 6.2 and see the function with 10 parameters. On the same page, there's a link to a sample usage of the process. I have confirmed that the problem lies with floating-point inaccuracy with this application.

The sample point I wish to rotate: ( 0, -9, 0) The point that intersects the axis of rotation: ( 0, -10, 0) The direction vector of the axis of rotation: < 1, 0, 0> The angle of rotation: 2 degrees

For each loop, the current point is rotated by 2 degrees, and is repeated 20 times. The correct Y value for the first iteration is -9.0006, while in my case it gives -8.9945 and continuously fluctuates until the 5th or 6th iteration, then the correct value shows after those iterations.

Here's the rotation function, take note that Point is the point you wish to rotate, AxisPoint is the point that intersects the axis of rotation, AxisDirection is the direction vector parallel to the axis of rotation, and Degrees is the amount of angle to rotate:

private static Vector3d RotateArbitrary(Vector3d Point, Vector3d AxisPoint, Vector3d AxisDirection, double Degrees)
    {
        return new Vector3d(
            ((AxisPoint.X * (Math.Pow(AxisDirection.Y, 2) + Math.Pow(AxisDirection.Z, 2)) - AxisDirection.X * (AxisPoint.Y * AxisDirection.Y + AxisPoint.Z * AxisDirection.Z - AxisDirection.X * Point.X - AxisDirection.Y * Point.Y - AxisDirection.Z * Point.Z)) * (1 - Math.Cos(MathHelper.DegreesToRadians(Degrees))) + Point.X * Math.Cos(MathHelper.DegreesToRadians(Degrees)) + (-AxisPoint.Z * AxisDirection.Y + AxisPoint.Y * AxisDirection.Z - AxisDirection.Z * Point.Y + AxisDirection.Y * Point.Z) * Math.Sin(MathHelper.DegreesToRadians(Degrees))),
            ((AxisPoint.X * (Math.Pow(AxisDirection.X, 2) + Math.Pow(AxisDirection.Z, 2)) - AxisDirection.Y * (AxisPoint.X * AxisDirection.X + AxisPoint.Z * AxisDirection.Z - AxisDirection.X * Point.X - AxisDirection.Y * Point.Y - AxisDirection.Z * Point.Z)) * (1 - Math.Cos(MathHelper.DegreesToRadians(Degrees))) + Point.Y * Math.Cos(MathHelper.DegreesToRadians(Degrees)) + (AxisPoint.Z * AxisDirection.X - AxisPoint.X * AxisDirection.Z + AxisDirection.Z * Point.X - AxisDirection.X * Point.Z) * Math.Sin(MathHelper.DegreesToRadians(Degrees))),
            ((AxisPoint.X * (Math.Pow(AxisDirection.X, 2) + Math.Pow(AxisDirection.Y, 2)) - AxisDirection.Z * (AxisPoint.X * AxisDirection.X + AxisPoint.Y * AxisDirection.Y - AxisDirection.X * Point.X - AxisDirection.Y * Point.Y - AxisDirection.Z * Point.Z)) * (1 - Math.Cos(MathHelper.DegreesToRadians(Degrees))) + Point.Z * Math.Cos(MathHelper.DegreesToRadians(Degrees)) + (-AxisPoint.Y * AxisDirection.X + AxisPoint.X * AxisDirection.Y - AxisDirection.Y * Point.X + AxisDirection.X * Point.Y) * Math.Sin(MathHelper.DegreesToRadians(Degrees)))
            );
    }

I have already tried both float and double data type, the result is still the same, any solutions you can offer, even a change in method of rotation, I'm all open. Please help.

Due to my lengthy web searches, I finally arrived to the conclusion that I should use Quaternions. Using the current method, I found out that excessive operations on floating point variables will increase round errors. Using Quaternions was simpler and cleaner.

Here's the code if anyone is interested:

private static Vector3 RotateArbitrary(Vector3 Point, Vector3 AxisPoint, Vector3 AxisDirection, float Radians)
        {
            return Vector3.Add(Vector3.Transform(Vector3.Subtract(Point, AxisPoint), Quaternion.FromAxisAngle(AxisDirection, Radians)), AxisPoint);
        }

Take note that the the Point was translated first such that the AxisPoint is at the origin, with that, rotations can be done. The result is then translated to its original position.