OpenGL weird cube rotation

Question

I'm learning modern OpenGL on my own. And I am currently running an issue while rendering 3D cube.

First of all I create a sprite, which is my class containing position, width, height and a texture of a single rectangular sprite. Then I pass it to my renderer class to create a cube which has every wall composed of the same sprite.

        void Renderer::constructWall(maths::vec3 bottomLeft, maths::vec3 topLeft, maths::vec3 topRight, maths::vec3 bottomRight, maths::vec4 color){
        m_Buffer->position = bottomLeft;
        m_Buffer->color = maths::vec4(color.x, color.y, color.z, 1.0f);
        m_Buffer->tex = maths::vec2(0.0f, 0.0f);
        m_Buffer++;

        m_Buffer->position = topLeft;
        m_Buffer->color = maths::vec4(color.x, color.y, color.z, 1.0f);
        m_Buffer->tex = maths::vec2(0.0f, 1.0f);
        m_Buffer++;

        m_Buffer->position = topRight;
        m_Buffer->color = maths::vec4(color.x, color.y, color.z, 1.0f);
        m_Buffer->tex = maths::vec2(1.0f, 1.0f);
        m_Buffer++;

        m_Buffer->position = bottomRight;
        m_Buffer->color = maths::vec4(color.x, color.y, color.z, 1.0f);
        m_Buffer->tex = maths::vec2(1.0f, 0.0f);
        m_Buffer++;

        m_IndexCount += 6;
    }

m_Buffer is a pointer of type VertexData, which is a struct:

        struct VertexData{
        maths::vec3 position;
        maths::vec4 color;
        maths::vec2 tex;
    };

m_IndexCount is a simple counter for how many indices I need to render when I commit a flush in my Renderer class. I create a cube by modelling vertices of my cube and use the constructWall function this way:

        maths::vec3 A(vertex.x, vertex.y, -depth/2);
        maths::vec3 B(vertex.x, vertex.y + size.y, -depth / 2);
        maths::vec3 C(vertex.x + size.x, vertex.y + size.y, -depth / 2);
        maths::vec3 D(vertex.x + size.x, vertex.y, -depth / 2);
        maths::vec3 E(vertex.x, vertex.y, depth/2);
        maths::vec3 F(vertex.x, vertex.y + size.y, depth/2);
        maths::vec3 G(vertex.x + size.x, vertex.y + size.y, depth/2);
        maths::vec3 H(vertex.x + size.x, vertex.y, depth/2);
        constructWall(A, B, C, D, color1);
        constructWall(E, F, G, H, color1);
        constructWall(A, B,F,E, color1);
        constructWall(D, C, G, H, color1);
        constructWall(A, E, H, D, color1);
        constructWall(B, F, G, C, color1);

When I create a cube it works perfectly, rotating it in one axis works too, but when I try to rotate it in more than one axis, my cube sometimes acts weirdly like that:

Source: http://imgur.com/fI3s9fC

There is also my matrix code:

        mat4::mat4()
    {
        for (int i = 0; i < 4 * 4; i++)
            elements[i] = 0.0f;
    }

    mat4::mat4(float diagonal)
    {
        for (int i = 0; i < 4 * 4; i++)
            elements[i] = 0.0f;

        elements[0 + 0 * 4] = diagonal;
        elements[1 + 1 * 4] = diagonal;
        elements[2 + 2 * 4] = diagonal;
        elements[3 + 3 * 4] = diagonal;
    }

    mat4 mat4::identity()
    {
        return mat4(1.0f);
    }

    mat4& mat4::multiply(const mat4& other)
    {
        float data[16];
        for (int y = 0; y < 4; y++)
        {
            for (int x = 0; x < 4; x++)
            {
                float sum = 0.0f;
                for (int e = 0; e < 4; e++)
                {
                    sum += elements[x + e * 4] * other.elements[e + y * 4];
                }
                data[x + y * 4] = sum;
            }
        }
        memcpy(elements, data, 4 * 4 * sizeof(float));

        return *this;
    }

    vec3 mat4::multiply(const vec3& other) const
    {
        return vec3(
            columns[0].x * other.x + columns[1].x * other.y + columns[2].x * other.z + columns[3].x,
            columns[0].y * other.x + columns[1].y * other.y + columns[2].y * other.z + columns[3].y,
            columns[0].z * other.x + columns[1].z * other.y + columns[2].z * other.z + columns[3].z
            );
    }

    vec4 mat4::multiply(const vec4& other) const
    {
        return vec4(
            columns[0].x * other.x + columns[1].x * other.y + columns[2].x * other.z + columns[3].x * other.w,
            columns[0].y * other.x + columns[1].y * other.y + columns[2].y * other.z + columns[3].y * other.w,
            columns[0].z * other.x + columns[1].z * other.y + columns[2].z * other.z + columns[3].z * other.w,
            columns[0].w * other.x + columns[1].w * other.y + columns[2].w * other.z + columns[3].w * other.w
            );
    }

    mat4 operator*(mat4 left, const mat4& right)
    {
        return left.multiply(right);
    }

    mat4& mat4::operator*=(const mat4& other)
    {
        return multiply(other);
    }

    vec3 operator*(const mat4& left, const vec3& right)
    {
        return left.multiply(right);
    }

    vec4 operator*(const mat4& left, const vec4& right)
    {
        return left.multiply(right);
    }

    mat4 mat4::orthographic(float left, float right, float bottom, float top, float near, float far)
    {
        mat4 result(1.0f);

        result.elements[0 + 0 * 4] = 2.0f / (right - left);

        result.elements[1 + 1 * 4] = 2.0f / (top - bottom);

        result.elements[2 + 2 * 4] = 2.0f / (near - far);

        result.elements[0 + 3 * 4] = (left + right) / (left - right);
        result.elements[1 + 3 * 4] = (bottom + top) / (bottom - top);
        result.elements[2 + 3 * 4] = (far + near) / (far - near);

        return result;
    }

    mat4 mat4::perspective(float fov, float aspectRatio, float near, float far)
    {
        mat4 result(1.0f);

        float q = 1.0f / tan(toRadians(0.5f * fov));
        float a = q / aspectRatio;

        float b = (near + far) / (near - far);
        float c = (2.0f * near * far) / (near - far);

        result.elements[0 + 0 * 4] = a;
        result.elements[1 + 1 * 4] = q;
        result.elements[2 + 2 * 4] = b;
        result.elements[3 + 2 * 4] = -1.0f;
        result.elements[2 + 3 * 4] = c;

        return result;
    }

    mat4 mat4::translation(const vec3& translation)
    {
        mat4 result(1.0f);

        result.elements[0 + 3 * 4] = translation.x;
        result.elements[1 + 3 * 4] = translation.y;
        result.elements[2 + 3 * 4] = translation.z;

        return result;
    }

    mat4 mat4::rotation(float angle, const vec3& axis)
    {
        mat4 result(1.0f);

        float r = toRadians(angle);
        float c = cos(r);
        float s = sin(r);
        float omc = 1.0f - c;

        float x = axis.x;
        float y = axis.y;
        float z = axis.z;

        result.elements[0 + 0 * 4] = x * omc + c;
        result.elements[1 + 0 * 4] = y * x * omc + z * s;
        result.elements[2 + 0 * 4] = x * z * omc - y * s;

        result.elements[0 + 1 * 4] = x * y * omc - z * s;
        result.elements[1 + 1 * 4] = y * omc + c;
        result.elements[2 + 1 * 4] = y * z * omc + x * s;

        result.elements[0 + 2 * 4] = x * z * omc + y * s;
        result.elements[1 + 2 * 4] = y * z * omc - x * s;
        result.elements[2 + 2 * 4] = z * omc + c;

        return result;
    }

    mat4 mat4::scale(const vec3& scale)
    {
        mat4 result(1.0f);

        result.elements[0 + 0 * 4] = scale.x;
        result.elements[1 + 1 * 4] = scale.y;
        result.elements[2 + 2 * 4] = scale.z;

        return result;
    }

Answer 1

The formulas seem wrong for calculating the diagonal elements in your code for rotating around an arbitrary axis. Try changing to this:

    result.elements[0 + 0 * 4] = x * x + (1.0f - (x * x)) * c;
    result.elements[1 + 0 * 4] = y * x * omc + z * s;
    result.elements[2 + 0 * 4] = x * z * omc - y * s;

    result.elements[0 + 1 * 4] = x * y * omc - z * s;
    result.elements[1 + 1 * 4] = y * y + (1.0f - (y * y)) * c;
    result.elements[2 + 1 * 4] = y * z * omc + x * s;

    result.elements[0 + 2 * 4] = x * z * omc + y * s;
    result.elements[1 + 2 * 4] = y * z * omc - x * s;
    result.elements[2 + 2 * 4] = z * z + (1.0f - (z * z)) * c;

Refer to section 5.2 of this mathematical explanation.