How does Blender figure out the normal of a face?
Yes, the order in which we define the vertex indices of the Face directly influence the Face's normal. 

See this thread on StackOverflow for some insight 

Usually when constructing geometry you can keep track of the order and never encounter this problem, just be consistent. If you suddenly change the order in which you specify the vertices, you will get an inverted normal on that face (which is what happens to you). There's no way for Blender to know if that's intentional or not.



Sometimes it can be more convenient to let an algorithm recalculate normals at the end, as an extra processing step.

bmesh.ops

There's a bmesh.ops to recalculate normals . 

bmesh.ops.recalc_face_normals(bm, faces=bm.faces)


See TextEditor > Templates > Python > Bmesh Simple for the template to get a bmesh representation from an object.

import bpy
import bmesh

verts = [(1.0, 1.0, -1.0),
         (1.0, -1.0, -1.0),
         (-1.0, -1.0, -1.0),
         (-1.0, 1.0, -1.0),
         (1.0, 1.0, 1.0),
         (1.0, -1.0, 1.0),
         (-1.0, -1.0, 1.0),
         (-1.0, 1.0, 1.0)]

faces = [(0, 1, 2, 3),
         (4, 5, 6, 7),         # (4, 7, 6, 5) --&gt; (4, 5, 6, 7)
         (0, 4, 5, 1),
         (1, 5, 6, 2),
         (2, 6, 7, 3),
         (4, 0, 3, 7)]

mesh = bpy.data.meshes.new('CubeMesh')
scn = bpy.context.scene
mesh.from_pydata(verts, [], faces) 
mesh.update()

bm = bmesh.new()
bm.from_mesh(mesh)
bmesh.ops.recalc_face_normals(bm, faces=bm.faces)
bm.to_mesh(mesh)
bm.free()

cube = bpy.data.objects.new('Cube', mesh)
scn.objects.link(cube)


You may want to work this into a function for clarity.

bpy.ops.mesh

You can use normals_make_consistent() if it's more convenient, there are however general reasons to avoid using bpy.ops within a scripting context. See the post on why to avoid bpy.ops

def make_object(name, verts, faces, normal_recalc=True):

    scn = bpy.context.scene

    mesh = bpy.data.meshes.new(name + "_Mesh")
    mesh.from_pydata(verts, [], faces) 
    mesh.update()

    ob = bpy.data.objects.new(name, mesh)
    scn.objects.link(ob)
    if normal_recalc:
        scn.objects.active = ob
        bpy.ops.object.mode_set(mode = 'EDIT')
        bpy.ops.mesh.normals_make_consistent()
        bpy.ops.object.mode_set(mode = 'OBJECT')

    return ob

make_object('my_cube', verts, faces)
# make_object(name, verts, faces, normal_recalc=False)