bugGNU Octave - Bugs: bug #29432, Sin(1e22) is wrong

 
 

bug #29432: Sin(1e22) is wrong

Submitted by:  Judd Storrs <judd>
Submitted on:  Sat 03 Apr 2010 04:37:14 PM UTC  
 
Category: NoneSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Incorrect Result
Status: InvalidAssigned to: None
Originator Name: Derek O'ConnorOpen/Closed: Closed
Release: 3.2.4Operating System: Microsoft Windows

Add a New Comment(Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

(Jump to the original submission Jump to the original submission)

Wed 07 Apr 2010 10:08:18 PM UTC, comment #8:

I think the originator choose 1e22 based on knowledge of floating point representation and possibly with knowledge of trig-function algorithm implementations--other than that I don't know why getting large sin()s correct would matter? Frankly, in terms of correctness, any number between -1 and +1 is plausible as the result of sin(1e22) with ieee doubles since:

I just assume this is a special diagnostic feature. Or I guess it could be an algorithmic fingerprint. I just thought I'd send in the test in the chance it was useful. It's not so that's fine.

Judd Storrs <judd>
Wed 07 Apr 2010 08:36:53 PM UTC, comment #7:

Should we really add this one test? I mean, do you really want to start adding tests for all C library functions in the Octave test suite? I suppose that would be OK, but maybe we should use a more systematic approach and create a lot of tests for library functions and clearly mark them as tests for system library functions, so any failures there should be reported to whoever is maintaining the library function, not to us. But even then, the results of those tests are only going to be seen by the people who build Octave and run the tests. The Octave users will still see the errors, without any clue that the problem is in a library function.

Or, how about just reporting this particular bug?

Also, 1e22 seems like a rather large argument for sin. Does this really come up a lot? I'm just curious about this point. I'm not trying to say that it isn't important that the function is correct, just wondering wtf you are doing if you are computing the sine of 1e22...

Also, you can simply write

%!assert (...)

instead of

%!test
%! assert (...)

John W. Eaton <jwe>
Project Administrator
Wed 07 Apr 2010 08:03:09 PM UTC, comment #6:

Here's a changeset against tip to add a test. Maybe it will pique interest in solving this problem on systems that might be affected.

(file #20152)

Judd Storrs <judd>
Wed 07 Apr 2010 05:58:22 PM UTC, comment #5:

I'm changing the status to "invalid" because this is not a bug in Octave, but a problem with a system library. It is outside the scope of the Octave project to try to fix problems that occur in system library functions.

If this is a problem in the mingw math library, then the problem should be fixed there.

Another possible place for a fix is gnulib, which does have a math module, but I'm not sure whether it is intended to fix problems like the one reported here.

John W. Eaton <jwe>
Project Administrator
Sun 04 Apr 2010 01:20:44 AM UTC, comment #4:

I'm pretty sure it's a i386-pc-mingw32 problem.

R ver 2.11.0 alpha under mingw32 has the same problem.

R x64 2.11.0 (dev.) under ming64 does not have the problem.

Octave 3.2.4 under "i686-pc-cygwin" does not have the problem.

Derek O'Connor

Derek O'Connor <derekroconnor>
Sun 04 Apr 2010 12:23:21 AM UTC, comment #3:

A recent hg-tip gets the correct answer on a second Linux machine. Is this bug still valid? If yes, is it really an Octave bug or a mingw bug?

Rik <rik5>
Project Administrator
Sat 03 Apr 2010 04:50:31 PM UTC, comment #2:

Sorry, I goofed. I can NOT reproduce this on linux. I don't know what I did the first time.

$ octave-3.3.51+ --quiet --eval 'sin(1e22)'
ans = -0.85220

$ octave-3.2.3 --quiet --eval 'sin(1e22)'
ans = -0.85220

$ uname -a
Linux cuneus 2.6.31-20-generic #58-Ubuntu SMP Fri Mar 12 04:38:19 UTC 2010 x86_64 GNU/Linux

Judd Storrs <judd>
Sat 03 Apr 2010 04:43:58 PM UTC, comment #1:

I can confirm this problem on 64-bit linux using a recently compiled hg-tip. Mathematica, matlab and libm agree that octave's result is incorrect:

< M A T L A B > Version 7.5.0.338 (R2007b)

>> sin(1e22)

ans =
-0.8522

Mathematica 7.0 for Linux x86 (64-bit)
In[1]:= N[Sin[10^22],20]
Out[1]= -0.85220084976718880177

A simple C program using libm:

#include <math.h>
#include <stdio.h>

void
main ( void )
{
double x = 1e22 ;
printf("%22.16g\n", sin(x)) ;
}

$ c99 test.c -lm
$ ./a.out
-0.8522008497671888

Judd Storrs <judd>
Sat 03 Apr 2010 04:37:14 PM UTC, original submission:

Bug report for Octave 3.2.4 configured for i686-pc-mingw32

Description:
-----------

octave-3.2.4.exe:9> sin(1e22)
ans =3D 4.12143367107085e-001

The correct answer is sin(10^22) =3D=20
-8.52200849767188801772705893753029368261762150410044e-001

Repeat-By:
---------
=20
sin(10^22). =20

Please note that 10^22 is exactly representable as an IEEE double =
precision floating point number, whose binary representation is

10^22 =3D 1000011110000110011110000011001001101110101011001001 x 2^22

Hence fl(10^22) =3D 10^22, whereas fl(10^23) ~=3D 10^23,=20

Fix:
---

CRlibm, an efficient and proven correctly-rounded mathematical library

http://lipforge.ens-lyon.fr/www/crlibm/

Configuration (please do not edit this section):
-----------------------------------------------

uname output: Windows
configure opts:
Fortran compiler: mingw32-Fortran-4.4.0-dw2
FFLAGS: -O -mieee-fp
FLIBS: -lgfortran
CPPFLAGS: -I/usr/local/octmgw32_gcc-4.4.0-dw2/include
INCFLAGS: -I. -I/octmgw32/octave/octave-3.2.4 -I. -I./liboctave -I./src -I./libcruft/misc -I/octmgw32/octave/octave-3.2.4 -I/octmgw32/octave/octave-3.2.4/liboctave -I/octmgw32/octave/octave-3.2.4/src -I/octmgw32/octave/octave-3.2.4/libcruft/misc
C compiler: mingw32-gcc-4.4.0-dw2, version4.4.0 (GCC)
CFLAGS: -march=i686 -mtune=generic -O3 -Wall
CPICFLAG:
C++ compiler: mingw32-g++-4.4.0-dw2, version4.4.0
CXXFLAGS: -D_GLIBCXX_DLL -march=i686 -mtune=generic -O3 -Wall
CXXPICFLAG:
LD_CXX: mingw32-g++-4.4.0-dw2
LDFLAGS: -shared-libgcc -Wl,--allow-multiple-definition
LIBFLAGS: -L.
RLD_FLAG:
BLAS_LIBS: -llapack -lblas
FFTW_LIBS: -lfftw3 -lfftw3f
LIBS: -lreadline -lncurses -liberty -lhdf5 -lz -lm -lgdi32 -lws2_32 -luser32 -lkernel32
LEXLIB:
LIBGLOB: -lglob
SED: /usr/bin/sed
DEFS:
-DPACKAGE_NAME="" -DPACKAGE_TARNAME="" -DPACKAGE_VERSION="" -DPACKAGE_STRING="" -DPACKAGE_BUGREPORT="" -DOCTAVE_SOURCE=1 -DSTDC_HEADERS=1 -DHAVE_SYS_TYPES_H=1 -DHAVE_SYS_STAT_H=1 -DHAVE_STDLIB_H=1 -DHAVE_STRING_H=1 -DHAVE_MEMORY_H=1 -DHAVE_STRINGS_H=1 -DHAVE_INTTYPES_H=1 -DHAVE_STDINT_H=1 -DHAVE_UNISTD_H=1 -D__EXTENSIONS__=1 -D_ALL_SOURCE=1 -D_GNU_SOURCE=1 -D_POSIX_PTHREAD_SEMANTICS=1 -D_TANDEM_SOURCE=1 -DSEPCHAR=';' -DSEPCHAR_STR=";" -D__NO_MATH_INLINES=1 -DCXX_NEW_FRIEND_TEMPLATE_DECL=1 -DCXX_ISO_COMPLIANT_LIBRARY=1 -DHAVE_PTHREAD=1 -DHAVE_X_WINDOWS=1 -DHAVE_LIBM=1 -DHAVE_QHULL=1 -DHAVE_PCRE_COMPILE=1 -DHAVE_PCRE=1 -DHAVE_REGEX=1 -DHAVE_ZLIB_H=1 -DHAVE_ZLIB=1 -DHAVE_HDF5_H=1 -DHAVE_HDF5=1 -DHAVE_H5GGET_NUM_OBJS=1 -DHAVE_FFTW3=1 -DHAVE_GLPK_GLPK_H=1 -DHAVE_GLPK=1 -DHAVE_CURL_CURL_H=1 -DHAVE_CURL=1 -DHAVE_MAGICK=1 -DHAVE_WINDOWS_H=1 -DHAVE_GL_GL_H=1 -DHAVE_GL_GLU_H=1 -DHAVE_OPENGL=1 -DHAVE_FTGL_FTGL_H=1 -DHAVE_FTGL=1 -DHAVE_FLTK=1 -DHAVE_IEEE754_DATA_FORMAT=1 -DF77_FUNC(name,NAME)=name ## -DF77_FUNC(name,NAME)=name ## _ -DHAVE_BLAS=1 -DHAVE_QRUPDATE=1 -DHAVE_SUITESPARSE_AMD_H=1 -DHAVE_AMD=1 -DHAVE_SUITESPARSE_UMFPACK_H=1 -DHAVE_UMFPACK=1 -DUMFPACK_SEPARATE_SPLIT=1 -DHAVE_SUITESPARSE_COLAMD_H=1 -DHAVE_COLAMD=1 -DHAVE_SUITESPARSE_CCOLAMD_H=1 -DHAVE_CCOLAMD=1 -DHAVE_SUITESPARSE_CHOLMOD_H=1 -DHAVE_CHOLMOD=1 -DHAVE_SUITESPARSE_CS_H=1 -DHAVE_CXSPARSE=1 -DHAVE_ARPACK=1 -Duid_t=int -Dgid_t=int -DHAVE_DEV_T=1 -DHAVE_INO_T=1 -DHAVE_LONG_LONG_INT=1 -DHAVE_UNSIGNED_LONG_LONG_INT=1 -DHAVE_SIGSET_T=1 -DHAVE_SIG_ATOMIC_T=1 -DSIZEOF_SHORT=2 -DSIZEOF_INT=4 -DSIZEOF_LONG=4 -DSIZEOF_LONG_LONG=8 -DHAVE_ALLOCA=1 -DHAVE_PLACEMENT_DELETE=1 -DHAVE_DYNAMIC_AUTO_ARRAYS=1 -DHAVE_FAST_INT_OPS=1 -DSIZEOF_LONG_DOUBLE=12 -DSTDC_HEADERS=1 -DHAVE_DIRENT_H=1 -DTIME_WITH_SYS_TIME=1 -DHAVE_ASSERT_H=1 -DHAVE_CURSES_H=1 -DHAVE_DIRECT_H=1 -DHAVE_FCNTL_H=1 -DHAVE_FLOAT_H=1 -DHAVE_INTTYPES_H=1 -DHAVE_LIMITS_H=1 -DHAVE_LOCALE_H=1 -DHAVE_MEMORY_H=1 -DHAVE_NCURSES_H=1 -DHAVE_PTHREAD_H=1 -DHAVE_STDINT_H=1 -DHAVE_STDLIB_H=1 -DHAVE_STRING_H=1 -DHAVE_SYS_PARAM_H=1 -DHAVE_SYS_STAT_H=1 -DHAVE_SYS_TIME_H=1 -DHAVE_SYS_TYPES_H=1 -DHAVE_SYS_UTIME_H=1 -DHAVE_TERMCAP_H=1 -DHAVE_UNISTD_H=1 -DHAVE_UTIME_H=1 -DHAVE_SSTREAM=1 -DHAVE_GLOB_H=1 -DHAVE_FNMATCH_H=1 -DHAVE_CONIO_H=1 -DHAVE_ATEXIT=1 -DHAVE_BASENAME=1 -DHAVE_CHMOD=1 -DHAVE_DUP2=1 -DHAVE_EXECVP=1 -DHAVE_EXPM1=1 -DHAVE_EXPM1F=1 -DHAVE_FSTAT=1 -DHAVE_GETCWD=1 -DHAVE_GETPID=1 -DHAVE_GETTIMEOFDAY=1 -DHAVE__KBHIT=1 -DHAVE_LGAMMA=1 -DHAVE_LGAMMAF=1 -DHAVE_LOG1P=1 -DHAVE_LOG1PF=1 -DHAVE_MEMMOVE=1 -DHAVE_MKDIR=1 -DHAVE_PUTENV=1 -DHAVE_RAISE=1 -DHAVE_RENAME=1 -DHAVE_RMDIR=1 -DHAVE_ROUNDL=1 -DHAVE_SETLOCALE=1 -DHAVE_SETVBUF=1 -DHAVE_SNPRINTF=1 -DHAVE_STAT=1 -DHAVE_STRCASECMP=1 -DHAVE_STRDUP=1 -DHAVE_STRERROR=1 -DHAVE_STRICMP=1 -DHAVE_STRNCASECMP=1 -DHAVE_STRNICMP=1 -DHAVE_TEMPNAM=1 -DHAVE_TGAMMAF=1 -DHAVE_TRUNC=1 -DHAVE_UMASK=1 -DHAVE_UNLINK=1 -DHAVE_USLEEP=1 -DHAVE_UTIME=1 -DHAVE_VFPRINTF=1 -DHAVE_VSPRINTF=1 -DHAVE_VSNPRINTF=1 -DHAVE__CHMOD=1 -DHAVE__SNPRINTF=1 -DHAVE_DECL_EXP2=1 -DHAVE_DECL_ROUND=1 -DHAVE_DECL_TGAMMA=1 -DHAVE_EXP2=1 -DHAVE_ROUND=1 -DHAVE_TGAMMA=1 -DHAVE_MKSTEMPS=1 -DHAVE_C99_VSNPRINTF=1 -DOCTAVE_HAVE_BROKEN_STRPTIME=1 -D_WIN32_WINNT=0x0403 -DHAVE_LOADLIBRARY_API=1 -DENABLE_DYNAMIC_LINKING=1 -DHAVE_TIMEVAL=1 -DHAVE_CMATH_ISNAN=1 -DHAVE_CMATH_ISNANF=1 -DHAVE_CMATH_ISINF=1 -DHAVE_CMATH_ISINFF=1 -DHAVE_CMATH_ISFINITE=1 -DHAVE_CMATH_ISFINITEF=1 -DHAVE_FINITE=1 -DHAVE_ISNAN=1 -DHAVE_COPYSIGN=1 -DHAVE_SIGNBIT=1 -DHAVE__FINITE=1 -DHAVE__ISNAN=1 -DHAVE__COPYSIGN=1 -DHAVE_DECL_SIGNBIT=1 -DHAVE_ACOSH=1 -DHAVE_ACOSHF=1 -DHAVE_ASINH=1 -DHAVE_ASINHF=1 -DHAVE_ATANH=1 -DHAVE_ATANHF=1 -DHAVE_ERF=1 -DHAVE_ERFF=1 -DHAVE_ERFC=1 -DHAVE_ERFCF=1 -DHAVE_EXP2F=1 -DHAVE_LOG2=1 -DHAVE_LOG2F=1 -DHAVE_HYPOTF=1 -DHAVE_STRUCT_STAT_ST_RDEV=1 -DHAVE_DECL_TZNAME=1 -DHAVE_TZNAME=1 -DMKDIR_TAKES_ONE_ARG=1 -DUSE_READLINE=1 -DRETSIGTYPE=void -DHAVE_DECL_SYS_SIGLIST=0 -DMUST_REINSTALL_SIGHANDLERS=1 -DRETSIGTYPE_IS_VOID=1

User-preferences (please do not edit this section):

EDITOR = D:\MathSoftware\Octave\3.2.4_gcc-4.4.0\tools\notepad++\notepad++.exe
EXEC_PATH = D:\MathSoftware\Octave\3.2.4_gcc-4.4.0\MINGW32\bin;D:\MathSoftware\Octave\3.2.4_gcc-4.4.0\MSYS\bin;D:\MathSoftware\Octave\3.2.4_gcc-4.4.0\libexec\octave\3.2.4\site\exec\i686-pc-mingw32;D:\MathSoftware\Octave\3.2.4_gcc-4.4.0\libexec\octave\api-v37\site\exec\i686-pc-mingw32;D:\MathSoftware\Octave\3.2.4_gcc-4.4.0\libexec\octave\site\exec\i686-pc-mingw32;D:\MathSoftware\Octave\3.2.4_gcc-4.4.0\libexec\octave\3.2.4\exec\i686-pc-mingw32;D:\MathSoftware\Octave\3.2.4_gcc-4.4.0\bin;D:\MathSoftware\MSF\Plugins;D:\Mendeley Desktop;D:\Converters\ImageMagick-6.5.8-Q16;C:\Program Files (x86)\Intel\Compiler\11.1\048\lib\Intel64;C:\Program Files (x86)\Intel\Compiler\11.1\048\lib\ia32;C:\Program Files (x86)\Intel\Compiler\11.1\048\mkl\em64t\bin;C:\Windows\system32;C:\Windows;C:\Windows\System32\Wbem;C:\Windows\System32\WindowsPowerShell\v1.0\;C:\Program Files (x86)\Common Files\DivX Shared\;C:\Program Files (x86)\Microsoft SQL Server\90\Tools\binn\;D:\Matlab R2008a\bin;D:\Matlab R2008a\bin\win64;D:\Matlab-R2007b\bin;D:\Converters\ImageConverter Plus;D:\Converters\ImageConverter Plus\Microsoft.VC80.CRT;D:\Converters\ImageConverter Plus\Microsoft.VC80.MFC;D:\Compilers\Lua5.1;D:\Compilers\Lua5.1\clibs;D:\Utilities\Media\QuickTime\QTSystem\;D:\MiKTeX2.8W7\miktex\bin;D:\Converters\ImageConverter Plus;D:\Converters\ImageConverter Plus\Microsoft.VC80.CRT;D:\Converters\ImageConverter Plus\Microsoft.VC80.MFC;D:\MathSoftware\FreeMat\bin
IMAGE_PATH = .;D:\MathSoftware\Octave\3.2.4_gcc-4.4.0\share\octave\3.2.4\imagelib
PAGER = D:\MathSoftware\Octave\3.2.4_gcc-4.4.0\bin\less.exe
PS1 = \s:\#>
PS2 = >
PS4 = +
beep_on_error = 0
completion_append_char =
crash_dumps_octave_core = 1
echo_executing_commands = 0
fixed_point_format = 0
gnuplot_binary = D:\MathSoftware\Octave\3.2.4_gcc-4.4.0\bin\gnuplot.exe
# gnuplot_command_end = <no value or error in displaying it>
# gnuplot_command_plot = <no value or error in displaying it>
# gnuplot_command_replot = <no value or error in displaying it>
# gnuplot_command_splot = <no value or error in displaying it>
# gnuplot_command_title = <no value or error in displaying it>
# gnuplot_command_using = <no value or error in displaying it>
# gnuplot_command_with = <no value or error in displaying it>
history_file = C:\Users\Derek O'Connor\.octave_hist
history_size = 1024
ignore_function_time_stamp = system
info_file = D:\MathSoftware\Octave\3.2.4_gcc-4.4.0\share\info\octave.info
info_program = info
makeinfo_program = makeinfo
max_recursion_depth = 256
output_max_field_width = 15
output_precision = 15
page_output_immediately = 0
page_screen_output = 1
# print_answer_id_name = <no value or error in displaying it>
print_empty_dimensions = 1
save_precision = 16
saving_history = 1
sighup_dumps_octave_core = 1
sigterm_dumps_octave_core = 1
silent_functions = 0
split_long_rows = 1
string_fill_char =
struct_levels_to_print = 2
suppress_verbose_help_message = 0

Judd Storrs <judd>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach File(s):
   
   
Comment:
   

Attached Files
file #20152:  bug29432.patch added by judd (1003B - text/x-diff - hg changeset against tip)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -unavailable- added by jwe (Posted a comment)
  • -unavailable- added by derekroconnor (Posted a comment)
  • -unavailable- added by rik5 (Posted a comment)
  • -unavailable- added by judd (Submitted the item)
  • -unavailable- added by judd
  •  

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 5 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Wed 07 Apr 2010 08:03:09 PM UTCjuddAttached File-=>Added bug29432.patch, #20152
    Wed 07 Apr 2010 05:58:22 PM UTCjweStatusWorks For Me=>Invalid
    Sun 04 Apr 2010 05:58:03 PM UTCrik5Open/ClosedOpen=>Closed
    Sun 04 Apr 2010 12:23:21 AM UTCrik5StatusNone=>Works For Me
    Sat 03 Apr 2010 04:37:14 PM UTCjuddCarbon-Copy-=>Added derek o'connor

    Back to the top


    Powered by Savane 3.1-cleanup1