APL (NARS), 2 bytes

√

NARS has built-in support for quaternions. ¯_(⍨)_/¯

Python 2, 72 bytes

def f(a,b,c,d):s=((a+(a*a+b*b+c*c+d*d)**.5)*2)**.5;print s/2,b/s,c/s,d/s

Try it online!

More or less a raw formula. I thought I could use list comprehensions to loop over b,c,d, but this seems to be longer. Python is really hurt here by a lack of vector operations, in particular scaling and norm.

Python 3, 77 bytes

def f(a,*l):r=a+sum(x*x for x in[a,*l])**.5;return[x/(r*2)**.5for x in[r,*l]]

Try it online!

Solving the quadratic directly was also shorter than using Python's complex-number square root to solve it like in the problem statement.

Wolfram Language (Mathematica), 19 bytes

Sqrt
<<Quaternions`

Try it online!

Mathematica has Quaternion built-in too, but is more verbose.

_{Although built-ins look cool, do upvote solutions that don't use built-ins too! I don't want votes on questions reaching HNQ be skewed.}

JavaScript (ES7), 55 53 bytes

Based on the direct formula used by xnor.

Takes input as an array.

q=>q.map(v=>1/q?v/2/q:q=((v+Math.hypot(...q))/2)**.5)

Try it online!

How?

Given an array $q=[a,b,c,d]$, this computes:

$$x=sqrtfraca+sqrta^2+b^2+c^2+d^22$$

And returns:

$$left[x,fracb2x,fracc2x,fracd2xright]$$

q => // q = input array
 q.map(v => // for each value v in q:
 1 / q ? // if q is numeric (2nd to 4th iteration):
 v / 2 / q // yield v / 2q
 : // else (1st iteration, with v = a):
 q = ( // compute x (as defined above) and store it in q
 (v + Math.hypot(...q)) // we use Math.hypot(...q) to compute:
 / 2 // (q[0]**2 + q[1]**2 + q[2]**2 + q[3]**2) ** 0.5
 ) ** .5 // yield x
 ) // end of map()

Haskell, 51 bytes

f(a:l)|r<-a+sqrt(sum$(^2)<$>a:l)=(/sqrt(r*2))<$>r:l

Try it online!

A direct formula. The main trick to express the real part of the output as r/sqrt(r*2) to parallel the imaginary part expression, which saves a few bytes over:

54 bytes

f(a:l)|s<-sqrt$2*(a+sqrt(sum$(^2)<$>a:l))=s/2:map(/s)l

Try it online!

Charcoal, 32 bytes

≔Ｘ⊗⁺§θ⁰ＸΣＥθ×ιι·⁵¦·⁵η≧∕ηθ§≔θ⁰⊘ηＩθ

Try it online! Link is to verbose version of code. Port of @xnor's Python answer. Explanation:

≔Ｘ⊗⁺§θ⁰ＸΣＥθ×ιι·⁵¦·⁵η

Square all of the elements of the input and take the sum, then take the square root. This calculates $ | x + yvecu | = sqrt x^2 + y^2 = sqrt (a^2 - b^2)^2 + (2ab)^2 = a^2 + b^2 $. Adding $ x $ gives $ 2a^2 $ which is then doubled and square rooted to give $ 2a $.

≧∕ηθ

Because $ y = 2ab $, calculate $ b $ by dividing by $ 2a $.

§≔θ⁰⊘η

Set the first element of the array (i.e. the real part) to half of $ 2a $.

Ｉθ

Cast the values to string and implicitly print.

Java 8, 84 bytes

(a,b,c,d)->(a=Math.sqrt(2*(a+Math.sqrt(a*a+b*b+c*c+d*d))))/2+" "+b/a+" "+c/a+" "+d/a

Port of @xnor's Python 2 answer.

Try it online.

Explanation:

(a,b,c,d)-> // Method with four double parameters and String return-type
 (a= // Change `a` to:
 Math.sqrt( // The square root of:
 2* // Two times:
 (a+ // `a` plus,
 Math.sqrt( // the square-root of:
 a*a // `a` squared,
 +b*b // `b` squared,
 +c*c // `c` squared,
 +d*d)))) // And `d` squared summed together
 /2 // Then return this modified `a` divided by 2
 +" "+b/a // `b` divided by the modified `a`
 +" "+c/a // `c` divided by the modified `a`
 +" "+d/a // And `d` divided by the modified `a`, with space delimiters

05AB1E, 14 bytes

nOtsн+·t©/¦®;š

Port of @xnor's Python 2 answer.

Try it online or verify all test cases.

Explanation:

n # Square each number in the (implicit) input-list
 O # Sum them
 t # Take the square-root of that
 sн+ # Add the first item of the input-list
 · # Double it
 t # Take the square-root of it
 © # Store it in the register (without popping)
 / # Divide each value in the (implicit) input with it
 ¦ # Remove the first item
 ®; # Push the value from the register again, and halve it
 š # Prepend it to the list (and output implicitly)

Wolfram Language (Mathematica), 28 bytes

s=#+Norm@##,##2/(2s)^.5&

Port of @xnor's Python 2 answer.

Try it online!

C# .NET, 88 bytes

(a,b,c,d)=>((a=System.Math.Sqrt(2*(a+System.Math.Sqrt(a*a+b*b+c*c+d*d))))/2,b/a,c/a,d/a)

Port of my Java 8 answer, but returns a Tuple instead of a String. I thought that would have been shorter, but unfortunately the Math.Sqrt require a System-import in C# .NET, ending up at 4 bytes longer instead of 10 bytes shorter.. >.>

The lambda declaration looks pretty funny, though:

System.Func<double, double, double, double, (double, double, double, double)> f =

Try it online.

Perl 6, 49 bytes

;(*+@^b>>².sum**.5*i).sqrt.&.re,(@b X/2*.re)

Try it online!

Curried function taking input as f(b,c,d)(a). Returns quaternion as a,(b,c,d).

Explanation

; # Block returning WhateverCode
 @^b>>².sum**.5 # Compute B of quaternion written as q = a + B*u
 # (length of vector (b,c,d))
 (*+ *i) # Complex number a + B*i
 .sqrt # Square root of complex number
 .& # Return
 .re, # Real part of square root
 (@b X/2*.re) # b,c,d divided by 2* real part