Bug - function roots numerical problem

[Marcelo-60]

Please see attached file.

To: bug@...
Cc: Marcelo
Subject: roots function
--------
Bug report for Octave 3.0.1 configured for i486-pc-linux-gnu

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

  * The roots function have numerical problems.

Repeat-By:
---------

  * For example: Find the root of (x+1)^n = binomial expansion.
    The answer should be x = -1 with multiplicity n.
    e.g., c = [1 4 6 4 1]; roots(c) should result -1. -1. -1. -1.
    The octave doesn't give the exact answer. It returns:
    -1.00006 + 0.00006i   -1.00006 - 0.00006i
    -0.99994 + 0.00006i   -0.99994 - 0.00006i
    The SciLab and Matlab do the job perfectly.

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 TEPLOTAXL 2.6.28-15-generic #49-Ubuntu SMP Tue Aug 18 18:40:08 UTC 2009 i686 GNU/Linux
configure opts:   '--host=i486-linux-gnu' '--build=i486-linux-gnu' '--prefix=/usr' '--datadir=/usr/share' '--libdir=/usr/lib' '--libexecdir=/usr/lib' '--infodir=/usr/share/info' '--mandir=/usr/share/man' '--with-blas=-lblas-3gf' '--with-lapack=-llapackgf-3' '--with-hdf5' '--with-fftw' '--enable-shared' '--enable-rpath' '--disable-static' 'build_alias=i486-linux-gnu' 'host_alias=i486-linux-gnu' 'CC=gcc' 'CFLAGS=-g -O2' 'LDFLAGS=-Wl,-Bsymbolic-functions' 'CPPFLAGS=' 'CXX=g++' 'CXXFLAGS=-g -O2' 'F77=gfortran' 'FFLAGS=-g -O2'
Fortran compiler: gfortran
FFLAGS:           -O2 -g
F2C:              @F2C@
F2CFLAGS:         @F2CFLAGS@
FLIBS:            -L/usr/lib/gcc/i486-linux-gnu/4.3.3 -L/usr/lib/gcc/i486-linux-gnu/4.3.3/../../../../lib -L/lib/../lib -L/usr/lib/../lib -L/usr/lib/gcc/i486-linux-gnu/4.3.3/../../.. -lhdf5 -lz -lgfortranbegin -lgfortran -lm
CPPFLAGS:        
INCFLAGS:         -I. -I. -I./liboctave -I./src -I./libcruft/misc
C compiler:       gcc, version 4.3.3 (Ubuntu 4.3.3-5ubuntu4)
CFLAGS:           -O2 -g
CPICFLAG:         -fPIC
C++ compiler:     g++, version 4.3.3
CXXFLAGS:         -O2 -g
CXXPICFLAG:       -fPIC
LD_CXX:           g++
LDFLAGS:          
LIBFLAGS:         -L.
RLD_FLAG:         -Wl,-rpath -Wl,/usr/lib/octave-3.0.1
BLAS_LIBS:        -llapackgf-3 -lblas-3gf
FFTW_LIBS:        -lfftw3
LIBS:             -lreadline  -lncurses -ldl -lhdf5 -lz -lm
LEXLIB:          
LIBGLOB:          
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 -D__EXTENSIONS__=1 -D_ALL_SOURCE=1 -D_GNU_SOURCE=1
  -D_POSIX_PTHREAD_SEMANTICS=1 -D_TANDEM_SOURCE=1 -D__EXTENSIONS__=1
  -D_ALL_SOURCE=1 -D_GNU_SOURCE=1 -D_POSIX_PTHREAD_SEMANTICS=1
  -D_TANDEM_SOURCE=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 -DCXX_ABI=unknown -DHAVE_LIBM=1
  -DHAVE_QHULL=1 -DHAVE_PCRE=1 -DHAVE_REGEXEC=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_H=1 -DHAVE_GLPK=1
  -DHAVE_CURL_CURL_H=1 -DHAVE_CURL=1 -DHAVE_IEEE754_DATA_FORMAT=1
  -DF77_FUNC(name,NAME)=name ## _ -DF77_FUNC_(name,NAME)=name ## _
  -DHAVE_BLAS=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_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_INTTYPES_H=1 -DHAVE_LIMITS_H=1
  -DHAVE_LOCALE_H=1 -DHAVE_MEMORY_H=1 -DHAVE_NCURSES_H=1 -DHAVE_POLL_H=1
  -DHAVE_PWD_H=1 -DHAVE_STDINT_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_UTIME_H=1 -DHAVE_SSTREAM=1 -DHAVE_TERMIO_H=1 -DHAVE_SGTTY_H=1
  -DHAVE_GLOB_H=1 -DHAVE_FNMATCH_H=1 -DHAVE_FNMATCH=1 -DHAVE_GLOB=1
  -DHAVE_ATEXIT=1 -DHAVE_BASENAME=1 -DHAVE_BCOPY=1 -DHAVE_BZERO=1
  -DHAVE_CANONICALIZE_FILE_NAME=1 -DHAVE_CHMOD=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_LGAMMA=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_REALPATH=1
  -DHAVE_RENAME=1 -DHAVE_RINDEX=1 -DHAVE_RMDIR=1 -DHAVE_ROUND=1
  -DHAVE_SELECT=1 -DHAVE_SETGRENT=1 -DHAVE_SETLOCALE=1
  -DHAVE_SETPWENT=1 -DHAVE_SETVBUF=1 -DHAVE_SIGACTION=1
  -DHAVE_SIGLONGJMP=1 -DHAVE_SIGPENDING=1 -DHAVE_SIGPROCMASK=1
  -DHAVE_SIGSUSPEND=1 -DHAVE_SNPRINTF=1 -DHAVE_STAT=1 -DHAVE_STRCASECMP=1
  -DHAVE_STRDUP=1 -DHAVE_STRERROR=1 -DHAVE_STRNCASECMP=1 -DHAVE_STRPTIME=1
  -DHAVE_STRSIGNAL=1 -DHAVE_SYMLINK=1 -DHAVE_TEMPNAM=1 -DHAVE_TGAMMA=1
  -DHAVE_UMASK=1 -DHAVE_UNAME=1 -DHAVE_UNLINK=1 -DHAVE_USLEEP=1
  -DHAVE_UTIME=1 -DHAVE_VFPRINTF=1 -DHAVE_VSPRINTF=1 -DHAVE_VSNPRINTF=1
  -DHAVE_WAITPID=1 -DHAVE_STRFTIME=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_DECL_SIGNBIT=1
  -DHAVE_ACOSH=1 -DHAVE_ASINH=1 -DHAVE_ATANH=1 -DHAVE_ERF=1 -DHAVE_ERFC=1
  -DHAVE_EXP2=1 -DHAVE_LOG2=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 -DHAVE_DECL_SYS_SIGLIST=1
  -DHAVE_POSIX_SIGNALS=1 -DRETSIGTYPE_IS_VOID=1 -DHAVE_GETRUSAGE=1
  -DHAVE_TIMES=1 -DYYTEXT_POINTER=1 -DGNUPLOT_BINARY="gnuplot"

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

  EDITOR = emacs
  EXEC_PATH = /usr/lib/octave/3.0.1/site/exec/i486-pc-linux-gnu:/usr/lib/octave/api-v32/site/exec/i486-pc-linux-gnu:/usr/lib/octave/site/exec/i486-pc-linux-gnu:/usr/lib/octave/3.0.1/exec/i486-pc-linux-gnu:/usr/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games
  IMAGE_PATH = .:/usr/share/octave/3.0.1/imagelib
  PAGER = pager
  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 = gnuplot
  gnuplot_command_end =

  gnuplot_command_plot = pl
  gnuplot_command_replot = rep
  gnuplot_command_splot = sp
  gnuplot_command_title = t
  gnuplot_command_using = u
  gnuplot_command_with = w
  history_file = /home/vesper/.octave_hist
  history_size = 1024
  ignore_function_time_stamp = system
  info_file = /usr/share/info/octave3.0.info
  info_program = info
  makeinfo_program = makeinfo
  max_recursion_depth = 256
  output_max_field_width = 5
  output_precision = 5
  page_output_immediately = 0
  page_screen_output = 1
  print_answer_id_name = 1
  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

_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave

[Dupuis]

Marcelo-60 wrote:

Please see attached file.

To: [e-mail]
Cc: Marcelo
Subject: roots function
--------
Bug report for Octave 3.0.1 configured for i486-pc-linux-gnu

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

  * The roots function have numerical problems.

Repeat-By:
---------

  * For example: Find the root of (x+1)^n = binomial expansion.
    The answer should be x = -1 with multiplicity n.
    e.g., c = [1 4 6 4 1]; roots(c) should result -1. -1. -1. -1.
    The octave doesn't give the exact answer. It returns:
    -1.00006 + 0.00006i   -1.00006 - 0.00006i
    -0.99994 + 0.00006i   -0.99994 - 0.00006i
    The SciLab and Matlab do the job perfectly.

It's not a bug, it's a feature. Roots extracting is performed through matrix computation, and in the case of a root with a multiplicity greater than 1, the matrix becomes singular. The fix is to implement an algorithm which computes the GCD between the polynomial and its first derivative: multiple roots are common factor of both. The division of the original polynomial by this GCD is a polynom where each root is of multiplicity 1, and this doesn't suffer for numerical inaccuracies.

I suggest you to have a look at Z. Zeng's homepage, http://www.neiu.edu/~zzeng/, then have a look at 'multroot', which implements this solution. It's Matlab code, easy to translate to Octave, but the license is very restricitive. See if you can cope with it.

Regards

Pascal

[Jaroslav Hajek-2]

On Sat, Sep 19, 2009 at 5:04 AM, Marcelo <marventur@...> wrote:
> Please see attached file.
>
> _______________________________________________
> Bug-octave mailing list
> Bug-octave@...
> https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave
>
>

The algorithm used by roots is fairly simplistic - it simply computes
the eigenvalues of the companion matrix. That is good enough for most
purposes, but is not suited for multiple roots. Contributing a better
algorithm (probably hybrid) would be welcome.

--
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave

[John W. Eaton-3]

On 21-Sep-2009, Jaroslav Hajek wrote:

| On Sat, Sep 19, 2009 at 5:04 AM, Marcelo <marventur@...> wrote:
| > Please see attached file.
| >
| > _______________________________________________
| > Bug-octave mailing list
| > Bug-octave@...
| > https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave
| >
| >
|
| The algorithm used by roots is fairly simplistic - it simply computes
| the eigenvalues of the companion matrix. That is good enough for most
| purposes, but is not suited for multiple roots. Contributing a better
| algorithm (probably hybrid) would be welcome.

If you do contribute an improved function, please remember that it
must not be derived from proprietary software.  For example, if
there happens to be a roots.m file in Matlab that you know about, you
must not use that as a basis for a function that will become a part of
Octave.

jwe
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave

[John W. Eaton-3]

On 21-Sep-2009, John W. Eaton wrote:

| If you do contribute an improved function, please remember that it
| must not be derived from proprietary software.  For example, if
| there happens to be a roots.m file in Matlab that you know about, you
| must not use that as a basis for a function that will become a part of
| Octave.

I should have also mentioned that some care is also needed even when
borrowing code from other free software.  For example, I recall that
someone once translated code from an RLaB function so that it would
run in Octave.  The resulting function was distributed as part of a
collection of Octave functions (this was long before Octave had the
concept of packages).  All of this should have been OK as RLaB was
distributed under the terms of the GPL as was the Octave package.
Unfortunatley, the author of the function later found that the code in
RLaB had simply been translated from Matlab to RLaB syntax by copying
it from the corresponding Matlab function.  In the end, the code had
to be removed from the collection of Octave functions.  It's important
that we avoid such problems in the future.

jwe
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/bug-octave