Generating random numbers a, b, c such that a^2 + b^2 + c^2 = 1

Question

To do some simulations in Python, I'm trying to generate numbers a,b,c such that a^2 + b^2 + c^2 = 1. I think generating some a between 0 and 1, then some b between 0 and sqrt(1 - a^2), and then c = sqrt(1 - a^2 - b^2) would work.

Floating point values are fine, the sum of squares should be close to 1. I want to keep generating them for some iterations.

Being new to Python, I'm not really sure how to do this. Negatives are allowed.

Edit: Thanks a lot for the answers!

Answer 1

According to this answer at stats.stackexchange.com, you should use normally distributed values to get uniformly distributed values on a sphere. That would mean, you could do:

import numpy as np
abc = np.random.normal(size=3)
a,b,c = abc/np.sqrt(sum(abc**2))

Answer 2

Just in case your interested in the probability densities I decided to do a comparison between the different approaches:

import numpy as np
import random
import math

def MSeifert():
    a = 1
    b = 1
    while a**2 + b**2 > 1:  # discard any a and b whose sum of squares already exceeds 1
        a = random.random()
        b = random.random()
    c = math.sqrt(1 - a**2 - b**2)  # fixed c
    return a, b, c

def VBB():
    x = np.random.uniform(0,1,3) # random numbers in [0, 1)
    x /= np.sqrt(x[0] ** 2 + x[1] ** 2 + x[2] ** 2)
    return x[0], x[1], x[2]

def user3684792():
    theta = random.uniform(0, 0.5*np.pi)
    phi = random.uniform(0, 0.5*np.pi)
    return np.sin(theta)* np.cos(phi), np.sin(theta)*np.sin(phi), np.cos(theta)

def JohanL():
    abc = np.random.normal(size=3)
    a,b,c = abc/np.sqrt(sum(abc**2))
    return a, b, c

def SeverinPappadeux():
    cos_th = 2.0*random.uniform(0, 1.0) - 1.0
    sin_th = math.sqrt(1.0 - cos_th*cos_th)
    phi = random.uniform(0, 2.0*math.pi)
    return sin_th * math.cos(phi), sin_th * math.sin(phi), cos_th

And plotting the distributions:

%matplotlib notebook

import matplotlib.pyplot as plt

f, axes = plt.subplots(3, 4)

for func_idx, func in enumerate([MSeifert, JohanL, user3684792, VBB]):
    axes[0, func_idx].set_title(str(func.__name__))
    res = [func() for _ in range(50000)]
    for idx in range(3):
        axes[idx, func_idx].hist([i[idx] for i in res], bins='auto')

axes[0, 0].set_ylabel('a')
axes[1, 0].set_ylabel('b')
axes[2, 0].set_ylabel('c')

plt.tight_layout()

With the result:

Explanation: The rows show the distributions for a, b and c respectively while the columns show the histograms (distributions) of the different approaches.

The only approaches that give a uniformly random distribution in the range (-1, 1) are JohanLs and Severin Pappadeux's approach. All other approaches have some features like spikes or a functional behavior in the range [0, 1). Note that these two solution currently gives values between -1 and 1 while all other approaches give values between 0 and 1.