My primary task to rotate around arbitrary axis, but I can’t use glRotatef because some other calculations after transformations to be done. So I need manually transform each vector by rotation matrix. The problem is that my manual matrix taken from here and here is not equal to matrix that glRotatef makes.

You can never expect to get binary-identical floating-point values like that. Whatever the driver is doing to compute glRotatef, it’s almost certainly not the exact same floating-point code you are using. As such, it will likely not provide the exact same, binary-identical floating-point results.

Just look at the two matrices, in text form. If they look reasonably close enough (within ~4-5 significant digits), then your rotate function is probably working.

glRotate can go around an arbitrary axis, and you can retrieve the resulting matrix (for further calculation) with glGetFloatv. This won’t stall the pipeline as none of this is performed in hardware - only the final multiplication of position by MVP is.

Are you normalizing your axis vector before doing your own rotation? The spec for glRotate requires this, so if you want your code to match as closely as possible you should also be doing it.