function rng(integer modifier)
    while true do
        atom r1 := rnd()
        if rnd() < modifier(r1) then
            return r1
        end if
    end while
end function

function modifier(atom x)
    return iff(x < 0.5 ? 2 * (0.5 - x)
                       : 2 * (x - 0.5))
end function

constant N              = 100000,
         NUM_BINS       = 20,
         HIST_CHAR_SIZE = 125,
         BIN_SIZE       = 1/NUM_BINS,
         LO = sq_mul(tagset(NUM_BINS-1,0),BIN_SIZE),
         HI = sq_mul(tagset(NUM_BINS),BIN_SIZE),
         LBLS = apply(true,sprintf,{{"[%4.2f,%4.2f)"},columnize({LO,HI})})

sequence bins := repeat(0, NUM_BINS)
for i=1 to N do
    bins[floor(rng(modifier) / BIN_SIZE)+1] += 1
end for

printf(1,"Modified random distribution with %,d samples in range [0, 1):\n\n",N)
for i=1 to NUM_BINS do
    sequence hist := repeat('#', round(bins[i]/HIST_CHAR_SIZE))
    printf(1,"%s  %s %,d\n", {LBLS[i], hist, bins[i]})
end for