How to vectorise obtaining random deviates from a Bernoulli random variable?

Given a sequence of independent but not identically distributed Bernoulli trials with success probabilities given by a vector, e.g.:

x <- seq(0, 50, 0.1)
prob <- - x*(x - 50)/1000 # trial probabilities for trials 1 to 501

What is the most efficient way to obtain a random variate from each trial? I am assuming that vectorisation is the way to go.

I know of two functions that give Bernoulli random variates:

1. rbernoulli from the package purr, which does not accept a vector of success probabilities as an input. In this case it may be possible to wrap the function in an apply type operation.

2. rbinom with arguments size = 1 gives Bernoulli random variates. It also accepts a vector of probabilities, so that:

   rbinom(n = length(prob), size = 1, prob = prob)

gives an output with the right length. However, I am not entirely sure that this is actually what I want. The bits in the helpfile ?rbinom that seem relevant are:

The length of the result is determined by n for rbinom, and is the maximum of the lengths of the numerical arguments for the other functions.

The numerical arguments other than n are recycled to the length of the result. Only the first elements of the logical arguments are used.

However, n is a parameter with no default, so I am not sure what the first sentence means. I presume the second sentence means that I get what I want, since only size = 1 should be recycled. However this thread seems to suggest that this method does not work.

This blog post gives some other methods as well.