Fix an `sqrt` crash, make it faster and more precise

[jpellegrini]

Hello @egallesio !

Currently, STklos crashes on this input: (sqrt 1-0.0i).

Looking at the algorithm used, it seems that we can both fix this crash and also make it a bit more precise.

With this PR, we calculate the square root of coplexes as this (same algorithm that Clisp uses):

sqrt(a+bi) = A+Bi

   if a < 0:
       B := sqrt( (|z|-a) /2) * sign(b)
       A := b/(2*B).

   if a >= 0:
       A := sqrt( (|z|+a) /2)
       if A != 0:
           B := (b / (2*A))
       else:
           B = 0.0

So we don't use ex2inex, angle (and consequently atan) anymore, only two square roots and one inversion (and it seems to be more precise - see below).

A tricky part is to test "a < 0" for negative zero (because -0.0 < 0 won't work, we need to get the REAL_PART and use signbit() to test it)

Looking at the limiting process, we see that

(sqrt -1) = i.

Then, (sqrt -1 + tiny*i) should get be very close to i.

Current algorithm:

(import (srfi 144))
(sqrt (+ -1 (* fl-least +i))) => 6.12323399573677e-17+1.0i

With this patch:

(sqrt (+ -1 (* fl-least +i))) => 0.0+1.0i

So we go from 6e-17 to 0.0 (great!)

We are also closer to the values of Clisp, Gambit, Guile, SBCL for complex square roots...

And... It seems faster!

(time
 (let ((x 2000-3000i))
   (dotimes (i 3000000)
     (when (= (modulo i 100)) (set! x 2000-3000i))
     (set! x (sqrt x)))
   x))

Old code: 3623.08 ms
New code: 3002.49 ms

[jpellegrini]

By the way, the line

     (when (= (modulo i 100)) (set! x 2000-3000i))

in the above micro-benchmark is necessary because otherwise the current system would try to do (sqrt 1.0-0.0i) and would crash.