How does B3D set the quaternion rotation?

public static class QuaternionExtensions
{        
    public static Vector3 GetPitchYawRoll(this Quaternion q)
    {
        return new Vector3(q.GetPitch(), q.GetYaw(), q.GetRoll());
    }

    public static float GetPitch(this Quaternion q)
    {
        return q.GetK().GetPitch();
    }

    public static float GetYaw(this Quaternion q)
    {
        return q.GetK().GetYaw();
    }

    public static float GetRoll(this Quaternion q)
    {
        // This is M12 * M22 of rotation matrix
        float xx = q.X * q.X;
        float xy = q.X * q.Y;
        float zz = q.Z * q.Z;
        float wz = q.W * q.Z;
        return (float)Math.Atan2(2f * (xy - wz), 1f - 2f * (xx + zz));
    }

    public static Vector3 GetK(this Quaternion q)
    {
        float xz = q.X * q.Z;
        float wy = q.W * q.Y;
        float yz = q.Y * q.Z;
        float wx = q.W * q.X;
        float xx = q.X * q.X;
        float yy = q.Y * q.Y;
        return new Vector3(
            2f * (xz - wy),
            2f * (yz + wx),
            1f - 2f * (xx + yy));
    }

}

public static class Vector3Extensions
{
    public static float GetPitch(this Vector3 v)
    {
        return (float)-Math.Atan2(v.Y, Math.Sqrt(v.X * v.X + v.Z * v.Z));
    }

    public static float GetYaw(this Vector3 v)
    {
        return (float)-Math.Atan2(v.X, v.Z);
    }
}

public Vector3 Rotation
{
    get
    {
        return _rotation.GetPitchYawRoll();
    }
    set
    {
        // Whoever gives me the solution here wins some PayPal
        UpdateViewMatrix();
    }
}

private void UpdateViewMatrix()
{
    Matrix rotation = Matrix.RotationQuaternion(_rotation);
    Matrix translation = Matrix.Translation(-_position.X, -_position.Y, -_position.Z);
    _view = translation * rotation;
}

public static void RotationQuaternion(ref Quaternion rotation, out Matrix result)
{
    float xx = rotation.X * rotation.X;
    float yy = rotation.Y * rotation.Y;
    float zz = rotation.Z * rotation.Z;
    float xy = rotation.X * rotation.Y;
    float zw = rotation.Z * rotation.W;
    float zx = rotation.Z * rotation.X;
    float yw = rotation.Y * rotation.W;
    float yz = rotation.Y * rotation.Z;
    float xw = rotation.X * rotation.W;
    result = Matrix.Identity;
    result.M11 = 1f - 2f * (yy + zz);
    result.M12 = 2f * (xy + zw);
    result.M13 = 2f * (zx - yw);
    result.M21 = 2f * (xy - zw);
    result.M22 = 1f - 2f * (zz + xx);
    result.M23 = 2f * (yz + xw);
    result.M31 = 2f * (zx + yw);
    result.M32 = 2f * (yz - xw);
    result.M33 = 1f - 2f * (yy + xx);
}

	Matrix( const Quat &q ){
		float xx=q.v.x*q.v.x,yy=q.v.y*q.v.y,zz=q.v.z*q.v.z;
		float xy=q.v.x*q.v.y,xz=q.v.x*q.v.z,yz=q.v.y*q.v.z;
		float wx=q.w*q.v.x,wy=q.w*q.v.y,wz=q.w*q.v.z;
		i=Vector( 1-2*(yy+zz),2*(xy-wz),2*(xz+wy) ),
		j=Vector( 2*(xy+wz),1-2*(xx+zz),2*(yz-wx) ),
		k=Vector( 2*(xz-wy),2*(yz+wx),1-2*(xx+yy) );
	}

	float Pitch() {
		return -asinf( 2*q.z*q.y+2*q.x*q.w );
	}

	float Yaw() {
		float x2=q.x*q.x, y2=q.y*q.y;
		return atan2f( 2*q.y*q.w-2*q.z*q.x,1-2*y2-2*x2 );
	}

	float Roll() {
		float x2=q.x*q.x, z2=q.z*q.z;
		return -atan2f( 2*q.z*q.w-2*q.y*q.x,1-2*z2-2*x2 );
	}

; Euler to/from Quaternion, adapted from
; www.euclideanspace.com/maths/geometry/rotations/conversions/eulerToQuaternion/index.htm
; That site uses heading to mean yaw, attitude for pitch, bank for roll
; It also interchanges the meaning of x and z and uses a right handed system.

Type Rotation
	Field pitch#, yaw#, roll#
End Type

Type Quat
	Field w#, x#, y#, z#
End Type


Function RotationToQuaternion( q.Quat, r.Rotation )

	Local c1#,c2#,c3#
	Local s1#,s2#,s3#
	Local w#,x#,y#,z#
	
	c1 = Cos(r\yaw/2)
	c2 = Cos(r\pitch/2)
	c3 = Cos(r\roll/2)
	
	s1 = Sin(r\yaw/2)
	s2 = Sin(r\pitch/2)
	s3 = Sin(r\roll/2)
	
	; x,z swapped and z made negative.
	
	q\w = c1*c2*c3 - s1*s2*s3
	q\z = s1*s2*c3 + c1*c2*s3  :  z = -z
	q\y = s1*c2*c3 + c1*s2*s3
	q\x = c1*s2*c3 - s1*c2*s3
	
End Function

Function QuaternionToRotation( r.Rotation, q.Quat )

	Local test#, sqx#,sqy#,sqz#

	test = q\z*q\y + q\x*q\w
	If test > 0.499	; singularity at north pole, 0.499 means about 86.3 degrees
		r\yaw = 2 * ATan2(q\z,q\w)
		r\pitch = 180
		r\roll = 0
		Return
	End If
	If (test < -0.499) ; singularity at south pole
		r\yaw = -2 * ATan2(q\z,q\w)
		r\pitch = -180
		r\roll = 0
		Return
	End If

	sqz = q\z*q\z
	sqy = q\y*q\y
	sqx = q\x*q\x
	r\yaw = ATan2(2*q\y*q\w-2*q\z*q\x , 1 - 2*sqy - 2*sqx)
	r\pitch = ASin(2*test)
	r\roll = ATan2(2*q\z*q\w-2*q\y*q\x , 1 - 2*sqz - 2*sqx)
	
End Function

internal static Vector3 GetPitchYawRoll(this Quaternion q)
{
    Vector3 rotation = new Vector3();

    double test = q.Z * q.Y + q.X * q.W;
    if (test > 0.499)
    {
        // Singularity at north pole, 0.499 means about 86.3 degrees
        rotation.Y = (float)(2.0 * Math.Atan2(q.Z, q.W));
        rotation.X = 180f;
    }
    else if (test < -0.499)
    {
        // Singularity at south pole
        rotation.Y = (float)(-2.0 * Math.Atan2(q.Z, q.W));
        rotation.X = -180f;
    }
    else
    {
        double zSquare = q.Z * q.Z;
        double ySquare = q.Y * q.Y;
        double xSquare = q.X * q.X;
        rotation.Y = (float)Math.Atan2(2.0 * q.Y * q.W - 2.0 * q.Z * q.X,
            1.0 - 2.0 * ySquare - 2f * xSquare);
        rotation.X = (float)Math.Asin(2.0 * test);
        rotation.Z = (float)Math.Atan2(2.0 * q.Z * q.W - 2.0 * q.Y * q.X,
            1.0 - 2.0 * zSquare - 2.0 * xSquare);
    }

    return rotation;
}

internal static Quaternion RotationQuaternion(this Vector3 v)
{
float c1 = (float)Math.Cos(v.Y / 2f);
float c2 = (float)Math.Cos(v.X / 2f);
float c3 = (float)Math.Cos(v.Z / 2f);
float s1 = (float)Math.Sin(v.Y / 2f);
float s2 = (float)Math.Sin(v.X / 2f);
float s3 = (float)Math.Sin(v.Z / 2f);

return new Quaternion(
    c1 * s2 * c3 - s1 * c2 * s3, //X
    s1 * c2 * c3 + c1 * s2 * s3, //Y
    -(s1 * s2 * c3 + c1 * c2 * s3), //Z
    c1 * c2 * c3 - s1 * s2 * s3); //W
}

	internal static Quaternion PitchYawRollQuat(this Vector3 v)
	{
		Quaternion r;

		float pitch = v.X;
		v.X = v.Y; // Swapping the pitch and yaw.
		v.Y = pitch;

		float c1=cosf(-v.Z/2),s1=sinf(-v.Z/2);
		float c2=cosf(-v.Y/2),s2=sinf(-v.Y/2);
		float c3=cosf( v.X/2),s3=sinf( v.X/2);
		float c1_c2=c1*c2,s1_s2=s1*s2;
		r.X=c1*s2*c3-s1*c2*s3;
		r.Y=c1_c2*s3+s1_s2*c3;
		r.Z.z=s1*c2*c3+c1*s2*s3;
		r.W=c1_c2*c3-s1_s2*s3;
		return r;
	}