text:
I see your BIDMAS and raise you a BADMIS
I see your BIDMAS and raise you a BADMIS
Challenge
Given a set of numbers with operators between them: "5 + 4 * 9 / 3 - 8", return all the possible results of the expression for every permutation of the order of basic operations: [/, *, +, -].
Rules

Standard loopholes forbidden
I/O

Input must be ordered with infix operations, but however that is easiest (string or array)
You are not required to support unary operators (e.g. "-3 * 8 / +2")
Integers can be replaced by floats for languages that implicitly parse type (e.g. 45 ⟶ 45.0)
Output must be all of the possible results of the expression, no specified format or order

All of the inputs are valid (e.g. do not need to deal with "7 / 3 + *"). This also means that you will never need to divide by zero.
Operators are all left-associative so "20 / 4 / 2" = "(20 / 4) / 2"
This is Code Golf so fewest number of bytes wins

Test Cases (With explanation)

"2 + 3 * 4" = [14, 20]

2 + (3 * 4)  ⟶  2 + (12) ⟶ 14
(2 + 3) * 4 ⟶ (5) * 4 ⟶ 20

"18 / 3 * 2 - 1" = [11, 2, 6]

((18 / 3) * 2) - 1 ⟶ ((6) * 2) - 1 ⟶ (12) - 1 ⟶ 11
(18 / 3) * (2 - 1) ⟶ (6) * (1) ⟶ 6
(18 / (3 * 2)) - 1 ⟶ (18 / (6)) - 1 ⟶ (3) - 1 ⟶ 2
18 / (3 * (2 - 1)) ⟶ 18 / (3 * (1)) ⟶ 6
18 / ((3 * 2) - 1) ⟶ 18 / 5 ⟶ 3.6

Test Cases (Without explanation)

"45 / 8 + 19 / 45 * 3" = [6.891666666666667, 18.141666666666666, 0.11111111111111113, 0.01234567901234568, 0.01234567901234568, 5.765740740740741]
"2 + 6 * 7 * 2 + 6 / 4" = [112 196 23 87.5]

Nice first challenge, though, by the way.

Closely related

Suggested test case 2 - 3 + 4 => [-5, 3]

Suggested test case: 2*3-6+2-9/6*8+5/2-9, giving 24 distinct results.

JavaScript (V8),  118  112 bytes
Prints the results.
f=(s,x='',y=x,o='+-*/')=>[...o].map(v=>f(s.split(v).join(x+v+y),x+')',y+'(',o.replace(v,'')))|print(eval(y+s+x))

Try it online!
Or see the deduplicated results.

C# (Visual C# Interactive Compiler), 285 bytes
x=>{int c=0,j,t=1,i;for(;c++<25;t=c){var r="*+-/".ToList();for(i=j=1;j++<4;t=t/j+1)(r[j-1],r[t%j])=(r[t%j],r[j-1]);float k(float z,int p=4){char d;int l;float m;return i<x.Count&&(l=r.IndexOf(d=x[i][0]))<p?k((m=k(x[(i+=2)-1],l))*0+d<43?z*m:d<44?z+m:d<46?z-m:z/m,p):z;}Print(k(x[0]));}}

Try it online!
x=>{                                          //Lambda taking in a List<dynamic>
  int c=0,j,t=1,i;                            //A bunch of delcarations jammed together to save bytes
  for(;c++<25;t=c){                           //Loop 24 times (amount of permutations a set of length 4 can have)
    var r="/+*-".ToList();                    //Initialize r as list of operators
    for(i=j=1;j++<4;t=t/j+1)                    //Create the Tth permutation, saving result in r, also reset i to 1
      (r[j-1],r[t%j])=(r[t%j],r[j-1]);
    float k(float z,int p=4) {                //Define local function 'k', with z as current value accumalated and p as current precedence
      char d;int l;float m;                   //Some helper variables
      return i<x.Count                        //If this is not the last number
        &&(l=r.IndexOf(d=x[i][0]))<p?         //  And the current operator's precedence is higher than the current precedence
      k(                                      //  Recursive call with the accumalative value as
        (m=k(x[(i+=2)-1],l))                  //    Another recursive call with the next number following the current operator as seed value,
                                              //    And the next operator's precedence as the precedence value, and store that in variable 'm'
        *0+d<43?z*m:d<44?z+m:d<46?z-m:z/m,    //    And doing the appropriate operation to m and current value ('z')
        p)                                    //  Passing in the current precedence
    :z;                                       //Else just return the current number
    }
    Print(k(x[0]));                           //Print the result of calling k with the first number as starting value
  }
}

I have fixed it so you do not need to omit duplicates as it is not a fundamental part of the problem as pointed out.

@Arnauld Fixed at the cost of 4 bytes, it was because my permutations algorithm was a bit wrong

Suggest for(i=j=1;j<4;t=t/j+1)(r[j],r[t%++j])=(r[t%j],r[j-1]); instead of for(i=j=1;j++<4;t=t/j+1)(r[j-1],r[t%j])=(r[t%j],r[j-1]);

JavaScript (Node.js), 132 bytes
a=>(w=[],F=(b,a)=>b?[...b].map(q=>F(b.replace(q,""),a.replace(eval(`/[\\d.-]+( \\${q} [\\d.-]+)+/g`),eval))):w.push(a))("+-*/",a)&&w

Try it online!
This allows duplicated outputs.
JavaScript (Node.js), 165 161 155 153 152  137 bytes
a=>Object.keys((F=(b,a)=>b?[...b].map(q=>F(b.replace(q,""),a.replace(eval(`/[\\d.-]+( \\${q} [\\d.-]+)+/g`),eval))):F[a]=1)("+-*/",a)&&F)

Try it online!
Takes a string with spaces between operators and numbers.
a=>                                             // Main function:
 Object.keys(                                   //  Return the keys of the -
  (
   F=(                                          //   Index container (helper function):
    b,                                          //    Operators
    a                                           //    The expression
   )=>
    b                                           //    If there are operators left:
    ?[...b].map(                                //     For each operator:
     q=>
      F(                                        //      Recur the helper function -
       b.replace(q,""),                         //       With the operator deleted
       a.replace(                               //       And all -
        eval(`/[\\d.-]+( \\${q} [\\d.-]+)+/g`), //        Expressions using the operator
        eval                                    //        Replaced with the evaluated result
       )
      )
    )
    :F[a]=1                                     //     Otherwise - set the result flag.
  )(
   "+-*/",                                      //    Starting with the four operators
   a                                            //    And the expression
  )
  &&F
 )

@JoKing Implemented the fix I've stated before, should output [3, -5] now.

Perl 6, 92 90 88 bytes
{map {[o](@_)($_)},<* / + ->>>.&{$^a;&{S:g{[\-?<[\d.]>+]+%"$a "}=$/.EVAL}}.permutations}

Try it online!
Takes a string with a space after any operators and returns a set of numbers. This mostly works by substituting all instances of n op n with the eval'ed result for all permutations of the operators.
Explanation:
{                                                                                }  # Anonymous code block
                    <* / + ->>>.&{                                    } # Map the operators to:
                                  $^a;&{                             }  # Functions that:
                                        S:g{                }      # Substitute all matches of:
                                            \-?<[\d.]>+]+        # Numbers
                                                         %$a     # Joined by the operator
                                                              =$/.EVAL   # With the match EVAL'd
 map {           },                                                    .permutations   # Map each of the permutations of these operators
      [o](@_)        # Join the functions
             ($_)    # And apply it to the input

You can remove the set, as the condition to eliminate duplicates has been removed. Nice code.

Python 3, 108 bytes
f=lambda e,s={*"+-*/"}:[str(eval(p.join(g)))for p in s for g in zip(*map(f,e.split(p),[s-{p}]*len(e)))]or[e]

Try it online!
The function takes a single string as input and returns a list of possible results.
Ungolfed
def get_all_eval_results(expr, operators={*"+-*/"}):
    results = []
    for operator in operators:
        remaining_operators = operators - {operator}

        # Split expression with the current operator and recursively evaluate each subexpression with remaining operators
        sub_expr_results = (get_all_eval_results(sub_expr, remaining_operators) for sub_expr in expr.split(operator))

        for result_group in zip(*sub_expr_results):   # Iterate over each group of subexpression evaluation outcomes
            expr_to_eval = operator.join(result_group)  # Join subexpression outcomes with current operator
            results.append(str(eval(expr_to_eval)))   # Evaluate and append outcome to result list of expr
    return results or [expr]  # If results is empty (no operators), return [expr]

Try it online!

Jelly, 30 bytes
œṡ⁹¹jṪḢƭ€jŒVɗßʋFL’$?
Ḋm2QŒ!çƒ€

Try it online!
A pair of links. The second one is the main link, and takes as its argument a Jelly list of floats/integers interspersed with the operators as characters. This is a flattened version of the way Jelly takes its input when run as a full program with command line arguments. The return value of the link is a list of list of single member lists, each of which is a possible value for the expression.
Explanation
Helper link
Takes a list of floats/integers alternating with operators (as characters) as its left argument and an operator as a character as its right argument; returns the input list after evaluating numbers separated by the relevant operator, working left to right.
œṡ⁹                  | Split once by the right argument (the operator currently being processed)
                   ? | If:
                  $  | - Following as a monad
                L    |   - Length
                 ’   |   - Decremented by 1
              ʋ      | Then, following as a dyad:
   ¹                 | - Identity function (used because of Jelly’s ordering of dyadic links at the start of a dyadic chain)
    j       ɗ        | - Join with the following as a dyad, using the original left and right arguments for this chain:
     ṪḢƭ€            |   - Tail of first item (popping from list) and head from second item (again popping from list); extracts the numbers that were either side of the operator, while removing them from the split list
         j           |   - Joined with the operator
          ŒV         |   - Evaluate as Python (rather than V because of Jelly’s handling of decimals with a leading zero)
            ß        | - Recursive call to this helper link (in case there are further of the same operator)
               F     | Else: Flatten

Main link
Takes a list of floats/integers alternating with operators (as characters)
Ḋ         | Remove first item (which will be a number)
 m2       | Every 2nd item, starting with the first (i.e. the operators)
   Q      | Uniquify
    Œ!    | Permutations
      çƒ€ | For each permuted list of operators, reduce using the helper link and with the input list as the starting point

Python 2, 182 172 bytes
import re
def f(s,P=set('+-/*')):
 S=[eval(s)]
 for p in P:
	t=s
	while p+' 'in t:t=re.sub(r'[-\d.]+ \%s [-\d.]+'%p,lambda m:`eval(m.group())`,t,1)
	S+=f(t,P-{p})
 return S

Try it online!
Takes input with ints formatted as floats, as per "Integers can be replaced by floats for languages that implicitly parse type".

Julia 1.2, 88 (82) bytes
f(t)=get(t,(),[f.([t[1:i-1];t[i+1](t[i],t[i+2]);t[i+3:end]] for i=1:2:length(t)-2)...;])

julia> f([2, +, 3, *, 4])
2-element Array{Int64,1}:
 20
 14

julia> f([18, /, 3, *, 2, -, 1])
6-element Array{Float64,1}:
 11.0
  6.0
  2.0
  3.6
  6.0
  6.0

Takes a tape in the form of a vector of numbers and infix functions, evaluates each single function call and recursively passes each resulting tape back to itself until only a single number is left. Unfortunately, get(t, (), ...) doesn't work properly in Julia 1.0, so a newer version is needed.
Six bytes can be saved, if a bunch of nested arrays is acceptable as an output:
f(t)=get(t,(),f.([t[1:i-1];t[i+1](t[i],t[i+2]);t[i+3:end]] for i=1:2:length(t)-2))

Output:

julia> f([18, /, 3, *, 2, -, 1])
3-element Array{Array{Array{Float64,1},1},1}:
 [[11.0], [6.0]]
 [[2.0], [3.6]]
 [[6.0], [6.0]]

Perl 5 (-alp), 89 bytes
my$x;map{$x.=$`.(eval$&.$1).$2.$"while/\d+[-+*\/](?=(\d+)(.*))/g}@F;$_=$x;/[-+*\/]/&&redo

TIO
or unique values, 99 bytes
my%H;map{$H{$`.(eval$&.$1).$2}++while/\d+[-+*\/](?=(\d+)(.*))/g}@F;$_=join$",keys%H;/[-+*\/]/&&redo


text:
Random point on a sphere
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

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:

Generate three random numbers \$x\$, \$y\$ and \$z\$ from a normal (Gaussian) distribution around \$0\$ and normalize them.

Implementation example

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

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

Provide in your answer a brief description of the algorithm that you are using.
Read more on sphere point picking on MathWorld.

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, so the answer using the fewest bytes in each language wins.
Standard rules, I/O rules and loophole rules apply.
Please include a Try it Online-link or equivalent to demonstrate your code working.
Please motivate your answer with an explanation of your code.

Is it okay to pick 3 reals uniformly in [-1, 1], then reject them (and repeat) if the sum of their squares isn't 1?

@Grimy I like that loophole. No, it is not allowed, because there is a theoretically zero chance of any output.

Isn't @Grimy's suggestion similar to the second example implementation mentioned by you? That solution also has theoretically zero chance of producing any output

@SaswatPadhi No, that has a chance pi/6 ≈ 0.5236 of producing an output. That's the area of the sphere inscribed in the unit-area cube

@LuisMendo I see, right. The probability is ~0.5 in that case, as you mentioned. For Grimy's proposal, it's ~0.

Are we allowed to assume that uniform distributions over $[0,1]$, $[0, 1)$ and $(0, 1)$ are effectively the same?

@ar4093 Yes, that's fine.

Wolfram Language (Mathematica), 20 bytes
RandomPoint@Sphere[]

Try it online!
Does exactly what it says on the tin.

R, 23 bytes
x=rnorm(3)
x/(x%*%x)^.5

Try it online!
Generates 3 realizations of the \$\mathcal N(0,1)\$ distribution and normalizes the resulting vector.
Plot of 1000 realizations:

Can you justify 3 axis distributed normally resulting in uniform distribution over the sphere ? (I don't see it)

@Jeffrey It's pretty well-known in probability/statistics; but the proof for 2D (which extends neatly to 3 dimensions) is approximately: let $X,Y\sim N(0,1)$ and independent. Then $f_X(x)=Ke^{-\frac{1}{2}x^2}$ and $f_Y(y)=Ke^{-\frac{1}{2}y^2}$, so by independence $f_{XY}(x,y)=K^2e^{-\frac{1}{2}(x^2+y^2)}=f_Z(z)=K^2e^{-\frac{1}{2}\lVert z\rVert^2}$ where $z=(x,y)$, so it's clear that the distribution of $z$ depends solely on the magnitude of $z$, and thus the direction is uniformly distributed.

So, the normal distribution gives us uniformly distributed points around the circle, and dividing by the magnitude ensures the points lie on the circle

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.

+1 Using aesenc to produce 128 "random" bits is just beautiful.

Python 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!
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 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, 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!
Generates three normals and scales the result. Thanks to ovs for saving 1 byte with eval in the second line.

Python 2 with numpy, 57 bytes
from numpy import*
a=random.randn(3)
print a/sum(a*a)**.5

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 @.

I like your first approach. I'm having trouble understanding how a bias towards the equator is avoided if z, from a uniform distribution, is left unmodified.

@Jitse The spherical distribution is in factor uniform over each coordinate. This is something special for dimension 3. See for instance this proof that the surface area of a slice of a sphere is proportional to its height. Regarding the intuition that this is biased to the equator, note that while slices near the equator have a bigger radius, ones near the pole are titled inward more which gives more area, and it turns out these two effects exactly cancel.

Very nice! Thanks for the clarification and the reference.

@Jitse Thanks, I added it to the body. I realized though that I was only sampling positive z though, and fixed that for a few bytes.

@Jitse Indeed, the surface area of a sphere equals the lateral surface area of the enclosing cylinder!

@Jitse Substituting $\phi$ into $x,y,z$ in the third provided algorithm results in the same approach.

@Joel You're absolutely right, not sure how I missed that.

Your second approach can be shortened using eval: 85 bytes

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!

For Octave only (i.e. not MATLAB), you can save a byte with this

@TomCarpenter Thanks! In this case as it is just one expression we can even omit the disp:)

Unity C#, 34 bytes
f=>UnityEngine.Random.onUnitSphere

Unity has a builtin for unit sphere random values, so I thought I'd post it.

Good use of a built in +1, You could just submit a function to be a bit shorter f=>Random.onUnitSphere

@LiefdeWen I knew about lambdas, I just wasn't sure if that was enough (in terms of validity on Code Golf) because it isn't declaring f's Type; using var only works inside a method and System.Func<Vector3> was longer.

In codegolf returning a function is perfectly fine, and you don't have to count the declaration either meaning you can do sneaky things with dynamic parameters. You also don't count the last semi-colon. You do however count all using statements you add. so your byte count needs to include the using. But f=>Random.onUnitSphere is a perfectly valid submission

@LiefdeWen Yeah, I just wasn't sure how the declaration was handled and wasn't really feeling up to "searching meta."

f=>UnityEngine.Random.onUnitSphere saves you the using

@Orace Doh. I forgot about doing that. Originally it was needed (because of the Vector3; 2 fully qualifieds or 1 using).

Ruby, 34 50 49 bytes
->{[z=rand*2-1]+((1-z*z)**0.5*1i**(rand*4)).rect}

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

@Jitse I forgot to scale back the x and y at high values of z. Effectively the points defined a cylinder. It's fixed now but it cost a lot of bytes! I could save a few if the the output can be expressed with a complex number [z, x+yi] I'll leave it as it is unless you say that's OK.

Looks good! I really like this approach. For consistency, the required output is three floats, so I suggest leaving it like this.

Why not use z*z instead of z**2?

@ValueInk yeah thanks I realised I'd missed that z*z. I've edited it in now. The other thing I could do is replace rand*4 with something like z*99 or x*9E9 (effectively limiting the possible values to a very fine spiral on the sphere) but I think that reduces the quality of the random.

MATL, 10 bytes
1&3Xrt2&|/

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

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, the rest are 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.

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!
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!
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

05AB1E, 23 22 bytes
[тε5°x<Ýs/<Ω}DnOtDî#}/

Implements the 2nd algorithm.
Try it online or get a few more random outputs.
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)

Using $l<1$ is equally as valid as $l\leq1$. The only criterion for $l\leq x$ is that $0< x\leq 1$. You may just as well accept vectors with $l<0.5$ if it would save any bytes. Any value smaller than or equal to $1$ removes the bias.

@Jitse Ok, implemented the normalization in both my Java and 05AB1E answers. I hope everything is correct now.

@Jitse Actually saved a byte by checking $v\leq1$ as $\lceil v\rceil==1$, instead of $v\lt1$. But thanks for the clarification that only $0\lt x\leq1$ is a requirement, and there isn't a strict $\leq$ requirement on $l$, as long as it's $\leq1$.

Processing 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;
  }

You might want to make it clearer that the implementation is not part of your byte count. I missed it on first reading, then did a double-take.

I like that the implementation uses essentially the same approach as me, though

Python 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!
Implements the first algorithm.

Python 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!
Implements the second algorithm.

@RobinRyder This implementation rejects vectors with an initial length >1, which is valid as specified in the challenge.

@Jitse Right, sorry. I misread the code.

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, 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.
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.
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.
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)

Using sqrt(1-k*k) actually saves more bytes in Java than it does in JS. :)

@Arnauld Yep. Instead of 3x M.sin, 1x M.cos and 1x M.acos, your approach uses 2x M.sin and 1x M.sqrt, which is where the additional saved bytes mostly come from. :)

108 bytes Uses a modified 2nd algorithm where I only allow values where s == 1 (instead of s<=1 and then normalizing). It sometimes gives an answer but mostly doesn't because of the timeout. Edit: Oops, I forgot to Math.sqrt the result

Actually, no, no need to sqrt because sqrt(1)==1. So I stand with my golf suggestion.

Aaah, it's been banned as a loophole -_-' My bad...

109 bytes (you can use your string output instead of double[] as that doesn't change the byte-count.)

@OlivierGrégoire Thanks! Nice switch of variable values. :)

Haskell, 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!
Does three uniforms randoms and rejection sampling.

It looks like your randoms are from the distribution (0,1) instead of (-1,1), so that only 1/8 of the sphere is covered.

@Jitse gotcha, thanks for noticing.

MathGolf, 21 19 18 bytes
{╘3Éƒ∞(ß_²Σ√_1>}▲/

Implementation of the 2nd algorithm.
Try it online or see a few more outputs at the same time.
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)

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.

Julia 1.0, 24 bytes
x=randn(3)
x/hypot(x...)

Try it online!
Draws a vector of 3 values, drawn from a normal distribution around 0 with standard deviation 1. Then just normalizes them.

randn(), from a couple of quick tests, doesn't seem to be bound to the required range. Also, this doesn't include a check for hypot() returning a value >1, which should be rejected.

@Shaggy it would appear randn simulates from a standard normal distribution rather than a uniform(0,1) one, so this approach is identical to the R one.

@Giuseppe Yes, exactly!

@Giuseppe, I think that I may not have a proper grasp on the maths behind this challenge but, if I'm understanding you correctly, you're saying that if any of the floats are outside the bounds of [-1,1) then dividing by them by the hypotenuse, which will be >1, offsets that? That leads me to wonder if the ternary in my solution is necessary ...

@Shaggy no, the normal/Gaussian distribution has some properties (specifically, rotational invariance) that the uniform doesn't have, see this comment, for example

Vyxal, 20 19 bytes
{≥|3ʁƛ∆Ṙd⌐;:²∑√:1}/

Try it Online!
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]

Pyth, 24 bytes
W<1Ks^R2JmtO2.0 3;cR@K2J

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)

Japt, 20 bytes
Port of Arnauld's implementation of the 2nd algorithm.
MhV=3ÆMrJ1
>1?ß:V®/U

Test it
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

OCaml, 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

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 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.

I'm not quite familiar with this language, but it looks like you use the random floats between 0 and 1 directly as spherical coordinates. This is incorrect, as shown in the challenge remarks 3 and 4, since you end up with a bias towards the poles of the sphere. You can correct this by applying the method shown in remark 4.

Thank you! Totally missed that. Fixed the bug and updated my answer

Looks good! Very nice first answer!

Thank you :) I was able to reduce it to under-100 bytes!

Nim, 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!
Uses the third method with a bit of substitution:
$$u = \cos(\phi) = 2j - 1$$
$$v = \sin(\phi) = \sqrt{1-u^2}$$

Raku, 65 bytes
{(&cos,&sin,{1})».(τ.rand)Z*(&sin,&sin,&cos)».(acos 2.rand-1)}

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.

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!

APL (Dyalog Classic), 22 bytes
⎕←{⍵÷((+/⍵*2)*0.5)}?⍳3

Try it online!

Hi, and welcome to Code Golf! Can you please provide a link so others can test it, for example Try It Online!?

Hej. Are you sure this solves the problem? As far as I can tell, it can only give 3 unique answers.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
Dividing Divisive Divisors
Given a positive integer \$n\$ you can always find a tuple \$(k_1,k_2,...,k_m)\$ of integers \$k_i \geqslant 2\$ such that \$k_1 \cdot k_2 \cdot ... \cdot k_m = n\$ and $$k_1 | k_2 \text{ , } k_2 | k_3 \text{ , } \ldots \text{ , }k_{m-1}|k_m.$$
Here \$a|b\$ means \$b\$ is a multiple of \$a\$, say "a divides b". If \$n>1\$ all entries \$k_i\$ must be at least \$2\$. For \$n=1\$ we have no such factor and therefore we get an empty tuple.
In case you're curious where this comes from: This decomposition is known as invariant factor decomposition in number theory and it is used in the classification of finitely generated
Abelian groups.
Challenge
Given \$n\$ output all such tuples \$(k_1,k_2,...,k_m)\$ for the given \$n\$ exactly once, in whatever order you like. The standard sequence output formats are allowed.
Examples
  1: () (empty tuple)
  2: (2)
  3: (3)
  4: (2,2), (4)
  5: (5)
  6: (6)
  7: (7)
  8: (2,2,2), (2,4), (8)
  9: (3,3), (9)
 10: (10)
 11: (11)
 12: (2,6), (12)
108: (2,54), (3,3,12), (3,6,6), (3,36), (6,18), (108)

Related: http://oeis.org/A000688, List all multiplicative partitions of n

May we output each tuple in reverse order? (e.g. 12,3,3)

@Arnauld Yes, I think as long as it is sorted in ascending or descending order it should be ok!

Can we limit input to integers >= 2? If not this would invalidate some of the existing answers?

No, the specs say clearly that any positive integer can be given as input which includes $n=1$. If I change it now everyone who actually adheres to the specs would have to change their answer.

Haskell, 66 62 60 bytes
f n=[n|n>1]:[k:l:m|k<-[2..n],l:m<-f$div n k,mod(gcd n l)k<1]

Try it online!

05AB1E, 13 bytes
Òœ€.œP€`êʒüÖP

Try it online!
Ò                      # prime factorization of the input
 œ€.œ                  # all partitions
     P                 # product of each sublist
      €`               # flatten
        ê              # sorted uniquified
         ʒ             # filter by:
          üÖ           #  pairwise divisible-by (yields list of 0s or 1s)
            P          #  product (will be 1 iff the list is all 1s)

Nice way of using Òœ€.œP to get the sublists. I indeed had trouble finding something shorter as well.. If only there was a builtin similar to Åœ but for product instead of sum. ;)

Fails for n=1 (see comments on question)

Jelly, 17 bytes
ÆfŒ!Œb€ẎP€€QḍƝẠ$Ƈ

Try it online!

JavaScript (V8),  73  70 bytes
Prints the tuples in descending order \$(k_m,k_{m-1},...,k_1)\$.
f=(n,d=2,a=[])=>n>1?d>n||f(n,d+1,a,d%a[0]||f(n/d,d,[d,...a])):print(a)

Try it online!
Commented
f = (             // f is a recursive function taking:
  n,              //   n   = input
  d = 2,          //   d   = current divisor
  a = []          //   a[] = list of divisors
) =>              //
  n > 1 ?         // if n is greater than 1:
    d > n ||      //   unless d is greater than n,
    f(            //   do a recursive call with:
      n,          //     -> n unchanged
      d + 1,      //     -> d + 1
      a,          //     -> a[] unchanged
      d % a[0] || //     unless the previous divisor does not divide the current one,
      f(          //     do another recursive call with:
        n / d,    //       -> n / d
        d,        //       -> d unchanged
        [d, ...a] //       -> d preprended to a[]
      )           //     end of inner recursive call
    )             //   end of outer recursive call
  :               // else:
    print(a)      //   this is a valid list of divisors: print it

05AB1E, 17 15 14 bytes
Ñ¦IиæʒPQ}êʒüÖP

Very slow for larger test cases.
-1 byte thanks to @Grimy.
Try it online.
Explanation:
Ñ               # Get all divisors of the (implicit) input-integer
 ¦              # Remove the first value (the 1)
  Iи            # Repeat this list (flattened) the input amount of times
                #  i.e. with input 4 we now have [2,4,2,4,2,4,2,4]
    æ           # Take the powerset of this list
     ʒ  }       # Filter it by:
      PQ        #  Where the product is equal to the (implicit) input
         ê      # Then sort and uniquify the filtered lists
          ʒ     # And filter it further by:
           ü    #  Loop over each overlapping pair of values
            Ö   #   And check if the first value is divisible by the second value
             P  #  Check if this is truthy for all pairs

                # (after which the result is output implicitly)

@Grimy Thanks. And good call on the divisors. It's still very slow for $n=8$, but all bits help, and if it doesn't cost any additional bytes to improve the performance, then why not use it. :)

13 and faster. Feels like it can be shorter still.

Wolfram Language (Mathematica), 78 76 ... 61 51 bytes
Do[i∣##&&##~#0~i,{i,2,-#/-##2/. 1:>Print@{##2}}]&

Try it online!
Recursive search tree.
Prints the list of decompositions.

JavaScript, 115 bytes
f=(n,a=[],i=1)=>{for(;i++<n;)n%i||(a=a.concat(f(n/i).filter(e=>!(e[0]%i)).map(e=>[i].concat(e))));return n>1?a:[a]}

I will write an explanation later

Japt, 22 19 bytes
ÆâÃcÅà â f@¶XcXäv)×

Try it
ÆâÃcÅà â f@¶XcXäv)×     :Implicit input of integer U
Æ                       :Map the range [0,U)
 â                      :  Divisors of U
  Ã                     :End map
   c                    :Flatten after
    Å                   :  Slicing the first element off each
     à                  :Combinations
       â                :Deduplicate
         f              :Filter by
          @             :Passing each X through the following function
           ¶            :  Test U for equality with
            Xc          :  Concatenate to X
              Xä        :    Consecutive pairs of X
                v       :    Reduced by testing divisibility
                 )      :  End concat
                  ×     :  Reduce by multiplication

Husk, 18 15 bytes
fΛ¦ufo=¹ΠṖmΠtṖp

Try it online!
Similar to Grimmy's 05AB1E answer.
-3 bytes from Zgarb.

oΛIẊ¦ can be just Λ¦. ΛV can take binary functions and check the condition on adjacent pairs.

Clean. I'll change it.

Jelly, 14 bytes
ḊœċRẎḍƝẠ×P=¥ʋƇ

Try it online!
How it works
ḊœċRẎḍƝẠ×P=¥ʋƇ - Main link. Takes n on the left
Ḋ              - Dequeue; Yield [2, 3, ..., n]
   R           - Range, Yield [1, 2, ..., n]
 œċ            - For each integer, 1 ≤ i ≤ n, yield combinations with
                  replacement of length i from the dequeued range
    Ẏ          - Tighten into a single list of lists
            ʋƇ - Keep those for which the following is true:
        ×  ¥   -   Both of the following are true:
      Ɲ        -     Over overlapping pairs:
     ḍ         -       Is the left divisible by the right?
       Ạ       -     Is all true?
        ×      -   And:
         P     -     The product
          =    -     Equals n?

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/