From bug-request at octave dot org Mon Jun 21 16:30:04 2004 Subject: Re: besselj bug From: David Bateman To: william turney Cc: bug-octave at bevo dot che dot wisc dot edu Date: Mon, 21 Jun 2004 23:29:12 +0200 According to william turney (on 06/18/04): > > To: bug-octave at bevo dot che dot wisc dot edu > Cc: turney > Subject: besselj generates incorrect values for negative integer indicies when the argument is very small. > -------- > Bug report for Octave 2.1.50 configured for i686-pc-linux-gnu > > Description: > ----------- > > It is well known that besselj(-n,x)=(-1)^n*besselj(n,x) for n an integer. > This relationship hold in octave until x is sufficiently small. For X=0, > and negative n, octave outputs a value of infinity (should be zero). > > Although I did not investigate the Hankle function, I suspect it has the > same issue. Hum the problem is that if negative values of n, besselj is implemented as besselj(-n,x) = cos(n*pi)*besselj(n,x) - sin(n*pi)*bessely(n,x) As x->0 bessely(n,x) goes to infinity and causes a divergent solution to the above as sin(n*pi) doesn't quite equal zero. This problem exists for all integer negative n values to besselj. A similar problem exists for bessely where bessely(-n,x) = cos(n*pi)*bessely(n,x) + sin(n*pi)*besselj(n,x) but in this case the problem is for negative n values where n+0.5 is an integer. The only way I see around the problem is special casing this situation. Give me a few minutes, as I'm working on a patch now, but 1) want to check that I got the implementation right, and 2) want to get all of the places where this occurs.. Cheers David > > Repeat-By: > --------- > > Simply calling the function with the above configuration consistently reproduces the problem. > This problem also occurs on another i686 based pc running the same version of octave. > > Fix: > --- > > * If possible, replace this item with a description of how to > fix the problem (if you don't have a fix for the problem, don't > include this section, but please do submit your report anyway). > > > > Configuration (please do not edit this section): > ----------------------------------------------- > > uname output: Linux gretna 2.6.5-1.358 #1 Sat May 8 09:04:50 EDT 2004 i686 i686 i386 GNU/Linux > configure opts: '--enable-dl' '--enable-shared=yes' '--enable-rpath' '--enable-lite-kernel' '--enable-picky-flags' '--enable-static=no' '--with-g77' '--prefix=/usr' '--infodir=/usr/share/info' 'CFLAGS=-fPIC' > Fortran compiler: g77 > FFLAGS: -O > F2C: > F2CFLAGS: > FLIBS: -L/usr/lib/gcc-lib/i386-redhat-linux/3.3.3 -L/usr/lib/gcc-lib/i386-redhat-linux/3.3.3/../../.. -lfrtbegin -lg2c -lm -lgcc_s > CPPFLAGS: > INCFLAGS: -I. -I. -I./liboctave -I./src -I./libcruft/misc -I./glob -I./glob > C compiler: gcc, version 3.3.3 20040311 (Red Hat Linux 3.3.3-3) > CFLAGS: -fPIC > CPICFLAG: -fPIC > C++ compiler: g++, version 3.3.3 > CXXFLAGS: -g -O2 -Wall -Wcast-align -Wcast-qual -Wid-clash-31 -Winline -Wmissing-prototypes -Wnested-externs -Wpointer-arith -Wstrict-prototypes -Wwrite-strings -Weffc++ -fno-nonnull-objects > CXXPICFLAG: -fPIC > LD_CXX: g++ > LDFLAGS: > LIBFLAGS: -L. > RLD_FLAG: -Wl,-rpath -Wl,/usr/lib/octave-2.1.50 > BLAS_LIBS: -llapack -lblas > FFTW_LIBS: > LIBS: -lreadline -lncurses -ldl -lm > LEXLIB: > LIBPLPLOT: > LIBDLFCN: > LIBGLOB: ./glob/glob.o ./glob/fnmatch.o > SED: /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 -DSEPCHAR=1 > -DSEPCHAR_STR=":" -D__NO_MATH_INLINES=1 -DCXX_NEW_FRIEND_TEMPLATE_DECL=1 > -DCXX_ISO_COMPLIANT_LIBRARY=1 -DCXX_ABI=gnu_v3 -DHAVE_LIBM=1 > -DF77_FUNC(name,NAME)=name ## _ -DF77_FUNC_(name,NAME)=name ## __ > -DHAVE_BLAS=1 -DHAVE_GETHOSTNAME=1 -DHAVE_GETPWNAM=1 -DHAVE_DEV_T=1 > -DHAVE_INO_T=1 -DHAVE_NLINK_T=1 -DHAVE_NLINK_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_H=1 -DHAVE_ALLOCA=1 -DNPOS=std::string::npos > -DHAVE_PLACEMENT_DELETE=1 -DHAVE_DYNAMIC_AUTO_ARRAYS=1 -DSTDC_HEADERS=1 > -DHAVE_DIRENT_H=1 -DTIME_WITH_SYS_TIME=1 -DHAVE_SYS_WAIT_H=1 > -DHAVE_ASSERT_H=1 -DHAVE_CURSES_H=1 -DHAVE_DLFCN_H=1 -DHAVE_FCNTL_H=1 > -DHAVE_FLOAT_H=1 -DHAVE_GRP_H=1 -DHAVE_LIMITS_H=1 -DHAVE_MEMORY_H=1 > -DHAVE_NCURSES_H=1 -DHAVE_POLL_H=1 -DHAVE_PWD_H=1 -DHAVE_STDLIB_H=1 > -DHAVE_STRING_H=1 -DHAVE_SYS_IOCTL_H=1 -DHAVE_SYS_PARAM_H=1 > -DHAVE_SYS_POLL_H=1 -DHAVE_SYS_RESOURCE_H=1 -DHAVE_SYS_SELECT_H=1 > -DHAVE_SYS_STAT_H=1 -DHAVE_SYS_TIME_H=1 -DHAVE_SYS_TIMES_H=1 > -DHAVE_SYS_TYPES_H=1 -DHAVE_SYS_UTSNAME_H=1 -DHAVE_TERMCAP_H=1 > -DHAVE_UNISTD_H=1 -DHAVE_SSTREAM=1 -DHAVE_TERMIO_H=1 -DHAVE_SGTTY_H=1 > -DHAVE_GLOB_H=1 -DHAVE_FNMATCH_H=1 -DHAVE_ATEXIT=1 -DHAVE_BASENAME=1 > -DHAVE_BCOPY=1 -DHAVE_BZERO=1 -DHAVE_DUP2=1 -DHAVE_ENDGRENT=1 > -DHAVE_ENDPWENT=1 -DHAVE_EXECVP=1 -DHAVE_FCNTL=1 -DHAVE_FORK=1 > -DHAVE_GETCWD=1 -DHAVE_GETEGID=1 -DHAVE_GETEUID=1 -DHAVE_GETGID=1 > -DHAVE_GETGRENT=1 -DHAVE_GETGRGID=1 -DHAVE_GETGRNAM=1 -DHAVE_GETPGRP=1 > -DHAVE_GETPID=1 -DHAVE_GETPPID=1 -DHAVE_GETPWENT=1 -DHAVE_GETPWUID=1 > -DHAVE_GETTIMEOFDAY=1 -DHAVE_GETUID=1 -DHAVE_GETWD=1 -DHAVE_KILL=1 > -DHAVE_LINK=1 -DHAVE_LOCALTIME_R=1 -DHAVE_LSTAT=1 -DHAVE_MEMMOVE=1 > -DHAVE_MKDIR=1 -DHAVE_MKFIFO=1 -DHAVE_MKSTEMP=1 -DHAVE_ON_EXIT=1 > -DHAVE_PIPE=1 -DHAVE_POLL=1 -DHAVE_PUTENV=1 -DHAVE_RAISE=1 > -DHAVE_READLINK=1 -DHAVE_RENAME=1 -DHAVE_RINDEX=1 -DHAVE_RMDIR=1 > -DHAVE_SELECT=1 -DHAVE_SETGRENT=1 -DHAVE_SETPWENT=1 -DHAVE_SETVBUF=1 > -DHAVE_SIGACTION=1 -DHAVE_SIGLONGJMP=1 -DHAVE_SIGPENDING=1 > -DHAVE_SIGPROCMASK=1 -DHAVE_SIGSUSPEND=1 -DHAVE_STAT=1 -DHAVE_STRCASECMP=1 > -DHAVE_STRDUP=1 -DHAVE_STRERROR=1 -DHAVE_STRFTIME=1 -DHAVE_STRNCASECMP=1 > -DHAVE_STRPTIME=1 -DHAVE_SYMLINK=1 -DHAVE_TEMPNAM=1 -DHAVE_UMASK=1 > -DHAVE_UNLINK=1 -DHAVE_USLEEP=1 -DHAVE_VFPRINTF=1 -DHAVE_VSPRINTF=1 > -DHAVE_VSNPRINTF=1 -DHAVE_WAITPID=1 -DHAVE_LIBDL=1 -DHAVE_DLOPEN=1 > -DHAVE_DLSYM=1 -DHAVE_DLERROR=1 -DHAVE_DLCLOSE=1 -DHAVE_DLOPEN_API=1 > -DENABLE_DYNAMIC_LINKING=1 -DHAVE_TIMEVAL=1 -DHAVE_FINITE=1 -DHAVE_ISNAN=1 > -DHAVE_ISINF=1 -DHAVE_COPYSIGN=1 -DHAVE_ACOSH=1 -DHAVE_ASINH=1 > -DHAVE_ATANH=1 -DHAVE_ERF=1 -DHAVE_ERFC=1 -DHAVE_STRUCT_STAT_ST_BLKSIZE=1 > -DHAVE_STRUCT_STAT_ST_BLOCKS=1 -DHAVE_STRUCT_STAT_ST_RDEV=1 > -DHAVE_STRUCT_TM_TM_ZONE=1 -DHAVE_TM_ZONE=1 -DUSE_READLINE=1 > -DEXCEPTION_IN_MATH=1 -DRETSIGTYPE=void -DSYS_SIGLIST_DECLARED=1 > -DHAVE_SYS_SIGLIST=1 -DHAVE_POSIX_SIGNALS=1 -DHAVE_GETRUSAGE=1 > -DHAVE_TIMES=1 -DYYTEXT_POINTER=1 -DGNUPLOT_BINARY="gnuplot" > -DGNUPLOT_HAS_FRAMES=1 > > User-preferences (please do not edit this section): > -------------------------------------------------- > > EDITOR = "emacs" > EXEC_PATH = ":/usr/kerberos/bin:/usr/local/bin:/bin:/usr/bin:/usr/X11R6/bin:/home/turney/bin:/sbin:/usr/sbin:/usr/user/bin" > IMAGEPATH = ".:/usr/share/octave/2.1.50/imagelib//" > INFO_FILE = "/usr/share/info/octave.info" > INFO_PROGRAM = "info" > LOADPATH = ":" > PAGER = "less -e -P'-- less ?pB(%pB\\%):--. (f)orward, (b)ack, (q)uit$'" > PS1 = "\\s:\\#> " > PS2 = "> " > PS4 = "+ " > automatic_replot = 0 > beep_on_error = 0 > completion_append_char = " " > default_eval_print_flag = 1 > # default_global_variable_value = > default_return_value = [] > default_save_format = "ascii" > define_all_return_values = 0 > do_fortran_indexing = 0 > echo_executing_commands = 0 > empty_list_elements_ok = 1 > fixed_point_format = 0 > gnuplot_binary = "gnuplot" > gnuplot_command_end = "\n" > gnuplot_command_plot = "pl" > gnuplot_command_replot = "rep" > gnuplot_command_splot = "sp" > gnuplot_command_title = "t" > gnuplot_command_using = "u" > gnuplot_command_with = "w" > gnuplot_has_frames = 1 > history_file = "/home/turney/.octave_hist" > history_size = 1024 > ignore_function_time_stamp = "system" > implicit_num_to_str_ok = 0 > implicit_str_to_num_ok = 0 > initialize_global_variables = 0 > max_recursion_depth = 256 > ok_to_lose_imaginary_part = "warn" > output_max_field_width = 10 > output_precision = 5 > page_output_immediately = 0 > page_screen_output = 1 > prefer_column_vectors = 1 > print_answer_id_name = 1 > print_empty_dimensions = 1 > print_rhs_assign_val = 0 > propagate_empty_matrices = 1 > resize_on_range_error = 1 > return_last_computed_value = 0 > save_precision = 15 > saving_history = 1 > silent_functions = 0 > split_long_rows = 1 > string_fill_char = " " > struct_levels_to_print = 2 > suppress_verbose_help_message = 0 > treat_neg_dim_as_zero = 0 > warn_assign_as_truth_value = 1 > warn_divide_by_zero = 1 > warn_function_name_clash = 1 > warn_future_time_stamp = 1 > warn_missing_semicolon = 0 > warn_variable_switch_label = 0 > whitespace_in_literal_matrix = -- David Bateman David dot Bateman at motorola dot com Motorola CRM +33 1 69 35 48 04 (Ph) Parc Les Algorithmes, Commune de St Aubin +33 1 69 35 77 01 (Fax) 91193 Gif-Sur-Yvette FRANCE The information contained in this communication has been classified as: [x] General Business Information [ ] Motorola Internal Use Only [ ] Motorola Confidential Proprietary ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web: http://www.octave.org How to fund new projects: http://www.octave.org/funding.html Subscription information: http://www.octave.org/archive.html -------------------------------------------------------------