Ticket #358: tt358-complex-nan.patch

File tt358-complex-nan.patch, 10.1 KB (added by rurban, 5 years ago)

unfinished

  • src/pmc/complex.pmc

    old new  
    7777        while (*t >= '0' && *t <= '9') 
    7878            t++; 
    7979    } 
     80    /* handle real NaN */ 
     81    if (*t == 'N' && *(t+1) == 'a' && *(t+2) == 'N') { 
     82        t++; 
     83        t++; 
     84        t++; 
     85    } 
    8086 
    8187    /* save the length of the real part */ 
    8288    first_num_length = t - first_num_offset; 
     
    128134                while (*t >= '0' && *t <= '9') 
    129135                    t++; 
    130136            } 
     137            /* handle imag NaN */ 
     138            if (*t == 'N' && *(t+1) == 'a' && *(t+2) == 'N') { 
     139                t++; 
     140                t++; 
     141                t++; 
     142            } 
    131143 
    132144            /* save the length of the imaginary part */ 
    133145            second_num_length = t - second_num_offset; 
     
    436448        FLOATVAL re, im; 
    437449        GET_ATTR_re(INTERP, SELF, re); 
    438450        GET_ATTR_im(INTERP, SELF, im); 
    439         return Parrot_sprintf_c(INTERP, "%vg%+vgi", re, im); 
     451        /* TT #358 NaN+i => NaN (not NaNi), 1+NaNi => NaN (not 1+NaNi) */ 
     452        if (isnan(im) || isnan(re)) 
     453            return Parrot_sprintf_c(INTERP, "NaN"); 
     454        else 
     455            return Parrot_sprintf_c(INTERP, "%vg%+vgi", re, im); 
    440456    } 
    441457 
    442458    VTABLE INTVAL get_bool() { 
     
    727743 
    728744        GET_ATTR_re(INTERP, SELF, re); 
    729745        GET_ATTR_im(INTERP, SELF, im); 
     746 
    730747        SET_ATTR_re(INTERP, dest, re + VTABLE_get_number(INTERP, value)); 
    731748        SET_ATTR_im(INTERP, dest, im); 
    732749 
     
    743760 
    744761        GET_ATTR_re(INTERP, SELF, re); 
    745762        GET_ATTR_im(INTERP, SELF, im); 
     763 
    746764        SET_ATTR_re(INTERP, dest, re + value); 
    747765        SET_ATTR_im(INTERP, dest, im); 
    748766 
     
    816834 
    817835        GET_ATTR_re(INTERP, SELF, re); 
    818836        GET_ATTR_im(INTERP, SELF, im); 
     837 
    819838        SET_ATTR_re(INTERP, dest, re - VTABLE_get_number(INTERP, value)); 
    820839        SET_ATTR_im(INTERP, dest, im); 
    821840 
     
    832851 
    833852        GET_ATTR_re(INTERP, SELF, re); 
    834853        GET_ATTR_im(INTERP, SELF, im); 
     854 
    835855        SET_ATTR_re(INTERP, dest, re - value); 
    836856        SET_ATTR_im(INTERP, dest, im); 
    837857 
     
    9991019 
    10001020=cut 
    10011021 
    1002 TODO: for better fp precision 
     1022For better fp precision: 
    10031023http://docs.sun.com/source/806-3568/ncg_goldberg.html 
    10041024(a+ib)/(c+id) = 
    10051025    (a + b(d/c)) / (c + d(d/c)) + i(b - a(d/c)) / (c + d(d/c)) if |d|<|c| 
     
    10081028*/ 
    10091029 
    10101030    MULTI PMC *divide(Complex value, PMC *dest) { 
    1011         FLOATVAL mod, re, im; 
     1031        FLOATVAL re, im; 
    10121032        FLOATVAL self_re, self_im, val_re, val_im; 
    10131033 
    10141034        complex_check_divide_zero(INTERP, value); 
     
    10191039        GET_ATTR_re(INTERP, value, val_re); 
    10201040        GET_ATTR_im(INTERP, value, val_im); 
    10211041 
    1022         /* a little speed optimisation: cache an intermediate number; 
    1023             I'm not sure the compiler does this */ 
    1024  
    10251042        if (self_im == 0.0 && val_im == 0.0) { 
    10261043            re = self_re / val_re; 
    10271044            im = 0.0; 
    10281045        } 
     1046        else if (abs(val_im) < abs(val_re)) { 
     1047            /* a: self_re, b: self_im, c: val_re, d: val_im */ 
     1048            /* (a + b(d/c)) / (c + d(d/c)) + i(b - a(d/c)) / (c + d(d/c)) if |d|<|c| */ 
     1049            FLOATVAL mod; 
     1050            /* A little speed optimisation: cache intermediate numbers; 
     1051               I'm not sure the compiler does this */ 
     1052            FLOATVAL d_c = val_im / val_re; 
     1053            mod = val_re + val_im * d_c; 
     1054            re  = (self_re + self_im * d_c) / mod; 
     1055            im  = (self_im - self_re * d_c) / mod; 
     1056        } 
    10291057        else { 
    1030             mod = (val_re * val_re + val_im * val_im); 
    1031             re  = (self_re * val_re + self_im * val_im) / mod; 
    1032             im  = (self_im * val_re - self_re * val_im) / mod; 
     1058            /* (b + a(c/d)) / (d + c(c/d)) + i(-a + b(c/d)) / (d + c(c/d)) if |d|>=|c| */ 
     1059            FLOATVAL mod; 
     1060            FLOATVAL c_d = val_re / val_im; 
     1061            mod = val_im + val_re * c_d; 
     1062            re  = (self_im + self_re * c_d) / mod; 
     1063            im  = (-self_re + self_im * c_d) / mod; 
    10331064        } 
    10341065 
    10351066        SET_ATTR_re(INTERP, dest, re); 
     
    10961127            re = self_re / val_re; 
    10971128            im = 0.0; 
    10981129        } 
     1130        /* For better fp precision: */ 
     1131        /* http://docs.sun.com/source/806-3568/ncg_goldberg.html */ 
     1132        else if (abs(val_im) < abs(val_re)) { 
     1133            FLOATVAL mod; 
     1134            FLOATVAL d_c = val_im / val_re; 
     1135            mod = val_re + val_im * d_c; 
     1136            re  = (self_re + self_im * d_c) / mod; 
     1137            im  = (self_im - self_re * d_c) / mod; 
     1138        } 
    10991139        else { 
    1100             /* a little speed optimisation: cache an intermediate number; 
    1101                I'm not sure the compiler does this */ 
    1102             mod = (val_re * val_re + val_im * val_im); 
    1103             re  = (self_re * val_re + self_im * val_im) / mod; 
    1104             im  = (self_im * val_re - self_re * val_im) / mod; 
     1140            FLOATVAL mod; 
     1141            FLOATVAL c_d = val_re / val_im; 
     1142            mod = val_im + val_re * c_d; 
     1143            re  = (self_im + self_re * c_d) / mod; 
     1144            im  = (-self_re + self_im * c_d) / mod; 
    11051145        } 
    11061146 
    11071147        SET_ATTR_re(INTERP, SELF, re); 
     
    12191259*/ 
    12201260 
    12211261/* 
     1262Straightforward: 
     1263  d = sqrt(re*re + im*im); 
    12221264 
    1223   TODO for better precision: hinted by vaxman according to "Numerical Recipes 
    1224   in Fortran 77", 2nd edition, Press, Vetterling, Teukolsky, Flannery, 
    1225   Cambridge University Press, 2001, pp. 171ff: 
    1226  
     1265For better precision: hinted by vaxman according to "Numerical Recipes 
     1266in Fortran 77", 2nd edition, Press, Vetterling, Teukolsky, Flannery, 
     1267Cambridge University Press, 2001, pp. 171ff: 
    12271268 
    12281269|a+ib|=|a|*sqrt(1+(b/a)**2), if |a|>=|b|, 
    12291270       |b|*sqrt(1+(a/b)**2)  else. 
     
    12341275        FLOATVAL re, im, d; 
    12351276        GET_ATTR_re(INTERP, SELF, re); 
    12361277        GET_ATTR_im(INTERP, SELF, im); 
    1237         d = sqrt(re*re + im*im); 
     1278        if (abs(re) >= abs(im)) { 
     1279            FLOATVAL b_a = im/re; 
     1280            d = abs(re)*sqrt(1+b_a*b_a); 
     1281        } 
     1282        else { 
     1283            FLOATVAL a_b = re/im; 
     1284            d = abs(im)*sqrt(1+a_b*a_b); 
     1285        } 
    12381286 
    12391287        dest = pmc_new(INTERP, 
    12401288            Parrot_get_ctx_HLL_type(INTERP, enum_class_Float)); 
     
    12471295        FLOATVAL re, im, d; 
    12481296        GET_ATTR_re(INTERP, SELF, re); 
    12491297        GET_ATTR_im(INTERP, SELF, im); 
    1250         d = sqrt(re*re + im*im); 
     1298        if (abs(re) >= abs(im)) { 
     1299            FLOATVAL b_a = im/re; 
     1300            d = abs(re)*sqrt(1+b_a*b_a); 
     1301        } 
     1302        else { 
     1303            FLOATVAL a_b = re/im; 
     1304            d = abs(im)*sqrt(1+a_b*a_b); 
     1305        } 
    12511306        pmc_reuse(INTERP, SELF, enum_class_Float, 0); 
    12521307        VTABLE_set_number_native(INTERP, SELF, d); 
    12531308    } 
  • t/pmc/complex.t

    old new  
    2121    .include 'fp_equality.pasm' 
    2222    .include "iglobals.pasm" 
    2323 
    24     plan(467) 
     24    plan(483) 
    2525 
    2626    string_parsing() 
    2727    exception_malformed_string__real_part() 
     
    7777    sech_of_complex_numbers() 
    7878    csch_of_complex_numbers() 
    7979    add_using_subclass_of_complex_bug_59630() 
     80    test_nan() 
    8081 
    8182    # END_OF_TESTS 
    8283 
     
    730731 
    731732    #XXX: can't do $P1.'$S2'() 
    732733    $P2 = $P1. $S2() 
    733     $S3 = sprintf "%f%+fi", $P2 
     734    $S4 = sprintf "%f%+fi", $P2 
    734735 
    735736    concat $S5, $S2, " of " 
    736737    concat $S5, $S5, $S4 
     738    concat $S5, $S5, " " 
     739    concat $S5, $S5, $S3 
    737740 
    738     $I0 = cmp_str $S1, $S3 
     741    $I0 = cmp_str $S1, $S4 
    739742    $I0 = not $I0 
    740743 
    741     todo( $I0, $S4 ) 
     744    todo( $I0, $S5 ) 
    742745.endm 
    743746 
    744747.sub ln_of_complex_numbers 
     
    11731176    .complex_op_is("-2-3i", "0.272549+0.040301i", 'csch' ) 
    11741177.end 
    11751178 
     1179.macro op_float_test_is(op, c_arg, f_arg, result) 
     1180    c  = .c_arg 
     1181    f  = .f_arg 
     1182    set $S2, .op 
     1183    c2 = c. $S2(f, $P1) 
     1184    $S0 = sprintf "%f%+fi", c2 
     1185    $S1 = .result 
     1186 
     1187    #is( $S0, $S1, $S1 ) 
     1188    concat $S3, $S2, " of " 
     1189    concat $S3, $S3, $S0 
     1190    concat $S3, $S3, " TT #358" 
     1191 
     1192    $I0 = cmp_str $S1, $S0 
     1193    $I0 = not $I0 
     1194 
     1195    todo( $I0, $S3 ) 
     1196.endm 
     1197 
     1198.macro op_c_test_is(op, arg1, arg2, result) 
     1199    c  = .arg1 
     1200    c1 = .arg2 
     1201    set $S2, .op 
     1202    c2 = c. $S2(c1, $P1) 
     1203    $S0 = sprintf "%f%+fi", c2 
     1204    $S1 = .result 
     1205    #is( $S0, $S1, $S1 ) 
     1206    concat $S3, $S2, " of " 
     1207    concat $S3, $S3, $S0 
     1208    concat $S3, $S3, " TT #358" 
     1209 
     1210    $I0 = cmp_str $S1, $S0 
     1211    $I0 = not $I0 
     1212 
     1213    todo( $I0, $S3 ) 
     1214.endm 
     1215 
     1216# TT #358 
     1217.sub test_nan 
     1218    .local pmc c, c1, c2, f 
     1219    c  = new ['Complex'] 
     1220    c1 = new ['Complex'] 
     1221    c2 = new ['Complex'] 
     1222    f = new ['Float'] 
     1223    set c, "NaN + i" 
     1224    is(c, 'NaN', "parse NaN + i") 
     1225    set c2, "1.0-NaNi" 
     1226    is(c2, 'NaN', "parse 1.0-NaNi") 
     1227    c2 = c.'multiply'($N1, 1.0) 
     1228    is(c2, 'NaN', "mul NaN+i, 1.0") 
     1229 
     1230    .op_float_test_is( "multiply", "NaN", 1.0,   "NaN" ) 
     1231    .op_float_test_is( "multiply", "NaN+i", "NaN", "NaN" ) 
     1232    .op_float_test_is( "multiply", "NaNi",  1.0,   "NaN" ) 
     1233    .op_c_test_is( "multiply", "NaN+i", "0+NaNi",   "NaN" ) 
     1234    .op_c_test_is( "add", "0+NaNi", "i", "NaN" ) 
     1235    .op_c_test_is( "add", "i", "0+NaNi", "NaN" ) 
     1236    .op_c_test_is( "subtract", "NaN", "i", "NaN" ) 
     1237    .op_c_test_is( "subtract", "i", "NaN", "NaN" ) 
     1238    .op_c_test_is( "pow", "NaN", "i", "NaN" ) 
     1239    .op_c_test_is( "pow", "i", "NaN", "NaN" ) 
     1240 
     1241    .complex_op_todo("NaN", "NaN", 'ln',   "TT #358" ) 
     1242    .complex_op_todo("NaN", "NaN", 'sin',  "TT #358" ) 
     1243    .complex_op_todo("NaN", "NaN", 'csch', "TT #358" ) 
     1244.end 
     1245 
    11761246.sub add_using_subclass_of_complex_bug_59630 
    11771247    skip( 3, 'add using subclass of Complex - RT #59630' ) 
    11781248    .return()