[Question]
  [
## The Challenge
Write a program or function that takes no input and outputs a 3-dimensional vector of length \$1\$ in a **theoretically uniform** random direction.
This is equivalent to a random point on the sphere described by $$x^2+y^2+z^2=1$$
resulting in a distribution like such
[![Random distribution of points on a sphere with radius 1.](https://i.stack.imgur.com/xM4hO.png)](https://i.stack.imgur.com/xM4hO.png)
### Output
Three floats from a theoretically uniform random distribution for which the equation \$x^2+y^2+z^2=1\$ holds true to precision limits.
## Challenge remarks
* The random distribution needs to be **theoretically uniform**. That is, if the pseudo-random number generator were to be replaced with a true RNG from the **real** numbers, it would result in a uniform random distribution of points on the sphere.
* Generating three random numbers from a uniform distribution and normalizing them is invalid: there will be a bias towards the corners of the three-dimensional space.
* Similarly, generating two random numbers from a uniform distribution and using them as spherical coordinates is invalid: there will be a bias towards the poles of the sphere.
* Proper uniformity can be achieved by algorithms including but not limited to:
	1. Generate three random numbers \$x\$, \$y\$ and \$z\$ from a *normal* (Gaussian) distribution around \$0\$ and normalize them.
		+ [Implementation example](https://tio.run/##ZVDLasMwELzrKxZyiBSKW5pboLdALz2WXkoJcrS2F@yV0CPU/XlXlt1XugehGWY1o3Fj7Czvp4kGZ30EToMbQQdgJ4TBBrxmY4fTBc/ReqkOAvIsCH7NA7yyqxZxxdYPupcKmiw6AfH8Sotyr94ANvDcecTMGnSYD46ryWxeow/Q@Aw0POoUAmkGQyF6qlMky6C9TWzgrgRZnOgDzZowB1kvt/kLVU@s@7YEkguvviNv4EgXMgiozx3UI8QOoUduYwe2KWhZKU4eY/L831BM@K4H1@NPgKvKhHCeOMrtSyEO2xv4u6K@BE/FfBaENEja7e5LgzQ3eLWixPQJ)
	2. Generate three random numbers \$x\$, \$y\$ and \$z\$ from a *uniform* distribution in the range \$(-1,1)\$. Calculate the length of the vector by \$l=\sqrt{x^2+y^2+z^2}\$. Then, if \$l>1\$, reject the vector and generate a new set of numbers. Else, if \$l \leq 1\$, normalize the vector and return the result.
		+ [Implementation example](https://tio.run/##jZLBTsMwDIbveQpLO6ydtkG3W0UvcEJCXJi4TunqrhFpEtIEGC9fkjYtWwGJHFrV8f/5t2t1MpUU27ZltZLagLC1OgFtQChCCixBU1HIev@GByN1FKcE3HmvGEfYaYsp/OfM4F6UTDCDwKVUZIjntMGADpHMFV73NbtXtI1Hxq7SiMBEgQrdQ5hgzpvOUTdQavdBwQpWSl1DwRqjWW4Nk8LJwFToFUeE6HoJSTzauHAQbGxgcWFvBUmw8VSx0nSwcGXkOXqV/MLecxRHU323yJmg/LgWzmfUp8SefUf5wXJqRnavG2msnABvMkjSv4Z@q5G@9FPxBv3oPSFIOTaNi1MBrg6@Wsp9Jwk5Z@Se0EW8U8rZJxbnPywbfF5NG53B46AIVx1Go7Fa/KSRFj9orfg472y6eYQozYSJ5s9dIJ0v4VISDwkPnQWf0Ng6YovFBtw@APM7MJHEpP0C)
	3. Generate two random numbers \$i\$ and \$j\$ from a *uniform* distribution in the range \$(0,1)\$ and convert them to spherical coordinates like so:\begin{align}\theta &= 2 \times \pi \times i\\\\\phi &= \cos^{-1}(2\times j -1)\end{align}so that \$x\$, \$y\$ and \$z\$ can be calculated by \begin{align}x &= \cos(\theta) \times \sin(\phi)\\\\y &= \sin(\theta) \times \sin(\phi)\\\\z &= \cos(\phi)\end{align}
		+ [Implementation example](https://tio.run/##lZJNT8MwDIbv/RWWOCydCoJym8SJKzcQF4SmrPVWT2sSOSm0@/MlH5voBkMih1ixXj@v3doMrtHqfhypNZodqK41A0gLymRZjWtgqWrdLj@wcppFvsjAn5BcUpHiFh68@iYJYxBlDpNzBS@fGkjVaNBfyh2gwWyFbGHN/iGhU7TW3EJN1jGtOkda@TJwDYaKDYK4LeAujy34pJPfHg9Qwjy0YcjH1N9JC8@mQaZK7qDSmmtS0mGCRJxppvI0keSq0laUR@AWrr37HzgPKfyb2X8srMEPczrLiqSNbj3AD7dgFfvJ0yCWlPDAPLg9SnZoSaqpWx9Rwy@oUPov1BBR@wtdpdLzcwG1j6i0LhPUW1/AUMD@Pe0Puo7VQZWN2MvW7PCwZF59tnRZZpiUE7PXmFjMCjgtyY@CJ1Qb1wSB7VpB83kZ/wKFNTorybPxCw)
* Provide in your answer a brief description of the algorithm that you are using.
* Read more on sphere point picking on [MathWorld](http://mathworld.wolfram.com/SpherePointPicking.html).
## Output examples
```
[ 0.72422852 -0.58643067  0.36275628]
[-0.79158628 -0.17595886  0.58517488]
[-0.16428481 -0.90804027  0.38532243]
[ 0.61238768  0.75123833 -0.24621596]
[-0.81111161 -0.46269121  0.35779156]
```
## General remarks
* This is [code-golf](/questions/tagged/code-golf "show questions tagged 'code-golf'"), so the answer using the fewest bytes in each language wins.
* [Standard rules](https://codegolf.meta.stackexchange.com/questions/2419/default-for-code-golf-program-function-or-snippet/), [I/O rules](https://codegolf.meta.stackexchange.com/questions/2447/default-for-code-golf-input-output-methods) and [loophole rules](https://codegolf.meta.stackexchange.com/questions/1061/loopholes-that-are-forbidden-by-default) apply.
* Please include a [Try it Online](https://tio.run/#)-link or equivalent to demonstrate your code working.
* Please motivate your answer with an explanation of your code.
  
[Answer]
# [Wolfram Language (Mathematica)](https://www.wolfram.com/wolframscript/), 20 bytes
```
RandomPoint@Sphere[]
```
[Try it online!](https://tio.run/##y00syUjNTSzJTE78H1CUmVcS/T8oMS8lPzcgH8hxCC7ISC1KjY79H/sfAA "Wolfram Language (Mathematica) – Try It Online")
Does exactly what it says on the tin.
[Answer]
# [R](https://www.r-project.org/), 23 bytes
```
x=rnorm(3)
x/(x%*%x)^.5
```
[Try it online!](https://tio.run/##K/r/v8K2KC@/KFfDWJOrQl@jQlVLtUIzTs/0/38A "R – Try It Online")
Generates 3 realizations of the \$\mathcal N(0,1)\$ distribution and normalizes the resulting vector.
Plot of 1000 realizations:
[![enter image description here](https://i.stack.imgur.com/emKP6.png)](https://i.stack.imgur.com/emKP6.png)
[Answer]
# x86-64 Machine Code - ~~63 62 55~~ 49 bytes
```
6A 4F                push        4Fh  
68 00 00 80 3F       push        3F800000h  
C4 E2 79 18 4C 24 05 vbroadcastss xmm1,dword ptr [rsp+5]  
rand:
0F C7 F0             rdrand      eax  
73 FB                jnc         rand  
66 0F 6E C0          movd        xmm0,eax  
greaterThanOne:
66 0F 38 DC C0       aesenc      xmm0,xmm0  
0F 5B C0             cvtdq2ps    xmm0,xmm0  
0F 5E C1             divps       xmm0,xmm1  
C4 E3 79 40 D0 7F    vdpps       xmm2,xmm0,xmm0,7Fh  
0F 2F 14 24          comiss      xmm2,dword ptr [rsp]  
75 E9                jne         greaterThanOne
58                   pop         rax  
58                   pop         rax  
C3                   ret  
```
Uses the second algorithm, modified. Returns vector of `[x, y, z, 0]` in xmm0.
Explanation:
```
push 4Fh
push 3f800000h
```
Pushes the value for 1 and 2^31 as a float to the stack. The data overlap due to the sign extension, saving a few bytes.
`vbroadcastss xmm1,dword ptr [rsp+5]` Loads the value for 2^31 into 4 positions of xmm1.
```
rdrand      eax  
jnc         rand  
movd        xmm0,eax
```
Generates random 32-bit integer and loads it to bottom of xmm0.
```
aesenc      xmm0,xmm0  
cvtdq2ps    xmm0,xmm0  
divps       xmm0,xmm1 
```
Generates a random 32 bit integer, convert it to float (signed) and divide by 2^31 to get numbers between -1 and 1.
`vdpps xmm2,xmm0,xmm0,7Fh` adds the squares of the lower 3 floats using a dot product by itself, masking out the top float. This gives the length
```
comiss      xmm2,dword ptr [rsp]  
jne          rand+9h (07FF7A1DE1C9Eh)
```
Compares the length squared with 1 and rejects the values if it is not equal to 1. If the length squared is one, then the length is also one. That means the vector is already normalised and saves a square root and divide.
```
pop         rax  
pop         rax 
```
Restore the stack.
`ret` returns value in xmm0
[Try it Online](https://tio.run/##ZZHbToQwEIavmaeYmHhl3QAewwN4643xRokpbXfFQFs7dcWnX3sAZbUJZfjmnzL/lBOpsRu@zndCHA6bXns1vNCX9nxCbaxT236CzW4w3YAj7zVs1OSV02hd0G6BvGtgw0n0PZ48nW4Z5qd9Lk8gFjRQCD4MuFP@UQlvHBT2g17RdRaK0exjwNDR/EV2YFhCQR9dhAyr65TglvApgJbhNI5lYkQxrlhKnF22K1jP8DZAsfdENclUuJSvYZVgdQzrBGsoBsXRyZ5hsJp7VHxieDH7mueQrcSGkzFrbA4m47I@bFA45QF@JtFAmsS8yunyDvEPuri7LeMKiX3nDJeCk198y0/jJFrvsterNqgc17IBJ@M7HxJ@HPibFsupmHMQOpYLSqPJyp1TPNzwwyvX91o1wBWppTjJ4hZ0YVjyvQ7X8ofLfp/hilexf2nXvGY/RaycbqJ1YcaeaKU4NthGF2rp@LhNiAP/9Rd9/Cdh@mE/HL4B).
[Answer]
# [Python 2](https://docs.python.org/2/), 86 bytes
```
from random import*;R=random
z=R()*2-1
a=(1-z*z)**.5*1j**(4*R())
print a.real,a.imag,z
```
[Try it online!](https://tio.run/##JclLCoAgFEbhuatoqD8lKDUKN9EOLvQy8sHFSW7egkYHvpOfcqZoW9s5hY4prl98yIkL5sX9IKpbpIIdjCAnzVBRFaAnmAuQI76pRGYfS0eaN7p70j7Q0dfWXg "Python 2 – Try It Online")
Generates the z-coordinate uniformly from -1 to 1. Then the x and y coordinates are sampled uniformly on a circle of radius `(1-z*z)**.5`.
It might not be obvious that the spherical distribution is in factor uniform over the z coordinate (and so over every coordinate). This is something special for dimension 3. See [this proof](http://mathworld.wolfram.com/Zone.html) that the surface area of a horizontal slice of a sphere is proportional to its height. Although slices near the equator have a bigger radius, slices near the pole are titled inward more, and it turns out these two effects exactly cancel.
To generate a random angle on this circle, we raise the imaginary unit `1j` to a uniformly random power between 0 and 4, which saves us from needing trig functions, pi, or e, any of which would need an import. We then extract the real imaginary part. If we can output a complex number for two of the coordinates, the last line could just be `print a,z`.
---
# [Python 2](https://docs.python.org/2/), 85 bytes
```
from random import*
a,b,c=eval('gauss(0,1),'*3)
R=(a*a+b*b+c*c)**.5
print a/R,b/R,c/R
```
[Try it online!](https://tio.run/##DcXJDYAgEADAP1XwE9aNeMSnTdDBgmeiQABNrB59TCa8efeuL2WN/uKR3Px3XMHHDIzQoJ2Wh05RbXSnJFrsJFYwSKYnQUC1AVNbsBKgGVmIh8uclEbzs0qX8gE "Python 2 – Try It Online")
Generates three normals and scales the result. Thanks to ovs for saving 1 byte with `eval` in the second line.
---
**[Python 2](https://docs.python.org/2/) with numpy, 57 bytes**
```
from numpy import*
a=random.randn(3)
print a/sum(a*a)**.5
```
[Try it online!](https://tio.run/##K6gsycjPM/r/P60oP1chrzS3oFIhM7cgv6hEiyvRtigxLyU/Vw9E5WkYa3IVFGXmlSgk6heX5mokaiVqamnpmf7/DwA "Python 2 – Try It Online")
`sum(a*a)**.5` is shorter than `linalg.norm(a)`. We could also do `dot(a,a)` for the same length as `sum(a*a)`. In Python 3, this can be shortened to `a@a` using the new operator `@`.
[Answer]
# [Octave](https://www.gnu.org/software/octave/), ~~40 33~~ 22 bytes
We sample form a 3d standard normal distribution and normalize the vector:
```
(x=randn(1,3))/norm(x)
```
[Try it online!](https://tio.run/##y08uSSxL/f9fo8K2KDEvJU/DUMdYU1M/L78oV6NC8/9/AA "Octave – Try It Online")
[Answer]
# [Unity C#](https://docs.unity3d.com/ScriptReference/Random-onUnitSphere.html), 34 bytes
```
f=>UnityEngine.Random.onUnitSphere
```
Unity has a builtin for unit sphere random values, so I thought I'd post it.
[Answer]
# [MATL](https://github.com/lmendo/MATL), 10 bytes
```
1&3Xrt2&|/
```
[Try it online!](https://tio.run/##y00syfn/31DNOKKoxEitRv//fwA "MATL – Try It Online")
### Explanation
This uses the first approach described in the challenge.
```
1&3Xr  % Generate a 1×3 vector of i.i.d standard Gaussian variables
t      % Duplicate
2&|    % Compute the 2-norm
/      % Divide, element-wise. Implicitly display
```
[Answer]
# [Ruby](https://www.ruby-lang.org/), ~~34 50~~ 49 bytes
```
->{[z=rand*2-1]+((1-z*z)**0.5*1i**(rand*4)).rect}
```
[Try it online!](https://tio.run/##KypNqvyfZvtf1646usq2KDEvRctI1zBWW0PDULdKq0pTS8tAz1TLMFNLSwMsaaKpqVeUmlxS@9/QQK8kMze1uLpAIS06tvY/AA "Ruby – Try It Online")
Returns an array of 3 numbers `[z,y,x]`.
`x` and `y` are generated by raising `i` (square root of -1) to a random power between 0 and 4. This complex number needs to be scaled appropriately according to the `z` value in accordance with Pythagoras theorem: `(x**2 + y**2) + z**2 = 1.` 
The `z` coordinate (which is generated first) is simply a uniformly distributed number between -1 and 1. Though not immediately obvious, dA/dz for a slice through a sphere is constant (and equal to the perimeter of a circle of the same radius as the whole sphere.) . 
This was apparently discovered by Archimedes who described it in a very non-calculus-like way, and it is known as Archimedes Hat-Box theorem. See <https://brilliant.org/wiki/surface-area-sphere/>
Another reference from comments on xnor's answer. A surprisingly short URL, describing a surprisingly simple formula: <http://mathworld.wolfram.com/Zone.html> 
[Answer]
# TI-BASIC, 15 bytes \*
```
:randNorm(0,1,3
:Ans/√(sum(Ans²
```
Using the algorithm "generate 3 normally distributed values and normalize that vector".
Ending a program with an expression automatically prints the result on the Homescreen after the program terminates, so the result is actually shown, not just generated and blackholed.
\*: `randNorm(` is a [two-byte token](http://tibasicdev.wikidot.com/miscellaneous-tokens), the rest are [one-byte tokens](http://tibasicdev.wikidot.com/one-byte-tokens). I've counted the initial (unavoidable) `:`, without that it would be 14 bytes. Saved as a program with a one-letter name, it takes 24 bytes of memory, which includes the 9 bytes of file-system overhead.
[Answer]
# [05AB1E](https://github.com/Adriandmen/05AB1E/wiki/Commands), ~~23~~ 22 [bytes](https://github.com/Adriandmen/05AB1E/wiki/Codepage)
```
[тε5°x<Ýs/<Ω}DnOtDî#}/
```
Implements the 2nd algorithm.
[Try it online](https://tio.run/##ATAAz/9vc2FiaWX//1vRgs61NcKweDzDnXMvPM6pfURuT3REw64jfS////8tLW5vLWxhenk) or [get a few more random outputs](https://tio.run/##ATMAzP9vc2FiaWX/NUb/W9GCzrU1wrB4PMOdcy88zql9RG5PdETDriN9L/8s//8tLW5vLWxhenk).
**Explanation:**
NOTE: 05AB1E doesn't have a builtin to get a random decimal value in the range \$[0,1)\$. Instead, I create a list in increments of \$0.00001\$, and pick random values from that list. This increment could be changed to \$0.000000001\$ by changing the `5` to `9` in the code (although it would become rather slow..).
```
[            # Start an infinite loop:
 тε          #  Push 100, and map (basically, create a list with 3 values):
   5°        #   Push 100,000 (10**5)
     x       #   Double it to 200,000 (without popping)
      <      #   Decrease it by 1 to 199,999
       Ý     #   Create a list in the range [0, 199,999]
        s/   #   Swap to get 100,000 again, and divide each value in the list by this
          <  #   And then decrease by 1 to change the range [0,2) to [-1,1)
           Ω #   And pop and push a random value from this list
  }          #  After the map, we have our three random values
   D         #   Duplicate this list
    n        #   Square each inner value
     O       #   Take the sum of these squares
      t      #   Take the square-root of that
       D     #   Duplicate that as well
        î    #   Ceil it, and if it's now exactly 1:
         #   #    Stop the infinite loop
}/           # After the infinite loop: normalize by dividing
             # (after which the result is output implicitly)
```
[Answer]
# JavaScript (ES7), ~~77 76~~ 75 bytes
Implements the 3rd algorithm, using \$\sin(\phi)=\sin(\cos^{-1}(z))=\sqrt{1-z^2}\$.
```
with(Math)f=_=>[z=2*(r=random)()-1,cos(t=2*PI*r(q=(1-z*z)**.5))*q,sin(t)*q]
```
[Try it online!](https://tio.run/##FYvLCsIwEADv/Yq9dTdtgxU82Xj3IHgXkdCHTalZmwSFiN8e62mGgZn0S/vWmWeoLHd9Sm8TRjzpMNKgbupwiWor0CmnbccPQqrqsmWPYc3no3C4KKyrKCIJIXdEYim9sRhWuaaBHVpQsNmDhQbqP4uC4JMBtGw9z72c@Y4Dkpx43fIScqLsm34 "JavaScript (Node.js) – Try It Online")
### Commented
```
with(Math)                       // use Math
f = _ =>                         //
  [ z = 2 * (r = random)() - 1,  // z = 2 * j - 1
    cos(                         //
      t =                        // θ =
        2 * PI *                 //   2 * π * i
        r(q = (1 - z * z) ** .5) // q = sin(ɸ) = sin(arccos(z)) = √(1 - z²)
                                 // NB: it is safe to compute q here because
                                 //     Math.random ignores its parameter(s)
    ) * q,                       // x = cos(θ) * sin(ɸ)
    sin(t) * q                   // y = sin(θ) * sin(ɸ)
  ]                              //
```
---
# JavaScript (ES6), 79 bytes
Implements the 2nd algorithm.
```
f=_=>(n=Math.hypot(...v=[0,0,0].map(_=>Math.random()*2-1)))>1?f():v.map(x=>x/n)
```
[Try it online!](https://tio.run/##HYtLDoJAEAX3nKJ3dAuM4FIcPIEnMMZM@AgEuglMCMZ49nEkb1GLV9Wb1Szl3E02Yalq5xr91AWyvhnbqvY9iUWl1Krvaez3UKOZ0Bv7PRuuZEQ6nJKMiIrs2iCd193ZdLEdmVwjMw61BQYNae5xgezPKCL4BACl8CJDrQZ5oa9VLx1jGENIFHzdDw "JavaScript (Node.js) – Try It Online")
### Commented
```
f = _ =>                         // f is a recursive function taking no parameter
  ( n = Math.hypot(...           // n is the Euclidean norm of
      v =                        // the vector v consisting of:
        [0, 0, 0].map(_ =>       //
          Math.random() * 2 - 1  //   3 uniform random values in [-1, 1]
        )                        //
  )) > 1 ?                       // if n is greater than 1:
    f()                          //   try again until it's not
  :                              // else:
    v.map(x => x / n)            //   return the normalized vector
```
[Answer]
# [Processing](https://processing.org/) 26 bytes
Full program
```
print(PVector.random3D());
```
This is the implementation <https://github.com/processing/processing/blob/master/core/src/processing/core/PVector.java>
```
  static public PVector random3D(PVector target, PApplet parent) {
    float angle;
    float vz;
    if (parent == null) {
      angle = (float) (Math.random()*Math.PI*2);
      vz    = (float) (Math.random()*2-1);
    } else {
      angle = parent.random(PConstants.TWO_PI);
      vz    = parent.random(-1,1);
    }
    float vx = (float) (Math.sqrt(1-vz*vz)*Math.cos(angle));
    float vy = (float) (Math.sqrt(1-vz*vz)*Math.sin(angle));
    if (target == null) {
      target = new PVector(vx, vy, vz);
      //target.normalize(); // Should be unnecessary
    } else {
      target.set(vx,vy,vz);
    }
    return target;
  }
```
[Answer]
# [Python 2](https://docs.python.org/2/), 86 bytes
```
from random import*
x,y,z=map(gauss,[0]*3,[1]*3);l=(x*x+y*y+z*z)**.5
print x/l,y/l,z/l
```
[Try it online!](https://tio.run/##DcVLCoAgFEDReatoqM9H9qFRtJJoIEQf8IcZ@Ny8Objneoq3s2MpZ3CmDcoedY/xLkRoEhLm1SjPLvW9L279DhNuQ5UvemUJkiAgkSFzgG5ufHhsbJPUSLUsdSk/ "Python 2 – Try It Online")
Implements the first algorithm.
---
# [Python 2](https://docs.python.org/2/), ~~107~~ 103 bytes
```
from random import*
l=2
while l>1:x,y,z=map(uniform,[-1]*3,[1]*3);l=(x*x+y*y+z*z)**.5
print x/l,y/l,z/l
```
[Try it online!](https://tio.run/##FcVNCoMwEEDhfU6RpY7TipZuKulFpAvBioHJD0OkmVw@1cX7XpS0Bz/WunFwmhe/nrMuBk6gyIzqt1v6anoPr4yCxbglNoe3W2CH8234wAPny3Yi02TInYB0BUoLcH@qyNYnnXtCOSs91foH "Python 2 – Try It Online")
Implements the second algorithm.
[Answer]
# [Haskell](https://www.haskell.org/), ~~125~~ ~~123~~ ~~119~~ 118 bytes
```
import System.Random
f=mapM(\_->randomRIO(-1,1))"lol">>= \a->last$f:[pure$(/n)<$>a|n<-[sqrt.sum$map(^2)a::Double],n<1]
```
[Try it online!](https://tio.run/##Hcy7CsIwFADQ3a8IJUMCTaWOJcnk4iBCHdsqV0xoMS@TdBD89/hYz3BmSA9lTCmLDT5mdH6lrGzTg7t7u9HCQjiS8cpk/Et/OBHW1i2llfGmklKgEZg0kDLW3RDWqDDZOsqxhLfjbEjPmJu0WvyNyGVHoev2fr0ZNdWOt1OxsDgkkEa/KsTFZVQ@ "Haskell – Try It Online")
Does three uniforms randoms and rejection sampling.
[Answer]
# JavaScript, 95 bytes
```
f=(a=[x,y,z]=[0,0,0].map(e=>Math.random()*2-1))=>(s=Math.sqrt(x*x+y*y+z*z))>1?f():a.map(e=>e/s)
```
You ~~don't~~ need *not* to input `a`.
[Answer]
# [Julia 1.0](http://julialang.org/), 24 bytes
```
x=randn(3)
x/hypot(x...)
```
[Try it online!](https://tio.run/##yyrNyUw0rPhfUJSZV5KTp5GUmp6Z97/CtigxLyVPw1iTq0I/o7Igv0SjQk9PT/N/al6K5n8A "Julia 1.0 – Try It Online")
Draws a vector of 3 values, drawn from a normal distribution around 0 with standard deviation 1. Then just normalizes them.
[Answer]
# [MathGolf](https://github.com/maxbergmark/mathgolf/blob/master/math_golf.txt), ~~21~~ ~~19~~ 18 [bytes](https://github.com/maxbergmark/mathgolf/blob/master/code_page.py)
```
{╘3Éƒ∞(ß_²Σ√_1>}▲/
```
Implementation of the 2nd algorithm.
[Try it online](https://tio.run/##ASsA1P9tYXRoZ29sZv//e@KVmDPDicaS4oieKMOfX8KyzqPiiJpfMT594payL///) or [see a few more outputs at the same time](https://tio.run/##AS4A0f9tYXRoZ29sZv9Fe/974pWYM8OJxpLiiJ4ow59fwrLOo@KIml8xPn3ilrIv/3D/).
**Explanation:**
```
{              }▲   # Do-while true by popping the value:
 ╘                  #  Discard everything on the stack to clean up previous iterations
  3É                #  Loop 3 times, executing the following three operations:
    ƒ               #   Push a random value in the range [0,1]
     ∞              #   Double it to make the range [0,2]
      (             #   Decrease it by 1 to make the range [-1,1]
       ß            #  Wrap these three values into a list
        _           #  Duplicate the list of random values
         ²          #  Square each value in the list
          Σ         #  Sum them
           √        #  And take the square-root of that
            _       #  Duplicate it as well
             1>     #  And check if it's larger than 1
                 /  # After the do-while, divide to normalize
                    # (after which the entire stack joined together is output implicitly,
                    #  which is why we need the `╘` to cleanup after every iteration)
```
[Answer]
# Java 8 (*@Arnauld*'s modified 3rd algorithm), ~~131~~ ~~126~~ ~~119~~ ~~111~~ 109 bytes
```
v->{double k=2*M.random()-1,t=M.sqrt(1-k*k),r[]={k,M.cos(k=2*M.PI*M.random())*t,M.sin(k)*t};return r;}
```
Port of [*@Arnauld*'s JavaScript answer](https://codegolf.stackexchange.com/a/191514/52210), so make sure to upvote him!  
-2 bytes thanks to *@OlivierGrégoire*.
This is implemented as:
\$k = N\cap[-1,1)\$  
\$t=\sqrt{1-k^2}\$  
\$u=2\pi×(N\cap[0,1))\$  
\$x,y,z = \{k, \cos(u)×t, \sin(u)×t\}\$
[Try it online.](https://tio.run/##TVA9b4MwEN3zK260aeyULh0QkTp2IKoUqUuUwTFOYz7s1D6oEOK30yN06HLf7@7eq1SvRFXWs25UjFAo68YNgHVowlVpA4clBSh9d2nM6QyafXpbQs8zqk8bMoXCGxQZwG4HT6@AHvBm4DKgEdp3DmnkAA5ymHuxH9dFUOcvSSGDcqVvGRfpFvNCxu@ALBV1UvNtOJ3zsd4WUvvI1umP938QniA1o3WspnDKgsEuOAjZNGfLV3c6YzVEVEiuX35uiRw7YrDui4govjK7@sCILtg8fc7ACrEnxx8tgOMQ0bTSdyjvhMPGsYoUkx3aRr6FoIYo0a87mTM/DwEZl05q5rqm4fxPp2n@BQ)
**Previous 3rd algorithm implementation (~~131~~ ~~126~~ 119 bytes):**
```
Math M;v->{double k=2*M.random()-1,t=2*M.PI*M.random();return k+","+M.cos(t)*M.sin(k=M.acos(k))+","+M.sin(t)*M.sin(k);}
```
Implemented as:
\$k = N\cap[-1,1)\$  
\$t=2\pi×(N\cap[0,1))\$  
\$x,y,z = \{k, \cos(t)×\sin(\arccos(k)), \sin(t)×\sin(\arccos(k))\}\$
[Try it online.](https://tio.run/##TZCxbsMgEIb3PMUpE9SBNF06IGfvQFSpUpeqA8EkIbaPCM6uosjP7mInUrtw4vsPcd@dTW/EuapH25iUQBuPtwWAR3LxYKyD3XQF@KDo8QiWfQZfQc9VpsMiH9rQCbQCWK@heAUKQCcH@ys5YUOHlFt2gFDC2IvtrQrdvnFQly9PWkaDVWgZF5sVzeD97R9V0VEXEepiuVoWWtqQGPHckDyyutTSTKTm/JFP@C/nahjVNN8lf@gtJDKUSz9N32ZJdhf6@gbD74aHEFnWBl9unhV4Iba58DnK@tdErpWhI3nJ76hBhu5nXhfjEqVl2DUNf6xlGH8B)
### Explanation:
```
Math M;                         // Math on class-level to use for static calls to save bytes
v->{                            // Method with empty unused parameter & double-array return
  double k=2*M.random()-1,      //  Get a random value in the range [-1,1)
         t=M.sqrt(1-k*k),       //  Calculate the square-root of 1-k^2
    r[]={                       //  Create the result-array, containing:
         k,                     //   X: the random value `k`
         M.cos(k=2*M.PI         //   Y: first change `k` to TAU (2*PI)
                     *M.random()//       multiplied by a random [0,1) value
                )               //      Take the cosine of that
                 *t,            //      and multiply it by `t`
         M.sin(k)               //   Z: Also take the sine of the new `k` (TAU * random)
                  *t};          //      And multiply it by `t` as well
  return r;}                    //  Return this array as result
```
---
# Java 8 (2nd algorithm), ~~153~~ 143 bytes
```
v->{double x=2,y=2,z=2,l;for(;(l=Math.sqrt(x*x+y*y+z*z))>1;y=m(),z=m())x=m();return x/l+","+y/l+","+z/l;};double m(){return Math.random()*2-1;}
```
[Try it online.](https://tio.run/##LVDBcoMgEL3nK3ZyAo2m5srgHySXzOTS6YEqaUkRUkArOn67XY2HfW9Y3sJ7@xCdyB71z1xp4T2chTLjDkCZIN1dVBIuyxHgGpwyX1CRm1U1dJRhd9ohXMAAh7nLyrG27aeW0PPTIWINWJrdrSOMaH4W4Tv3vy6QPunTmMR0SAZKy4JF3hCKakTaL8icDK0z0B91uj/s07jxcNRsYtsvqBs33fq0E6a22ExOWcGmGZ09Uacq8EEEpG7x3WA88ory/gGCvrItFjEwKF68MVBZViLR9QqDRx9kk9s25E@cC9oQI//WRRGam7wiptWabguZ5n8)
*2nd algorithm:*
```
v->{                              // Method with empty unused parameter & String return-type
  double x=2,y=2,z=2,l;           //  Start results a,b,c all at 2
  for(;(l=Math.sqrt(x*x+y*y+z*z)) //  Loop as long as the hypotenuse of x,y,z
       >1;                        //  is larger than 1
    y=m(),z=m())x=m();            //   Calculate a new x, y, and z
  return x/l+","+y/l+","+z/l;}    //  And return the normalized x,y,z as result
double m(){                       // Separated method to reduce bytes, which will:
  return Math.random()*2-1;}      //  Return a random value in the range [-1,1)
```
[Answer]
# [Vyxal](https://github.com/Lyxal/Vyxal), ~~20~~ 19 bytes
```
{≥|3ʁƛ∆Ṙd⌐;:²∑√:1}/
```
[Try it Online!](http://lyxal.pythonanywhere.com?flags=&code=%7B%E2%89%A5%7C3%CA%81%C6%9B%E2%88%86%E1%B9%98d%E2%8C%90%3B%3A%C2%B2%E2%88%91%E2%88%9A%3A1%7D%2F&inputs=&header=&footer=)
Uses the second method.
## Explanation
```
{≥|3ʁƛ∆Ṙd⌐;:²∑√:1}/
{                }  While
 ≥                    Greater or equal
  |                 do:
   3ʁ                 [0,1,2]
     ƛ    ;           Map:
      ∆Ṙ                Random number from 0 to 1
        d               Double
         ⌐              Subtract from 1
           :          Duplicate
            ²         Square
             ∑        Sum
              √       Square root
               :      Duplicate
                1     Push 1
                  / Divide [the vector by the length]
```
[Answer]
# [Japt](https://github.com/ETHproductions/japt), 20 [bytes](https://en.wikipedia.org/wiki/ISO/IEC_8859-1)
Port of [Arnauld's implementation](https://codegolf.stackexchange.com/a/191514/58974) of the 2nd algorithm.
```
MhV=3ÆMrJ1
>1?ß:V®/U
```
[Test it](https://petershaggynoble.github.io/Japt-Interpreter/?v=1.4.6&code=TWhWPTPGTXJKMQo%2bMT/fOlauL1U)
```
MhV=3ÆMrJ1
Mh             :Get the hypotenuse of
  V=           :  Assign to V
    3Æ         :  Map the range [0,3)
      Mr       :    Random float
        J1     :    In range [-1,1)
>1?ß:V®/U      :Assign result to U
>1?            :If U is greater than 1
   ß           :  Run the programme again
    :V®/U      :Else map V, dividing all elements by U
```
[Answer]
# [Pyth](https://github.com/isaacg1/pyth), 24 bytes
```
W<1Ks^R2JmtO2.0 3;cR@K2J
```
[Try it online!](https://tio.run/##K6gsyfj/P9zG0Ls4LsjIK7fE30jPQMHYOjnIwdvI6/9/AA "Pyth – Try It Online")
Uses algorithm #2
```
W                         # while 
 <1                       #   1 < 
   Ks                     #       K := sum(
     ^R2                  #               map(lambda x:x**2,
        Jm      3         #                    J := map(                            , range(3))
          tO2.0           #                             lambda x: random(0, 2.0) - 1           )):
                 ;        #   pass
                   R   J  # [return] map(lambda x:            , J)
                  c @K2   #                        x / sqrt(K)
```
[Answer]
# [OCaml](https://ocaml.org), ~~110~~ ~~99~~ 95 bytes
```
(fun f a c s->let t,p=f 4.*.a 0.,a(f 2.-.1.)in[c t*.s p;s t*.s p;c p])Random.float acos cos sin
```
EDIT: Shaved off some bytes by inlining \$ i \$ and \$ j \$, replacing the first `let ... in` with a `fun`, and taking advantage of operator associativity to avoid some parens `()`.
[**Try it online**](https://tio.run/##NY7BasMwEETv@YrBJ8kkS1p6C@0X9JQeSzCLvEoF8kpIm/6@m5iUYWCGB8OUwEte3YjhfNMBVnAVlcYmmFOM0kQNjXUuC34lWGkdwuEHlhbB6He78wapS45T0mRwHifUltQm0TknFbgvu/crhaKBDcNd7jN1o4UrXN/gVOIUc2HzcKuLN0UEI6AfPrIYbF/fI95oJMaR9uwiXulAL@STfgfYSB311P9DQL3457VtFRxKx8M96eq9X/8A)
---
Original solution:
```
Random.(let a,c,s,i,j=acos,cos,sin,float 4.,float 2. in let t,p=i*.(a 0.),a (j-.1.) in[c t*.s p;s t*.s p;c p])
```
First I define:
$$
a = \arccos,\ \ c = \cos,\ \ s = \sin \\
i \sim \textsf{unif}(0,4),\ \ j \sim \textsf{unif}(0,2)
$$
OCaml's [`Random.float` function](https://caml.inria.fr/pub/docs/manual-ocaml/libref/Random.html#VALfloat) includes the bounds. Then,
$$
t = i \cdot a(0) = \frac{i\pi}{2},\ \ p = a (j-1)
$$
This is very similar to the 3rd example implementation (with \$ \phi = p \$ and \$ \theta = t \$) \$ - \$ except that I pick \$ i \$ and \$ j \$ within larger intervals to avoid multiplication (with 2) later on.
[Answer]
# [Nim](http://nim-lang.org/), 105 101 100 99 bytes
```
import math,random
randomize()
let
 t=TAU.rand
 u=2.0.rand-1
 v=sqrt 1-u*u
echo [t.cos*v,t.sin*v,u]
```
[Try it online!](https://tio.run/##y8vM/f8/M7cgv6hEITexJEOnKDEvJT@XC0JlVqVqaHLlpJZwKZTYhjiG6oGEuRRKbY30DMBsXUMuhTLb4kKgbkPdUq1SrtTkjHyF6BK95PxirTKdEr3izDwgXRr7/z8A "Nim – Try It Online")
Uses the third method with a bit of substitution:
$$u = \cos(\phi) = 2j - 1$$
$$v = \sin(\phi) = \sqrt{1-u^2}$$
[Answer]
# [Raku](http://raku.org/), 65 bytes
```
{(&cos,&sin,{1})».(τ.rand)Z*(&sin,&sin,&cos)».(acos 2.rand-1)}
```
[Try it online!](https://tio.run/##K0gtyjH7n1upoJamYPu/WkMtOb9YR604M0@n2rBW89BuPY3zLXpFiXkpmlFaGmBxCAFUBpZNBDIUjMAqdA01a/@n6RUnViqk5RcpxBka/AcA "Perl 6 – Try It Online")
This function uses the "convert two random numbers to spherical coordinates" method. `τ.rand` generates the first number in the range of zero to τ (tau, that is, twice pi) and `acos 2.rand - 1` generates the second with the acos argument in the range of -1 to 1.
`(&cos, &sin, {1})».(...)` applies the cosine function, the sine function, and an anonymous function that always returns `1` to the first random parameter, producing a list of the results. Similarly, `(&sin, &sin, &cos)».(...)` applies the sine, sine, and cosine functions to the second random parameter, producing a second list. Finally, `Z*` zips those two lists together using multiplication.
[Answer]
# [Pxem](https://esolangs.org/wiki/Pxem), 58 bytes (filename).
Because Pxem does not handle decimal values nor negative values (negative values can be supplied via input but cannot be generated in other ways), I am representing them as integers.
Unprintables are represented as a backslash followed by three-digit octal value of its code point.
```
x\003\144.w\002.r\002.!+.+.o.v21.y.c\013.r.c.t.c.!.a.m.n.m.m.!.-.v\001.-.c.a
```
## Output format
It outputs 10x, 10y, and 10z rather than x, y, and z.
```
"%c%d%c%d%c%d\n", sign, value, sign value, sign, value
```
Sign is as either `+` or `-`; always output. Value is multiplied by 10. Thus the precision of the program is up to 0.1.
## With comments
```
XX.z
# memory layout: counter, l
# until counter is zero; do
.ax\003\144.wXX.z
  # output sign
  .a\002.r\002.!+.+.oXX.z
  # while :; do
  .a.v21.yXX.z
    # generate a random number 0<=x<=10
    .a.c\013.rXX.z
    # continue if x^2>l; otherwise break
    .a.c.t.c.!XX.z
  # done
  .a.aXX.z
  # output x
  .a.m.nXX.z
  # update l to l minus x^2
  .a.m.m.!.-XX.z
# decrement counter; done
.a.v\001.-.c.a
```
[Try it online!](https://tio.run/##fVdrV@JIE/7ur2iykXQTcumIisToMOjrsKuMOq7H0aAnQFAwJGwSFA3sX5@3upMImZ0znoPpp6q67n17jp5@/IEM3aBI30G0rqK93a1th@7s7BgbfoAJLiG3/xQgMaFK9LQ0Am9wUDaIyVjBAClNHykv8Ywu4hAprQjZd7qy10XSne1XutLiMXSnSCVm5MZIcWfmbOJEz0jXDcN059MgjNFp66F5emq1TBy/TV306Mb9wB@Wyxz1g8nE8QfkQBu4L5o/8zxkHJRpuZxNPm9efbEEEWdySJnmCjiLiAn7yI1lOlCWQm7270775uHb1ZFl6PpWTjz/@q19c/r9ofX18vK4dWVR8w4JYvK5@e3Lw/Xx5bf2104DdKBuucwDCtA0iEbzFMmfU@nb/wq7k5nnxC6KnrhoMI1h@BqEg2jqjeJy2fFGToSURyRBmmVB/CQsJUtiX6lc7vx9eto6O7IaXD1mpgZm9DQaxjyrA8TkMkJaqj@IgCxEmWU/MN0Xx0NCCQ0dL3KFxQJIXBGXnS8id4CkSJtrb5oZwb93TeLT31GXySLgm30oNheXhkGIPafnepYARk2SgJc4WGCuRbuv7Yha0miaYs/smJrtc3wOeu9V1fY1zew1zWUKVUoNUZtqZpzhul6AlFJRC1Z4d68AKQVxf41dL0Cqb4vaaI29VYB7wH1Yg6C7v4I7uwVIKUhHa45uFSClYPpljb1TgFSHMIdrnugFSHUI012zvVeAlNZELVxTbhQgpSD@usbeLUBqgLH5GrtegNQA229r7L0CpAYYe1/Lkl6EYMtZc3y7ACmFNMRrumsFSHWwNVnLym4BUh1sDdaU1wuwBhWQ1yCYVlZwC7ilNQiOiGsQLG0yuJTS5TIg6RrQ@M6lJSese@fmABr1hM2xsWoTG9s@@69WbCJq9pZt2FQDqUhTK9DW83ynS1U1PFMsJdrd6/zt3Rl0tVLSMXseKFyab7AqNKSZ2t19xmXNJGpooBV0KD6SMvN395@qXbALHuC7dAguGJ8Qc@Epk0npSK0Ai3uWWod4ptxNUasyag9cgDIlgEpJBKOmNjV7F0BNxaqfQMGH4OBDcKAd5YIFsZyT5aERmmKSmvuEmGrYC7g0auaUuMlTK5Y6wLH9O6S8dysDpi0O1ygQEs9PNw1W4ywHxKa8kmlwTjprCd42Lj5Uqty3@KJYkzwfXyAzR03IyzwrLs/sKm9pfk4yTfyXupJ6klX9VxOlfC9dq4c2Gmb7pSXYVCAJDwTfTQN/9NCPXoZuGE/kkripdJm0TTExtcdUhrWjndCqvWSsc2zTjOesexTNelhg9RCqbPYhjIUGbM5VbpYsM4NF@d8JH4laL3SdZw6gbE0seEEwRX4Qo4kT959G/qMAoh/hgnaZuY@J1gMnMSWEu1niBD0FIoB/crDJAQUgkaVANn78r9M8O7akuQ2HsU1rNfUVRoYa8v8lWZXVQH0xqPqm9m2dbqmh2ldj@JVUR52oPvwmMFbUF5hA4dtXHWmj9bVzddy5sjY2@k7kwpHHrQho5CNc0SqLim3r27sVkl5wJI0xuIgENxzTdCOnvzENR348/Ji7eO0jpQ/sH0NrOPP78SjwbWRGYd@SPh@ftDtJFMI9BI5Dd0CWcDYOBekcH5Hk291z9S9MutaRCcOuLGfMY0yS0I1noY9KMghk5D9XZCAe0Iz814rMtGTUqwKVmVFozrvBlyS5tECE21UUM5O8zARmuEWSNErBhr/NQKi2Mt4JvqzefuiGi5U38t0D/VDUG0ru0pfUgEhN3lj3vGtl6C9BIAUSSkk/mf@bmZ/5Vm6yic8yd9DZQug7MeSa3dlGMcotTiHc16eR5@ISZI8EH0mDUQKUxWKGb4CRUX3crl5Xj6oeSWAt8inJtQUCpmd5rv8YP@FrYrJ7jdne90wyw/jIAsejOMTXVVluV6Gf/xV4EMLhkVyrN2rbZJlpH4HNc3yyMveAz6qt6m31O0nOLEF9aKAOrBsHDdz@aAIXMYjMfXRD@E5nsWC2LJhr3sJNk3tgcidru/utcrm1v10nfCWaQGv9K9zX7ra2uyJsIrdWraa0ssnpYgURXE@n0RrkoGVZWwZJUhG4EMcjf@YuWXqZaGu/VgeZg@1dwkip8ZXZVDP5buHvBLWUWt1kdTLP8W3lex5oP8/2Ob5ahR/l1JuPsrxAAcbVa5K0LZ0bGlu8RSHdY9aU1bGidK1rVhWG2l1ObMtyl@NxN891uLIIWcR8oZEKK3Vu6nW1En5yYo4vqtxECZYWy2cJ0AVvg2vrYp8JZhOvsxlvv51x8IsZ7x/W2RRuXVl5EOfeP8E6tFb0CdCf2dYDb4NzzJg5p4cdSFzat3@ytj3HziFMlB0WdOMmC95cLfqKRRWjgp19nVRYVeCbJ@@fn5WNszVwk0WHxxWPQDeADUHdhFNBhU4TGujWDeGONHoZRbDXCawJWPKdw/GBdzje9Bre5rjBxxofQy2W6U4Ibll6F2rOngg/76N9IkiZZ0OcviI@RLJ9W1gE/ISx0nMmPdvh6gTnDnurKaH/TBf9GYwGEpKQMjTyAwmuIue4TNhDJjfi8iZMT7kEGo434rMsm7I83sdsa2TNBilM0znOIAyI@WvnWAQD13PhWQdFgxbmKtsWNDZsGvswXaZpK/MNv801yhT6O92K5RTDRm0GC@f1GZ7RiB0b7DkLV6A@4k/eGevyYMKeckAM/0NcPcEHNaR0agWBxQJP4VnZhGfqaFB149HErU7701n1JXo3B/AgJYu1J3ztt0/4ReyMPLiV1tayzG7GSZ1dT/gdSFpg0GDbPpLTm7Vi0Npurb61U6uTRQ8SBvmDo1LIw4Xt1xIxRiLOH@3Ny5OHs@bNYsHn1/S9HaIZCjUQIUjKduue1UOiTg4AZidFz@xZgrA87hzlhKWEFnMnfITo23OUFW/@4/8 "ksh – Try It Online")
[Answer]
# [APL (Dyalog Classic)](https://www.dyalog.com/), 22 bytes
```
⎕←{⍵÷((+/⍵*2)*0.5)}?⍳3
```
[Try it online!](https://tio.run/##SyzI0U2pTMzJT9dNzkksLs5M/v//Ud/UR20Tqh/1bj28XUNDWx/I0DLS1DLQM9WstX/Uu9n4/38A "APL (Dyalog Classic) – Try It Online")
]