var q1 = new Quaternion(2, 3, 4, 5);
var q2 = new Quaternion(5, 4, 3, 2);
var q3 = Quaternion.Multiply(q1, q2);

public void main()
{
    Vector3 result;
    Quaternion left, right, multiplied;

    left = Quaternion.RotationAxis(Vector3.UnitY, (float)(Math.PI / 2.0));
    right = Quaternion.RotationAxis(Vector3.UnitZ, (float)(Math.PI / 2.0));

    // Case1 : Matrix(Quaternion * Quaternion)
    multiplied = left * right;

    result = Vector3.TransformCoordinate(Vector3.UnitZ, Matrix.RotationQuaternion(multiplied));
    System.Diagnostics.Debug.WriteLine(result.ToString());
    //result X:0.9999999 Y:0 Z:1.192093E-07

    // Case2 : Matrix(Quaternion) * Matrix(Quaternion)
    result = Vector3.TransformCoordinate(Vector3.UnitZ, Matrix.RotationQuaternion(left) * Matrix.RotationQuaternion(right));
    System.Diagnostics.Debug.WriteLine(result.ToString());
    //result X:5.960464E-08 Y:0.9999999 Z:5.960464E-08

    // Case3 : Matrix(Multiply(Quaternion, Quaternion)) by using revised method
    multiplied = Multiply(left, right);

    result = Vector3.TransformCoordinate(Vector3.UnitZ, Matrix.RotationQuaternion(multiplied));
    System.Diagnostics.Debug.WriteLine(result.ToString());
    //result X:0 Y:0.9999999 Z:1.192093E-07
}

public Quaternion Multiply(Quaternion left, Quaternion right)
{
    Quaternion result;

    float lx = left.X;
    float ly = left.Y;
    float lz = left.Z;
    float lw = left.W;
    float rx = right.X;
    float ry = right.Y;
    float rz = right.Z;
    float rw = right.W;
    float a = (lz * ry - ly * rz);
    float b = (lx * rz - lz * rx);
    float c = (ly * rx - lx * ry);
    float d = (lx * rx + ly * ry + lz * rz);
    result.X = (lx * rw + rx * lw) + a;
    result.Y = (ly * rw + ry * lw) + b;
    result.Z = (lz * rw + rz * lw) + c;
    result.W = lw * rw - d;

    return result;
}