I found many source for rotate x-z plane using quaternion, but if I want to write a space shooting game, the problem is how can I also rotate Z-axis. Anyone can help? thanks!

Use a 4x4 matrix to handle you position and rotation information. You will then be free to rotate about whatever axis you want as well as translate in the same unit. Otherwise you have to cascade your transform matricies. AFAIK.

As an addition to my above post, I suggest you check out this wonderful resource provided by this site.

http://www.gamedev.net/reference/articles/article1095.asp

The use of a 4x4 or quaternion and post multiplication also allows you to avoid gimbal lock.

You will also probably want to implement some sort of SLERP algorithm to smooth things out a bit.

SLERP general definition: http://en.wikipedia.org/wiki/Slerp

A good resource with some example code:
http://number-none.com/product/Understanding%20Slerp,%20Then%20Not%20Using%20It/

The code from the above:

 Quaternion slerp(Quaternion const &v0, Quaternion const &v1, double t) {    // v0 and v1 should be unit length or else    // something broken will happen.    // Compute the cosine of the angle between the two vectors.    double dot = dot_product(v0, v1);    const double DOT_THRESHOLD = 0.9995;    if (dot > DOT_THRESHOLD) {        // If the inputs are too close for comfort, linearly interpolate        // and normalize the result.        Quaternion result = v0 + t*(v1 – v0);        result.normalize();        return result;    }    Clamp(dot, -1, 1);           // Robustness: Stay within domain of acos()    double theta_0 = acos(dot);  // theta_0 = angle between input vectors    double theta = theta_0*t;    // theta = angle between v0 and result     Quaternion v2 = v1 – v0*dot;    v2.normalize();              // { v0, v2 } is now an orthonormal basis    return v0*cos(theta) + v2*sin(theta);}