Merge pull request #533 from justinethier/issue-530-2

Resolve Issue 530

CHANGELOG.md

@@ -5,6 +5,7 @@
 Bug Fixes
 
 - Added a fix from Yorick Hardy to define `*ai-v4mapped*` to zero on platforms where `AI_V4MAPPED` is undefined.
+- Updated `sqrt` to properly handle complex results given non-complex input. EG: `(sqrt -1) ==> 1i`. And updated the parser to properly handle `+i` and `-i`. Thanks to Christopher Hebert for the bug reports!
 
 ## 0.36.0 - February 14, 2024
 

runtime.c

@@ -7931,9 +7931,6 @@ static void _read_return_number(void *data, port_type * p, int base, int exact)
  */
 static void _read_return_complex_number(void *data, port_type * p, int len)
 {
-//      TODO: return complex num, see _read_return_number for possible template
-//      probably want to have that function extract/identify the real/imaginary components.
-//      can just scan the buffer and read out start/end index of each number.
   int i;
   make_empty_vector(vec);
   make_string(str, p->tok_buf);

scheme/base.sld

@@ -1437,10 +1437,10 @@
       (error "exact non-negative integer required" k))
     (let* ((s (if (bignum? k)
                   (bignum-sqrt k)
-                  (exact (truncate (sqrt k)))))
+                  (exact (truncate (_sqrt k)))))
            (r (- k (* s s))))
       (values s r)))
-  (define-c sqrt
+  (define-c _sqrt
     "(void *data, int argc, closure _, object k, object z)"
     " return_inexact_double_op(data, k, sqrt, z);"
     "(void *data, object ptr, object z)"

scheme/inexact.sld

@@ -69,7 +69,6 @@
             (/ (c-log z1) (c-log z2*)))))
     (define-inexact-op c-log "log"   "clog")
     (define-inexact-op exp   "exp"   "cexp")
-    (define-inexact-op sqrt  "sqrt"  "csqrt")
     (define-inexact-op sin   "sin"   "csin")
     (define-inexact-op cos   "cos"   "ccos")
     (define-inexact-op tan   "tan"   "ctan")
@@ -93,4 +92,58 @@
                         (* (if (eqv? y -0.0) -1 1)
                            (if (eqv? x -0.0) 3.141592653589793 x))
                         (atan1 (/ y x))))))))
+
+ (define-c
+   sqrt
+   "(void *data, int argc, closure _, object k, object z)"
+   " double complex result;
+     Cyc_check_num(data, z);
+     if (obj_is_int(z)) {
+       result = csqrt(obj_obj2int(z));
+     } else if (type_of(z) == integer_tag) {
+       result = csqrt(((integer_type *)z)->value);
+     } else if (type_of(z) == bignum_tag) {
+       result = csqrt(mp_get_double(&bignum_value(z)));
+     } else if (type_of(z) == complex_num_tag) {
+       result = csqrt(complex_num_value(z));
+     } else {
+       result = csqrt(((double_type *)z)->value);
+     }
+   
+     if (cimag(result) == 0.0) {
+       if (obj_is_int(z) && creal(result) == round(creal(result))) {
+         return_closcall1(data, k, obj_int2obj(creal(result)));
+       }
+       make_double(d, creal(result));
+       return_closcall1(data, k, &d);
+     } else {
+       complex_num_type cn;
+       assign_complex_num((&cn), result);
+       return_closcall1(data, k, &cn);
+     } "
+   "(void *data, object ptr, object z)"
+   " double complex result;
+     Cyc_check_num(data, z);
+     if (obj_is_int(z)) {
+       result = csqrt(obj_obj2int(z));
+     } else if (type_of(z) == integer_tag) {
+       result = csqrt(((integer_type *)z)->value);
+     } else if (type_of(z) == bignum_tag) {
+       result = csqrt(mp_get_double(&bignum_value(z)));
+     } else if (type_of(z) == complex_num_tag) {
+       result = csqrt(complex_num_value(z));
+     } else {
+       result = csqrt(((double_type *)z)->value);
+     }
+  
+     if (cimag(result) == 0.0) {
+       if (obj_is_int(z) && creal(result) == round(creal(result))) {
+         return obj_int2obj(creal(result));
+       }
+       assign_double(ptr, creal(result));
+     } else {
+       assign_complex_num(ptr, result);
+     }
+     return ptr;
+     ")
 ))

scheme/read.sld

@@ -294,7 +294,10 @@
                                    (substring t 0 end)
                                    (substring t end (- len 1))))
                      (real (string->number real-str))
-                     (imag (string->number imag-str))
+                     (imag (cond
+                             ((equal? "+" imag-str) 1)  ;; Special case, +i w/no number
+                             ((equal? "-" imag-str) -1) ;; Special case, -i
+                             (else (string->number imag-str))))
                     )
                 (Cyc-make-rect real imag)))
              (else

tests/base.scm

@@ -9,6 +9,7 @@
 
 (import 
   (scheme base)
+  (scheme inexact)
   (cyclone test))
 
 
@@ -102,6 +103,18 @@
   (test 2.0 (denominator (inexact (/ 6 4))))
 )
 
+(test-group
+  "sqrt"
+  (test 1i (sqrt -1))
+  (test 1i (sqrt -1.0))
+  (test +i (sqrt -1.0))
+  (test 2 (sqrt 4))
+  (test 2.0 (sqrt 4.0))
+  (test 2i (sqrt -4.0))
+  (test #t (complex? (sqrt -1)))
+  (test #t (complex? (sqrt -i)))
+)
+
 (test-group
   "exact"
   (test -1 (exact -1))