Ticket #358: tt358-complex-nan.patch

src/pmc/complex.pmc

@@ -77 +77 @@
         while (*t >= '0' && *t <= '9')
             t++;
     }
+    /* handle real NaN */
+    if (*t == 'N' && *(t+1) == 'a' && *(t+2) == 'N') {
+        t++;
+        t++;
+        t++;
+    }
 
     /* save the length of the real part */
     first_num_length = t - first_num_offset;
@@ -128 +134 @@
                 while (*t >= '0' && *t <= '9')
                     t++;
             }
+            /* handle imag NaN */
+            if (*t == 'N' && *(t+1) == 'a' && *(t+2) == 'N') {
+                t++;
+                t++;
+                t++;
+            }
 
             /* save the length of the imaginary part */
             second_num_length = t - second_num_offset;
@@ -436 +448 @@
         FLOATVAL re, im;
         GET_ATTR_re(INTERP, SELF, re);
         GET_ATTR_im(INTERP, SELF, im);
-        return Parrot_sprintf_c(INTERP, "%vg%+vgi", re, im);
+        /* TT #358 NaN+i => NaN (not NaNi), 1+NaNi => NaN (not 1+NaNi) */
+        if (isnan(im) || isnan(re))
+            return Parrot_sprintf_c(INTERP, "NaN");
+        else
+            return Parrot_sprintf_c(INTERP, "%vg%+vgi", re, im);
     }
 
     VTABLE INTVAL get_bool() {
@@ -727 +743 @@
 
         GET_ATTR_re(INTERP, SELF, re);
         GET_ATTR_im(INTERP, SELF, im);
+
         SET_ATTR_re(INTERP, dest, re + VTABLE_get_number(INTERP, value));
         SET_ATTR_im(INTERP, dest, im);
 
@@ -743 +760 @@
 
         GET_ATTR_re(INTERP, SELF, re);
         GET_ATTR_im(INTERP, SELF, im);
+
         SET_ATTR_re(INTERP, dest, re + value);
         SET_ATTR_im(INTERP, dest, im);
 
@@ -816 +834 @@
 
         GET_ATTR_re(INTERP, SELF, re);
         GET_ATTR_im(INTERP, SELF, im);
+
         SET_ATTR_re(INTERP, dest, re - VTABLE_get_number(INTERP, value));
         SET_ATTR_im(INTERP, dest, im);
 
@@ -832 +851 @@
 
         GET_ATTR_re(INTERP, SELF, re);
         GET_ATTR_im(INTERP, SELF, im);
+
         SET_ATTR_re(INTERP, dest, re - value);
         SET_ATTR_im(INTERP, dest, im);
 
@@ -999 +1019 @@
 
 =cut
 
-TODO: for better fp precision
+For better fp precision:
 http://docs.sun.com/source/806-3568/ncg_goldberg.html
 (a+ib)/(c+id) =
     (a + b(d/c)) / (c + d(d/c)) + i(b - a(d/c)) / (c + d(d/c)) if |d|<|c|
@@ -1008 +1028 @@
 */
 
     MULTI PMC *divide(Complex value, PMC *dest) {
-        FLOATVAL mod, re, im;
+        FLOATVAL re, im;
         FLOATVAL self_re, self_im, val_re, val_im;
 
         complex_check_divide_zero(INTERP, value);
@@ -1019 +1039 @@
         GET_ATTR_re(INTERP, value, val_re);
         GET_ATTR_im(INTERP, value, val_im);
 
-        /* a little speed optimisation: cache an intermediate number;
-            I'm not sure the compiler does this */
-
         if (self_im == 0.0 && val_im == 0.0) {
             re = self_re / val_re;
             im = 0.0;
         }
+        else if (abs(val_im) < abs(val_re)) {
+            /* a: self_re, b: self_im, c: val_re, d: val_im */
+            /* (a + b(d/c)) / (c + d(d/c)) + i(b - a(d/c)) / (c + d(d/c)) if |d|<|c| */
+            FLOATVAL mod;
+            /* A little speed optimisation: cache intermediate numbers;
+               I'm not sure the compiler does this */
+            FLOATVAL d_c = val_im / val_re;
+            mod = val_re + val_im * d_c;
+            re  = (self_re + self_im * d_c) / mod;
+            im  = (self_im - self_re * d_c) / mod;
+        }
         else {
-            mod = (val_re * val_re + val_im * val_im);
-            re  = (self_re * val_re + self_im * val_im) / mod;
-            im  = (self_im * val_re - self_re * val_im) / mod;
+            /* (b + a(c/d)) / (d + c(c/d)) + i(-a + b(c/d)) / (d + c(c/d)) if |d|>=|c| */
+            FLOATVAL mod;
+            FLOATVAL c_d = val_re / val_im;
+            mod = val_im + val_re * c_d;
+            re  = (self_im + self_re * c_d) / mod;
+            im  = (-self_re + self_im * c_d) / mod;
         }
 
         SET_ATTR_re(INTERP, dest, re);
@@ -1096 +1127 @@
             re = self_re / val_re;
             im = 0.0;
         }
+        /* For better fp precision: */
+        /* http://docs.sun.com/source/806-3568/ncg_goldberg.html */
+        else if (abs(val_im) < abs(val_re)) {
+            FLOATVAL mod;
+            FLOATVAL d_c = val_im / val_re;
+            mod = val_re + val_im * d_c;
+            re  = (self_re + self_im * d_c) / mod;
+            im  = (self_im - self_re * d_c) / mod;
+        }
         else {
-            /* a little speed optimisation: cache an intermediate number;
-               I'm not sure the compiler does this */
-            mod = (val_re * val_re + val_im * val_im);
-            re  = (self_re * val_re + self_im * val_im) / mod;
-            im  = (self_im * val_re - self_re * val_im) / mod;
+            FLOATVAL mod;
+            FLOATVAL c_d = val_re / val_im;
+            mod = val_im + val_re * c_d;
+            re  = (self_im + self_re * c_d) / mod;
+            im  = (-self_re + self_im * c_d) / mod;
         }
 
         SET_ATTR_re(INTERP, SELF, re);
@@ -1219 +1259 @@
 */
 
 /*
+Straightforward:
+  d = sqrt(re*re + im*im);
 
-  TODO for better precision: hinted by vaxman according to "Numerical Recipes
-  in Fortran 77", 2nd edition, Press, Vetterling, Teukolsky, Flannery,
-  Cambridge University Press, 2001, pp. 171ff:
-
+For better precision: hinted by vaxman according to "Numerical Recipes
+in Fortran 77", 2nd edition, Press, Vetterling, Teukolsky, Flannery,
+Cambridge University Press, 2001, pp. 171ff:
 
 |a+ib|=|a|*sqrt(1+(b/a)**2), if |a|>=|b|,
        |b|*sqrt(1+(a/b)**2)  else.
@@ -1234 +1275 @@
         FLOATVAL re, im, d;
         GET_ATTR_re(INTERP, SELF, re);
         GET_ATTR_im(INTERP, SELF, im);
-        d = sqrt(re*re + im*im);
+        if (abs(re) >= abs(im)) {
+            FLOATVAL b_a = im/re;
+            d = abs(re)*sqrt(1+b_a*b_a);
+        }
+        else {
+            FLOATVAL a_b = re/im;
+            d = abs(im)*sqrt(1+a_b*a_b);
+        }
 
         dest = pmc_new(INTERP,
             Parrot_get_ctx_HLL_type(INTERP, enum_class_Float));
@@ -1247 +1295 @@
         FLOATVAL re, im, d;
         GET_ATTR_re(INTERP, SELF, re);
         GET_ATTR_im(INTERP, SELF, im);
-        d = sqrt(re*re + im*im);
+        if (abs(re) >= abs(im)) {
+            FLOATVAL b_a = im/re;
+            d = abs(re)*sqrt(1+b_a*b_a);
+        }
+        else {
+            FLOATVAL a_b = re/im;
+            d = abs(im)*sqrt(1+a_b*a_b);
+        }
         pmc_reuse(INTERP, SELF, enum_class_Float, 0);
         VTABLE_set_number_native(INTERP, SELF, d);
     }

t/pmc/complex.t

@@ -21 +21 @@
     .include 'fp_equality.pasm'
     .include "iglobals.pasm"
 
-    plan(467)
+    plan(483)
 
     string_parsing()
     exception_malformed_string__real_part()
@@ -77 +77 @@
     sech_of_complex_numbers()
     csch_of_complex_numbers()
     add_using_subclass_of_complex_bug_59630()
+    test_nan()
 
     # END_OF_TESTS
 
@@ -730 +731 @@
 
     #XXX: can't do $P1.'$S2'()
     $P2 = $P1. $S2()
-    $S3 = sprintf "%f%+fi", $P2
+    $S4 = sprintf "%f%+fi", $P2
 
     concat $S5, $S2, " of "
     concat $S5, $S5, $S4
+    concat $S5, $S5, " "
+    concat $S5, $S5, $S3
 
-    $I0 = cmp_str $S1, $S3
+    $I0 = cmp_str $S1, $S4
     $I0 = not $I0
 
-    todo( $I0, $S4 )
+    todo( $I0, $S5 )
 .endm
 
 .sub ln_of_complex_numbers
@@ -1173 +1176 @@
     .complex_op_is("-2-3i", "0.272549+0.040301i", 'csch' )
 .end
 
+.macro op_float_test_is(op, c_arg, f_arg, result)
+    c  = .c_arg
+    f  = .f_arg
+    set $S2, .op
+    c2 = c. $S2(f, $P1)
+    $S0 = sprintf "%f%+fi", c2
+    $S1 = .result
+
+    #is( $S0, $S1, $S1 )
+    concat $S3, $S2, " of "
+    concat $S3, $S3, $S0
+    concat $S3, $S3, " TT #358"
+
+    $I0 = cmp_str $S1, $S0
+    $I0 = not $I0
+
+    todo( $I0, $S3 )
+.endm
+
+.macro op_c_test_is(op, arg1, arg2, result)
+    c  = .arg1
+    c1 = .arg2
+    set $S2, .op
+    c2 = c. $S2(c1, $P1)
+    $S0 = sprintf "%f%+fi", c2
+    $S1 = .result
+    #is( $S0, $S1, $S1 )
+    concat $S3, $S2, " of "
+    concat $S3, $S3, $S0
+    concat $S3, $S3, " TT #358"
+
+    $I0 = cmp_str $S1, $S0
+    $I0 = not $I0
+
+    todo( $I0, $S3 )
+.endm
+
+# TT #358
+.sub test_nan
+    .local pmc c, c1, c2, f
+    c  = new ['Complex']
+    c1 = new ['Complex']
+    c2 = new ['Complex']
+    f = new ['Float']
+    set c, "NaN + i"
+    is(c, 'NaN', "parse NaN + i")
+    set c2, "1.0-NaNi"
+    is(c2, 'NaN', "parse 1.0-NaNi")
+    c2 = c.'multiply'($N1, 1.0)
+    is(c2, 'NaN', "mul NaN+i, 1.0")
+
+    .op_float_test_is( "multiply", "NaN", 1.0,   "NaN" )
+    .op_float_test_is( "multiply", "NaN+i", "NaN", "NaN" )
+    .op_float_test_is( "multiply", "NaNi",  1.0,   "NaN" )
+    .op_c_test_is( "multiply", "NaN+i", "0+NaNi",   "NaN" )
+    .op_c_test_is( "add", "0+NaNi", "i", "NaN" )
+    .op_c_test_is( "add", "i", "0+NaNi", "NaN" )
+    .op_c_test_is( "subtract", "NaN", "i", "NaN" )
+    .op_c_test_is( "subtract", "i", "NaN", "NaN" )
+    .op_c_test_is( "pow", "NaN", "i", "NaN" )
+    .op_c_test_is( "pow", "i", "NaN", "NaN" )
+
+    .complex_op_todo("NaN", "NaN", 'ln',   "TT #358" )
+    .complex_op_todo("NaN", "NaN", 'sin',  "TT #358" )
+    .complex_op_todo("NaN", "NaN", 'csch', "TT #358" )
+.end
+
 .sub add_using_subclass_of_complex_bug_59630
     skip( 3, 'add using subclass of Complex - RT #59630' )
     .return()