static package
bf_random(Var arglist, Byte next, void *vdata, Objid progr)
{
    int nargs = arglist.v.list[0].v.num;
    int num = (nargs >= 1 ? arglist.v.list[1].v.num : INTNUM_MAX);
    Var r;
    int e;
    int rnd;

    free_var(arglist);

    if (num <= 0)
	return make_error_pack(E_INVARG);

    const int range_l =
	((INTNUM_MAX > RAND_MAX ? RAND_MAX : (RAND_MAX - num)) + 1) % num;

    r.type = TYPE_INT;

#if ((RAND_MAX <= 0) || 0!=(RAND_MAX & (RAND_MAX+1)))
#   error RAND_MAX+1 is not a positive power of 2 ??
#endif

#if (INTNUM_MAX > RAND_MAX)
    /* num >= RAND_MAX possible; launch general algorithm */

#   define RANGE       (RAND_MAX+1)
#   define OR_ZERO(n)  (n)

    rnd = 0;
    e = 1;
#else
    /* num >= RAND_MAX not possible; unroll first loop iteration */

#   define OR_ZERO(n)  0

    rnd = RANDOM();
    e = range_l;

    if (rnd >= e) {
	r.v.num = 1 + rnd % num;
	return make_var_pack(r);
    }
#endif

    for (;;) {
	/* INVARIANT: rnd uniform over [0..e-1] */
	int rnd_next = RANDOM();

#if RAND_MAX < INTNUM_MAX
	/* compiler should turn [/%*]RANGE into bitwise ops */
	while (e < (num/RANGE) ||
	       ((e == num/RANGE) && (num%RANGE != 0))) {
	    rnd = rnd*RANGE + rnd_next;
	    e *= RANGE;
	    rnd_next = RANDOM();
	}
#endif
	/* INVARIANTS:
	 *   e*RANGE >= num
	 *   rnd uniform over [0..e-1]
	 *   rnd*RANGE + rnd_next uniform over [0..e*RANGE-1]
	 */
	if (rnd > OR_ZERO(num/RANGE)) {
	    /* rnd*RANGE > num */
	    r.v.num = 1 + muladdmod(rnd, range_l, rnd_next, num);
	    break;
	}
	rnd = OR_ZERO(rnd*RANGE) + rnd_next;
	e = muladdmod(e, range_l, 0, num);

	if (rnd >= e) {
	    r.v.num = 1 + rnd % num;
	    break;
	}
    }
    return make_var_pack(r);
#undef RANGE
#undef OR_ZERO
}