Return to Answer

EDIT: I see from the comments that you probably want an integer point -- that is, uniformly random from $$B_\mathbb{Z}(d,r)=\left\{ x \in \mathbb{Z}^d : \|x\|_2 = \sqrt{\sum_{i=1}^d x_i^2} \leq r \right\}.$$ Peter Shor's answer covers this case. Below is an alternative approach.

The approach is to sample one coordinate at a time.

We can partition $$B_\mathbb{Z}(d,r) = \bigcup_{a \in \mathbb{Z}} ~~ \{a\} \times B_\mathbb{Z}\left(d-1,\sqrt{r^2-a^2}\right).$$ Thus we can (i) sample $a \in [-r,r] \cap \mathbb{Z}$ from the appropriate marginal distribution and (ii) recursively sample from $B_\mathbb{Z}\left(d-1,\sqrt{r^2-a^2}\right)$. Combining the two gives a uniform sample from $B_\mathbb{Z}(d,r)$.

Define $f(d,r^2)=\left|B_\mathbb{Z}(d,r)\right|$. We first use $f$ to sample from $B_\mathbb{Z}(d,r)$ and then show how to compute $f$.

To sample $a$, we compute $f(d-1,r^2-a^2)$ for all $a \in [-r,r] \cap \mathbb{Z}$ and then scale these values to obtain a probability mass function for $a$.

All that remains is to determine how to compute $f$. We can use the above partition to obtain a recursive formula: $$f(d,r^2) = \sum_{a \in [-r,r] \cap \mathbb{Z}} f(d-1,r^2-a^2).$$ The base case of the recursion is $f(0,r^2)=1$. So we can compute all necessary values of $f$ using dynamic programming with a lookup table of size $O(r^2d)$.

Packing this into a single algorithm:

Algorithm 3: Recursion

• Input: $d,r^2 \in \mathbb{N}$.
• Output: $x \in B_\mathbb{Z}(d,r)$ uniformly random.

1. Allocate a lookup table for $f : \{0, \cdots, d\} \times \{0, \cdots, r^2\} \to \mathbb{N}$.
2. Initialize $f(0,s)=1$ for $s = 0 \cdots r^2$.
3. For $m = 1 \cdots d$:
4.    For $s = 0 \cdots r^2$:
5.       Set $f(m,s) = \sum_{a = \lceil - \sqrt{s} \rceil}^{\lfloor \sqrt{s} \rfloor} f(m-1, s-a^2 ).$
6. Set $s_1=r^2$ and allocate space for $x \in \{\lceil -r \rceil, \cdots ,\lfloor r \rfloor\}^d$.
7. For $i = 1 \cdots d$:
8.    Sample $U_i \in \{1, \cdots, ~ f(d+1-i,s_i)\}$ uniformly at random.
9.    Pick the largest $x_i \in \{\lceil - \sqrt{s_i} \rceil, \cdots, \lfloor \sqrt{s_i} \rfloor\}$ such that $$\sum_{a = \lceil - \sqrt{s_i} \rceil}^{x_i} f(d-i, s_i-a^2 ) \leq U_i.$$
10.    Set $s_{i+1}=s_i-x_i^2$.
11. Return $x \in B_\mathbb{Z}(d,r)$.

The running time of this algorithm is $O(r^3d)$. Assuming the $U_i$s are uniform and independent, it outputs a uniformly random point from $B_\mathbb{Z}(d,r)$.

---

Before EDIT: My understanding of your question is that you want to sample a point uniformly (i.e. Lebesgue measure) from the ball $$B(d,r)=\left\{ x \in \mathbb{R}^d : \|x\|_2 = \sqrt{\sum_{i=1}^d x_i^2} \leq r \right\}.$$ I will give two algorithms for this. For simplicity, I assume $r=1$, as you can always scale up.

---

Algorithm 1: Rejection Sampling

1. Sample $U_1, \cdots, U_d \in [-1,1]$ independently and uniformly at random.
2. Let $x = (U_1, \cdots, U_d)$.
3. If $\|x\|_2 \leq 1$, return $x$; else GOTO 1.

This algorithm is simple. We just use the fact that it is easy to sample from the "box" $[-1,1]^d$ that contains the unit ball. And, by rejecting samples that fall outside the ball, we are able to sample from the desired conditional distribution. The downside of this algorithm is that it takes time exponential in $d$, as the probability of a random point from the box being in the ball is decaying exponentially with $d$.

---

Algorithm 2: Scaled Spherical Sample

1. Sample $G_1, \cdots, G_d \in \mathbb{R}$ independently from a standard Gaussian distribution ($\mathcal{N}(0,1)$).
2. Let $x = (G_1, \cdots, G_d)$ and $y=x/\|x\|_2$.
3. Sample $U \in [0,1]$ uniformly at random.
4. Return $z=y \cdot U^{1/d}$.

Since the Gaussian distribution is spherically symmetric, $y$ is a uniformly random unit vector. Now we scale down $y$ by $U^{1/d}$ so that, instead of uniformly sampling from the sphere, we uniformly sample from the ball. To show that this is the right scaling, we just need to show that $$\mathbb{P}\left[\|z\|_2 \leq t\right] = \frac{\text{Volume of $d$-ball of radius $t$}}{\text{Volume of $d$-ball of radius $1$}}.$$ We have $$\frac{\text{Volume of $d$-ball of radius $t$}}{\text{Volume of $d$-ball of radius $1$}} = \frac{\frac{\pi^{d/2}}{\Gamma(d/2+1)} t^{d}}{\frac{\pi^{d/2}}{\Gamma(d/2+1)} 1^{d}} = t^{d}$$ and $$\mathbb{P}\left[\|z\|_2 \leq t\right] = \mathbb{P}\left[U^{1/d} \leq t\right] = \mathbb{P}\left[U \leq t^d\right] = t^d,$$ as required. Finally note that we can efficiently sample standard Gaussians using the Box-Muller transform.

EDIT: I see from the comments that you probably want an integer point -- that is, uniformly random from $$B_\mathbb{Z}(d,r)=\left\{ x \in \mathbb{Z}^d : \|x\|_2 = \sqrt{\sum_{i=1}^d x_i^2} \leq r \right\}.$$ Peter Shor's answer covers this case. Below is an alternative approach.

The approach is to sample one coordinate at a time.

We can partition $$B_\mathbb{Z}(d,r) = \bigcup_{a \in \mathbb{Z}} ~~ \{a\} \times B_\mathbb{Z}\left(d-1,\sqrt{r^2-a^2}\right).$$ Thus we can (i) sample $a \in [-r,r] \cap \mathbb{Z}$ from the appropriate marginal distribution and (ii) recursively sample from $B_\mathbb{Z}\left(d-1,\sqrt{r^2-a^2}\right)$. Combining the two gives a uniform sample from $B_\mathbb{Z}(d,r)$.

Define $f(d,r^2)=\left|B_\mathbb{Z}(d,r)\right|$. We first use $f$ to sample from $B_\mathbb{Z}(d,r)$ and then show how to compute $f$.

To sample $a$, we compute $f(d-1,r^2-a^2)$ for all $a \in [-r,r] \cap \mathbb{Z}$ and then scale these values to obtain a probability mass function for $a$.

All that remains is to determine how to compute $f$. We can use the above partition to obtain a recursive formula: $$f(d,r^2) = \sum_{a \in [-r,r] \cap \mathbb{Z}} f(d-1,r^2-a^2).$$ The base case of the recursion is $f(0,r^2)=1$. So we can compute all necessary values of $f$ using dynamic programming with a lookup table of size $O(r^2d)$.

Packing this into a single algorithm:

Algorithm 3: Recursion

• Input: $d,r^2 \in \mathbb{N}$.
• Output: $x \in B_\mathbb{Z}(d,r)$ uniformly random.

1. Allocate a lookup table for $f : \{0, \cdots, d\} \times \{0, \cdots, r^2\} \to \mathbb{N}$.
2. Initialize $f(0,s)=1$ for $s = 0 \cdots r^2$.
3. For $m = 1 \cdots d$:
4.    For $s = 0 \cdots r^2$:
5.       Set $f(m,s) = \sum_{a = \lceil - \sqrt{s} \rceil}^{\lfloor \sqrt{s} \rfloor} f(m-1, s-a^2 ).$
6. Set $s_1=r^2$ and allocate space for $x \in \{\lceil -r \rceil, \cdots ,\lfloor r \rfloor\}^d$.
7. For $i = 1 \cdots d$:
8.    Sample $U_i \in \{1, \cdots, ~ f(d+1-i,s_i)\}$ uniformly at random.
9.    Pick the largest $x_i \in \{\lceil - \sqrt{s_i} \rceil, \cdots, \lfloor \sqrt{s_i} \rfloor\}$ such that $$\sum_{a = \lceil - \sqrt{s_i} \rceil}^{x_i} f(d-i, s_i-a^2 ) \leq U_i.$$
10.    Set $s_{i+1}=s_i-x_i^2$.
11. Return $x \in B_\mathbb{Z}(d,r)$.

The running time of this algorithm is $O(r^3d)$. Assuming the $U_i$s are uniform and independent, it outputs a uniformly random point from $B_\mathbb{Z}(d,r)$.

---

Before EDIT: My understanding of your question is that you want to sample a point uniformly (i.e. Lebesgue measure) from the ball $$B(d,r)=\left\{ x \in \mathbb{R}^d : \|x\|_2 = \sqrt{\sum_{i=1}^d x_i^2} \leq r \right\}.$$ I will give two algorithms for this. For simplicity, I assume $r=1$, as you can always scale up.

---

Algorithm 1: Rejection Sampling

1. Sample $U_1, \cdots, U_d \in [-1,1]$ independently and uniformly at random.
2. Let $x = (U_1, \cdots, U_d)$.
3. If $\|x\|_2 \leq 1$, return $x$; else GOTO 1.

This algorithm is simple. We just use the fact that it is easy to sample from the "box" $[-1,1]^d$ that contains the unit ball. And, by rejecting samples that fall outside the ball, we are able to sample from the desired conditional distribution. The downside of this algorithm is that it takes time exponential in $d$, as the probability of a random point from the box being in the ball is decaying exponentially with $d$.

---

Algorithm 2: Scaled Spherical Sample

1. Sample $G_1, \cdots, G_d \in \mathbb{R}$ independently from a standard Gaussian distribution ($\mathcal{N}(0,1)$).
2. Let $x = (G_1, \cdots, G_d)$ and $y=x/\|x\|_2$.
3. Sample $U \in [0,1]$ uniformly at random.
4. Return $z=y \cdot U^{1/d}$.

Since the Gaussian distribution is spherically symmetric, $y$ is a uniformly random unit vector. Now we scale down $y$ by $U^{1/d}$ so that, instead of uniformly sampling from the sphere, we uniformly sample from the ball. To show that this is the right scaling, we just need to show that $$\mathbb{P}\left[\|z\|_2 \leq t\right] = \frac{\text{Volume of $d$-ball of radius $t$}}{\text{Volume of $d$-ball of radius $1$}}.$$ We have $$\frac{\text{Volume of $d$-ball of radius $t$}}{\text{Volume of $d$-ball of radius $1$}} = \frac{\frac{\pi^{d/2}}{\Gamma(d/2+1)} t^{d}}{\frac{\pi^{d/2}}{\Gamma(d/2+1)} 1^{d}} = t^{d}$$ and $$\mathbb{P}\left[\|z\|_2 \leq t\right] = \mathbb{P}\left[U^{1/d} \leq t\right] = \mathbb{P}\left[U \leq t^d\right] = t^d,$$ as required. Finally note that we can efficiently sample standard Gaussians using the Box-Muller transform.

EDIT: I see from the comments that you probably want an integer point -- that is, uniformly random from $$B_\mathbb{Z}(d,r)=\left\{ x \in \mathbb{Z}^d : \|x\|_2 = \sqrt{\sum_{i=1}^d x_i^2} \leq r \right\}.$$ Peter Shor's answer covers this case. Below is an alternative approach.

The approach is to sample one coordinate at a time.

We can partition $$B_\mathbb{Z}(d,r) = \bigcup_{a \in \mathbb{Z}} ~~ \{a\} \times B_\mathbb{Z}\left(d-1,\sqrt{r^2-a^2}\right).$$ Thus we can (i) sample $a \in [-r,r] \cap \mathbb{Z}$ from the appropriate marginal distribution and (ii) recursively sample from $B_\mathbb{Z}\left(d-1,\sqrt{r^2-a^2}\right)$. Combining the two gives a uniform sample from $B_\mathbb{Z}(d,r)$.

Define $f(d,r^2)=\left|B_\mathbb{Z}(d,r)\right|$. We first use $f$ to sample from $B_\mathbb{Z}(d,r)$ and then show how to compute $f$.

To sample $a$, we compute $f(d-1,r^2-a^2)$ for all $a \in [-r,r] \cap \mathbb{Z}$ and then scale these values to obtain a probability mass function for $a$.

All that remains is to determine how to compute $f$. We can use the above partition to obtain a recursive formula: $$f(d,r^2) = \sum_{a \in [-r,r] \cap \mathbb{Z}} f(d-1,r^2-a^2).$$ The base case of the recursion is $f(0,r^2)=1$. So we can compute all necessary values of $f$ using dynamic programming with a lookup table of size $O(r^2d)$.

Packing this into a single algorithm:

Algorithm 3: Recursion

1. Input: $d,r^2 \in \mathbb{N}$.
Allocate a lookup table for $f : \{0, \cdots, d\} \times \{0, \cdots, r^2\} \to \mathbb{Z}$.
3. Initialize $f(0,s)=1$ for $s = 0 \cdots r^2$.
4. For $m = 1 \cdots d$:
5.    For $s = 0 \cdots r^2$:
6.       Set $f(m,s) = \sum_{a = \lceil - \sqrt{s} \rceil}^{\lfloor \sqrt{s} \rfloor} f(m-1, s-a^2 ).$
7. Set $s_1=r^2$ and allocate space for $x \in \{\lceil -r \rceil, \cdots ,\lfloor r \rfloor\}^d$.
8. For $i = 1 \cdots d$:
9.    Sample $U_i \in \{1, \cdots, f(d+1-i,s_i)\}$ uniformly at random.
10.    Pick the largest $x_i \in \{\lceil - \sqrt{s_i} \rceil, \cdots, \lfloor \sqrt{s_i} \rfloor\}$ such that $$\sum_{a = \lceil - \sqrt{s_i} \rceil}^{x_i} f(d-i, s_i-a^2 ) \leq U_i.$$
11.    Set $s_{i+1}=s_i-x_i^2$.
12. Return $x \in B_\mathbb{Z}(d,r)$.

The running time of this algorithm is $O(r^3d)$. Assuming the $U_i$s are uniform and independent, it outputs a uniformly random point from $B_\mathbb{Z}(d,r)$.

---

Before EDIT: My understanding of your question is that you want to sample a point uniformly (i.e. Lebesgue measure) from the ball $$B(d,r)=\left\{ x \in \mathbb{R}^d : \|x\|_2 = \sqrt{\sum_{i=1}^d x_i^2} \leq r \right\}.$$ I will give two algorithms for this. For simplicity, I assume $r=1$, as you can always scale up.

---

Algorithm 1: Rejection Sampling

1. Sample $U_1, \cdots, U_d \in [-1,1]$ independently and uniformly at random.
2. Let $x = (U_1, \cdots, U_d)$.
3. If $\|x\|_2 \leq 1$, return $x$; else GOTO 1.

This algorithm is simple. We just use the fact that it is easy to sample from the "box" $[-1,1]^d$ that contains the unit ball. And, by rejecting samples that fall outside the ball, we are able to sample from the desired conditional distribution. The downside of this algorithm is that it takes time exponential in $d$, as the probability of a random point from the box being in the ball is decaying exponentially with $d$.

---

Algorithm 2: Scaled Spherical Sample

1. Sample $G_1, \cdots, G_d \in \mathbb{R}$ independently from a standard Gaussian distribution ($\mathcal{N}(0,1)$).
2. Let $x = (G_1, \cdots, G_d)$ and $y=x/\|x\|_2$.
3. Sample $U \in [0,1]$ uniformly at random.
4. Return $z=y \cdot U^{1/d}$.

Since the Gaussian distribution is spherically symmetric, $y$ is a uniformly random unit vector. Now we scale down $y$ by $U^{1/d}$ so that, instead of uniformly sampling from the sphere, we uniformly sample from the ball. To show that this is the right scaling, we just need to show that $$\mathbb{P}\left[\|z\|_2 \leq t\right] = \frac{\text{Volume of $d$-ball of radius $t$}}{\text{Volume of $d$-ball of radius $1$}}.$$ We have $$\frac{\text{Volume of $d$-ball of radius $t$}}{\text{Volume of $d$-ball of radius $1$}} = \frac{\frac{\pi^{d/2}}{\Gamma(d/2+1)} t^{d}}{\frac{\pi^{d/2}}{\Gamma(d/2+1)} 1^{d}} = t^{d}$$ and $$\mathbb{P}\left[\|z\|_2 \leq t\right] = \mathbb{P}\left[U^{1/d} \leq t\right] = \mathbb{P}\left[U \leq t^d\right] = t^d,$$ as required. Finally note that we can efficiently sample standard Gaussians using the Box-Muller transform.

EDIT: I see from the comments that you probably want an integer point -- that is, uniformly random from $$B_\mathbb{Z}(d,r)=\left\{ x \in \mathbb{Z}^d : \|x\|_2 = \sqrt{\sum_{i=1}^d x_i^2} \leq r \right\}.$$ Peter Shor's answer covers this case. Below is an alternative approach.

The approach is to sample one coordinate at a time.

We can partition $$B_\mathbb{Z}(d,r) = \bigcup_{a \in \mathbb{Z}} ~~ \{a\} \times B_\mathbb{Z}\left(d-1,\sqrt{r^2-a^2}\right).$$ Thus we can (i) sample $a \in [-r,r] \cap \mathbb{Z}$ from the appropriate marginal distribution and (ii) recursively sample from $B_\mathbb{Z}\left(d-1,\sqrt{r^2-a^2}\right)$. Combining the two gives a uniform sample from $B_\mathbb{Z}(d,r)$.

Define $f(d,r^2)=\left|B_\mathbb{Z}(d,r)\right|$. We first use $f$ to sample from $B_\mathbb{Z}(d,r)$ and then show how to compute $f$.

To sample $a$, we compute $f(d-1,r^2-a^2)$ for all $a \in [-r,r] \cap \mathbb{Z}$ and then scale these values to obtain a probability mass function for $a$.

All that remains is to determine how to compute $f$. We can use the above partition to obtain a recursive formula: $$f(d,r^2) = \sum_{a \in [-r,r] \cap \mathbb{Z}} f(d-1,r^2-a^2).$$ The base case of the recursion is $f(0,r^2)=1$. So we can compute all necessary values of $f$ using dynamic programming with a lookup table of size $O(r^2d)$.

Packing this into a single algorithm:

Algorithm 3: Recursion

• Input: $d,r^2 \in \mathbb{N}$.
• Output: $x \in B_\mathbb{Z}(d,r)$ uniformly random.

1. Allocate a lookup table for $f : \{0, \cdots, d\} \times \{0, \cdots, r^2\} \to \mathbb{N}$.
2. Initialize $f(0,s)=1$ for $s = 0 \cdots r^2$.
3. For $m = 1 \cdots d$:
4.    For $s = 0 \cdots r^2$:
5.       Set $f(m,s) = \sum_{a = \lceil - \sqrt{s} \rceil}^{\lfloor \sqrt{s} \rfloor} f(m-1, s-a^2 ).$
6. Set $s_1=r^2$ and allocate space for $x \in \{\lceil -r \rceil, \cdots ,\lfloor r \rfloor\}^d$.
7. For $i = 1 \cdots d$:
8.    Sample $U_i \in \{1, \cdots, ~ f(d+1-i,s_i)\}$ uniformly at random.
9.    Pick the largest $x_i \in \{\lceil - \sqrt{s_i} \rceil, \cdots, \lfloor \sqrt{s_i} \rfloor\}$ such that $$\sum_{a = \lceil - \sqrt{s_i} \rceil}^{x_i} f(d-i, s_i-a^2 ) \leq U_i.$$
10.    Set $s_{i+1}=s_i-x_i^2$.
11. Return $x \in B_\mathbb{Z}(d,r)$.

The running time of this algorithm is $O(r^3d)$. Assuming the $U_i$s are uniform and independent, it outputs a uniformly random point from $B_\mathbb{Z}(d,r)$.

---

Before EDIT: My understanding of your question is that you want to sample a point uniformly (i.e. Lebesgue measure) from the ball $$B(d,r)=\left\{ x \in \mathbb{R}^d : \|x\|_2 = \sqrt{\sum_{i=1}^d x_i^2} \leq r \right\}.$$ I will give two algorithms for this. For simplicity, I assume $r=1$, as you can always scale up.

---

Algorithm 1: Rejection Sampling

1. Sample $U_1, \cdots, U_d \in [-1,1]$ independently and uniformly at random.
2. Let $x = (U_1, \cdots, U_d)$.
3. If $\|x\|_2 \leq 1$, return $x$; else GOTO 1.

This algorithm is simple. We just use the fact that it is easy to sample from the "box" $[-1,1]^d$ that contains the unit ball. And, by rejecting samples that fall outside the ball, we are able to sample from the desired conditional distribution. The downside of this algorithm is that it takes time exponential in $d$, as the probability of a random point from the box being in the ball is decaying exponentially with $d$.

---

Algorithm 2: Scaled Spherical Sample

1. Sample $G_1, \cdots, G_d \in \mathbb{R}$ independently from a standard Gaussian distribution ($\mathcal{N}(0,1)$).
2. Let $x = (G_1, \cdots, G_d)$ and $y=x/\|x\|_2$.
3. Sample $U \in [0,1]$ uniformly at random.
4. Return $z=y \cdot U^{1/d}$.

Since the Gaussian distribution is spherically symmetric, $y$ is a uniformly random unit vector. Now we scale down $y$ by $U^{1/d}$ so that, instead of uniformly sampling from the sphere, we uniformly sample from the ball. To show that this is the right scaling, we just need to show that $$\mathbb{P}\left[\|z\|_2 \leq t\right] = \frac{\text{Volume of $d$-ball of radius $t$}}{\text{Volume of $d$-ball of radius $1$}}.$$ We have $$\frac{\text{Volume of $d$-ball of radius $t$}}{\text{Volume of $d$-ball of radius $1$}} = \frac{\frac{\pi^{d/2}}{\Gamma(d/2+1)} t^{d}}{\frac{\pi^{d/2}}{\Gamma(d/2+1)} 1^{d}} = t^{d}$$ and $$\mathbb{P}\left[\|z\|_2 \leq t\right] = \mathbb{P}\left[U^{1/d} \leq t\right] = \mathbb{P}\left[U \leq t^d\right] = t^d,$$ as required. Finally note that we can efficiently sample standard Gaussians using the Box-Muller transform.

EDIT: I see from the comments that you probably want an integer point -- that is, uniformly random from $$B_\mathbb{Z}(d,r)=\left\{ x \in \mathbb{Z}^d : \|x\|_2 = \sqrt{\sum_{i=1}^d x_i^2} \leq r \right\}.$$ Peter Shor's answer covers this case. Below is an alternative approach.

The approach is to sample one coordinate at a time.

We can partition $$B_\mathbb{Z}(d,r) = \bigcup_{a \in \mathbb{Z}} \{a\} \times B_\mathbb{Z}\left(d-1,\sqrt{r^2-a^2}\right).$$ Thus we can (i) sample $a \in [-r,r] \cap \mathbb{Z}$ from the appropriate marginal distribution and (ii) recursively sample from $B_\mathbb{Z}\left(d-1,\sqrt{r^2-a^2}\right)$. Combining the two gives a uniform sample from $B_\mathbb{Z}(d,r)$.

Define $f(d,r^2)=\left|B_\mathbb{Z}(d,r)\right|$. We first use $f$ to sample from $B_\mathbb{Z}(d,r)$ and then show how to compute $f$.

To sample $a$, we compute $f(d,r^2-a^2)$ for all $a \in [-r,r] \cap \mathbb{Z}$ and then scale these values to gain a probability mass function for $a$.

All that remains is to determine how to compute $f$. We can use the partition to obtain a recursive formula $$f(d,r^2) = \sum_{a \in [-r,r] \cap \mathbb{Z}} f(d-1,r^2-a^2).$$ The base case of the recursion is $f(0,r^2)=1$. So we can compute all necessary values of $f$ using dynamic programming with a lookup table of size $O(r^2d)$.

Packing this into a single algorithm:

Algorithm 3: Recursion

1. Input: $d,r^2 \in \mathbb{N}$.
2. Allocate a lookup table for $f : \{0, \cdots, d\} \times \{0, \cdots, r^2\} \to \mathbb{Z}$.
3. Initialize $f(0,s)=1$ for $s = 0 \cdots r^2$.
4. For $m = 1 \cdots d$:
5.    For $s = 0 \cdots r^2$:
6.       Set $f(m,s) = \sum_{a = \lceil - \sqrt{s} \rceil}^{\lfloor \sqrt{s} \rfloor} f(m-1, s-a^2 ).$
7. Set $s_1=r^2$ and allocate space for $x \in \{\lceil -r \rceil, \cdots ,\lfloor r \rfloor\}^d$.
8. For $i = 1 \cdots d$:
9.    Sample $U_i \in \{1, \cdots, f(d+1-i,s_i)\}$ uniformly at random.
10.    Pick the largest $x_i \in \{\lceil - \sqrt{s_i} \rceil, \cdots, \lfloor \sqrt{s_i} \rfloor\}$ such that $$\sum_{a = \lceil - \sqrt{s_i} \rceil}^{x_i} f(d-i, s_i-a^2 ) \leq U_i.$$
11.    Set $s_{i+1}=s_i-x_i^2$.
12. Return $x \in B_\mathbb{Z}(d,r)$.

The running time of this algorithm is $O(r^3d)$.

---

Before EDIT: My understanding of your question is that you want to sample a point uniformly (i.e. Lebesgue measure) from the ball $$B(d,r)=\left\{ x \in \mathbb{R}^d : \|x\|_2 = \sqrt{\sum_{i=1}^d x_i^2} \leq r \right\}.$$ I will give two algorithms for this. For simplicity, I assume $r=1$, as you can always scale up.

---

Algorithm 1: Rejection Sampling

1. Sample $U_1, \cdots, U_d \in [-1,1]$ independently and uniformly at random.
2. Let $x = (U_1, \cdots, U_d)$.
3. If $\|x\|_2 \leq 1$, return $x$; else GOTO 1.

This algorithm is simple. We just use the fact that it is easy to sample from the "box" $[-1,1]^d$ that contains the unit ball. And, by rejecting samples that fall outside the ball, we are able to sample from the desired conditional distribution. The downside of this algorithm is that it takes time exponential in $d$, as the probability of a random point from the box being in the ball is decaying exponentially with $d$.

---

Algorithm 2: Scaled Spherical Sample

1. Sample $G_1, \cdots, G_d \in \mathbb{R}$ independently from a standard Gaussian distribution ($\mathcal{N}(0,1)$).
2. Let $x = (G_1, \cdots, G_d)$ and $y=x/\|x\|_2$.
3. Sample $U \in [0,1]$ uniformly at random.
4. Return $z=y \cdot U^{1/d}$.

Since the Gaussian distribution is spherically symmetric, $y$ is a uniformly random unit vector. Now we scale down $y$ by $U^{1/d}$ so that, instead of uniformly sampling from the sphere, we uniformly sample from the ball. To show that this is the right scaling, we just need to show that $$\mathbb{P}\left[\|z\|_2 \leq t\right] = \frac{\text{Volume of $d$-ball of radius $t$}}{\text{Volume of $d$-ball of radius $1$}}.$$ We have $$\frac{\text{Volume of $d$-ball of radius $t$}}{\text{Volume of $d$-ball of radius $1$}} = \frac{\frac{\pi^{d/2}}{\Gamma(d/2+1)} t^{d}}{\frac{\pi^{d/2}}{\Gamma(d/2+1)} 1^{d}} = t^{d}$$ and $$\mathbb{P}\left[\|z\|_2 \leq t\right] = \mathbb{P}\left[U^{1/d} \leq t\right] = \mathbb{P}\left[U \leq t^d\right] = t^d,$$ as required. Finally note that we can efficiently sample standard Gaussians using the Box-Muller transform.

EDIT: I see from the comments that you probably want an integer point -- that is, uniformly random from $$B_\mathbb{Z}(d,r)=\left\{ x \in \mathbb{Z}^d : \|x\|_2 = \sqrt{\sum_{i=1}^d x_i^2} \leq r \right\}.$$ Peter Shor's answer covers this case. Below is an alternative approach.

The approach is to sample one coordinate at a time.

We can partition $$B_\mathbb{Z}(d,r) = \bigcup_{a \in \mathbb{Z}} ~~ \{a\} \times B_\mathbb{Z}\left(d-1,\sqrt{r^2-a^2}\right).$$ Thus we can (i) sample $a \in [-r,r] \cap \mathbb{Z}$ from the appropriate marginal distribution and (ii) recursively sample from $B_\mathbb{Z}\left(d-1,\sqrt{r^2-a^2}\right)$. Combining the two gives a uniform sample from $B_\mathbb{Z}(d,r)$.

Define $f(d,r^2)=\left|B_\mathbb{Z}(d,r)\right|$. We first use $f$ to sample from $B_\mathbb{Z}(d,r)$ and then show how to compute $f$.

To sample $a$, we compute $f(d-1,r^2-a^2)$ for all $a \in [-r,r] \cap \mathbb{Z}$ and then scale these values to obtain a probability mass function for $a$.

All that remains is to determine how to compute $f$. We can use the above partition to obtain a recursive formula: $$f(d,r^2) = \sum_{a \in [-r,r] \cap \mathbb{Z}} f(d-1,r^2-a^2).$$ The base case of the recursion is $f(0,r^2)=1$. So we can compute all necessary values of $f$ using dynamic programming with a lookup table of size $O(r^2d)$.

Packing this into a single algorithm:

Algorithm 3: Recursion

1. Input: $d,r^2 \in \mathbb{N}$.
2. Allocate a lookup table for $f : \{0, \cdots, d\} \times \{0, \cdots, r^2\} \to \mathbb{Z}$.
3. Initialize $f(0,s)=1$ for $s = 0 \cdots r^2$.
4. For $m = 1 \cdots d$:
5.    For $s = 0 \cdots r^2$:
6.       Set $f(m,s) = \sum_{a = \lceil - \sqrt{s} \rceil}^{\lfloor \sqrt{s} \rfloor} f(m-1, s-a^2 ).$
7. Set $s_1=r^2$ and allocate space for $x \in \{\lceil -r \rceil, \cdots ,\lfloor r \rfloor\}^d$.
8. For $i = 1 \cdots d$:
9.    Sample $U_i \in \{1, \cdots, f(d+1-i,s_i)\}$ uniformly at random.
10.    Pick the largest $x_i \in \{\lceil - \sqrt{s_i} \rceil, \cdots, \lfloor \sqrt{s_i} \rfloor\}$ such that $$\sum_{a = \lceil - \sqrt{s_i} \rceil}^{x_i} f(d-i, s_i-a^2 ) \leq U_i.$$
11.    Set $s_{i+1}=s_i-x_i^2$.
12. Return $x \in B_\mathbb{Z}(d,r)$.

The running time of this algorithm is $O(r^3d)$. Assuming the $U_i$s are uniform and independent, it outputs a uniformly random point from $B_\mathbb{Z}(d,r)$.

---

Before EDIT: My understanding of your question is that you want to sample a point uniformly (i.e. Lebesgue measure) from the ball $$B(d,r)=\left\{ x \in \mathbb{R}^d : \|x\|_2 = \sqrt{\sum_{i=1}^d x_i^2} \leq r \right\}.$$ I will give two algorithms for this. For simplicity, I assume $r=1$, as you can always scale up.

---

Algorithm 1: Rejection Sampling

1. Sample $U_1, \cdots, U_d \in [-1,1]$ independently and uniformly at random.
2. Let $x = (U_1, \cdots, U_d)$.
3. If $\|x\|_2 \leq 1$, return $x$; else GOTO 1.

This algorithm is simple. We just use the fact that it is easy to sample from the "box" $[-1,1]^d$ that contains the unit ball. And, by rejecting samples that fall outside the ball, we are able to sample from the desired conditional distribution. The downside of this algorithm is that it takes time exponential in $d$, as the probability of a random point from the box being in the ball is decaying exponentially with $d$.

---

Algorithm 2: Scaled Spherical Sample

1. Sample $G_1, \cdots, G_d \in \mathbb{R}$ independently from a standard Gaussian distribution ($\mathcal{N}(0,1)$).
2. Let $x = (G_1, \cdots, G_d)$ and $y=x/\|x\|_2$.
3. Sample $U \in [0,1]$ uniformly at random.
4. Return $z=y \cdot U^{1/d}$.

Since the Gaussian distribution is spherically symmetric, $y$ is a uniformly random unit vector. Now we scale down $y$ by $U^{1/d}$ so that, instead of uniformly sampling from the sphere, we uniformly sample from the ball. To show that this is the right scaling, we just need to show that $$\mathbb{P}\left[\|z\|_2 \leq t\right] = \frac{\text{Volume of $d$-ball of radius $t$}}{\text{Volume of $d$-ball of radius $1$}}.$$ We have $$\frac{\text{Volume of $d$-ball of radius $t$}}{\text{Volume of $d$-ball of radius $1$}} = \frac{\frac{\pi^{d/2}}{\Gamma(d/2+1)} t^{d}}{\frac{\pi^{d/2}}{\Gamma(d/2+1)} 1^{d}} = t^{d}$$ and $$\mathbb{P}\left[\|z\|_2 \leq t\right] = \mathbb{P}\left[U^{1/d} \leq t\right] = \mathbb{P}\left[U \leq t^d\right] = t^d,$$ as required. Finally note that we can efficiently sample standard Gaussians using the Box-Muller transform.