Ticket #358: tt358complexnan.patch
File tt358complexnan.patch, 10.1 KB (added by rurban, 6 years ago) 


src/pmc/complex.pmc
old new 77 77 while (*t >= '0' && *t <= '9') 78 78 t++; 79 79 } 80 /* handle real NaN */ 81 if (*t == 'N' && *(t+1) == 'a' && *(t+2) == 'N') { 82 t++; 83 t++; 84 t++; 85 } 80 86 81 87 /* save the length of the real part */ 82 88 first_num_length = t  first_num_offset; … … 128 134 while (*t >= '0' && *t <= '9') 129 135 t++; 130 136 } 137 /* handle imag NaN */ 138 if (*t == 'N' && *(t+1) == 'a' && *(t+2) == 'N') { 139 t++; 140 t++; 141 t++; 142 } 131 143 132 144 /* save the length of the imaginary part */ 133 145 second_num_length = t  second_num_offset; … … 436 448 FLOATVAL re, im; 437 449 GET_ATTR_re(INTERP, SELF, re); 438 450 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); 440 456 } 441 457 442 458 VTABLE INTVAL get_bool() { … … 727 743 728 744 GET_ATTR_re(INTERP, SELF, re); 729 745 GET_ATTR_im(INTERP, SELF, im); 746 730 747 SET_ATTR_re(INTERP, dest, re + VTABLE_get_number(INTERP, value)); 731 748 SET_ATTR_im(INTERP, dest, im); 732 749 … … 743 760 744 761 GET_ATTR_re(INTERP, SELF, re); 745 762 GET_ATTR_im(INTERP, SELF, im); 763 746 764 SET_ATTR_re(INTERP, dest, re + value); 747 765 SET_ATTR_im(INTERP, dest, im); 748 766 … … 816 834 817 835 GET_ATTR_re(INTERP, SELF, re); 818 836 GET_ATTR_im(INTERP, SELF, im); 837 819 838 SET_ATTR_re(INTERP, dest, re  VTABLE_get_number(INTERP, value)); 820 839 SET_ATTR_im(INTERP, dest, im); 821 840 … … 832 851 833 852 GET_ATTR_re(INTERP, SELF, re); 834 853 GET_ATTR_im(INTERP, SELF, im); 854 835 855 SET_ATTR_re(INTERP, dest, re  value); 836 856 SET_ATTR_im(INTERP, dest, im); 837 857 … … 999 1019 1000 1020 =cut 1001 1021 1002 TODO: for better fp precision 1022 For better fp precision: 1003 1023 http://docs.sun.com/source/8063568/ncg_goldberg.html 1004 1024 (a+ib)/(c+id) = 1005 1025 (a + b(d/c)) / (c + d(d/c)) + i(b  a(d/c)) / (c + d(d/c)) if d<c … … 1008 1028 */ 1009 1029 1010 1030 MULTI PMC *divide(Complex value, PMC *dest) { 1011 FLOATVAL mod,re, im;1031 FLOATVAL re, im; 1012 1032 FLOATVAL self_re, self_im, val_re, val_im; 1013 1033 1014 1034 complex_check_divide_zero(INTERP, value); … … 1019 1039 GET_ATTR_re(INTERP, value, val_re); 1020 1040 GET_ATTR_im(INTERP, value, val_im); 1021 1041 1022 /* a little speed optimisation: cache an intermediate number;1023 I'm not sure the compiler does this */1024 1025 1042 if (self_im == 0.0 && val_im == 0.0) { 1026 1043 re = self_re / val_re; 1027 1044 im = 0.0; 1028 1045 } 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 } 1029 1057 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; 1033 1064 } 1034 1065 1035 1066 SET_ATTR_re(INTERP, dest, re); … … 1096 1127 re = self_re / val_re; 1097 1128 im = 0.0; 1098 1129 } 1130 /* For better fp precision: */ 1131 /* http://docs.sun.com/source/8063568/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 } 1099 1139 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; 1105 1145 } 1106 1146 1107 1147 SET_ATTR_re(INTERP, SELF, re); … … 1219 1259 */ 1220 1260 1221 1261 /* 1262 Straightforward: 1263 d = sqrt(re*re + im*im); 1222 1264 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 1265 For better precision: hinted by vaxman according to "Numerical Recipes 1266 in Fortran 77", 2nd edition, Press, Vetterling, Teukolsky, Flannery, 1267 Cambridge University Press, 2001, pp. 171ff: 1227 1268 1228 1269 a+ib=a*sqrt(1+(b/a)**2), if a>=b, 1229 1270 b*sqrt(1+(a/b)**2) else. … … 1234 1275 FLOATVAL re, im, d; 1235 1276 GET_ATTR_re(INTERP, SELF, re); 1236 1277 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 } 1238 1286 1239 1287 dest = pmc_new(INTERP, 1240 1288 Parrot_get_ctx_HLL_type(INTERP, enum_class_Float)); … … 1247 1295 FLOATVAL re, im, d; 1248 1296 GET_ATTR_re(INTERP, SELF, re); 1249 1297 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 } 1251 1306 pmc_reuse(INTERP, SELF, enum_class_Float, 0); 1252 1307 VTABLE_set_number_native(INTERP, SELF, d); 1253 1308 } 
t/pmc/complex.t
old new 21 21 .include 'fp_equality.pasm' 22 22 .include "iglobals.pasm" 23 23 24 plan(4 67)24 plan(483) 25 25 26 26 string_parsing() 27 27 exception_malformed_string__real_part() … … 77 77 sech_of_complex_numbers() 78 78 csch_of_complex_numbers() 79 79 add_using_subclass_of_complex_bug_59630() 80 test_nan() 80 81 81 82 # END_OF_TESTS 82 83 … … 730 731 731 732 #XXX: can't do $P1.'$S2'() 732 733 $P2 = $P1. $S2() 733 $S 3= sprintf "%f%+fi", $P2734 $S4 = sprintf "%f%+fi", $P2 734 735 735 736 concat $S5, $S2, " of " 736 737 concat $S5, $S5, $S4 738 concat $S5, $S5, " " 739 concat $S5, $S5, $S3 737 740 738 $I0 = cmp_str $S1, $S 3741 $I0 = cmp_str $S1, $S4 739 742 $I0 = not $I0 740 743 741 todo( $I0, $S 4)744 todo( $I0, $S5 ) 742 745 .endm 743 746 744 747 .sub ln_of_complex_numbers … … 1173 1176 .complex_op_is("23i", "0.272549+0.040301i", 'csch' ) 1174 1177 .end 1175 1178 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.0NaNi" 1226 is(c2, 'NaN', "parse 1.0NaNi") 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 1176 1246 .sub add_using_subclass_of_complex_bug_59630 1177 1247 skip( 3, 'add using subclass of Complex  RT #59630' ) 1178 1248 .return()