How to model 14,16 and 18 sided dice?
Equal distribution on sphere
Using the script from This answer
Need to alter the parameters. Here is what I used for the 14, 16, 18 side die. (where n is number of sides.  Had to run the 16 side twice, as it didn't find solution the first time, kill the script by hitting ⎈ Ctrl C
 over the system console.  The rat value determines the size of the hole. Too small 0.8 and the holes overlap, too big 0.999and the holes will be too small
# parameters
n = 18  # number of points on sphere
rat = (n - 1) / n # how far along the radius to bisect
u_segments = 32  # UV sphere settings
v_segments = 32 
thickness = 0.2  # solidify thickness
TOL = 1e-7


14, 16, and 18 sided die, ret = (n - 1) / n
To make the die &quot;solid&quot; remove the solidify modifier, and grid fill each circular hole.
Solid Version
Here is an edit to create a solid version, where the holes can overlap. The holes are filled with an ngon. Poke it to make tris.

14, 16, 18 sided solid, where ret = (n - 2) / n.
Edit update for 2.8x for prior see revisions
import bpy
import bmesh
from mathutils import Vector
import random
from math import pi, asin, atan2, cos, sin, radians

# parmaaters
n = 23  # number of points on sphere
rat = (n - 2) / n # how far along the radius to bisect
u_segments = 32  # UV sphere settings
v_segments = 32 
thickness = 0.2  # solidify thickness
TOL = 1e-5

points = [Vector((0, 0, 1))]
for i in range(n - 1):
    theta = random.random() * radians(360)
    phi = 2 * asin(random.random() * 2 - 1)
    points.append(Vector((cos(theta) * cos(phi), 
           sin(theta) * cos(phi), 
           sin(phi))))

while True:
    # Determine the total force acting on each point.
    forces = []
    for i in range(len(points)):
        p = points[i]
        f = Vector()
        ftotal = 0
        for j in range(len(points)):
            if j == i: continue
            q = points[j]
            # Find the distance vector, and its length.
            dv = p - q
            dl = dv.length
            dl3 = dl * dl * dl
            fv = dv / dl3
            # Add to the total force on the point p.
            f = f + fv
        # Stick this in the forces array.
        forces.append(f)
        # Add to the running sum of the total forces/distances.
        ftotal = ftotal + f.length

    fscale = 1 if ftotal &lt;= 0.25 else 0.25 / ftotal

    # Move each point, and normalise. While we do this, also track
    # the distance each point ends up moving.
    dist = 0
    for i in range(len(points)):
        p = points[i]
        f = forces[i]
        p2 = (p + fscale * f).normalized()

        dv = p - p2
        dist = dist + dv.length
        points[i] = p2
    # Done. Check for convergence and finish.
    if dist &lt; TOL: # TOL
        break

context = bpy.context
view_layer = context.view_layer
scene = context.scene
coll = context.collection or scene.collection

# make one point north pole.
R = points[0].rotation_difference(Vector((0, 0, 1))).to_matrix()
points = [R @ p for p in points]

bm = bmesh.new()
#bmesh.ops.create_icosphere(bm, diameter=1, subdivisions=5 )
bmesh.ops.create_uvsphere(bm, 
        diameter=1, 
        u_segments=u_segments, 
        v_segments=v_segments)

for p in points:
    ret = bmesh.ops.bisect_plane(bm, 
            geom=bm.faces[:]+bm.edges[:]+bm.verts[:],
            plane_co= rat * p,
            plane_no=-p,
            clear_outer=False,
            clear_inner=True)

    # fill the holes
    bmesh.ops.contextual_create(
            bm, 
            geom=[e for e in ret[&quot;geom_cut&quot;] 
                  if isinstance(e, bmesh.types.BMEdge)]
            )

me = bpy.data.meshes.new(f&quot;{n} Sided Dice&quot;)
bm.to_mesh(me)
ob = bpy.data.objects.new(me.name, me)   
coll.objects.link(ob)
view_layer.objects.active = ob
ob.select_set(True)
ob.location = scene.cursor.location