Index: parrot-svn/src/pmc/complex.pmc =================================================================== --- parrot-svn.orig/src/pmc/complex.pmc +++ parrot-svn/src/pmc/complex.pmc @@ -77,6 +77,12 @@ complex_parse_string(PARROT_INTERP, FLOA 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,6 +134,12 @@ complex_parse_string(PARROT_INTERP, FLOA 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,7 +448,11 @@ Returns true if the complex number is no 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,6 +743,7 @@ Adds C to the complex number, pla 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,6 +760,7 @@ Adds C to the complex number, pla 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,6 +834,7 @@ Subtracts C from the complex numb 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,6 +851,7 @@ Subtracts C from the complex numb 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,7 +1019,7 @@ Throws divide by zero exception if divis =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,7 +1028,7 @@ http://docs.sun.com/source/806-3568/ncg_ */ 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,17 +1039,28 @@ http://docs.sun.com/source/806-3568/ncg_ 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,12 +1127,21 @@ http://docs.sun.com/source/806-3568/ncg_ 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,11 +1259,12 @@ Sets C to the absolute value of SE */ /* +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,7 +1275,14 @@ Sets C to the absolute value of SE 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,7 +1295,14 @@ Sets C to the absolute value of SE 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); } Index: parrot-svn/t/pmc/complex.t =================================================================== --- parrot-svn.orig/t/pmc/complex.t +++ parrot-svn/t/pmc/complex.t @@ -21,7 +21,7 @@ Tests the Complex PMC. .include 'fp_equality.pasm' .include "iglobals.pasm" - plan(467) + plan(483) string_parsing() exception_malformed_string__real_part() @@ -77,6 +77,7 @@ Tests the Complex PMC. sech_of_complex_numbers() csch_of_complex_numbers() add_using_subclass_of_complex_bug_59630() + test_nan() # END_OF_TESTS @@ -730,15 +731,17 @@ handler: #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,6 +1176,73 @@ todo: .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()