Getting smooth, big points in OpenGL

Question

Unlike what was said previously, this is possible with the fixed-function pipeline, even with the GL_POINTS primitive type, as long as you have support for OpenGL 1.4 or the GL_ARB_point_sprite extension. Consult this document, or the OpenGL core specification of your choice : http://www.opengl.org/registry/specs/ARB/point_sprite.txt

GL_ARB_point_sprite converts points into "quads", i.e a polygon with the form of a plane. The exact primitive type it gets converted to is not defined by the specification, though it is not important. What is important is that GL_COORD_REPLACE auto-generates texture coordinates for the surface when enabled, so you can texture-map them with a sphere-shaped RGBA-texture.

EDIT: It seems like you (the poster) is right. Anti-aliased points get rounded with respect to their radius. (I've used OpenGL since 2003, and I didn't know this. [/shame]) So enabling GL_POINT_SMOOTH while you have a multisample-able visual/pixelformat, you get rounded points. Still, multisampling can be slow, so I'd implement both. Textured quads are cheap.

To request a visual with multisampling with XLib, use these two attributes in the list to glXChooseFBConfig():

GLX_SAMPLE_BUFFERS - its value should be True. This is an on/off toggle.
GLX_SAMPLES - the number of samples.

To request a pixelformat with Win32, use these two attributes in the list to ChoosePixelFormat() or wglChoosePixelFormatARB():

WGL_SAMPLE_BUFFERS_ARB Same as above, a toggle.
WGL_SAMPLES_ARB Same as above, the number of samples.

It seem that you can OR in the flag GLUT_MULTISAMPLE to glutInitDisplayMode to get multisampling in GLUT, but you can't request the number of sample buffers.

Here is how alpha-blended quads could be implemented using your test case.

void onInitialization( ) 
{
    glEnable( GL_POINT_SPRITE ); // GL_POINT_SPRITE_ARB if you're
                                 // using the functionality as an extension.

    glEnable( GL_POINT_SMOOTH );
    glEnable( GL_BLEND );
    glBlendFunc( GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA );
    glPointSize( 6.0 );

    /* assuming you have setup a 32-bit RGBA texture with a legal name */
    glActiveTexture(GL_TEXTURE0);
    glEnable( GL_TEXTURE_2D );
    glTexEnv(GL_POINT_SPRITE, GL_COORD_REPLACE, GL_TRUE);
    glTexEnv(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
    glBindTexture(GL_TEXTURE_2D, texture_name);
}    

void onDisplay()
{
    glClearColor( 1.0f, 1.0f, 1.0f, 1.0f );
    glClear( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT );

    glBegin( GL_POINTS );
        glColor4f( 0.95f, 0.207, 0.031f, 1.0f );
    for ( int i = 0; i < g_numPoints; ++i )
    {
        glVertex2f( g_points[i].X, g_points[i].Y );
    }
    glEnd();
    glFinish();
    glutSwapBuffers();
}

Image of rounded points using per-fragment alpha blending + textures:
_{(source: mechcore.net)}
Image of rounded points by using GL_POINT_SMOOTH and multisampling:
_{(source: mechcore.net)}
A little sample I made which shows both techniques. Requires libSDL and libGLEW to compile:

#include <iostream>
#include <exception>
#include <memory>
#include <SDL/SDL.h> 
#include <cmath>
#include <GL/glew.h>
#include <GL/glu.h>

#define ENABLE_TEXTURE
#define ENABLE_MULTISAMPLE

int Width = 800;
int Height = 600;

void Draw(void);
void Init(void);

inline float maxf(float a, float b)
{
    if(a < b)
        return b;
    return a;
}

inline float minf(float a, float b)
{
    if(a > b)
        return b;
    return a;
}

GLuint texture_name;

int main(void)
{
    try {
        SDL_Init(SDL_INIT_VIDEO);
        SDL_GL_SetAttribute(SDL_GL_RED_SIZE, 8);
        SDL_GL_SetAttribute(SDL_GL_GREEN_SIZE, 8);
        SDL_GL_SetAttribute(SDL_GL_BLUE_SIZE, 8);
        SDL_GL_SetAttribute(SDL_GL_ALPHA_SIZE, 8);
        SDL_GL_SetAttribute(SDL_GL_STENCIL_SIZE, 8);
        SDL_GL_SetAttribute(SDL_GL_DEPTH_SIZE, 24);
        SDL_GL_SetAttribute(SDL_GL_DOUBLEBUFFER, 1);
        #ifdef ENABLE_MULTISAMPLE
            SDL_GL_SetAttribute(SDL_GL_MULTISAMPLEBUFFERS, 1);
            SDL_GL_SetAttribute(SDL_GL_MULTISAMPLESAMPLES, 4);
        #endif
        SDL_GL_SetAttribute(SDL_GL_SWAP_CONTROL, 1);
        SDL_SetVideoMode(Width, Height, 32, SDL_OPENGL);

        glewInit();
        Init();

        SDL_Event event;
        bool running = true;

        while(running){
            while(SDL_PollEvent(&event)){
                switch(event.type)
                {
                    case SDL_KEYDOWN:
                        if(event.key.keysym.sym == SDLK_ESCAPE)
                            running = false;
                    break;
                    case SDL_QUIT:
                        running = false;
                    break;
                }
            }
            Draw();
            SDL_GL_SwapBuffers();
        }
        SDL_Quit();
    }
    catch(std::bad_alloc& e)
    {
        std::cout << "Out of memory. " << e.what() << std::endl;
        exit(-1);
    }
    catch(std::exception& e)
    {
        std::cout << "Runtime exception: " << e.what() << std::endl;
        exit(-1);
    }
    catch(...)
    {
        std::cout << "Runtime exception of unknown type." << std::endl;
        exit(-1);
    }
    return 0;
}

void Init(void)
{
    const GLint texWidth = 256;
    const GLint texHeight = 256;
    const float texHalfWidth = 128.0f;
    const float texHalfHeight = 128.0f;
    printf("INIT: \n");

    unsigned char* pData = new unsigned char[texWidth*texHeight*4];
    for(int y=0; y<texHeight; ++y){
        for(int x=0; x<texWidth; ++x){
            int offs = (x + y*texWidth) * 4;
            float xoffs = ((float)x - texHalfWidth) / texHalfWidth;
            float yoffs = ((float)y - texHalfWidth) / texHalfHeight;
            float alpha = 1.0f - std::sqrt(xoffs*xoffs + yoffs*yoffs);
            if(alpha < 0.0f)
                alpha = 0.0f;
            pData[offs + 0] = 255; //r
            pData[offs + 1] = 0; //g
            pData[offs + 2] = 0; //b
            pData[offs + 3] = 255.0f * alpha; // * 
            //printf("alpha: %f\n", pData[x + y*texWidth + 3]);
        }
    }

    #ifdef ENABLE_TEXTURE
    glGenTextures(1, &texture_name);
    glActiveTexture(GL_TEXTURE0);
    glEnable(GL_TEXTURE_2D);
    glBindTexture(GL_TEXTURE_2D, texture_name);
    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, texWidth, texHeight, 0, GL_RGBA, GL_UNSIGNED_BYTE, pData);
    glEnable(GL_POINT_SPRITE);
    glTexEnvi(GL_POINT_SPRITE, GL_COORD_REPLACE, GL_TRUE);
    glTexEnvi(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
    #endif

    glPointSize(32.0f);

    glMatrixMode(GL_PROJECTION);
    glOrtho(0, Width, 0, Height, -1.0f, 1.0f);
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();

    glEnable(GL_BLEND);
    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
    glDisable(GL_DEPTH_TEST);

    #ifdef ENABLE_MULTISAMPLE
        glEnable(GL_POINT_SMOOTH);
    #endif

    GLenum e;
    do{
        e = glGetError();
        printf("%s\n",gluErrorString(e));
    } while(e != GL_NO_ERROR);

    delete [] pData;
}

void Draw(void)
{
    const int gridWidth = 1024;
    const int gridHeight = 1024;
    float t1, t2;

    t1 = t2 = (float)SDL_GetTicks() * 0.001f;
    t1 = fmod(t1, 10.0f) / 10.0f;
    t2 = fmod(t2, 4.0f) / 4.0f;
    float scale = 0.5f + (-sin(t2 * 2.0 * M_PI) + 1.0f) * 1.2f;
    //glColor4f(0.4f, 0.5f, 0.9f, 1.0f);
    glClear(GL_COLOR_BUFFER_BIT);
    glLoadIdentity();

    glTranslatef((Width>>1), (Height>>1), 0.0f);
    glScalef(scale,scale,scale);
    glRotatef(t1 * 360.0f, 0.0f, 0.0f, 1.0f);

    glBegin(GL_POINTS);
    for(int j=0; j<gridHeight; j+=64){
        for(int i=0; i<gridWidth; i+=64){ 
            glVertex2i(i-(gridWidth>>1),j-(gridHeight>>1));
        }
    }
    glEnd();
}

Mads' answer provides everything you need if you go for the fixed function pipeline. However, if you have a system that does not provide the ARB_point_sprite extension or with a broken implementation (some ATI drivers), you can solve this part also with geometry shaders. The ARB_geometry_shader4 extension allows you to convert a point primitive to two triangles, which can be used as the quad created by the ARB_point_sprite extension. On OpenGL 3.2, geometry shaders are already supported in core, no extension needed. The OpenGL wiki has two examples.