* bug #13469: Scilab contained non free code 67/18567/5
Clément DAVID [Mon, 3 Oct 2016 17:06:18 +0000 (19:06 +0200)]
fsultra removed as this code license is "commercial use restriction".
Its non-free (Violate DFSG 6)

rpoly.f removed as this code license is "commercial use restriction".
The BSD rpoly_plus_plus is used instead.

Change-Id: I6212105c7014dc3590ea044d6e40ef3bfcadeed9

41 files changed:
scilab/ACKNOWLEDGEMENTS
scilab/CHANGES.md
scilab/modules/development_tools/macros/testexamples.sce
scilab/modules/polynomials/Makefile.am
scilab/modules/polynomials/Makefile.in
scilab/modules/polynomials/help/en_US/roots.xml
scilab/modules/polynomials/help/fr_FR/roots.xml
scilab/modules/polynomials/help/ja_JP/roots.xml
scilab/modules/polynomials/help/pt_BR/roots.xml
scilab/modules/polynomials/license.txt
scilab/modules/polynomials/src/c/polynomials.vcxproj
scilab/modules/polynomials/src/c/polynomials.vcxproj.filters
scilab/modules/polynomials/src/cpp/find_polynomial_roots_jenkins_traub.cc [new file with mode: 0644]
scilab/modules/polynomials/src/cpp/find_polynomial_roots_jenkins_traub.h [new file with mode: 0644]
scilab/modules/polynomials/src/cpp/polynomial.cc [new file with mode: 0644]
scilab/modules/polynomials/src/cpp/polynomial.h [new file with mode: 0644]
scilab/modules/polynomials/src/cpp/rpoly.cpp [new file with mode: 0644]
scilab/modules/polynomials/src/fortran/polynomials_f.vfproj
scilab/modules/polynomials/src/fortran/polynomials_f2c.vcxproj
scilab/modules/polynomials/src/fortran/polynomials_f2c.vcxproj.filters
scilab/modules/polynomials/src/fortran/rpoly.f [deleted file]
scilab/modules/polynomials/tests/unit_tests/roots.dia.ref
scilab/modules/polynomials/tests/unit_tests/roots.tst
scilab/modules/polynomials/tests/unit_tests/roots_difficult.tst [new file with mode: 0644]
scilab/modules/randlib/Makefile.am
scilab/modules/randlib/Makefile.in
scilab/modules/randlib/help/en_US/grand.xml
scilab/modules/randlib/help/fr_FR/grand.xml
scilab/modules/randlib/help/ja_JP/grand.xml
scilab/modules/randlib/help/ru_RU/grand.xml
scilab/modules/randlib/includes/others_generators.h
scilab/modules/randlib/license.txt
scilab/modules/randlib/sci_gateway/cpp/sci_grand.cpp
scilab/modules/randlib/src/c/fsultra.c [deleted file]
scilab/modules/randlib/src/c/randlib.vcxproj
scilab/modules/randlib/src/c/randlib.vcxproj.filters
scilab/modules/randlib/src/cpp/grand.cpp
scilab/modules/randlib/tests/unit_tests/grand_generators.dia.ref
scilab/modules/randlib/tests/unit_tests/grand_generators.tst
scilab/modules/randlib/tests/unit_tests/grand_hypermat.dia.ref
scilab/modules/randlib/tests/unit_tests/grand_hypermat.tst

index bf70e7c..8a508ef 100644 (file)
@@ -201,9 +201,6 @@ calelm: low level routines (INRIA).
 
 control: LINPACK + EISPACK + INRIA routines.
          dsubsp and exchnqz: Paul van Dooren.
-         rpoly: copyrighted by the ACM (alg. 493), which grants
-                 general permission to distribute provided 
-                 the copies are not made for direct commercial advantage. 
          lybsc, lydsr, lybad,sydsr and sybad are adapted from SLICE 
                 (M. Denham).
          sszer: Emami-naeini, A. and van Dooren, P. (Automatica paper).
index 30e0b68..2d7aa83 100644 (file)
@@ -197,6 +197,8 @@ will zoom all axes in the current figure.
 * The graphics entity "Datatip" has a new property `detached_position` which accepts `[]`
 or a 3-components vector to set the position in axes coordinates to draw the datatip text box.
 * `MPI_Create_comm` create a new communicator from MPI_COMM_WORLD using MPI world ranks.
+* The `grand` non-free `fsultra` generator was removed.
+* The original `rpoly` algorithm was removed in favor of a C++11 implementation
 
 Help pages:
 -----------
@@ -327,6 +329,7 @@ Bug Fixes
 * [#11959](http://bugzilla.scilab.org/show_bug.cgi?id=11959): Allow "Zoom Area" to be clicked out of axes
 * [#12110](http://bugzilla.scilab.org/show_bug.cgi?id=12110): Unable to zoom multiple axes at once
 * [#13166](http://bugzilla.scilab.org/show_bug.cgi?id=13166): `l` and `b` endian flags used with `mget` and `mgeti` were sticky
+* [#13469](http://bugzilla.scilab.org/show_bug.cgi?id=13469): Scilab contained non free code
 * [#13470](http://bugzilla.scilab.org/show_bug.cgi?id=13470): `histplot(0,0,%t)` crashed
 * [#13583](http://bugzilla.scilab.org/show_bug.cgi?id=13583): `getd` loading a script including a `clear` instruction yielded an error
 * [#13597](http://bugzilla.scilab.org/show_bug.cgi?id=13597): `help format` claimed setting a number of digits instead of characters
index 2ff8cbb..070ab4a 100644 (file)
@@ -51,7 +51,6 @@ function reinit_for_test()
     grand("setgen","kiss");grand("setsd",362436069,521288629,123456789,380116160);
     grand("setgen","clcg2");grand("setsd",1234567890,123456789);
     grand("setgen","urand");grand("setsd",0);
-    grand("setgen","fsultra");grand("setsd",1234567,7654321);
     grand("setgen","mt");grand("setsd",5489);
     rand("seed",0);
     format("v",10);
index e844973..b4582fc 100644 (file)
@@ -36,7 +36,6 @@ POLYNOMIALS_FORTRAN_SOURCES = \
     src/fortran/fxshfr.f \
     src/fortran/mpdiag.f \
     src/fortran/dmpcle.f \
-    src/fortran/rpoly.f \
     src/fortran/wpodiv.f \
     src/fortran/wdmpmu.f \
     src/fortran/wmptra.f \
@@ -61,6 +60,12 @@ POLYNOMIALS_FORTRAN_SOURCES = \
     src/fortran/writebufsfact.f \
     src/fortran/chkvar.f
 
+POLYNOMIALS_CXX_SOURCES = \
+    src/cpp/rpoly.cpp \
+    src/cpp/find_polynomial_roots_jenkins_traub.cc \
+    src/cpp/polynomial.cc
+
+
 GATEWAY_CXX_SOURCES = \
     sci_gateway/cpp/polynomials_gw.cpp \
     sci_gateway/cpp/sci_poly.cpp \
@@ -93,6 +98,7 @@ libscipolynomials_la_CPPFLAGS = \
     -I$(top_srcdir)/modules/string/includes \
     -I$(top_srcdir)/modules/console/includes \
     -I$(top_srcdir)/modules/elementary_functions/includes/ \
+    $(EIGEN_CPPFLAGS) \
     $(AM_CPPFLAGS)
 
 if MAINTAINER_MODE
@@ -104,7 +110,7 @@ endif
 
 
 
-libscipolynomials_algo_la_SOURCES = $(POLYNOMIALS_FORTRAN_SOURCES)
+libscipolynomials_algo_la_SOURCES = $(POLYNOMIALS_FORTRAN_SOURCES) $(POLYNOMIALS_CXX_SOURCES)
 libscipolynomials_la_SOURCES = $(GATEWAY_C_SOURCES) $(GATEWAY_CXX_SOURCES)
 libscipolynomials_algo_la_CPPFLAGS = $(libscipolynomials_la_CPPFLAGS)
 
index 6e7248d..617c285 100644 (file)
@@ -191,19 +191,23 @@ am__objects_1 = src/fortran/residu.lo src/fortran/imptra.lo \
        src/fortran/mptri.lo src/fortran/horner.lo \
        src/fortran/idegre.lo src/fortran/fxshfr.lo \
        src/fortran/mpdiag.lo src/fortran/dmpcle.lo \
-       src/fortran/rpoly.lo src/fortran/wpodiv.lo \
-       src/fortran/wdmpmu.lo src/fortran/wmptra.lo \
-       src/fortran/wmpins.lo src/fortran/mpdegr.lo \
-       src/fortran/ddmpev.lo src/fortran/realit.lo \
-       src/fortran/dpmul.lo src/fortran/sfact2.lo \
-       src/fortran/dmpad.lo src/fortran/nextk.lo \
-       src/fortran/blktit.lo src/fortran/dimin.lo \
-       src/fortran/newest.lo src/fortran/dwmpmu.lo \
-       src/fortran/impcnc.lo src/fortran/wdmpad.lo \
-       src/fortran/wmpcle.lo src/fortran/quadit.lo \
-       src/fortran/quad.lo src/fortran/calcsc.lo \
-       src/fortran/writebufsfact.lo src/fortran/chkvar.lo
-am_libscipolynomials_algo_la_OBJECTS = $(am__objects_1)
+       src/fortran/wpodiv.lo src/fortran/wdmpmu.lo \
+       src/fortran/wmptra.lo src/fortran/wmpins.lo \
+       src/fortran/mpdegr.lo src/fortran/ddmpev.lo \
+       src/fortran/realit.lo src/fortran/dpmul.lo \
+       src/fortran/sfact2.lo src/fortran/dmpad.lo \
+       src/fortran/nextk.lo src/fortran/blktit.lo \
+       src/fortran/dimin.lo src/fortran/newest.lo \
+       src/fortran/dwmpmu.lo src/fortran/impcnc.lo \
+       src/fortran/wdmpad.lo src/fortran/wmpcle.lo \
+       src/fortran/quadit.lo src/fortran/quad.lo \
+       src/fortran/calcsc.lo src/fortran/writebufsfact.lo \
+       src/fortran/chkvar.lo
+am__objects_2 = src/cpp/libscipolynomials_algo_la-rpoly.lo \
+       src/cpp/libscipolynomials_algo_la-find_polynomial_roots_jenkins_traub.lo \
+       src/cpp/libscipolynomials_algo_la-polynomial.lo
+am_libscipolynomials_algo_la_OBJECTS = $(am__objects_1) \
+       $(am__objects_2)
 libscipolynomials_algo_la_OBJECTS =  \
        $(am_libscipolynomials_algo_la_OBJECTS)
 AM_V_lt = $(am__v_lt_@AM_V@)
@@ -213,7 +217,7 @@ am__v_lt_1 =
 @MAINTAINER_MODE_FALSE@am_libscipolynomials_algo_la_rpath =
 @MAINTAINER_MODE_TRUE@am_libscipolynomials_algo_la_rpath =
 libscipolynomials_la_DEPENDENCIES = libscipolynomials-algo.la
-am__objects_2 =  \
+am__objects_3 =  \
        sci_gateway/cpp/libscipolynomials_la-polynomials_gw.lo \
        sci_gateway/cpp/libscipolynomials_la-sci_poly.lo \
        sci_gateway/cpp/libscipolynomials_la-sci_varn.lo \
@@ -225,7 +229,7 @@ am__objects_2 =  \
        sci_gateway/cpp/libscipolynomials_la-sci_simp.lo \
        sci_gateway/cpp/libscipolynomials_la-sci_sfact.lo \
        sci_gateway/cpp/libscipolynomials_la-sci_bezout.lo
-am_libscipolynomials_la_OBJECTS = $(am__objects_2)
+am_libscipolynomials_la_OBJECTS = $(am__objects_3)
 libscipolynomials_la_OBJECTS = $(am_libscipolynomials_la_OBJECTS)
 @MAINTAINER_MODE_FALSE@am_libscipolynomials_la_rpath =
 @MAINTAINER_MODE_TRUE@am_libscipolynomials_la_rpath = -rpath \
@@ -653,7 +657,6 @@ POLYNOMIALS_FORTRAN_SOURCES = \
     src/fortran/fxshfr.f \
     src/fortran/mpdiag.f \
     src/fortran/dmpcle.f \
-    src/fortran/rpoly.f \
     src/fortran/wpodiv.f \
     src/fortran/wdmpmu.f \
     src/fortran/wmptra.f \
@@ -678,6 +681,11 @@ POLYNOMIALS_FORTRAN_SOURCES = \
     src/fortran/writebufsfact.f \
     src/fortran/chkvar.f
 
+POLYNOMIALS_CXX_SOURCES = \
+    src/cpp/rpoly.cpp \
+    src/cpp/find_polynomial_roots_jenkins_traub.cc \
+    src/cpp/polynomial.cc
+
 GATEWAY_CXX_SOURCES = \
     sci_gateway/cpp/polynomials_gw.cpp \
     sci_gateway/cpp/sci_poly.cpp \
@@ -709,12 +717,13 @@ libscipolynomials_la_CPPFLAGS = \
     -I$(top_srcdir)/modules/string/includes \
     -I$(top_srcdir)/modules/console/includes \
     -I$(top_srcdir)/modules/elementary_functions/includes/ \
+    $(EIGEN_CPPFLAGS) \
     $(AM_CPPFLAGS)
 
 @MAINTAINER_MODE_TRUE@pkglib_LTLIBRARIES = libscipolynomials.la
 @MAINTAINER_MODE_FALSE@noinst_LTLIBRARIES = libscipolynomials-algo.la libscipolynomials.la
 @MAINTAINER_MODE_TRUE@noinst_LTLIBRARIES = libscipolynomials-algo.la
-libscipolynomials_algo_la_SOURCES = $(POLYNOMIALS_FORTRAN_SOURCES)
+libscipolynomials_algo_la_SOURCES = $(POLYNOMIALS_FORTRAN_SOURCES) $(POLYNOMIALS_CXX_SOURCES)
 libscipolynomials_la_SOURCES = $(GATEWAY_C_SOURCES) $(GATEWAY_CXX_SOURCES)
 libscipolynomials_algo_la_CPPFLAGS = $(libscipolynomials_la_CPPFLAGS)
 
@@ -813,7 +822,7 @@ HELP_CHAPTERLANG = en_US fr_FR pt_BR
 all: all-am
 
 .SUFFIXES:
-.SUFFIXES: .sci .bin .cpp .f .lo .o .obj
+.SUFFIXES: .sci .bin .cc .cpp .f .lo .o .obj
 $(srcdir)/Makefile.in: @MAINTAINER_MODE_TRUE@ $(srcdir)/Makefile.am $(top_srcdir)/Makefile.incl.am $(am__configure_deps)
        @for dep in $?; do \
          case '$(am__configure_deps)' in \
@@ -960,8 +969,6 @@ src/fortran/mpdiag.lo: src/fortran/$(am__dirstamp) \
        src/fortran/$(DEPDIR)/$(am__dirstamp)
 src/fortran/dmpcle.lo: src/fortran/$(am__dirstamp) \
        src/fortran/$(DEPDIR)/$(am__dirstamp)
-src/fortran/rpoly.lo: src/fortran/$(am__dirstamp) \
-       src/fortran/$(DEPDIR)/$(am__dirstamp)
 src/fortran/wpodiv.lo: src/fortran/$(am__dirstamp) \
        src/fortran/$(DEPDIR)/$(am__dirstamp)
 src/fortran/wdmpmu.lo: src/fortran/$(am__dirstamp) \
@@ -1008,9 +1015,21 @@ src/fortran/writebufsfact.lo: src/fortran/$(am__dirstamp) \
        src/fortran/$(DEPDIR)/$(am__dirstamp)
 src/fortran/chkvar.lo: src/fortran/$(am__dirstamp) \
        src/fortran/$(DEPDIR)/$(am__dirstamp)
+src/cpp/$(am__dirstamp):
+       @$(MKDIR_P) src/cpp
+       @: > src/cpp/$(am__dirstamp)
+src/cpp/$(DEPDIR)/$(am__dirstamp):
+       @$(MKDIR_P) src/cpp/$(DEPDIR)
+       @: > src/cpp/$(DEPDIR)/$(am__dirstamp)
+src/cpp/libscipolynomials_algo_la-rpoly.lo: src/cpp/$(am__dirstamp) \
+       src/cpp/$(DEPDIR)/$(am__dirstamp)
+src/cpp/libscipolynomials_algo_la-find_polynomial_roots_jenkins_traub.lo:  \
+       src/cpp/$(am__dirstamp) src/cpp/$(DEPDIR)/$(am__dirstamp)
+src/cpp/libscipolynomials_algo_la-polynomial.lo:  \
+       src/cpp/$(am__dirstamp) src/cpp/$(DEPDIR)/$(am__dirstamp)
 
 libscipolynomials-algo.la: $(libscipolynomials_algo_la_OBJECTS) $(libscipolynomials_algo_la_DEPENDENCIES) $(EXTRA_libscipolynomials_algo_la_DEPENDENCIES) 
-       $(AM_V_F77LD)$(F77LINK) $(am_libscipolynomials_algo_la_rpath) $(libscipolynomials_algo_la_OBJECTS) $(libscipolynomials_algo_la_LIBADD) $(LIBS)
+       $(AM_V_CXXLD)$(CXXLINK) $(am_libscipolynomials_algo_la_rpath) $(libscipolynomials_algo_la_OBJECTS) $(libscipolynomials_algo_la_LIBADD) $(LIBS)
 sci_gateway/cpp/$(am__dirstamp):
        @$(MKDIR_P) sci_gateway/cpp
        @: > sci_gateway/cpp/$(am__dirstamp)
@@ -1058,6 +1077,8 @@ mostlyclean-compile:
        -rm -f *.$(OBJEXT)
        -rm -f sci_gateway/cpp/*.$(OBJEXT)
        -rm -f sci_gateway/cpp/*.lo
+       -rm -f src/cpp/*.$(OBJEXT)
+       -rm -f src/cpp/*.lo
        -rm -f src/fortran/*.$(OBJEXT)
        -rm -f src/fortran/*.lo
 
@@ -1075,8 +1096,11 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@sci_gateway/cpp/$(DEPDIR)/libscipolynomials_la-sci_simp.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@sci_gateway/cpp/$(DEPDIR)/libscipolynomials_la-sci_simpMode.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@sci_gateway/cpp/$(DEPDIR)/libscipolynomials_la-sci_varn.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@src/cpp/$(DEPDIR)/libscipolynomials_algo_la-find_polynomial_roots_jenkins_traub.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@src/cpp/$(DEPDIR)/libscipolynomials_algo_la-polynomial.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@src/cpp/$(DEPDIR)/libscipolynomials_algo_la-rpoly.Plo@am__quote@
 
-.cpp.o:
+.cc.o:
 @am__fastdepCXX_TRUE@  $(AM_V_CXX)depbase=`echo $@ | sed 's|[^/]*$$|$(DEPDIR)/&|;s|\.o$$||'`;\
 @am__fastdepCXX_TRUE@  $(CXXCOMPILE) -MT $@ -MD -MP -MF $$depbase.Tpo -c -o $@ $< &&\
 @am__fastdepCXX_TRUE@  $(am__mv) $$depbase.Tpo $$depbase.Po
@@ -1084,7 +1108,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@     DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXXCOMPILE) -c -o $@ $<
 
-.cpp.obj:
+.cc.obj:
 @am__fastdepCXX_TRUE@  $(AM_V_CXX)depbase=`echo $@ | sed 's|[^/]*$$|$(DEPDIR)/&|;s|\.obj$$||'`;\
 @am__fastdepCXX_TRUE@  $(CXXCOMPILE) -MT $@ -MD -MP -MF $$depbase.Tpo -c -o $@ `$(CYGPATH_W) '$<'` &&\
 @am__fastdepCXX_TRUE@  $(am__mv) $$depbase.Tpo $$depbase.Po
@@ -1092,7 +1116,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@     DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXXCOMPILE) -c -o $@ `$(CYGPATH_W) '$<'`
 
-.cpp.lo:
+.cc.lo:
 @am__fastdepCXX_TRUE@  $(AM_V_CXX)depbase=`echo $@ | sed 's|[^/]*$$|$(DEPDIR)/&|;s|\.lo$$||'`;\
 @am__fastdepCXX_TRUE@  $(LTCXXCOMPILE) -MT $@ -MD -MP -MF $$depbase.Tpo -c -o $@ $< &&\
 @am__fastdepCXX_TRUE@  $(am__mv) $$depbase.Tpo $$depbase.Plo
@@ -1100,6 +1124,27 @@ distclean-compile:
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@     DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LTCXXCOMPILE) -c -o $@ $<
 
+src/cpp/libscipolynomials_algo_la-rpoly.lo: src/cpp/rpoly.cpp
+@am__fastdepCXX_TRUE@  $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libscipolynomials_algo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT src/cpp/libscipolynomials_algo_la-rpoly.lo -MD -MP -MF src/cpp/$(DEPDIR)/libscipolynomials_algo_la-rpoly.Tpo -c -o src/cpp/libscipolynomials_algo_la-rpoly.lo `test -f 'src/cpp/rpoly.cpp' || echo '$(srcdir)/'`src/cpp/rpoly.cpp
+@am__fastdepCXX_TRUE@  $(AM_V_at)$(am__mv) src/cpp/$(DEPDIR)/libscipolynomials_algo_la-rpoly.Tpo src/cpp/$(DEPDIR)/libscipolynomials_algo_la-rpoly.Plo
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     $(AM_V_CXX)source='src/cpp/rpoly.cpp' object='src/cpp/libscipolynomials_algo_la-rpoly.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libscipolynomials_algo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o src/cpp/libscipolynomials_algo_la-rpoly.lo `test -f 'src/cpp/rpoly.cpp' || echo '$(srcdir)/'`src/cpp/rpoly.cpp
+
+src/cpp/libscipolynomials_algo_la-find_polynomial_roots_jenkins_traub.lo: src/cpp/find_polynomial_roots_jenkins_traub.cc
+@am__fastdepCXX_TRUE@  $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libscipolynomials_algo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT src/cpp/libscipolynomials_algo_la-find_polynomial_roots_jenkins_traub.lo -MD -MP -MF src/cpp/$(DEPDIR)/libscipolynomials_algo_la-find_polynomial_roots_jenkins_traub.Tpo -c -o src/cpp/libscipolynomials_algo_la-find_polynomial_roots_jenkins_traub.lo `test -f 'src/cpp/find_polynomial_roots_jenkins_traub.cc' || echo '$(srcdir)/'`src/cpp/find_polynomial_roots_jenkins_traub.cc
+@am__fastdepCXX_TRUE@  $(AM_V_at)$(am__mv) src/cpp/$(DEPDIR)/libscipolynomials_algo_la-find_polynomial_roots_jenkins_traub.Tpo src/cpp/$(DEPDIR)/libscipolynomials_algo_la-find_polynomial_roots_jenkins_traub.Plo
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     $(AM_V_CXX)source='src/cpp/find_polynomial_roots_jenkins_traub.cc' object='src/cpp/libscipolynomials_algo_la-find_polynomial_roots_jenkins_traub.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libscipolynomials_algo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o src/cpp/libscipolynomials_algo_la-find_polynomial_roots_jenkins_traub.lo `test -f 'src/cpp/find_polynomial_roots_jenkins_traub.cc' || echo '$(srcdir)/'`src/cpp/find_polynomial_roots_jenkins_traub.cc
+
+src/cpp/libscipolynomials_algo_la-polynomial.lo: src/cpp/polynomial.cc
+@am__fastdepCXX_TRUE@  $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libscipolynomials_algo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT src/cpp/libscipolynomials_algo_la-polynomial.lo -MD -MP -MF src/cpp/$(DEPDIR)/libscipolynomials_algo_la-polynomial.Tpo -c -o src/cpp/libscipolynomials_algo_la-polynomial.lo `test -f 'src/cpp/polynomial.cc' || echo '$(srcdir)/'`src/cpp/polynomial.cc
+@am__fastdepCXX_TRUE@  $(AM_V_at)$(am__mv) src/cpp/$(DEPDIR)/libscipolynomials_algo_la-polynomial.Tpo src/cpp/$(DEPDIR)/libscipolynomials_algo_la-polynomial.Plo
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     $(AM_V_CXX)source='src/cpp/polynomial.cc' object='src/cpp/libscipolynomials_algo_la-polynomial.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libscipolynomials_algo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o src/cpp/libscipolynomials_algo_la-polynomial.lo `test -f 'src/cpp/polynomial.cc' || echo '$(srcdir)/'`src/cpp/polynomial.cc
+
 sci_gateway/cpp/libscipolynomials_la-polynomials_gw.lo: sci_gateway/cpp/polynomials_gw.cpp
 @am__fastdepCXX_TRUE@  $(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libscipolynomials_la_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT sci_gateway/cpp/libscipolynomials_la-polynomials_gw.lo -MD -MP -MF sci_gateway/cpp/$(DEPDIR)/libscipolynomials_la-polynomials_gw.Tpo -c -o sci_gateway/cpp/libscipolynomials_la-polynomials_gw.lo `test -f 'sci_gateway/cpp/polynomials_gw.cpp' || echo '$(srcdir)/'`sci_gateway/cpp/polynomials_gw.cpp
 @am__fastdepCXX_TRUE@  $(AM_V_at)$(am__mv) sci_gateway/cpp/$(DEPDIR)/libscipolynomials_la-polynomials_gw.Tpo sci_gateway/cpp/$(DEPDIR)/libscipolynomials_la-polynomials_gw.Plo
@@ -1177,6 +1222,30 @@ sci_gateway/cpp/libscipolynomials_la-sci_bezout.lo: sci_gateway/cpp/sci_bezout.c
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@     DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libscipolynomials_la_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o sci_gateway/cpp/libscipolynomials_la-sci_bezout.lo `test -f 'sci_gateway/cpp/sci_bezout.cpp' || echo '$(srcdir)/'`sci_gateway/cpp/sci_bezout.cpp
 
+.cpp.o:
+@am__fastdepCXX_TRUE@  $(AM_V_CXX)depbase=`echo $@ | sed 's|[^/]*$$|$(DEPDIR)/&|;s|\.o$$||'`;\
+@am__fastdepCXX_TRUE@  $(CXXCOMPILE) -MT $@ -MD -MP -MF $$depbase.Tpo -c -o $@ $< &&\
+@am__fastdepCXX_TRUE@  $(am__mv) $$depbase.Tpo $$depbase.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     $(AM_V_CXX)source='$<' object='$@' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXXCOMPILE) -c -o $@ $<
+
+.cpp.obj:
+@am__fastdepCXX_TRUE@  $(AM_V_CXX)depbase=`echo $@ | sed 's|[^/]*$$|$(DEPDIR)/&|;s|\.obj$$||'`;\
+@am__fastdepCXX_TRUE@  $(CXXCOMPILE) -MT $@ -MD -MP -MF $$depbase.Tpo -c -o $@ `$(CYGPATH_W) '$<'` &&\
+@am__fastdepCXX_TRUE@  $(am__mv) $$depbase.Tpo $$depbase.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     $(AM_V_CXX)source='$<' object='$@' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXXCOMPILE) -c -o $@ `$(CYGPATH_W) '$<'`
+
+.cpp.lo:
+@am__fastdepCXX_TRUE@  $(AM_V_CXX)depbase=`echo $@ | sed 's|[^/]*$$|$(DEPDIR)/&|;s|\.lo$$||'`;\
+@am__fastdepCXX_TRUE@  $(LTCXXCOMPILE) -MT $@ -MD -MP -MF $$depbase.Tpo -c -o $@ $< &&\
+@am__fastdepCXX_TRUE@  $(am__mv) $$depbase.Tpo $$depbase.Plo
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     $(AM_V_CXX)source='$<' object='$@' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@     DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(LTCXXCOMPILE) -c -o $@ $<
+
 .f.o:
        $(AM_V_F77)$(F77COMPILE) -c -o $@ $<
 
@@ -1192,6 +1261,7 @@ mostlyclean-libtool:
 clean-libtool:
        -rm -rf .libs _libs
        -rm -rf sci_gateway/cpp/.libs sci_gateway/cpp/_libs
+       -rm -rf src/cpp/.libs src/cpp/_libs
        -rm -rf src/fortran/.libs src/fortran/_libs
 install-libscipolynomials_la_etcDATA: $(libscipolynomials_la_etc_DATA)
        @$(NORMAL_INSTALL)
@@ -1375,6 +1445,8 @@ distclean-generic:
        -test . = "$(srcdir)" || test -z "$(CONFIG_CLEAN_VPATH_FILES)" || rm -f $(CONFIG_CLEAN_VPATH_FILES)
        -rm -f sci_gateway/cpp/$(DEPDIR)/$(am__dirstamp)
        -rm -f sci_gateway/cpp/$(am__dirstamp)
+       -rm -f src/cpp/$(DEPDIR)/$(am__dirstamp)
+       -rm -f src/cpp/$(am__dirstamp)
        -rm -f src/fortran/$(DEPDIR)/$(am__dirstamp)
        -rm -f src/fortran/$(am__dirstamp)
 
@@ -1387,7 +1459,7 @@ clean-am: clean-generic clean-libtool clean-local \
        clean-noinstLTLIBRARIES clean-pkglibLTLIBRARIES mostlyclean-am
 
 distclean: distclean-am
-       -rm -rf sci_gateway/cpp/$(DEPDIR)
+       -rm -rf sci_gateway/cpp/$(DEPDIR) src/cpp/$(DEPDIR)
        -rm -f Makefile
 distclean-am: clean-am distclean-compile distclean-generic \
        distclean-local distclean-tags
@@ -1436,7 +1508,7 @@ install-ps-am:
 installcheck-am:
 
 maintainer-clean: maintainer-clean-am
-       -rm -rf sci_gateway/cpp/$(DEPDIR)
+       -rm -rf sci_gateway/cpp/$(DEPDIR) src/cpp/$(DEPDIR)
        -rm -f Makefile
 maintainer-clean-am: distclean-am maintainer-clean-generic
 
index ea6d89e..54336db 100644 (file)
@@ -170,12 +170,4 @@ max(abs(r1-r2))
             ACM TOMS 1, 1 (March 1975), pp. 26-34
         </para>
     </refsection>
-    <refsection>
-        <title>Used Functions</title>
-        <para>The rpoly.f source codes can be found in the directory
-            SCI/modules/polynomials/src/fortran of a Scilab source distribution. In the case where the
-            companion matrix is used, the eigenvalue computation is perfomed using
-            DGEEV and ZGEEV LAPACK codes.
-        </para>
-    </refsection>
 </refentry>
index 167738d..c0ece1d 100644 (file)
@@ -168,13 +168,4 @@ max(abs(r1-r2))
             Polynomial", ACM TOMS Volume 1, Issue 2 (June 1975), pp. 178-189
         </para>
     </refsection>
-    <refsection>
-        <title>Fonctions Utilisées</title>
-        <para>
-            Le code source de rpoly.f peut être trouvé dans le répertoire
-            SCI/modules/polynomials/src/fortran/ de la distribution source de Scilab. Dans le cas où la
-            matrix compagnon est utilisée, le calcul des valeurs propres est effectué
-            en utilisant les routines DGEEV et ZGEEV de LAPACK.
-        </para>
-    </refsection>
 </refentry>
index 27cb8c7..1e986d4 100644 (file)
@@ -87,12 +87,4 @@ spec(A)
             ACM TOMS 1, 1 (March 1975), pp. 26-34
         </para>
     </refsection>
-    <refsection>
-        <title>使用される関数</title>
-        <para>ソースコードrpoly.f はScilabの配布ソースのディレクトリ
-            SCI/modules/polynomials/src/fortran にある.
-            コンパニオン行列が使用される場合,固有値計算は
-            LapackコードであるDGEEVおよびZGEEVにより行われる.
-        </para>
-    </refsection>
 </refentry>
index 07681cc..4ecb47a 100644 (file)
@@ -86,11 +86,4 @@ spec(A)
             26-34
         </para>
     </refsection>
-    <refsection>
-        <title> Funções Utilizadas </title>
-        <para>O código fonte de rpoly.f pode ser achado no diretório
-            SCI/modules/polynomials/src/fortran de uma distribuição fonte do Scilab. A computação de
-            autovalores é feita utilizando-se os códigos do LAPACK DGEEV e ZGEEV.
-        </para>
-    </refsection>
 </refentry>
index d7f8e2f..e26d1ff 100644 (file)
@@ -12,13 +12,6 @@ and continues to be available under such terms.
 For more information, see the COPYING file which you should have received
 along with this program.
 
-Non-free source files:
-======================
-
-With "commercial use restriction"
-
-File: modules/polynomials/src/fortran/rpoly.f
-
 Eispack:
 ========
 
@@ -29,3 +22,16 @@ Copyright:
 
 License:
 Public domain
+
+RPolyPlusPlus:
+==============
+
+Files: src/cpp/find_polynomial_roots_jenkins_traub.cc src/cpp/polynomial.cc
+       src/cpp/find_polynomial_roots_jenkins_traub.h  src/cpp/polynomial.h
+
+Copyright:
+2015 - Chris Sweeney
+
+License:
+BSD
+
index 92ee44c..c2d6c05 100644 (file)
@@ -80,7 +80,7 @@
   <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">
     <ClCompile>
       <Optimization>Disabled</Optimization>
-      <AdditionalIncludeDirectories>../../includes;../../../../libs/intl;../../../localization/includes;../../../core/includes;../../../output_stream/includes;../../../api_scilab/includes;../../../threads/includes;../../../string/includes;../../../console/includes;../../../elementary_functions/includes;../../../dynamic_link/includes;../../../ast/includes/types;../../../ast/includes/ast;../../../ast/includes/analysis;../../../ast/includes/exps;../../../ast/includes/operations;../../../ast/includes/symbol;../../../ast/includes/system_env;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
+      <AdditionalIncludeDirectories>../../src/cpp;../../includes;../../../../libs/Eigen/includes;../../../../libs/intl;../../../localization/includes;../../../core/includes;../../../output_stream/includes;../../../api_scilab/includes;../../../threads/includes;../../../string/includes;../../../console/includes;../../../elementary_functions/includes;../../../dynamic_link/includes;../../../ast/includes/types;../../../ast/includes/ast;../../../ast/includes/analysis;../../../ast/includes/exps;../../../ast/includes/operations;../../../ast/includes/symbol;../../../ast/includes/system_env;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
       <PreprocessorDefinitions>_CRT_SECURE_NO_DEPRECATE;FORDLL;_DEBUG;_WINDOWS;_USRDLL;POLYNOMIALS_EXPORTS;%(PreprocessorDefinitions)</PreprocessorDefinitions>
       <RuntimeLibrary>MultiThreadedDebugDLL</RuntimeLibrary>
       <WarningLevel>Level3</WarningLevel>
@@ -108,7 +108,7 @@ lib /DEF:"$(ProjectDir)elementary_functions_f_Import.def" /SUBSYSTEM:WINDOWS /MA
     </Midl>
     <ClCompile>
       <Optimization>Disabled</Optimization>
-      <AdditionalIncludeDirectories>../../includes;../../../../libs/intl;../../../localization/includes;../../../core/includes;../../../output_stream/includes;../../../api_scilab/includes;../../../threads/includes;../../../string/includes;../../../console/includes;../../../elementary_functions/includes;../../../dynamic_link/includes;../../../ast/includes/types;../../../ast/includes/ast;../../../ast/includes/analysis;../../../ast/includes/exps;../../../ast/includes/operations;../../../ast/includes/symbol;../../../ast/includes/system_env;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
+      <AdditionalIncludeDirectories>../../src/cpp;../../includes;../../../../libs/Eigen/includes;../../../../libs/intl;../../../localization/includes;../../../core/includes;../../../output_stream/includes;../../../api_scilab/includes;../../../threads/includes;../../../string/includes;../../../console/includes;../../../elementary_functions/includes;../../../dynamic_link/includes;../../../ast/includes/types;../../../ast/includes/ast;../../../ast/includes/analysis;../../../ast/includes/exps;../../../ast/includes/operations;../../../ast/includes/symbol;../../../ast/includes/system_env;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
       <PreprocessorDefinitions>_CRT_SECURE_NO_DEPRECATE;FORDLL;_DEBUG;_WINDOWS;_USRDLL;POLYNOMIALS_EXPORTS;%(PreprocessorDefinitions)</PreprocessorDefinitions>
       <RuntimeLibrary>MultiThreadedDebugDLL</RuntimeLibrary>
       <WarningLevel>Level3</WarningLevel>
@@ -134,7 +134,7 @@ lib /DEF:"$(ProjectDir)elementary_functions_f_Import.def" /SUBSYSTEM:WINDOWS /MA
     <ClCompile>
       <FavorSizeOrSpeed>Speed</FavorSizeOrSpeed>
       <WholeProgramOptimization>false</WholeProgramOptimization>
-      <AdditionalIncludeDirectories>../../includes;../../../../libs/intl;../../../localization/includes;../../../core/includes;../../../output_stream/includes;../../../api_scilab/includes;../../../threads/includes;../../../string/includes;../../../console/includes;../../../elementary_functions/includes;../../../dynamic_link/includes;../../../ast/includes/types;../../../ast/includes/ast;../../../ast/includes/analysis;../../../ast/includes/exps;../../../ast/includes/operations;../../../ast/includes/symbol;../../../ast/includes/system_env;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
+      <AdditionalIncludeDirectories>../../src/cpp;../../includes;../../../../libs/Eigen/includes;../../../../libs/intl;../../../localization/includes;../../../core/includes;../../../output_stream/includes;../../../api_scilab/includes;../../../threads/includes;../../../string/includes;../../../console/includes;../../../elementary_functions/includes;../../../dynamic_link/includes;../../../ast/includes/types;../../../ast/includes/ast;../../../ast/includes/analysis;../../../ast/includes/exps;../../../ast/includes/operations;../../../ast/includes/symbol;../../../ast/includes/system_env;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
       <PreprocessorDefinitions>_CRT_SECURE_NO_DEPRECATE;FORDLL;NDEBUG;_WINDOWS;_USRDLL;POLYNOMIALS_EXPORTS;%(PreprocessorDefinitions)</PreprocessorDefinitions>
       <StringPooling>true</StringPooling>
       <RuntimeLibrary>MultiThreadedDLL</RuntimeLibrary>
@@ -167,7 +167,7 @@ lib /DEF:"$(ProjectDir)elementary_functions_f_Import.def" /SUBSYSTEM:WINDOWS /MA
     <ClCompile>
       <FavorSizeOrSpeed>Speed</FavorSizeOrSpeed>
       <WholeProgramOptimization>false</WholeProgramOptimization>
-      <AdditionalIncludeDirectories>../../includes;../../../../libs/intl;../../../localization/includes;../../../core/includes;../../../output_stream/includes;../../../api_scilab/includes;../../../threads/includes;../../../string/includes;../../../console/includes;../../../elementary_functions/includes;../../../dynamic_link/includes;../../../ast/includes/types;../../../ast/includes/ast;../../../ast/includes/analysis;../../../ast/includes/exps;../../../ast/includes/operations;../../../ast/includes/symbol;../../../ast/includes/system_env;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
+      <AdditionalIncludeDirectories>../../src/cpp;../../includes;../../../../libs/Eigen/includes;../../../../libs/intl;../../../localization/includes;../../../core/includes;../../../output_stream/includes;../../../api_scilab/includes;../../../threads/includes;../../../string/includes;../../../console/includes;../../../elementary_functions/includes;../../../dynamic_link/includes;../../../ast/includes/types;../../../ast/includes/ast;../../../ast/includes/analysis;../../../ast/includes/exps;../../../ast/includes/operations;../../../ast/includes/symbol;../../../ast/includes/system_env;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
       <PreprocessorDefinitions>_CRT_SECURE_NO_DEPRECATE;FORDLL;NDEBUG;_WINDOWS;_USRDLL;POLYNOMIALS_EXPORTS;%(PreprocessorDefinitions)</PreprocessorDefinitions>
       <StringPooling>true</StringPooling>
       <RuntimeLibrary>MultiThreadedDLL</RuntimeLibrary>
@@ -205,12 +205,17 @@ lib /DEF:"$(ProjectDir)elementary_functions_f_Import.def" /SUBSYSTEM:WINDOWS /MA
     <ClCompile Include="..\..\sci_gateway\cpp\sci_simp.cpp" />
     <ClCompile Include="..\..\sci_gateway\cpp\sci_simpMode.cpp" />
     <ClCompile Include="..\..\sci_gateway\cpp\sci_varn.cpp" />
+    <ClCompile Include="..\cpp\find_polynomial_roots_jenkins_traub.cc" />
+    <ClCompile Include="..\cpp\polynomial.cc" />
+    <ClCompile Include="..\cpp\rpoly.cpp" />
     <ClCompile Include="DllmainPolynomial.c" />
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="..\..\includes\dynlib_polynomials.h" />
     <ClInclude Include="..\..\includes\dynlib_polynomials_gw.h" />
     <ClInclude Include="..\..\includes\polynomials_gw.hxx" />
+    <ClInclude Include="..\cpp\find_polynomial_roots_jenkins_traub.h" />
+    <ClInclude Include="..\cpp\polynomial.h" />
   </ItemGroup>
   <ItemGroup>
     <None Include="..\..\locales\polynomials.pot" />
index 7efb540..5fa7867 100644 (file)
     <ClCompile Include="..\..\sci_gateway\cpp\sci_simpMode.cpp">
       <Filter>Source Files</Filter>
     </ClCompile>
+    <ClCompile Include="..\cpp\polynomial.cc">
+      <Filter>Source Files</Filter>
+    </ClCompile>
+    <ClCompile Include="..\cpp\rpoly.cpp">
+      <Filter>Source Files</Filter>
+    </ClCompile>
+    <ClCompile Include="..\cpp\find_polynomial_roots_jenkins_traub.cc">
+      <Filter>Source Files</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="..\..\includes\dynlib_polynomials.h">
     <ClInclude Include="..\..\includes\dynlib_polynomials_gw.h">
       <Filter>Header Files</Filter>
     </ClInclude>
+    <ClInclude Include="..\cpp\polynomial.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
+    <ClInclude Include="..\cpp\find_polynomial_roots_jenkins_traub.h">
+      <Filter>Header Files</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <None Include="polynomials_f_Import.def">
diff --git a/scilab/modules/polynomials/src/cpp/find_polynomial_roots_jenkins_traub.cc b/scilab/modules/polynomials/src/cpp/find_polynomial_roots_jenkins_traub.cc
new file mode 100644 (file)
index 0000000..66eff60
--- /dev/null
@@ -0,0 +1,844 @@
+// Copyright (C) 2015 Chris Sweeney. All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are
+// met:
+//
+//     * Redistributions of source code must retain the above copyright
+//       notice, this list of conditions and the following disclaimer.
+//
+//     * Redistributions in binary form must reproduce the above
+//       copyright notice, this list of conditions and the following
+//       disclaimer in the documentation and/or other materials provided
+//       with the distribution.
+//
+//     * Neither the name of Chris Sweeney nor the names of its contributors may
+//       be used to endorse or promote products derived from this software
+//       without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Please contact the author of this library if you have any questions.
+// Author: Chris Sweeney (cmsweeney@cs.ucsb.edu)
+
+#include "find_polynomial_roots_jenkins_traub.h"
+
+#include <Eigen/Core>
+
+#include <cmath>
+#include <complex>
+#include <iostream>
+#include <limits>
+#include <vector>
+
+#include "polynomial.h"
+
+namespace rpoly_plus_plus {
+
+using Eigen::MatrixXd;
+using Eigen::Vector3d;
+using Eigen::VectorXd;
+using Eigen::Vector3cd;
+using Eigen::VectorXcd;
+
+namespace {
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846264338327950288
+#endif
+
+// Machine precision constants.
+static const double mult_eps = std::numeric_limits<double>::epsilon();
+static const double sum_eps = std::numeric_limits<double>::epsilon();
+
+enum ConvergenceType{
+  NO_CONVERGENCE = 0,
+  LINEAR_CONVERGENCE = 1,
+  QUADRATIC_CONVERGENCE = 2
+};
+
+// Solves for the root of the equation ax + b = 0.
+double FindLinearPolynomialRoots(const double a, const double b) {
+  return -b / a;
+}
+
+// Stable quadratic roots according to BKP Horn.
+// http://people.csail.mit.edu/bkph/articles/Quadratics.pdf
+void FindQuadraticPolynomialRoots(const double a,
+                                  const double b,
+                                  const double c,
+                                  std::complex<double>* roots) {
+  const double D = b * b - 4 * a * c;
+  const double sqrt_D = std::sqrt(std::abs(D));
+
+  // Real roots.
+  if (D >= 0) {
+    if (b >= 0) {
+      roots[0] = std::complex<double>((-b - sqrt_D) / (2.0 * a), 0);
+      roots[1] = std::complex<double>((2.0 * c) / (-b - sqrt_D), 0);
+    } else {
+      roots[0] = std::complex<double>((2.0 * c) / (-b + sqrt_D), 0);
+      roots[1] = std::complex<double>((-b + sqrt_D) / (2.0 * a), 0);
+    }
+    return;
+  }
+
+  // Use the normal quadratic formula for the complex case.
+  roots[0] = std::complex<double>(-b / (2.0 * a), sqrt_D / (2.0 * a));
+  roots[1] = std::complex<double>(-b / (2.0 * a), -sqrt_D / (2.0 * a));
+}
+
+// Perform division by a linear term of the form (z - x) and evaluate P at x.
+void SyntheticDivisionAndEvaluate(const VectorXd& polynomial,
+                                  const double x,
+                                  VectorXd* quotient,
+                                  double* eval) {
+  quotient->setZero(polynomial.size() - 1);
+  (*quotient)(0) = polynomial(0);
+  for (int i = 1; i < polynomial.size() - 1; i++) {
+    (*quotient)(i) = polynomial(i) + (*quotient)(i - 1) * x;
+  }
+  const VectorXd::ReverseReturnType& creverse_quotient = quotient->reverse();
+  *eval = polynomial.reverse()(0) + creverse_quotient(0) * x;
+}
+
+// Perform division of a polynomial by a quadratic factor. The quadratic divisor
+// should have leading 1s.
+void QuadraticSyntheticDivision(const VectorXd& polynomial,
+                                const VectorXd& quadratic_divisor,
+                                VectorXd* quotient,
+                                VectorXd* remainder) {
+  quotient->setZero(polynomial.size() - 2);
+  remainder->setZero(2);
+
+  (*quotient)(0) = polynomial(0);
+
+  // If the quotient is a constant then polynomial is degree 2 and the math is
+  // simple.
+  if (quotient->size() == 1) {
+    *remainder =
+        polynomial.tail<2>() - polynomial(0) * quadratic_divisor.tail<2>();
+    return;
+  }
+
+  (*quotient)(1) = polynomial(1) - polynomial(0) * quadratic_divisor(1);
+  for (int i = 2; i < polynomial.size() - 2; i++) {
+    (*quotient)(i) = polynomial(i) - (*quotient)(i - 2) * quadratic_divisor(2) -
+                     (*quotient)(i - 1) * quadratic_divisor(1);
+  }
+
+  const VectorXd::ReverseReturnType &creverse_quotient = quotient->reverse();
+  (*remainder)(0) = polynomial.reverse()(1) -
+                    quadratic_divisor(1) * creverse_quotient(0) -
+                    quadratic_divisor(2) * creverse_quotient(1);
+  (*remainder)(1) =
+      polynomial.reverse()(0) - quadratic_divisor(2) * creverse_quotient(0);
+}
+
+// Determines whether the iteration has converged by examining the three most
+// recent values for convergence.
+template<typename T>
+bool HasConverged(const T& sequence) {
+  const bool convergence_condition_1 =
+      std::abs(sequence(1) - sequence(0)) < std::abs(sequence(0)) / 2.0;
+  const bool convergence_condition_2 =
+      std::abs(sequence(2) - sequence(1)) < std::abs(sequence(1)) / 2.0;
+
+  // If the sequence has converged then return true.
+  return convergence_condition_1 && convergence_condition_2;
+}
+
+// Determines if the root has converged by measuring the relative and absolute
+// change in the root value. This stopping criterion is a simple measurement
+// that proves to work well. It is referred to as "Ward's method" in the
+// following reference:
+//
+// Nikolajsen, Jorgen L. "New stopping criteria for iterative root finding."
+// Royal Society open science (2014)
+template <typename T>
+bool HasRootConverged(const std::vector<T>& roots) {
+  static const double kRootMagnitudeTolerance = 1e-8;
+  static const double kAbsoluteTolerance = 1e-14;
+  static const double kRelativeTolerance = 1e-10;
+
+  if (roots.size() != 3) {
+    return false;
+  }
+
+  const double e_i = std::abs(roots[2] - roots[1]);
+  const double e_i_minus_1 = std::abs(roots[1] - roots[0]);
+  const double mag_root = std::abs(roots[1]);
+  if (e_i <= e_i_minus_1) {
+    if (mag_root < kRootMagnitudeTolerance) {
+      return e_i < kAbsoluteTolerance;
+    } else {
+      return e_i / mag_root <= kRelativeTolerance;
+    }
+  }
+
+  return false;
+}
+
+// Implementation closely follows the three-stage algorithm for finding roots of
+// polynomials with real coefficients as outlined in: "A Three-Stage Algorithm
+// for Real Polynomaials Using Quadratic Iteration" by Jenkins and Traub, SIAM
+// 1970. Please note that this variant is different than the complex-coefficient
+// version, and is estimated to be up to 4 times faster.
+class JenkinsTraubSolver {
+ public:
+  JenkinsTraubSolver(const VectorXd& coeffs,
+                     VectorXd* real_roots,
+                     VectorXd* complex_roots)
+      : polynomial_(coeffs),
+        real_roots_(real_roots),
+        complex_roots_(complex_roots),
+        num_solved_roots_(0) {}
+
+  // Extracts the roots using the Jenkins Traub method.
+  bool ExtractRoots();
+
+ private:
+  // Removes any zero roots and divides polynomial by z.
+  void RemoveZeroRoots();
+
+  // Computes the magnitude of the roots to provide and initial search radius
+  // for the iterative solver.
+  double ComputeRootRadius();
+
+  // Computes the zero-shift applied to the K-Polynomial.
+  void ComputeZeroShiftKPolynomial();
+
+  // Stage 1 of the Jenkins-Traub method. This stage is not technically
+  // necessary, but helps separate roots that are close to zero.
+  void ApplyZeroShiftToKPolynomial(const int num_iterations);
+
+  // Computes and returns the update of sigma(z) based on the current
+  // K-polynomial.
+  //
+  // NOTE: This function is used by the fixed shift iterations (which hold sigma
+  // constant) so sigma is *not* modified internally by this function. If you
+  // want to change sigma, simply call
+  //    sigma = ComputeNextSigma();
+  VectorXd ComputeNextSigma();
+
+  // Updates the K-polynomial based on the current value of sigma for the fixed
+  // or variable shift stage.
+  void UpdateKPolynomialWithQuadraticShift(
+      const VectorXd& polynomial_quotient,
+      const VectorXd& k_polynomial_quotient);
+
+  // Apply fixed-shift iterations to the K-polynomial to separate the
+  // roots. Based on the convergence of the K-polynomial, we apply a
+  // variable-shift linear or quadratic iteration to determine a real root or
+  // complex conjugate pair of roots respectively.
+  ConvergenceType ApplyFixedShiftToKPolynomial(const std::complex<double>& root,
+                                               const int max_iterations);
+
+  // Applies one of the variable shifts to the K-Polynomial. Returns true upon
+  // successful convergence to a good root, and false otherwise.
+  bool ApplyVariableShiftToKPolynomial(
+      const ConvergenceType& fixed_shift_convergence,
+      const std::complex<double>& root);
+
+  // Applies a quadratic shift to the K-polynomial to determine a pair of roots
+  // that are complex conjugates. Return true if a root was successfully found.
+  bool ApplyQuadraticShiftToKPolynomial(const std::complex<double>& root,
+                                        const int max_iterations);
+
+  // Applies a linear shift to the K-polynomial to determine a single real root.
+  // Return true if a root was successfully found.
+  bool ApplyLinearShiftToKPolynomial(const std::complex<double>& root,
+                                     const int max_iterations);
+
+  // Adds the root to the output variables.
+  void AddRootToOutput(const double real, const double imag);
+
+  // Solves polynomials of degree <= 2.
+  bool SolveClosedFormPolynomial();
+
+  // Helper variables to manage the polynomials as they are being manipulated
+  // and deflated.
+  VectorXd polynomial_;
+  VectorXd k_polynomial_;
+  // Sigma is the quadratic factor the divides the K-polynomial.
+  Vector3d sigma_;
+
+  // Let us define a, b, c, and d such that:
+  //   P(z) = Q_P * sigma(z) + b * (z + u) + a
+  //   K(z) = Q_K * sigma(z) + d * (z + u ) + c
+  //
+  // where Q_P and Q_K are the quotients from polynomial division of
+  // sigma(z). Note that this means for a given a root s of sigma:
+  //
+  //   P(s)      = a - b * s_conj
+  //   P(s_conj) = a - b * s
+  //   K(s)      = c - d * s_conj
+  //   K(s_conj) = c - d * s
+  double a_, b_, c_, d_;
+
+  // Output reference variables.
+  VectorXd* real_roots_;
+  VectorXd* complex_roots_;
+  int num_solved_roots_;
+
+  // Keeps track of whether the linear and quadratic shifts have been attempted
+  // yet so that we do not attempt the same shift twice.
+  bool attempted_linear_shift_;
+  bool attempted_quadratic_shift_;
+
+  // Number of zero-shift iterations to perform.
+  static const int kNumZeroShiftIterations = 20;
+
+  // The number of fixed shift iterations is computed as
+  //   # roots found * this multiplier.
+  static const int kFixedShiftIterationMultiplier = 20;
+
+  // If the fixed shift iterations fail to converge, we restart this many times
+  // before considering the solve attempt as a failure.
+  static const int kMaxFixedShiftRestarts = 20;
+
+  // The maximum number of linear shift iterations to perform before considering
+  // the shift as a failure.
+  static const int kMaxLinearShiftIterations = 20;
+
+  // The maximum number of quadratic shift iterations to perform before
+  // considering the shift as a failure.
+  static const int kMaxQuadraticShiftIterations = 20;
+
+  // When quadratic shift iterations are stalling, we attempt a few fixed shift
+  // iterations to help convergence.
+  static const int kInnerFixedShiftIterations = 5;
+
+  // During quadratic iterations, the real values of the root pairs should be
+  // nearly equal since the root pairs are complex conjugates. This tolerance
+  // measures how much the real values may diverge before consider the quadratic
+  // shift to be failed.
+  static const double kRootPairTolerance;;
+};
+
+const double JenkinsTraubSolver::kRootPairTolerance = 0.01;
+
+bool JenkinsTraubSolver::ExtractRoots() {
+  if (polynomial_.size() == 0) {
+    // std::cout << "Invalid polynomial of size 0 passed to "
+    //              "FindPolynomialRootsJenkinsTraub" << std::endl;
+    return false;
+  }
+
+  // Remove any leading zeros of the polynomial.
+  polynomial_ = RemoveLeadingZeros(polynomial_);
+  // Normalize the polynomial.
+  polynomial_ /= polynomial_(0);
+  const int degree = polynomial_.size() - 1;
+
+  // Allocate the output roots.
+  if (real_roots_ != NULL) {
+    real_roots_->setZero(degree);
+  }
+  if (complex_roots_ != NULL) {
+    complex_roots_->setZero(degree);
+  }
+
+  // Remove any zero roots.
+  RemoveZeroRoots();
+
+  // Choose the initial starting value for the root-finding on the complex
+  // plane.
+  const double kDegToRad = M_PI / 180.0;
+  double phi = 49.0 * kDegToRad;
+
+  // Iterate until the polynomial has been completely deflated.
+  for (int i = 0; i < degree; i++) {
+    // Compute the root radius.
+    const double root_radius = ComputeRootRadius();
+
+    // Solve in closed form if the polynomial is small enough.
+    if (polynomial_.size() <= 3) {
+      break;
+    }
+
+    // Stage 1: Apply zero-shifts to the K-polynomial to separate the small
+    // zeros of the polynomial.
+    ApplyZeroShiftToKPolynomial(kNumZeroShiftIterations);
+
+    // Stage 2: Apply fixed shift iterations to the K-polynomial to separate the
+    // roots further.
+    std::complex<double> root;
+    ConvergenceType convergence = NO_CONVERGENCE;
+    for (int j = 0; j < kMaxFixedShiftRestarts; j++) {
+      root = root_radius * std::complex<double>(std::cos(phi), std::sin(phi));
+      convergence = ApplyFixedShiftToKPolynomial(
+          root, kFixedShiftIterationMultiplier * (i + 1));
+      if (convergence != NO_CONVERGENCE) {
+        break;
+      }
+
+      // Rotate the initial root value on the complex plane and try again.
+      phi += 94.0 * kDegToRad;
+    }
+
+    // Stage 3: Find the root(s) with variable shift iterations on the
+    // K-polynomial. If this stage was not successful then we return a failure.
+    if (!ApplyVariableShiftToKPolynomial(convergence, root)) {
+      return false;
+    }
+  }
+  return SolveClosedFormPolynomial();
+}
+
+// Stage 1: Generate K-polynomials with no shifts (i.e. zero-shifts).
+void JenkinsTraubSolver::ApplyZeroShiftToKPolynomial(
+    const int num_iterations) {
+  // K0 is the first order derivative of polynomial.
+  k_polynomial_ = DifferentiatePolynomial(polynomial_) / polynomial_.size();
+  for (int i = 1; i < num_iterations; i++) {
+    ComputeZeroShiftKPolynomial();
+  }
+}
+
+ConvergenceType JenkinsTraubSolver::ApplyFixedShiftToKPolynomial(
+    const std::complex<double>& root, const int max_iterations) {
+  // Compute the fixed-shift quadratic:
+  // sigma(z) = (x - m - n * i) * (x - m + n * i) = x^2 - 2 * m + m^2 + n^2.
+  sigma_(0) = 1.0;
+  sigma_(1) = -2.0 * root.real();
+  sigma_(2) = root.real() * root.real() + root.imag() * root.imag();
+
+  // Compute the quotient and remainder for divinding P by the quadratic
+  // divisor. Since this iteration involves a fixed-shift sigma these may be
+  // computed once prior to any iterations.
+  VectorXd polynomial_quotient, polynomial_remainder;
+  QuadraticSyntheticDivision(
+      polynomial_, sigma_, &polynomial_quotient, &polynomial_remainder);
+
+  // Compute a and b from the above equations.
+  b_ = polynomial_remainder(0);
+  a_ = polynomial_remainder(1) - b_ * sigma_(1);
+
+  // Precompute P(s) for later using the equation above.
+  const std::complex<double> p_at_root = a_ - b_ * std::conj(root);
+
+  // These two containers hold values that we test for convergence such that the
+  // zero index is the convergence value from 2 iterations ago, the first
+  // index is from one iteration ago, and the second index is the current value.
+  Vector3cd t_lambda = Vector3cd::Zero();
+  Vector3d sigma_lambda = Vector3d::Zero();
+  VectorXd k_polynomial_quotient, k_polynomial_remainder;
+  for (int i = 0; i < max_iterations; i++) {
+    k_polynomial_ /= k_polynomial_(0);
+
+    // Divide the shifted polynomial by the quadratic polynomial.
+    QuadraticSyntheticDivision(
+        k_polynomial_, sigma_, &k_polynomial_quotient, &k_polynomial_remainder);
+    d_ = k_polynomial_remainder(0);
+    c_ = k_polynomial_remainder(1) - d_ * sigma_(1);
+
+    // Test for convergence.
+    const VectorXd variable_shift_sigma = ComputeNextSigma();
+    const std::complex<double> k_at_root = c_ - d_ * std::conj(root);
+
+    t_lambda.head<2>() = t_lambda.tail<2>().eval();
+    sigma_lambda.head<2>() = sigma_lambda.tail<2>().eval();
+    t_lambda(2) = root - p_at_root / k_at_root;
+    sigma_lambda(2) = variable_shift_sigma(2);
+
+    // Return with the convergence code if the sequence has converged.
+    if (HasConverged(sigma_lambda)) {
+      return QUADRATIC_CONVERGENCE;
+    } else if (HasConverged(t_lambda)) {
+      return LINEAR_CONVERGENCE;
+    }
+
+    // Compute K_next using the formula above.
+    UpdateKPolynomialWithQuadraticShift(polynomial_quotient,
+                                        k_polynomial_quotient);
+  }
+
+  return NO_CONVERGENCE;
+}
+
+bool JenkinsTraubSolver::ApplyVariableShiftToKPolynomial(
+    const ConvergenceType& fixed_shift_convergence,
+    const std::complex<double>& root) {
+  attempted_linear_shift_ = false;
+  attempted_quadratic_shift_ = false;
+
+  if (fixed_shift_convergence == LINEAR_CONVERGENCE) {
+    return ApplyLinearShiftToKPolynomial(root, kMaxLinearShiftIterations);
+  } else if (fixed_shift_convergence == QUADRATIC_CONVERGENCE) {
+    return ApplyQuadraticShiftToKPolynomial(root, kMaxQuadraticShiftIterations);
+  }
+  return false;
+}
+
+// Generate K-polynomials with variable-shifts. During variable shifts, the
+// quadratic shift is computed as:
+//                | K0(s1)  K0(s2)  z^2 |
+//                | K1(s1)  K1(s2)    z |
+//                | K2(s1)  K2(s2)    1 |
+//    sigma(z) = __________________________
+//                  | K1(s1)  K2(s1) |
+//                  | K2(s1)  K2(s2) |
+// Where K0, K1, and K2 are successive zero-shifts of the K-polynomial.
+//
+// The K-polynomial shifts are otherwise exactly the same as Stage 2 after
+// accounting for a variable-shift sigma.
+bool JenkinsTraubSolver::ApplyQuadraticShiftToKPolynomial(
+    const std::complex<double>& root, const int max_iterations) {
+  // Only proceed if we have not already tried a quadratic shift.
+  if (attempted_quadratic_shift_) {
+    return false;
+  }
+
+  const double kTinyRelativeStep = 0.01;
+
+  // Compute the fixed-shift quadratic:
+  // sigma(z) = (x - m - n * i) * (x - m + n * i) = x^2 - 2 * m + m^2 + n^2.
+  sigma_(0) = 1.0;
+  sigma_(1) = -2.0 * root.real();
+  sigma_(2) = root.real() * root.real() + root.imag() * root.imag();
+
+  // These two containers hold values that we test for convergence such that the
+  // zero index is the convergence value from 2 iterations ago, the first
+  // index is from one iteration ago, and the second index is the current value.
+  VectorXd polynomial_quotient, polynomial_remainder, k_polynomial_quotient,
+      k_polynomial_remainder;
+  double poly_at_root(0), prev_poly_at_root(0), prev_v(0);
+  bool tried_fixed_shifts = false;
+
+  // These containers maintain a history of the predicted roots. The convergence
+  // of the algorithm is determined by the convergence of the root value.
+  std::vector<std::complex<double> > roots1, roots2;
+  roots1.push_back(root);
+  roots2.push_back(std::conj(root));
+  for (int i = 0; i < max_iterations; i++) {
+    // Terminate if the root evaluation is within our tolerance. This will
+    // return false if we do not have enough samples.
+    if (HasRootConverged(roots1) && HasRootConverged(roots2)) {
+      AddRootToOutput(roots1[1].real(), roots1[1].imag());
+      AddRootToOutput(roots2[1].real(), roots2[1].imag());
+      polynomial_ = polynomial_quotient;
+      return true;
+    }
+
+    QuadraticSyntheticDivision(
+        polynomial_, sigma_, &polynomial_quotient, &polynomial_remainder);
+
+    // Compute a and b from the above equations.
+    b_ = polynomial_remainder(0);
+    a_ = polynomial_remainder(1) - b_ * sigma_(1);
+
+    std::complex<double> roots[2];
+    FindQuadraticPolynomialRoots(sigma_(0), sigma_(1), sigma_(2), roots);
+
+    // Check that the roots are close. If not, then try a linear shift.
+    if (std::abs(std::abs(roots[0].real()) - std::abs(roots[1].real())) >
+        kRootPairTolerance * std::abs(roots[1].real())) {
+      return ApplyLinearShiftToKPolynomial(root, kMaxLinearShiftIterations);
+    }
+
+    // If the iteration is stalling at a root pair then apply a few fixed shift
+    // iterations to help convergence.
+    poly_at_root =
+        std::abs(a_ - roots[0].real() * b_) + std::abs(roots[0].imag() * b_);
+    const double rel_step = std::abs((sigma_(2) - prev_v) / sigma_(2));
+    if (!tried_fixed_shifts && rel_step < kTinyRelativeStep &&
+        prev_poly_at_root > poly_at_root) {
+      tried_fixed_shifts = true;
+      ApplyFixedShiftToKPolynomial(roots[0], kInnerFixedShiftIterations);
+    }
+
+    // Divide the shifted polynomial by the quadratic polynomial.
+    QuadraticSyntheticDivision(
+        k_polynomial_, sigma_, &k_polynomial_quotient, &k_polynomial_remainder);
+    d_ = k_polynomial_remainder(0);
+    c_ = k_polynomial_remainder(1) - d_ * sigma_(1);
+
+    prev_v = sigma_(2);
+    sigma_ = ComputeNextSigma();
+
+    // Compute K_next using the formula above.
+    UpdateKPolynomialWithQuadraticShift(polynomial_quotient,
+                                        k_polynomial_quotient);
+    k_polynomial_ /= k_polynomial_(0);
+    prev_poly_at_root = poly_at_root;
+
+    // Save the roots for convergence testing.
+    roots1.push_back(roots[0]);
+    roots2.push_back(roots[1]);
+    if (roots1.size() > 3) {
+      roots1.erase(roots1.begin());
+      roots2.erase(roots2.begin());
+    }
+  }
+
+  attempted_quadratic_shift_ = true;
+  return ApplyLinearShiftToKPolynomial(root, kMaxLinearShiftIterations);
+}
+
+// Generate K-Polynomials with variable-shifts that are linear. The shift is
+// computed as:
+//   K_next(z) = 1 / (z - s) * (K(z) - K(s) / P(s) * P(z))
+//   s_next = s - P(s) / K_next(s)
+bool JenkinsTraubSolver::ApplyLinearShiftToKPolynomial(
+    const std::complex<double>& root, const int max_iterations) {
+  if (attempted_linear_shift_) {
+    return false;
+  }
+
+  // Compute an initial guess for the root.
+  double real_root = (root -
+                      EvaluatePolynomial(polynomial_, root) /
+                          EvaluatePolynomial(k_polynomial_, root)).real();
+
+  VectorXd deflated_polynomial, deflated_k_polynomial;
+  double polynomial_at_root, k_polynomial_at_root;
+
+  // This container maintains a history of the predicted roots. The convergence
+  // of the algorithm is determined by the convergence of the root value.
+  std::vector<double> roots;
+  roots.push_back(real_root);;
+  for (int i = 0; i < max_iterations; i++) {
+    // Terminate if the root evaluation is within our tolerance. This will
+    // return false if we do not have enough samples.
+    if (HasRootConverged(roots)) {
+      AddRootToOutput(roots[1], 0);
+      polynomial_ = deflated_polynomial;
+      return true;
+    }
+
+    const double prev_polynomial_at_root = polynomial_at_root;
+    SyntheticDivisionAndEvaluate(
+        polynomial_, real_root, &deflated_polynomial, &polynomial_at_root);
+
+    // If the root is exactly the root then end early. Otherwise, the k
+    // polynomial will be filled with inf or nans.
+    if (polynomial_at_root == 0) {
+      AddRootToOutput(real_root, 0);
+      polynomial_ = deflated_polynomial;
+      return true;
+    }
+
+    // Update the K-Polynomial.
+    SyntheticDivisionAndEvaluate(k_polynomial_, real_root,
+                                 &deflated_k_polynomial, &k_polynomial_at_root);
+    k_polynomial_ = AddPolynomials(
+        deflated_k_polynomial,
+        -k_polynomial_at_root / polynomial_at_root * deflated_polynomial);
+
+    k_polynomial_ /= k_polynomial_(0);
+
+    // Compute the update for the root estimation.
+    k_polynomial_at_root = EvaluatePolynomial(k_polynomial_, real_root);
+    const double delta_root = polynomial_at_root / k_polynomial_at_root;
+    real_root -= polynomial_at_root / k_polynomial_at_root;
+
+    // Save the root so that convergence can be measured. Only the 3 most
+    // recently root values are needed.
+    roots.push_back(real_root);
+    if (roots.size() > 3) {
+      roots.erase(roots.begin());
+    }
+
+    // If the linear iterations appear to be stalling then we may have found a
+    // double real root of the form (z - x^2). Attempt a quadratic variable
+    // shift from the current estimate of the root.
+    if (i >= 2 &&
+        std::abs(delta_root) < 0.001 * std::abs(real_root) &&
+        std::abs(prev_polynomial_at_root) < std::abs(polynomial_at_root)) {
+      const std::complex<double> new_root(real_root, 0);
+      return ApplyQuadraticShiftToKPolynomial(new_root,
+                                              kMaxQuadraticShiftIterations);
+    }
+  }
+
+  attempted_linear_shift_ = true;
+  return ApplyQuadraticShiftToKPolynomial(root, kMaxQuadraticShiftIterations);
+}
+
+void JenkinsTraubSolver::AddRootToOutput(const double real, const double imag) {
+  if (real_roots_ != NULL) {
+    (*real_roots_)(num_solved_roots_) = real;
+  }
+  if (complex_roots_ != NULL) {
+    (*complex_roots_)(num_solved_roots_) = imag;
+  }
+  ++num_solved_roots_;
+}
+
+void JenkinsTraubSolver::RemoveZeroRoots() {
+  int num_zero_roots = 0;
+
+  const VectorXd::ReverseReturnType& creverse_polynomial =
+      polynomial_.reverse();
+  while (creverse_polynomial(num_zero_roots) == 0) {
+    ++num_zero_roots;
+  }
+
+  // The output roots have 0 as the default value so there is no need to
+  // explicitly add the zero roots.
+  polynomial_ = polynomial_.head(polynomial_.size() - num_zero_roots).eval();
+}
+
+bool JenkinsTraubSolver::SolveClosedFormPolynomial() {
+  const int degree = polynomial_.size() - 1;
+
+  // Is the polynomial constant?
+  if (degree == 0) {
+    // std::cout << "Trying to extract roots from a constant "
+    //          << "polynomial in FindPolynomialRoots" << std::endl;
+    // We return true with no roots, not false, as if the polynomial is constant
+    // it is correct that there are no roots. It is not the case that they were
+    // there, but that we have failed to extract them.
+    return true;
+  }
+  
+  // Linear
+  if (degree == 1) {
+    AddRootToOutput(FindLinearPolynomialRoots(polynomial_(0), polynomial_(1)),
+                    0);
+    return true;
+  }
+
+  // Quadratic
+  if (degree == 2) {
+    std::complex<double> roots[2];
+    FindQuadraticPolynomialRoots(polynomial_(0), polynomial_(1), polynomial_(2),
+                                 roots);
+    AddRootToOutput(roots[0].real(), roots[0].imag());
+    AddRootToOutput(roots[1].real(), roots[1].imag());
+    return true;
+  }
+
+  return false;
+}
+
+// Computes a lower bound on the radius of the roots of polynomial by examining
+// the Cauchy sequence:
+//
+//    z^n + |a_1| * z^{n - 1} + ... + |a_{n-1}| * z - |a_n|
+//
+// The unique positive zero of this polynomial is an approximate lower bound of
+// the radius of zeros of the original polynomial.
+double JenkinsTraubSolver::ComputeRootRadius() {
+  static const double kEpsilon = 1e-2;
+  static const int kMaxIterations = 100;
+
+  VectorXd poly = polynomial_;
+  // Take the absolute value of all coefficients.
+  poly = poly.array().abs();
+  // Negate the last coefficient.
+  poly.reverse()(0) *= -1.0;
+
+  // Find the unique positive zero using Newton-Raphson iterations.
+  double x0 = 1.0;
+  return FindRootIterativeNewton(poly, x0, kEpsilon, kMaxIterations);
+}
+
+// The k polynomial with a zero-shift is
+//  (K(x) - K(0) / P(0) * P(x)) / x.
+//
+// This is equivalent to:
+//    K(x) - K(0)      K(0)     P(x) - P(0)
+//    ___________   -  ____  *  ___________
+//         x           P(0)          x
+//
+// Note that removing the constant term and dividing by x is equivalent to
+// shifting the polynomial to one degree lower in our representation.
+void JenkinsTraubSolver::ComputeZeroShiftKPolynomial() {
+  // Evaluating the polynomial at zero is equivalent to the constant term
+  // (i.e. the last coefficient). Note that reverse() is an expression and does
+  // not actually reverse the vector elements.
+  const double polynomial_at_zero = polynomial_(polynomial_.size() - 1);
+  const double k_at_zero = k_polynomial_(k_polynomial_.size() - 1);
+
+  k_polynomial_ = AddPolynomials(k_polynomial_.head(k_polynomial_.size() - 1),
+                                 -k_at_zero / polynomial_at_zero *
+                                     polynomial_.head(polynomial_.size() - 1));
+}
+
+// The iterations are computed with the following equation:
+//              a^2 + u * a * b + v * b^2
+//   K_next =  ___________________________ * Q_K
+//                    b * c - a * d
+//
+//                      a * c + u * a * d + v * b * d
+//             +  (z - _______________________________) * Q_P + b.
+//                              b * c - a * d
+//
+// This is done using *only* realy arithmetic so it can be done very fast!
+void JenkinsTraubSolver::UpdateKPolynomialWithQuadraticShift(
+    const VectorXd& polynomial_quotient,
+    const VectorXd& k_polynomial_quotient) {
+  const double coefficient_q_k =
+      (a_ * a_ + sigma_(1) * a_ * b_ + sigma_(2) * b_ * b_) /
+      (b_ * c_ - a_ * d_);
+  VectorXd linear_polynomial(2);
+  linear_polynomial(0) = 1.0;
+  linear_polynomial(1) =
+      -(a_ * c_ + sigma_(1) * a_ * d_ + sigma_(2) * b_ * d_) /
+      (b_ * c_ - a_ * d_);
+  k_polynomial_ = AddPolynomials(
+      coefficient_q_k * k_polynomial_quotient,
+      MultiplyPolynomials(linear_polynomial, polynomial_quotient));
+  k_polynomial_(k_polynomial_.size() - 1) += b_;
+}
+
+// Using a bit of algebra, the update of sigma(z) can be computed from the
+// previous value along with a, b, c, and d defined above. The details of this
+// simplification can be found in "Three Stage Variable-Shift Iterations for the
+// Solution of Polynomial Equations With a Posteriori Error Bounds for the
+// Zeros" by M.A. Jenkins, Doctoral Thesis, Stanford Univeristy, 1969.
+//
+// NOTE: we assume the leading term of quadratic_sigma is 1.0.
+VectorXd JenkinsTraubSolver::ComputeNextSigma() {
+  const double u = sigma_(1);
+  const double v = sigma_(2);
+
+  const VectorXd::ReverseReturnType& creverse_k_polynomial =
+      k_polynomial_.reverse();
+  const VectorXd::ReverseReturnType& creverse_polynomial =
+      polynomial_.reverse();
+
+  const double b1 = -creverse_k_polynomial(0) / creverse_polynomial(0);
+  const double b2 = -(creverse_k_polynomial(1) + b1 * creverse_polynomial(1)) /
+                    creverse_polynomial(0);
+
+  const double a1 = b_* c_ - a_ * d_;
+  const double a2 = a_ * c_ + u * a_ * d_ + v * b_* d_;
+  const double c2 = b1 * a2;
+  const double c3 = b1 * b1 * (a_ * a_ + u * a_ * b_ + v * b_ * b_);
+  const double c4 = v * b2 * a1 - c2 - c3;
+  const double c1 = c_ * c_ + u * c_ * d_ + v * d_ * d_ +
+                    b1 * (a_ * c_ + u * b_ * c_ + v * b_ * d_) - c4;
+  const double delta_u = -(u * (c2 + c3) + v * (b1 * a1 + b2 * a2)) / c1;
+  const double delta_v = v * c4 / c1;
+
+  // Update u and v in the quadratic sigma.
+  VectorXd new_quadratic_sigma(3);
+  new_quadratic_sigma(0) = 1.0;
+  new_quadratic_sigma(1) = u + delta_u;
+  new_quadratic_sigma(2) = v + delta_v;
+  return new_quadratic_sigma;
+}
+
+}  // namespace
+
+bool FindPolynomialRootsJenkinsTraub(const VectorXd& polynomial,
+                                     VectorXd* real_roots,
+                                     VectorXd* complex_roots) {
+  JenkinsTraubSolver solver(polynomial, real_roots, complex_roots);
+  return solver.ExtractRoots();
+}
+
+}  // namespace rpoly_plus_plus
diff --git a/scilab/modules/polynomials/src/cpp/find_polynomial_roots_jenkins_traub.h b/scilab/modules/polynomials/src/cpp/find_polynomial_roots_jenkins_traub.h
new file mode 100644 (file)
index 0000000..f098910
--- /dev/null
@@ -0,0 +1,58 @@
+// Copyright (C) 2015 Chris Sweeney. All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are
+// met:
+//
+//     * Redistributions of source code must retain the above copyright
+//       notice, this list of conditions and the following disclaimer.
+//
+//     * Redistributions in binary form must reproduce the above
+//       copyright notice, this list of conditions and the following
+//       disclaimer in the documentation and/or other materials provided
+//       with the distribution.
+//
+//     * Neither the name of Chris Sweeney nor the names of its contributors may
+//       be used to endorse or promote products derived from this software
+//       without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Please contact the author of this library if you have any questions.
+// Author: Chris Sweeney (cmsweeney@cs.ucsb.edu)
+
+#ifndef RPOLY_PLUS_PLUS_FIND_POLYNOMIAL_ROOTS_JENKINS_TRAUB_H_
+#define RPOLY_PLUS_PLUS_FIND_POLYNOMIAL_ROOTS_JENKINS_TRAUB_H_
+
+#include <Eigen/Core>
+
+namespace rpoly_plus_plus
+{
+
+// A three-stage algorithm for finding roots of polynomials with real
+// coefficients as outlined in: "A Three-Stage Algorithm for Real Polynomaials
+// Using Quadratic Iteration" by Jenkins and Traub, SIAM 1970. Please note that
+// this variant is different than the complex-coefficient version, and is
+// estimated to be up to 4 times faster.
+//
+// The algorithm works by computing shifts in so-called "K-polynomials" that
+// exaggerate the roots. Once a root is found (or in the real-polynomial case, a
+// pair of roots) then it is divided from the polynomial and the process is
+// repeated.
+bool FindPolynomialRootsJenkinsTraub(const Eigen::VectorXd& polynomial,
+                                     Eigen::VectorXd* real_roots,
+                                     Eigen::VectorXd* complex_roots);
+
+}  // namespace rpoly_plus_plus
+
+#endif  // RPOLY_PLUS_PLUS_FIND_POLYNOMIAL_ROOTS_JENKINS_TRAUB_H_
diff --git a/scilab/modules/polynomials/src/cpp/polynomial.cc b/scilab/modules/polynomials/src/cpp/polynomial.cc
new file mode 100644 (file)
index 0000000..98b0153
--- /dev/null
@@ -0,0 +1,113 @@
+// Copyright (C) 2015 Chris Sweeney. All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are
+// met:
+//
+//     * Redistributions of source code must retain the above copyright
+//       notice, this list of conditions and the following disclaimer.
+//
+//     * Redistributions in binary form must reproduce the above
+//       copyright notice, this list of conditions and the following
+//       disclaimer in the documentation and/or other materials provided
+//       with the distribution.
+//
+//     * Neither the name of Chris Sweeney nor the names of its contributors may
+//       be used to endorse or promote products derived from this software
+//       without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Please contact the author of this library if you have any questions.
+// Author: Chris Sweeney (cmsweeney@cs.ucsb.edu)
+
+#include "polynomial.h"
+
+#include <Eigen/Core>
+#include <cmath>
+
+namespace rpoly_plus_plus {
+
+using Eigen::MatrixXd;
+using Eigen::VectorXd;
+using Eigen::VectorXcd;
+
+// Remove leading terms with zero coefficients.
+VectorXd RemoveLeadingZeros(const VectorXd& polynomial_in) {
+  int i = 0;
+  while (i < (polynomial_in.size() - 1) && polynomial_in(i) == 0) {
+    ++i;
+  }
+  return polynomial_in.tail(polynomial_in.size() - i);
+}
+
+VectorXd DifferentiatePolynomial(const VectorXd& polynomial) {
+  const int degree = polynomial.rows() - 1;
+
+  // Degree zero polynomials are constants, and their derivative does
+  // not result in a smaller degree polynomial, just a degree zero
+  // polynomial with value zero.
+  if (degree == 0) {
+    return VectorXd::Zero(1);
+  }
+
+  VectorXd derivative(degree);
+  for (int i = 0; i < degree; ++i) {
+    derivative(i) = (degree - i) * polynomial(i);
+  }
+
+  return derivative;
+}
+
+VectorXd MultiplyPolynomials(const VectorXd& poly1, const VectorXd& poly2) {
+  VectorXd multiplied_poly = VectorXd::Zero(poly1.size() + poly2.size() - 1);;
+  for (int i = 0; i < poly1.size(); i++) {
+    for (int j = 0; j < poly2.size(); j++) {
+      multiplied_poly.reverse()(i + j) +=
+          poly1.reverse()(i) * poly2.reverse()(j);
+    }
+  }
+  return multiplied_poly;
+}
+
+VectorXd AddPolynomials(const VectorXd& poly1, const VectorXd& poly2) {
+  if (poly1.size() > poly2.size()) {
+    VectorXd sum = poly1;
+    sum.tail(poly2.size()) += poly2;
+    return sum;
+  } else {
+    VectorXd sum = poly2;
+    sum.tail(poly1.size()) += poly1;
+    return sum;
+  }
+}
+
+double FindRootIterativeNewton(const Eigen::VectorXd& polynomial,
+                               const double x0,
+                               const double epsilon,
+                               const int max_iterations) {
+  double root = x0;
+  const Eigen::VectorXd derivative = DifferentiatePolynomial(polynomial);
+  double prev;
+  for (int i = 0; i < max_iterations; i++) {
+    prev = root;
+    root -= EvaluatePolynomial(polynomial, root) /
+            EvaluatePolynomial(derivative, root);
+    if (std::abs(prev - root) < epsilon) {
+      break;
+    }
+  }
+  return root;
+}
+
+}  // namespace rpoly_plus_plus
diff --git a/scilab/modules/polynomials/src/cpp/polynomial.h b/scilab/modules/polynomials/src/cpp/polynomial.h
new file mode 100644 (file)
index 0000000..fcb3e97
--- /dev/null
@@ -0,0 +1,84 @@
+// Copyright (C) 2015 Chris Sweeney
+// All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are
+// met:
+//
+//     * Redistributions of source code must retain the above copyright
+//       notice, this list of conditions and the following disclaimer.
+//
+//     * Redistributions in binary form must reproduce the above
+//       copyright notice, this list of conditions and the following
+//       disclaimer in the documentation and/or other materials provided
+//       with the distribution.
+//
+//     * Neither the name of Chris Sweeney nor the names of its contributors may
+//       be used to endorse or promote products derived from this software
+//       without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Please contact the author of this library if you have any questions.
+// Author: Chris Sweeney (cmsweeney@cs.ucsb.edu)
+
+#ifndef RPOLY_PLUS_PLUS_POLYNOMIAL_H_
+#define RPOLY_PLUS_PLUS_POLYNOMIAL_H_
+
+#include <Eigen/Core>
+
+namespace rpoly_plus_plus
+{
+
+// All polynomials are assumed to be the form
+//
+//   sum_{i=0}^N polynomial(i) x^{N-i}.
+//
+// and are given by a vector of coefficients of size N + 1.
+
+// Remove leading terms with zero coefficients.
+Eigen::VectorXd RemoveLeadingZeros(const Eigen::VectorXd& polynomial_in);
+
+// Evaluate the polynomial at x using the Horner scheme.
+template <typename T>
+inline T EvaluatePolynomial(const Eigen::VectorXd& polynomial, const T& x)
+{
+    T v = 0.0;
+    for (int i = 0; i < polynomial.size(); ++i)
+    {
+        v = v * x + polynomial(i);
+    }
+    return v;
+}
+
+// Return the derivative of the given polynomial. It is assumed that
+// the input polynomial is at least of degree zero.
+Eigen::VectorXd DifferentiatePolynomial(const Eigen::VectorXd& polynomial);
+
+// Multiplies the two polynoimals together.
+Eigen::VectorXd MultiplyPolynomials(const Eigen::VectorXd& poly1,
+                                    const Eigen::VectorXd& poly2);
+
+// Adds two polynomials together.
+Eigen::VectorXd AddPolynomials(const Eigen::VectorXd& poly1,
+                               const Eigen::VectorXd& poly2);
+
+// Find a root from the starting guess using Newton's method.
+double FindRootIterativeNewton(const Eigen::VectorXd& polynomial,
+                               const double x0,
+                               const double epsilon,
+                               const int max_iterations);
+
+}  // namespace rpoly_plus_plus
+
+#endif  // RPOLY_PLUS_PLUS_POLYNOMIAL_H_
diff --git a/scilab/modules/polynomials/src/cpp/rpoly.cpp b/scilab/modules/polynomials/src/cpp/rpoly.cpp
new file mode 100644 (file)
index 0000000..0008947
--- /dev/null
@@ -0,0 +1,73 @@
+/*
+ * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+ * Copyright (C) 2016 - Scilab Enterprises - Clement DAVID
+ *
+ * Copyright (C) 2012 - 2016 - Scilab Enterprises
+ *
+ * This file is hereby licensed under the terms of the GNU GPL v2.0,
+ * pursuant to article 5.3.4 of the CeCILL v.2.1.
+ * This file was originally licensed under the terms of the CeCILL v2.1,
+ * and continues to be available under such terms.
+ * For more information, see the COPYING file which you should have received
+ * along with this program.
+ *
+ */
+
+extern "C" {
+#include "dynlib_polynomials.h"
+#include "machine.h"            /* C2F */
+
+    POLYNOMIALS_IMPEXP int C2F(rpoly)(double* op, int* degree, double* zeror, double* zeroi, int* fail);
+}
+
+#include "find_polynomial_roots_jenkins_traub.h"
+#include <Eigen/Core>
+
+/**
+ * finds the zeros of a real polynomial (compatible with the original rpoly.f)
+ *
+ * \param op  double precision vector of coefficients in
+ *            order of decreasing powers.
+ * \param degree integer degree of polynomial.
+ * \param zeror real part of the zeros
+ * \param zeroi imaginary part of the zeros
+ * \param fail output parameter,
+ *             2 if  leading coefficient is zero
+ *             1 for non convergence or if rpoly has found fewer than degree
+ *               zeros. In the latter case degree is reset to the number of
+ *               zeros found.
+ *             3 if degree>100
+ */
+POLYNOMIALS_IMPEXP int C2F(rpoly)(double* op, int* degree, double* zeror, double* zeroi, int* fail)
+{
+    if (*degree > 100)
+    {
+        *fail = 3;
+        return 0;
+    }
+
+    // let's copy there as Map<VectorXd> is not supported yet on rpoly_plus_plus
+    Eigen::Map<Eigen::VectorXd> mapped_op(op, *degree + 1);
+    Eigen::Map<Eigen::VectorXd> mapped_zeror(zeror, *degree);
+    Eigen::Map<Eigen::VectorXd> mapped_zeroi(zeroi, *degree);
+
+    // reverse as the polynomials are used in different ordering
+    Eigen::VectorXd polynomial(mapped_op);
+    Eigen::VectorXd real_roots(*degree);
+    Eigen::VectorXd complex_roots(*degree);
+
+    bool valid = rpoly_plus_plus::FindPolynomialRootsJenkinsTraub(polynomial, &real_roots, &complex_roots);
+    if (!valid)
+    {
+        *fail = 1;
+        return 0;
+    }
+
+    mapped_zeror = real_roots;
+    mapped_zeroi = complex_roots;
+
+    *fail = 0;
+    return 0;
+}
+
+
index d532e37..22d4a9f 100644 (file)
@@ -95,7 +95,6 @@
                <File RelativePath=".\realit.f"/>
                <File RelativePath=".\recbez.f"/>
                <File RelativePath=".\residu.f"/>
-               <File RelativePath=".\rpoly.f"/>
                <File RelativePath=".\sfact1.f"/>
                <File RelativePath=".\sfact2.f"/>
                <File RelativePath=".\wdmpad.f"/>
index 80297aa..3bf928a 100644 (file)
@@ -304,7 +304,6 @@ cd ..
     <ClCompile Include="realit.c" />
     <ClCompile Include="recbez.c" />
     <ClCompile Include="residu.c" />
-    <ClCompile Include="rpoly.c" />
     <ClCompile Include="sfact1.c" />
     <ClCompile Include="sfact2.c" />
     <ClCompile Include="wdmpad.c" />
@@ -362,7 +361,6 @@ cd ..
     <f2c_rule Include="realit.f" />
     <f2c_rule Include="recbez.f" />
     <f2c_rule Include="residu.f" />
-    <f2c_rule Include="rpoly.f" />
     <f2c_rule Include="sfact1.f" />
     <f2c_rule Include="sfact2.f" />
     <f2c_rule Include="wdmpad.f" />
@@ -400,4 +398,4 @@ cd ..
   <ImportGroup Label="ExtensionTargets">
     <Import Project="..\..\..\..\Visual-Studio-settings\f2c.targets" />
   </ImportGroup>
-</Project>
\ No newline at end of file
+</Project>
index 28a5d5a..8311adf 100644 (file)
     <ClCompile Include="residu.c">
       <Filter>Source Files</Filter>
     </ClCompile>
-    <ClCompile Include="rpoly.c">
-      <Filter>Source Files</Filter>
-    </ClCompile>
     <ClCompile Include="sfact1.c">
       <Filter>Source Files</Filter>
     </ClCompile>
     <f2c_rule Include="residu.f">
       <Filter>Fortran files</Filter>
     </f2c_rule>
-    <f2c_rule Include="rpoly.f">
-      <Filter>Fortran files</Filter>
-    </f2c_rule>
     <f2c_rule Include="sfact1.f">
       <Filter>Fortran files</Filter>
     </f2c_rule>
       <Filter>Libraries Dependencies</Filter>
     </None>
   </ItemGroup>
-</Project>
\ No newline at end of file
+</Project>
diff --git a/scilab/modules/polynomials/src/fortran/rpoly.f b/scilab/modules/polynomials/src/fortran/rpoly.f
deleted file mode 100644 (file)
index 241cb31..0000000
+++ /dev/null
@@ -1,276 +0,0 @@
-c
-c From NETLIB : TOMS/493
-c Based on Jenkins-Traub algorithm :
-c http://en.wikipedia.org/wiki/Jenkins-Traub_method
-c
-      subroutine rpoly(op, degree, zeror, zeroi, fail)
-c!purpose
-c finds the zeros of a real polynomial
-c!calling sequence
-c op  - double precision vector of coefficients in
-c       order of decreasing powers.
-c degree   - integer degree of polynomial.
-c zeror, zeroi - output double precision vectors of
-c                real and imaginary parts of the
-c                zeros.
-c fail  - output parameter,
-c        2 if  leading coefficient is zero
-c        1 for non convergence or if rpoly
-c         has found fewer than degree zeros.
-c         in the latter case degree is reset to
-c         the number of zeros found.
-c        3 if degree>100
-c!comments
-c to change the size of polynomials which can be
-c solved, reset the dimensions of the arrays in the
-c common area and in the following declarations.
-c the subroutine uses single precision calculations
-c for scaling, bounds and error calculations. all
-c calculations for the iterations are done in double
-c precision.
-c!
-      external dlamch,slamch
-      double precision dlamch
-      real slamch
-c
-      common /gloglo/ p, qp, k, qk, svk, sr, si, u,
-     * v, a, b, c, d, a1, a2, a3, a6, a7, e, f, g,
-     * h, szr, szi, lzr, lzi, eta, are, mre, n, nn
-      double precision p(101), qp(101), k(101),
-     * qk(101), svk(101), sr, si, u, v, a, b, c, d,
-     * a1, a2, a3, a6, a7, e, f, g, h, szr, szi,
-     * lzr, lzi
-      double precision op(*), temp(101),
-     * zeror(*), zeroi(*), t, aa, bb, cc,factor
-      real ptt(101), lo, maxi, mini, xx, yy, cosr,
-     * sinr, xxx, x, sc, bnd, xm, ff, df, dx, infin,
-     * smalno, base,eta,are,mre
-      integer degree, cnt, nz, i, j, jj, nm1, n, nn
-      logical  zerok
-      integer fail
-c
-      real ZERO
-      parameter (ZERO = 0.0E+0)
-
-      if(degree.gt.100) goto 300
-c the following statements set machine constants used
-c in various parts of the program. the meaning of the
-c four constants are...
-c eta     the maximum relative representation error
-c         which can be described as the smallest
-c         positive floating point number such that
-c         1.do+eta is greater than 1.
-c infiny  the largest floating-point number.
-c smalno  the smallest positive floating-point number
-c         if the exponent range differs in single and
-c         double precision then smalno and infin
-c         should indicate the smaller range.
-c base    the base of the floating-point number
-c         system used.
-
-c
-c Rely on compiler instead of Lapack
-c slamch function does not work under macosX
-c replace this by compiler stuffs working on each platform
-c http://www.netlib.org/lapack/util/slamch.f
-c
-c     slamch('u') <=> TINY(ZERO)
-      smalno=TINY(ZERO)
-c     slamch('o') <=> HUGE(ZERO)
-      infin=HUGE(ZERO)
-c     slamch('b') <=> RADIX(ZERO)
-      base=RADIX(ZERO)
-      eta=real(dlamch('p'))
-c are and mre refer to the unit error in + and *
-c respectively. they are assumed to be the same as
-c eta.
-      are = eta
-      mre = eta
-      lo = smalno/eta
-c initialization of constants for shift rotation
-      xx = 0.707106780d+0
-      yy = -xx
-      cosr = -0.0697564740d+0
-      sinr =  0.997564050d+0
-      fail = 0.
-      n = degree
-      nn = n + 1
-c algorithm fails if the leading coefficient is zero.
-      if (op(1).ne.0.0d+0) go to 10
-      fail = 2
-      degree = 0
-      return
-c make a copy of the coefficients
-   10 do 20 i=1,nn
-        p(i) = op(i)
-   20 continue
-c remove the zeros at the origin if any
-   30 if (p(nn).ne.0.0d+0) go to 40
-      j = degree - n + 1
-      zeror(j) = 0.0d+0
-      zeroi(j) = 0.0d+0
-      nn = nn - 1
-      n = n - 1
-      go to 30
-c start the algorithm for one zero
-   40 if (n.gt.2) go to 60
-      if (n.lt.1) return
-c calculate the final zero or pair zeros
-      if (n.eq.2) go to 50
-      zeror(degree) = - p(2)/p(1)
-      zeroi(degree) = 0.0d+0
-      return
-   50 call quad(p(1), p(2), p(3), zeror(degree-1),
-     * zeroi(degree-1), zeror(degree), zeroi(degree))
-      return
-c find largest and smallest moduli of coefficients.
-   60 maxi = 0.
-      mini = infin
-      do 70 i=1,nn
-        x = abs(real(p(i)))
-        if (x.gt.maxi) maxi = x
-        if (x.ne.0. .and. x.lt.mini) mini = x
-   70 continue
-C      maxi=min(infin,maxi) bug "f77 -mieee-with-inexact"
-      if (infin.lt.maxi) maxi=infin
-c scale if there are large or very small coefficients
-c computes a scale factor to multiply the
-c coefficients of the polynomial. the scaling is done
-c to avoid overflow and to avoid undetected underflow
-c interfering with the convergence criterion.
-c the factor is a power of the base
-      sc = lo/mini
-      if (sc.gt.1.0) go to 80
-      if (maxi.lt.10.) go to 110
-      if (sc.eq.0.) sc = smalno
-      go to 90
-   80 if (infin/sc.lt.maxi) go to 110
-   90 l = log(sc)/log(base) + .5
-      factor = (base*1.0d+0)**l
-      if (factor.eq.1.0d+0) go to 110
-      do 100 i=1,nn
-        p(i) = factor*p(i)
-  100 continue
-c compute lower bound on moduli of zeros.
-  110 do 120 i=1,nn
-c        ptt(i) = min(infin,abs(real(p(i)))) bug "f77 -mieee-with-inexact"
-         ptt(i) = abs(real(p(i)))
-         if (infin.lt.abs(real(p(i)))) ptt(i)=infin
-  120 continue
-      ptt(nn) = -ptt(nn)
-c compute upper estimate of bound
-      x = exp((log(-ptt(nn))-log(ptt(1)))/real(n))
-      if (ptt(n).eq.0.) go to 130
-c if newton step at the origin is better, use it.
-      xm = -ptt(nn)/ptt(n)
-      if (xm.lt.x) x = xm
-c chop the interval (0,x) until ff .le. 0
-  130 xm = x*.1
-      ff = ptt(1)
-      do 140 i=2,nn
-        ff = ff*xm + ptt(i)
-  140 continue
-      if (ff.le.0) go to 150
-      if(ff.gt.infin) goto 310
-      x = xm
-      go to 130
-  150 dx = x
-c do newton iteration until x converges to two
-c decimal places
-  160 if (abs(dx/x).le..005) go to 180
-      ff = ptt(1)
-      df = ff
-      do 170 i=2,n
-        ff = ff*x + ptt(i)
-        df = df*x + ff
-  170 continue
-      ff = ff*x + ptt(nn)
-      if(ff.gt.infin) goto 310
-      dx = ff/df
-      x = x - dx
-      go to 160
-  180 bnd = x
-c compute the derivative as the initial k polynomial
-c and do 5 steps with no shift
-      nm1 = n - 1
-      do 190 i=2,n
-        k(i) = real(nn-i)*p(i)/real(n)
-  190 continue
-      k(1) = p(1)
-      aa = p(nn)
-      bb = p(n)
-      zerok = k(n).eq.0.0d+0
-      do 230 jj=1,5
-        cc = k(n)
-        if (zerok) go to 210
-c use scaled form of recurrence if value of k at 0 is
-c nonzero
-        t = -aa/cc
-        do 200 i=1,nm1
-          j = nn - i
-          k(j) = t*k(j-1) + p(j)
-  200   continue
-        k(1) = p(1)
-        zerok = abs(k(n)).le.abs(bb)*eta*10.
-        go to 230
-c use unscaled form form of recurrence
-  210   do 220 i=1,nm1
-          j = nn - i
-          k(j) = k(j-1)
-  220   continue
-        k(1) = 0.0d+0
-        zerok = k(n).eq.0.0d+0
-  230 continue
-c save k for restarts with new shifts
-      do 240 i=1,n
-        temp(i) = k(i)
-  240 continue
-c loop to select the quadratic  corresponding to each
-c new shift
-      do 280 cnt=1,20
-c quadratic corresponds to a double shift to a
-c non-real point and its complex conjugate. the point
-c has modulus bnd and amplitude rotated by 94 degrees
-c from the previous shift
-        xxx = cosr*xx - sinr*yy
-        yy = sinr*xx + cosr*yy
-        xx = xxx
-        sr = bnd*xx
-        si = bnd*yy
-        u = -2.0d+0*sr
-        v = bnd
-c second stage calculation, fixed quadratic
-        call fxshfr(20*cnt, nz)
-        if (nz.eq.0) go to 260
-c the second stage jumps directly to one of the third
-c stage iterations and returns here if successful.
-c deflate the polynomial, store the zero or zeros and
-c return to the main algorithm.
-        j = degree - n + 1
-        zeror(j) = szr
-        zeroi(j) = szi
-        nn = nn - nz
-        n = nn - 1
-        do 250 i=1,nn
-          p(i) = qp(i)
-  250   continue
-        if (nz.eq.1) go to 40
-        zeror(j+1) = lzr
-        zeroi(j+1) = lzi
-        go to 40
-c if the iteration is unsuccessful another quadratic
-c is chosen after restoring k
-  260   do 270 i=1,n
-           k(i) = temp(i)
-  270   continue
-  280 continue
-c return with failure if no convergence with 20
-c shifts
-      fail = 1
-      degree = degree - n
-      return
-  300 fail=3
-      return
- 310  fail=1
-      return
-      end
index 1b768fc..0f8118c 100644 (file)
@@ -14,30 +14,30 @@ function sortedRoots=sortRoots(rootsToSort)
     sortedRoots = rootsToSort(kRoots);
 endfunction
 function checkroots(p,expectedroots,varargin)
-     // Checks the roots function against given roots.
-     //
-     // 1. Check default algorithm
-     myroots=roots(p);
-     computedroots = sortRoots(myroots);
-     expectedroots  = sortRoots(expectedroots);
-     assert_checkalmostequal(computedroots,expectedroots,varargin(:));
-     //
-     // 2. Check "e" algorithm
-     myroots=roots(p,"e");
-     computedroots = sortRoots(myroots);
-     expectedroots  = sortRoots(expectedroots);
-     assert_checkalmostequal(computedroots,expectedroots,varargin(:));
-     //
-     // 3. Check "f" algorithm
-     if ( isreal(p) ) then
-         myroots=roots(p,"f");
-         computedroots = sortRoots(myroots);
-         expectedroots  = sortRoots(expectedroots);
-         assert_checkalmostequal(computedroots,expectedroots,varargin(:));
-     end
+    // Checks the roots function against given roots.
+    //
+    // 1. Check default algorithm
+    myroots=roots(p);
+    computedroots = sortRoots(myroots);
+    expectedroots  = sortRoots(expectedroots);
+    assert_checkalmostequal(computedroots,expectedroots,varargin(:));
+    //
+    // 2. Check "e" algorithm
+    myroots=roots(p,"e");
+    computedroots = sortRoots(myroots);
+    expectedroots  = sortRoots(expectedroots);
+    assert_checkalmostequal(computedroots,expectedroots,varargin(:));
+    //
+    // 3. Check "f" algorithm
+    if ( isreal(p) ) then
+        myroots=roots(p,"f");
+        computedroots = sortRoots(myroots);
+        expectedroots  = sortRoots(expectedroots);
+        assert_checkalmostequal(computedroots,expectedroots,varargin(:));
+    end
 endfunction
 //   Check the computation of the roots of a polynomial
-//   with different kinds of polynomials and different 
+//   with different kinds of polynomials and different
 //   kinds of roots :
 //   - real poly,
 //   - complex poly,
@@ -79,42 +79,6 @@ checkroots(p,expectedroots,%eps);
 p=(%s-%pi)^2;
 expectedroots  = [%pi;%pi];
 checkroots(p,expectedroots,10*%eps);
-//
-// Caution !
-// The following are difficult root-finding problems
-// with expected precision problems.
-// See "Principles for testing polynomial 
-// zerofinding programs"
-// Jenkins, Traub
-// 1975
-// p.28
-// "The accuracy which one may expect to achieve in calculating
-// zeros is limited by the condition of these zeros. In particular,
-// for multiple zeros perturbations of size epsilon in the 
-// coefficients cause perturbations of size epsilon^(1/m)
-// in the zeros."
-// 
-//
-// 3 real roots with a zero derivate at the root
-// Really difficult problem : only simple precision computed, instead of double precision ***
-p=(%s-%pi)^3;
-expectedroots  = [%pi;%pi;%pi];
-checkroots(p,expectedroots,%eps^(1/3),5*%eps^(1/3));
-// 4 real roots with a zero derivate at the root
-// Really difficult problem : only simple precision
-p=(%s-%pi)^4;
-expectedroots  = [%pi;%pi;%pi;%pi];
-checkroots(p,expectedroots,%eps^(1/4),5*%eps^(1/4))
-// 10 real roots with a zero derivate at the root
-// Really difficult problem : only one correct digit
-p=(%s-%pi)^10;
-expectedroots  = [%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi];
-checkroots(p,expectedroots,%eps^(1/10),8*%eps^(1/10))
-// "Numerical computing with Matlab", Cleve Moler.
-A = diag(1:20);
-p = poly(A,'x');
-e = [1:20]';
-checkroots(p,e,%eps,0.2);
 // Tests from CPOLY
 // M. A. Jenkins and J. F. Traub. 1972.
 // Algorithm 419: zeros of a complex polynomial.
@@ -242,7 +206,7 @@ E = [
 E = sortRoots(E);
 R = sortRoots(R);
 assert_checkalmostequal(R, E, 1.e-3, 1.e-3);
-// EXAMPLE 5. 12 ZEROS EVENLY DISTRIBUTE ON A CIRCLE OF 
+// EXAMPLE 5. 12 ZEROS EVENLY DISTRIBUTE ON A CIRCLE OF
 // RADIUS 1. CENTERED AT 0+2I
 // Real part:
 // 4095 - 67584*x^2 + 126720*x^4 - 59136*x^6 + 7920*x^8 - 264*x^10 + x^12
@@ -296,5 +260,5 @@ E = [
 E = sortRoots(E);
 R = sortRoots(R);
 assert_checkalmostequal(R, E, 1.e-10, 1.e-8);
-assert_checkequal(roots([4 3 2 1]), roots(poly([1 2 3 4], 'x', 'coeff')));
-assert_checkequal(roots([4 3 2 1] + [1 2 3 4]*%i), roots(poly([1 2 3 4]+[4 3 2 1]*%i,'x','coeff')));
+assert_checkequal(roots([4 3 2 1]), roots(poly([1 2 3 4], "x", "coeff")));
+assert_checkequal(roots([4 3 2 1] + [1 2 3 4]*%i), roots(poly([1 2 3 4]+[4 3 2 1]*%i,"x","coeff")));
index 11a8441..b942015 100644 (file)
@@ -9,7 +9,6 @@
 
 // <-- CLI SHELL MODE -->
 
-
 function sortedRoots=sortRoots(rootsToSort)
     //Sort roots using rounded values to avoid rounding errors
     //Here 10000 is ok due to roots values
@@ -18,31 +17,31 @@ function sortedRoots=sortRoots(rootsToSort)
 endfunction
 
 function checkroots(p,expectedroots,varargin)
-     // Checks the roots function against given roots.
-     //
-     // 1. Check default algorithm
-     myroots=roots(p);
-     computedroots = sortRoots(myroots);
-     expectedroots  = sortRoots(expectedroots);
-     assert_checkalmostequal(computedroots,expectedroots,varargin(:));
-     //
-     // 2. Check "e" algorithm
-     myroots=roots(p,"e");
-     computedroots = sortRoots(myroots);
-     expectedroots  = sortRoots(expectedroots);
-     assert_checkalmostequal(computedroots,expectedroots,varargin(:));
-     //
-     // 3. Check "f" algorithm
-     if ( isreal(p) ) then
-         myroots=roots(p,"f");
-         computedroots = sortRoots(myroots);
-         expectedroots  = sortRoots(expectedroots);
-         assert_checkalmostequal(computedroots,expectedroots,varargin(:));
-     end
+    // Checks the roots function against given roots.
+    //
+    // 1. Check default algorithm
+    myroots=roots(p);
+    computedroots = sortRoots(myroots);
+    expectedroots  = sortRoots(expectedroots);
+    assert_checkalmostequal(computedroots,expectedroots,varargin(:));
+    //
+    // 2. Check "e" algorithm
+    myroots=roots(p,"e");
+    computedroots = sortRoots(myroots);
+    expectedroots  = sortRoots(expectedroots);
+    assert_checkalmostequal(computedroots,expectedroots,varargin(:));
+    //
+    // 3. Check "f" algorithm
+    if ( isreal(p) ) then
+        myroots=roots(p,"f");
+        computedroots = sortRoots(myroots);
+        expectedroots  = sortRoots(expectedroots);
+        assert_checkalmostequal(computedroots,expectedroots,varargin(:));
+    end
 endfunction
 
 //   Check the computation of the roots of a polynomial
-//   with different kinds of polynomials and different 
+//   with different kinds of polynomials and different
 //   kinds of roots :
 //   - real poly,
 //   - complex poly,
@@ -86,43 +85,6 @@ checkroots(p,expectedroots,%eps);
 p=(%s-%pi)^2;
 expectedroots  = [%pi;%pi];
 checkroots(p,expectedroots,10*%eps);
-//
-// Caution !
-// The following are difficult root-finding problems
-// with expected precision problems.
-// See "Principles for testing polynomial 
-// zerofinding programs"
-// Jenkins, Traub
-// 1975
-// p.28
-// "The accuracy which one may expect to achieve in calculating
-// zeros is limited by the condition of these zeros. In particular,
-// for multiple zeros perturbations of size epsilon in the 
-// coefficients cause perturbations of size epsilon^(1/m)
-// in the zeros."
-// 
-//
-// 3 real roots with a zero derivate at the root
-// Really difficult problem : only simple precision computed, instead of double precision ***
-p=(%s-%pi)^3;
-expectedroots  = [%pi;%pi;%pi];
-checkroots(p,expectedroots,%eps^(1/3),5*%eps^(1/3));
-// 4 real roots with a zero derivate at the root
-// Really difficult problem : only simple precision
-p=(%s-%pi)^4;
-expectedroots  = [%pi;%pi;%pi;%pi];
-checkroots(p,expectedroots,%eps^(1/4),5*%eps^(1/4))
-// 10 real roots with a zero derivate at the root
-// Really difficult problem : only one correct digit
-p=(%s-%pi)^10;
-expectedroots  = [%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi];
-checkroots(p,expectedroots,%eps^(1/10),8*%eps^(1/10))
-
-// "Numerical computing with Matlab", Cleve Moler.
-A = diag(1:20);
-p = poly(A,'x');
-e = [1:20]';
-checkroots(p,e,%eps,0.2);
 
 // Tests from CPOLY
 // M. A. Jenkins and J. F. Traub. 1972.
@@ -255,7 +217,7 @@ E = sortRoots(E);
 R = sortRoots(R);
 assert_checkalmostequal(R, E, 1.e-3, 1.e-3);
 
-// EXAMPLE 5. 12 ZEROS EVENLY DISTRIBUTE ON A CIRCLE OF 
+// EXAMPLE 5. 12 ZEROS EVENLY DISTRIBUTE ON A CIRCLE OF
 // RADIUS 1. CENTERED AT 0+2I
 // Real part:
 // 4095 - 67584*x^2 + 126720*x^4 - 59136*x^6 + 7920*x^8 - 264*x^10 + x^12
@@ -311,6 +273,6 @@ E = sortRoots(E);
 R = sortRoots(R);
 assert_checkalmostequal(R, E, 1.e-10, 1.e-8);
 
-assert_checkequal(roots([4 3 2 1]), roots(poly([1 2 3 4], 'x', 'coeff')));
-assert_checkequal(roots([4 3 2 1] + [1 2 3 4]*%i), roots(poly([1 2 3 4]+[4 3 2 1]*%i,'x','coeff')));
+assert_checkequal(roots([4 3 2 1]), roots(poly([1 2 3 4], "x", "coeff")));
+assert_checkequal(roots([4 3 2 1] + [1 2 3 4]*%i), roots(poly([1 2 3 4]+[4 3 2 1]*%i,"x","coeff")));
 
diff --git a/scilab/modules/polynomials/tests/unit_tests/roots_difficult.tst b/scilab/modules/polynomials/tests/unit_tests/roots_difficult.tst
new file mode 100644 (file)
index 0000000..d3d4a58
--- /dev/null
@@ -0,0 +1,50 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008 - 2009 - INRIA - Michael Baudin
+// Copyright (C) 2011 - DIGITEO - Michael Baudin
+// Copyright (C) 2012 - INRIA - Serge Steer
+// Copyright (C) 2016 - Scilab Enterprises - Clement David
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+
+// <-- CLI SHELL MODE -->
+// <-- NOT FIXED YET -->
+
+//
+// Caution !
+// The following are difficult root-finding problems
+// with expected precision problems.
+// See "Principles for testing polynomial
+// zerofinding programs"
+// Jenkins, Traub
+// 1975
+// p.28
+// "The accuracy which one may expect to achieve in calculating
+// zeros is limited by the condition of these zeros. In particular,
+// for multiple zeros perturbations of size epsilon in the
+// coefficients cause perturbations of size epsilon^(1/m)
+// in the zeros."
+//
+//
+// 3 real roots with a zero derivate at the root
+// Really difficult problem : only simple precision computed, instead of double precision ***
+p=(%s-%pi)^3;
+expectedroots  = [%pi;%pi;%pi];
+checkroots(p,expectedroots,%eps^(1/3),5*%eps^(1/3));
+// 4 real roots with a zero derivate at the root
+// Really difficult problem : only simple precision
+p=(%s-%pi)^4;
+expectedroots  = [%pi;%pi;%pi;%pi];
+checkroots(p,expectedroots,%eps^(1/4),5*%eps^(1/4))
+// 10 real roots with a zero derivate at the root
+// Really difficult problem : only one correct digit
+p=(%s-%pi)^10;
+expectedroots  = [%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi;%pi];
+checkroots(p,expectedroots,%eps^(1/10),8*%eps^(1/10))
+
+// "Numerical computing with Matlab", Cleve Moler.
+A = diag(1:20);
+p = poly(A,"x");
+e = [1:20]';
+checkroots(p,e,%eps,0.2);
index 7981fcf..e350858 100644 (file)
@@ -1,11 +1,12 @@
 # Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 # Copyright (C) 2006 - INRIA - Sylvestre LEDRU
+# Copyright (C) 2016 - Scilab Enterprises - Clement DAVID
+#
 #
 # This file is distributed under the same license as the Scilab package.
 
 
 RANDLIB_C_SOURCES = \
-    src/c/fsultra.c \
     src/c/mt.c \
     src/c/igngeom.c \
     src/c/kiss.c \
index 7132194..5f58d98 100644 (file)
@@ -16,6 +16,8 @@
 
 # Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 # Copyright (C) 2006 - INRIA - Sylvestre LEDRU
+# Copyright (C) 2016 - Scilab Enterprises - Clement DAVID
+#
 #
 # This file is distributed under the same license as the Scilab package.
 
@@ -176,8 +178,7 @@ LTLIBRARIES = $(noinst_LTLIBRARIES) $(pkglib_LTLIBRARIES)
 libscirandlib_algo_la_LIBADD =
 am__dirstamp = $(am__leading_dot)dirstamp
 am__objects_1 = src/cpp/libscirandlib_algo_la-grand.lo
-am__objects_2 = src/c/libscirandlib_algo_la-fsultra.lo \
-       src/c/libscirandlib_algo_la-mt.lo \
+am__objects_2 = src/c/libscirandlib_algo_la-mt.lo \
        src/c/libscirandlib_algo_la-igngeom.lo \
        src/c/libscirandlib_algo_la-kiss.lo \
        src/c/libscirandlib_algo_la-urand.lo \
@@ -596,7 +597,6 @@ top_builddir = @top_builddir@
 top_srcdir = @top_srcdir@
 yacc_present = @yacc_present@
 RANDLIB_C_SOURCES = \
-    src/c/fsultra.c \
     src/c/mt.c \
     src/c/igngeom.c \
     src/c/kiss.c \
@@ -844,8 +844,6 @@ src/c/$(am__dirstamp):
 src/c/$(DEPDIR)/$(am__dirstamp):
        @$(MKDIR_P) src/c/$(DEPDIR)
        @: > src/c/$(DEPDIR)/$(am__dirstamp)
-src/c/libscirandlib_algo_la-fsultra.lo: src/c/$(am__dirstamp) \
-       src/c/$(DEPDIR)/$(am__dirstamp)
 src/c/libscirandlib_algo_la-mt.lo: src/c/$(am__dirstamp) \
        src/c/$(DEPDIR)/$(am__dirstamp)
 src/c/libscirandlib_algo_la-igngeom.lo: src/c/$(am__dirstamp) \
@@ -937,7 +935,6 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@sci_gateway/cpp/$(DEPDIR)/libscirandlib_la-sci_grand.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/c/$(DEPDIR)/libscirandlib_algo_la-clcg2.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/c/$(DEPDIR)/libscirandlib_algo_la-clcg4.Plo@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@src/c/$(DEPDIR)/libscirandlib_algo_la-fsultra.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/c/$(DEPDIR)/libscirandlib_algo_la-genbet.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/c/$(DEPDIR)/libscirandlib_algo_la-ignbin.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@src/c/$(DEPDIR)/libscirandlib_algo_la-igngeom.Plo@am__quote@
@@ -973,13 +970,6 @@ distclean-compile:
 @AMDEP_TRUE@@am__fastdepCC_FALSE@      DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCC_FALSE@  $(AM_V_CC@am__nodep@)$(LTCOMPILE) -c -o $@ $<
 
-src/c/libscirandlib_algo_la-fsultra.lo: src/c/fsultra.c
-@am__fastdepCC_TRUE@   $(AM_V_CC)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libscirandlib_algo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT src/c/libscirandlib_algo_la-fsultra.lo -MD -MP -MF src/c/$(DEPDIR)/libscirandlib_algo_la-fsultra.Tpo -c -o src/c/libscirandlib_algo_la-fsultra.lo `test -f 'src/c/fsultra.c' || echo '$(srcdir)/'`src/c/fsultra.c
-@am__fastdepCC_TRUE@   $(AM_V_at)$(am__mv) src/c/$(DEPDIR)/libscirandlib_algo_la-fsultra.Tpo src/c/$(DEPDIR)/libscirandlib_algo_la-fsultra.Plo
-@AMDEP_TRUE@@am__fastdepCC_FALSE@      $(AM_V_CC)source='src/c/fsultra.c' object='src/c/libscirandlib_algo_la-fsultra.lo' libtool=yes @AMDEPBACKSLASH@
-@AMDEP_TRUE@@am__fastdepCC_FALSE@      DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
-@am__fastdepCC_FALSE@  $(AM_V_CC@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libscirandlib_algo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o src/c/libscirandlib_algo_la-fsultra.lo `test -f 'src/c/fsultra.c' || echo '$(srcdir)/'`src/c/fsultra.c
-
 src/c/libscirandlib_algo_la-mt.lo: src/c/mt.c
 @am__fastdepCC_TRUE@   $(AM_V_CC)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libscirandlib_algo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT src/c/libscirandlib_algo_la-mt.lo -MD -MP -MF src/c/$(DEPDIR)/libscirandlib_algo_la-mt.Tpo -c -o src/c/libscirandlib_algo_la-mt.lo `test -f 'src/c/mt.c' || echo '$(srcdir)/'`src/c/mt.c
 @am__fastdepCC_TRUE@   $(AM_V_at)$(am__mv) src/c/$(DEPDIR)/libscirandlib_algo_la-mt.Tpo src/c/$(DEPDIR)/libscirandlib_algo_la-mt.Plo
index b9e2531..ebe691c 100644 (file)
                     <itemizedlist>
                         <listitem>
                             <para>
-                                <literal>[0, 2^32 - 1]</literal> for mt, kiss and fsultra;
+                                <literal>[0, 2^32 - 1]</literal> for mt and kiss;
                             </para>
                         </listitem>
                         <listitem>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>fsultra</term>
-                <listitem>
-                    <para>
-                        A Subtract-with-Borrow generator mixing with a congruential
-                        generator of Arif Zaman and George Marsaglia, period more than <literal>10^356</literal>,
-                        state given by an array of <literal>37</literal> integers (plus an index onto this array, a flag (<literal>0</literal> or <literal>1</literal>)
-                        and another integer).
-                    </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
                 <term>urand</term>
                 <listitem>
                     <para>
                     <para>
                         <code>S = grand("getgen")</code> returns the current base generator.
                         In this case <varname>S</varname> is
-                        a string among <literal>"mt"</literal>, <literal>"kiss"</literal>, <literal>"clcg2"</literal>, <literal>"clcg4"</literal>, <literal>"urand"</literal>, <literal>"fsultra"</literal>.
+                        a string among <literal>"mt"</literal>, <literal>"kiss"</literal>, <literal>"clcg2"</literal>, <literal>"clcg4"</literal>, <literal>"urand"</literal>.
                     </para>
                 </listitem>
             </varlistentry>
                 <listitem>
                     <para>
                         <code>grand("setgen",gen)</code> sets the current base generator to be <varname>gen</varname>
-                        a string among <literal>"mt"</literal>, <literal>"kiss"</literal>, <literal>"clcg2"</literal>, <literal>"clcg4"</literal>, <literal>"urand"</literal>, <literal>"fsultra"</literal>.
+                        a string among <literal>"mt"</literal>, <literal>"kiss"</literal>, <literal>"clcg2"</literal>, <literal>"clcg4"</literal>, <literal>"urand"</literal>.
                         Notice that this call returns the new current generator, i.e. <varname>gen</varname>.
                     </para>
                 </listitem>
                         <code>S = grand("getsd")</code> gets the current state (the current seeds) of the current base
                         generator ; <varname>S</varname> is given as a column vector (of integers) of dimension <literal>625</literal>
                         for mt (the first being an index in <literal>[1,624]</literal>), <literal>4</literal> for kiss, <literal>2</literal>
-                        for clcg2,  <literal>40</literal> for fsultra, <literal>4</literal> for clcg4
+                        for clcg2,  <literal>4</literal> for clcg4
                         (for this last one you get the current state of the current virtual generator) and <literal>1</literal>
                         for urand.
                     </para>
                                 </para>
                             </listitem>
                         </varlistentry>
-                        <varlistentry>
-                            <term>for fsultra</term>
-                            <listitem>
-                                <para>
-                                    <varname>S</varname> is a vector of integers of dim <literal>40</literal> (the first component
-                                    is an index and must be in <literal>[0,37]</literal>, the 2nd component is a flag (0 or 1), the 3rd component is
-                                    an integer in <literal>[1,2^32[</literal> and the 37 others integers are in <literal>[0,2^32[</literal>) ; a simpler (and recommended)
-                                    initialization may be done with two integers <varname>s1</varname> and <varname>s2</varname> in
-                                    <literal>[0,2^32[</literal>.
-                                </para>
-                            </listitem>
-                        </varlistentry>
                     </variablelist>
                 </listitem>
             </varlistentry>
index 4100c07..53321de 100644 (file)
                     <itemizedlist>
                         <listitem>
                             <para>
-                                <literal>[0, 2^32 - 1]</literal> for mt, kiss and fsultra;
+                                <literal>[0, 2^32 - 1]</literal> for mt and kiss;
                             </para>
                         </listitem>
                         <listitem>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>fsultra</term>
-                <listitem>
-                    <para>
-                        Un générateur SWB (subtract-with-borrow) mixé avec un générator congruentiel
-                        concu par Arif Zaman et George Marsaglia. Sa période est supérieure à <literal>10^356</literal>,
-                        et son état interne est constitué d'un tableau de 37 entiers, d'un index sur
-                        ce tableau et d'un drapeau (0 ou 1) ainsi qu'un autre entier donnant l'état interne
-                        du générateur congruentiel.
-                    </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
                 <term>urand</term>
                 <listitem>
                     <para>
                 <listitem>
                     <para>
                         <literal>S = grand("getgen")</literal> retourne le nom du générateur de base actuel (<literal>S</literal> est
-                        l'une des chaînes de caractères "mt", "kiss", "clcg2", "clcg4", "urand", "fsultra").
+                        l'une des chaînes de caractères "mt", "kiss", "clcg2", "clcg4", "urand").
                     </para>
                 </listitem>
             </varlistentry>
                 <listitem>
                     <para>
                         <literal>grand("setgen", gen)</literal> permet de changer le générateur de base : <literal>gen</literal>
-                        doit être l'une des chaînes de caractères "mt", "kiss", "clcg2", "clcg4", "urand", "fsultra".
+                        doit être l'une des chaînes de caractères "mt", "kiss", "clcg2", "clcg4", "urand".
                         En cas de succès la fonction retourne cette même chaîne.
                     </para>
                 </listitem>
                         <literal>S</literal> est un vecteur colonne (d'entiers) de dimension <literal>625</literal>
                         pour mt (la première composante étant un 'index' sur l'état, c-a-d un entier de l'intervalle
                         <literal>[1,624]</literal>), <literal>4</literal>
-                        pour kiss, <literal>2</literal> pour clcg2 , <literal>40</literal>pour fsultra, <literal>4</literal> pour clcg4
+                        pour kiss, <literal>2</literal> pour clcg2 , <literal>4</literal> pour clcg4
                         (pour ce dernier vous obtenez l'état interne du générateur virtuel courant), et <literal>1</literal>
                         pour urand.
                     </para>
                                 </para>
                             </listitem>
                         </varlistentry>
-                        <varlistentry>
-                            <term>for fsultra</term>
-                            <listitem>
-                                <para>
-                                    <literal>S</literal> est un vecteur de <literal>40</literal> entiers (son premier élément doit être dans
-                                    l'intervalle<literal>[0, 37]</literal>, son deuxième (drapeau) doit être 0 ou 1, le troisième un
-                                    entier de <literal>[1, 2^32[</literal> et les 37 composantes suivantes, des entiers de <literal>[0, 2^32[</literal>) ; il est recommandé
-                                    d'utiliser l'autre procédure d'initialisation (plus simple) avec deux entiers <literal>s1</literal> et
-                                    <literal>s2</literal> de <literal>[0, 2^32[</literal>.
-                                </para>
-                            </listitem>
-                        </varlistentry>
                     </variablelist>
                 </listitem>
             </varlistentry>
index bd739be..c4e3e11 100644 (file)
                     <itemizedlist>
                         <listitem>
                             <para>
-                                mt, kiss および fsultra の場合は <literal>[0, 2^32 - 1]</literal>
+                                mt および kiss の場合は <literal>[0, 2^32 - 1]</literal>
                             </para>
                         </listitem>
                         <listitem>
                     </para>
                 </listitem>
             </varlistentry>
-            <varlistentry>
-                <term>fsultra</term>
-                <listitem>
-                    <para>
-                        Arif Zaman および George Marsagliaの合同生成器に基づく
-                        Subtract-with-Borrow 生成器,
-                        周期は <literal>10^356</literal>以上,
-                        状態は37個の整数の配列(およびこの配列へのインデックス,
-                        フラグ (0または1)および他の整数)により指定されます.
-                    </para>
-                </listitem>
-            </varlistentry>
         </variablelist>
         <para>全ての生成器に共通なアクションの動作を以下に示します:
         </para>
                     <para>
                         <literal>S=grand('getgen')</literal> はカレントの基底生成器を返します
                         ( <literal>S</literal> は以下の文字列のどれかです:
-                        'mt', 'kiss', 'clcg2', 'clcg4', 'urand', 'fsultra'.
+                        'mt', 'kiss', 'clcg2', 'clcg4', 'urand'.
                     </para>
                 </listitem>
             </varlistentry>
                     <para>
                         <literal>grand('setgen',gen)</literal>は
                         カレントの基底生成器を<literal>gen</literal>に設定します.
-                        この文字列には 'mt', 'kiss', 'clcg2', 'clcg4', 'urand', 'fsultra'
+                        この文字列には 'mt', 'kiss', 'clcg2', 'clcg4', 'urand'
                         のどれかを指定します(
                         このコールは新しいカレントの生成器を返すことに注意してください).
                     </para>
                         <literal>S</literal> には, mt の場合は<literal>625</literal>次
                         (先頭は<literal>[1,624]</literal>の範囲のインデックスです),
                         kiss の場合は <literal>4</literal> 次, clcg2の場合は <literal>2</literal>次,
-                        fsultraの場合は <literal>40</literal> , clcg4 の場合は <literal>4</literal>次
+                        clcg4 の場合は <literal>4</literal>次
                         (最後の1つについてはカレントの仮想生成器のカレントの状態が取得されます),
                         および urand の場合は <literal>1</literal>次
                         の(整数)列ベクトルが出力されます.
                                 </para>
                             </listitem>
                         </varlistentry>
-                        <varlistentry>
-                            <term>fsultraの場合</term>
-                            <listitem>
-                                <para>
-                                    <literal>S</literal> は<literal>40</literal>次の整数ベクトルで,
-                                    (最初の要素はインデックスで<literal>[0,37]</literal>の範囲とします,
-                                    2番目の要素はフラグ (0または1),3番目は[1,2^32[ の範囲の整数,
-                                    367個のその他の整数(範囲: [0,2^32[)));
-                                    <literal>[0,2^32[</literal>の範囲の
-                                    整数を2つだけ (<literal>s1</literal> および <literal>s2</literal>)
-                                    指定することでより簡単に(そして推奨される)初期化を行うことができます.
-                                </para>
-                            </listitem>
-                        </varlistentry>
                     </variablelist>
                 </listitem>
             </varlistentry>
index 3c76886..0f8456a 100644 (file)
                     <itemizedlist>
                         <listitem>
                             <para>
-                                <literal>[0, 2^32 - 1]</literal>  для mt, kiss и fsultra;
+                                <literal>[0, 2^32 - 1]</literal>  для mt и kiss;
                             </para>
                         </listitem>
                         <listitem>
                 </listitem>
             </varlistentry>
             <varlistentry>
-                <term>fsultra</term>
-                <listitem>
-                    <para>
-                        Генератор вычитания с займом (Subtract-with-Borrow), смешанный с
-                        конгруэнтным генератором Арифа Замана (Arif Zaman) и Джорджа
-                        Марсальи (George Marsaglia), с периодом свыше
-                        <literal>10^356</literal>, состояние задаётся массивом из
-                        <literal>37</literal> целых чисел (плюс индекс на этот массив,
-                        флаг (<literal>0</literal> или <literal>1</literal>) и другое
-                        целое число).
-                    </para>
-                </listitem>
-            </varlistentry>
-            <varlistentry>
                 <term>urand</term>
                 <listitem>
                     <para>
                         генератор. В данном случае <varname>S</varname> является одной
                         строкой из <literal>"mt"</literal>, <literal>"kiss"</literal>,
                         <literal>"clcg2"</literal>, <literal>"clcg4"</literal>,
-                        <literal>"urand"</literal>, <literal>"fsultra"</literal>.
+                        <literal>"urand"</literal>.
                     </para>
                 </listitem>
             </varlistentry>
                         генератор. В данном случае <varname>gen</varname> может быть одной
                         строкой из <literal>"mt"</literal>, <literal>"kiss"</literal>,
                         <literal>"clcg2"</literal>, <literal>"clcg4"</literal>,
-                        <literal>"urand"</literal>, <literal>"fsultra"</literal>.
+                        <literal>"urand"</literal>.
                         Заметьте, что этот вызов возвращает новый текущий генератор, т.е.
                         <varname>gen</varname>.
                     </para>
                         <literal>"mt"</literal> (первое значение является индексом в
                         интервале <literal>[1,624]</literal>),  <literal>4</literal> для
                         <literal>"kiss"</literal>,  <literal>2</literal> для
-                        <literal>"clcg2"</literal>, <literal>40</literal> для
-                        <literal>"fsultra"</literal>,  <literal>4</literal> для
+                        <literal>"clcg2"</literal>, <literal>4</literal> для
                         <literal>"clcg4"</literal> (для последнего вы получите текущее
                         состояние текущего виртуального генератора) и
                         <literal>1</literal> для <literal>"urand"</literal>.
                                 </para>
                             </listitem>
                         </varlistentry>
-                        <varlistentry>
-                            <term>для fsultra</term>
-                            <listitem>
-                                <para>
-                                    <varname>S</varname> является вектором целых чисел, состоящий из 40 элементов  (первый элемент является
-                                    индексом и должен быть на интервале <literal>[0,37]</literal>, второй элемент является
-                                    флагом (0 или 1), третий элемент является целым числом на интервале  <literal>[1,2^32[</literal>, а 37
-                                    остальных целых чисел должны быть на интервале <literal>[0,2^32[</literal>. Более простая (и
-                                    рекомендованная) инициализация может быть сделана с помощью двух целых чисел <varname>s1</varname> и
-                                    <varname>s2</varname> на интервале <literal>[0,2^32[</literal>.
-                                </para>
-                            </listitem>
-                        </varlistentry>
                     </variablelist>
                 </listitem>
             </varlistentry>
index b591388..c1ef35a 100644 (file)
@@ -40,12 +40,6 @@ RANDLIB_IMPEXP unsigned long int urandc(void);
 RANDLIB_IMPEXP int set_state_urand(double g);
 RANDLIB_IMPEXP void get_state_urand(double g[]);
 
-/* header for scilab fsultra */
-RANDLIB_IMPEXP unsigned long int fsultra(void);
-RANDLIB_IMPEXP int set_state_fsultra(double g[]);
-RANDLIB_IMPEXP int set_state_fsultra_simple(double g1, double g2);
-RANDLIB_IMPEXP void get_state_fsultra(double g[]);
-
 #endif /** SCI_OTHER_GEN   **/
 
 
index 030bfc8..ce74026 100644 (file)
@@ -14,13 +14,6 @@ and continues to be available under such terms.
 For more information, see the COPYING file which you should have received
 along with this program.
 
-Non-free source files:
-======================
-
-With "commercial use restriction"
-
-File: modules/randlib/src/c/fsultra.c
-
 mt:
 ===
 
index fc5f5f4..aef63c8 100644 (file)
@@ -43,11 +43,11 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
 {
     enum
     {
-        MT, KISS, CLCG4, CLCG2, URAND, FSULTRA
+        MT, KISS, CLCG4, CLCG2, URAND
     };
 
     //  names at the scilab level
-    const wchar_t* names_gen[6] = { L"mt", L"kiss", L"clcg4", L"clcg2", L"urand", L"fsultra" };
+    const wchar_t* names_gen[6] = { L"mt", L"kiss", L"clcg4", L"clcg2", L"urand" };
 
     types::String* pStrMethod = NULL;
     types::String* pStrGenOrPhr = NULL;
@@ -238,20 +238,6 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
             case URAND:
                 iNumInputArg = 2;
                 break;
-            case FSULTRA:
-            {
-                if (in.size() != 2 && in.size() != 3)
-                {
-                    char* pstMeth = wide_string_to_UTF8(wcsMeth);
-                    Scierror(77, _("%s: Wrong number of input argument(s) for method %s: %d or %d expected.\n"), "grand", pstMeth, 2, 3);
-                    FREE(pstMeth);
-                    delete[] itab;
-                    return types::Function::Error;
-                }
-
-                iNumInputArg = (int)in.size();
-                break;
-            }
         }
         meth = 24;
     }
@@ -1275,13 +1261,9 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
             {
                 ConfigVariable::setCurrentBaseGen(URAND);
             }
-            else if (wcscmp(wcsGen, L"fsultra") == 0)
-            {
-                ConfigVariable::setCurrentBaseGen(FSULTRA);
-            }
             else
             {
-                Scierror(999, _("%s: Wrong value for input argument #%d: '%s', '%s', '%s', '%s', '%s' or '%s' expected.\n"), "grand", 2, "mt", "kiss", "clcg4", "clcg2", "urand", "fsultra");
+                Scierror(999, _("%s: Wrong value for input argument #%d: '%s', '%s', '%s', '%s', '%s' or '%s' expected.\n"), "grand", 2, "mt", "kiss", "clcg4", "clcg2", "urand");
                 return types::Function::Error;
             }
 
@@ -1325,12 +1307,6 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
                     get_state_urand(pDblOut->get());
                     break;
                 }
-                case FSULTRA:
-                {
-                    pDblOut = new types::Double(40, 1);
-                    get_state_fsultra(pDblOut->get());
-                    break;
-                }
             }
 
             out.push_back(pDblOut);
@@ -1415,34 +1391,6 @@ types::Function::ReturnValue sci_grand(types::typed_list &in, int _iRetCount, ty
                     ierr = set_state_urand(vectpDblInput[0]->get(0));
                     break;
                 }
-                case FSULTRA:
-                {
-                    if (in.size() == 2)
-                    {
-                        if (vectpDblInput[0]->getRows() != 40 || vectpDblInput[0]->getCols() != 1)
-                        {
-                            Scierror(999, _("%s: Wrong size for input argument #%d : A vector of size %d x %d expected.\n"), "grand", 2, 40, 1);
-                            return types::Function::Error;
-                        }
-
-                        ierr = set_state_fsultra(vectpDblInput[0]->get());
-                    }
-                    else // in.size() == 3
-                    {
-                        for (int i = 0; i < 2; i++)
-                        {
-                            if (vectpDblInput[i]->isScalar() == false)
-                            {
-                                Scierror(999, _("%s: Wrong type for input argument #%d : A scalar expected.\n"), "grand", i + 2);
-                                return types::Function::Error;
-                            }
-                        }
-
-                        ierr = set_state_fsultra_simple(vectpDblInput[0]->get(0), vectpDblInput[1]->get(0));
-                    }
-
-                    break;
-                }
             }
 
             if (ierr == 0)
diff --git a/scilab/modules/randlib/src/c/fsultra.c b/scilab/modules/randlib/src/c/fsultra.c
deleted file mode 100644 (file)
index bf72e7a..0000000
+++ /dev/null
@@ -1,257 +0,0 @@
-/*
- *
- * Copyright (C) 2010 - DIGITEO - Michael Baudin
- * Copyright (C) 2004 - Bruno Pincon
- * Copyright (C) 1992 - Arif Zaman
- * Copyright (C) 1992 - George Marsaglia
- *
-FSU - ULTRA    The greatest random number generator that ever was
-               or ever will be.  Way beyond Super-Duper.
-               (Just kidding, but we think its a good one.)
-
-Authors:       Arif Zaman (arif@stat.fsu.edu) and
-               George Marsaglia (geo@stat.fsu.edu).
-
-Date:          27 May 1992
-
-Version:       1.05
-
-Copyright:     To obtain permission to incorporate this program into
-               any commercial product, please contact the authors at
-               the e-mail address given above or at
-
-               Department of Statistics and
-               Supercomputer Computations Research Institute
-               Florida State University
-               Tallahassee, FL 32306.
-
-See Also:      README          for a brief description
-               ULTRA.DOC       for a detailed description
-
------------------------------------------------------------------------
-*/
-/*
-   File: ULTRA.C
-
-   This is the ULTRA random number generator written entirely in C.
-
-   This may serve as a model for an assembler version of this routine.
-   The programmer should avoid simply duplicating and instead use the
-   usual assembler features to increase the speed of this routine.
-
-   Especially the subroutine SWB should be replaced by the one
-   machine instruction (usually called subtract-with-borrow) that
-   is available in almost every hardware.
-
-   For people not familiar with 8086 assembler, it may help to
-   consult this when reading the assembler code. This program should
-   be a dropin replacement for the assembler versions, but is about
-   half as fast.
-*/
-
-/* Slight modifications by Bruno Pincon (4 december 2004) for inclusion
-   in scilab and nsp:
-
-   1/ in scilab we use only i32bit output ( renamed here fsultra )
-      and  I have deleted the others;
-
-   2/ only one array is now used (swbseed which is renamed
-      swb_state) and the xor with the linear congruential generator
-      is done only just before the output.
-
-   3/ add a var is_init (to say if the generator is initialised)
-
-   4/ add routine to set/get the state
-
-*/
-#include <math.h>               /* to use floor    */
-#include "others_generators.h"
-#include "sciprint.h"
-#include "localization.h"
-
-#define N  37           /* size of table        */
-#define N2 24           /* The shorter lag      */
-
-static int is_init = 0;
-static int swb_state[N];          /* state of the swb generator */
-static int swb_index = N;          /* an index on the swb state */
-static int swb_flag;              /* the carry flag for the SWB generator */
-static unsigned int cong_state;   /* state of the congruential generator */
-
-/* for this generator the state seems completly defined by:
-swb_state, swb_index, swb_flag (which define the state of the swb generator)
-cong_state (which defines the state of the congruential generator)
-*/
-
-/* those are the default for the simple initialisation routine */
-static  double DEFAULT_SEED1 = 1234567.0, DEFAULT_SEED2 = 7654321.0;
-
-
-
-
-/* SWB is the subtract-with-borrow operation which should be one line
-in assembler code. This should be done by using the hardware s-w-b
-operation in the SWBfill routine (renamed advance_state_swb here).
-
-What has been done here is to look at the msb of x, y and z=x-y-c.
-Using these three bits, one can determine if a borrow bit is needed
-or not according to the following table:
-
-msbz=0  msby=0  msby=1          msbz=1  msby=0  msby=1
-
-msbx=0  0       1               msbx=0  1       1
-msbx=1  0       0               msbx=1  0       1
-
-PS: note that the definition is very carefully written because the
-calls to SWB have y and z as the same memory location, so y must
-be tested before z is assigned a value.
-*/
-#define SWB(c,x,y,z) c = (y<0) ? (((z=x-y-c) < 0) || (x>=0)) : (((z=x-y-c) < 0) && (x>=0));
-
-static int advance_state_swb(void)
-{
-    int i;
-    /*
-    *  The following are the heart of the system and should be
-    *  written is assembler to be as fast as possible. It may even make sense
-    *  to unravel the loop and simply wirte 37 consecutive SWB operations!
-    */
-    for (i = 0;  i < N2; i++)
-    {
-        SWB(swb_flag, swb_state[i + N - N2], swb_state[i], swb_state[i]);
-    }
-    for (i = N2; i < N;  i++)
-    {
-        SWB(swb_flag, swb_state[i  - N2], swb_state[i], swb_state[i]);
-    }
-    return swb_index = 0;
-}
-
-/* set_state_fsultra_simple initializes the state from 2 integers
-
-it defines the constants and fills the swb_state array one bit at
-a time by taking the leading bit of the xor of a shift register
-and a congruential sequence. The same congruential generator continues
-to be used as a mixing generator for the Subtract-with-borrow generator
-to produce the `ultra' random numbers
-
-Since this is called just once, speed doesn't matter much and it might
-be fine to leave this subroutine coded just as it is.
-
-PS:    there are quick and easy ways to fill this, but since random number
-generators are really "randomness amplifiers", it is important to
-start off on the right foot. This is why we take such care here.
-*/
-
-int set_state_fsultra_simple(double s1, double s2)
-{
-    unsigned int shrgx, tidbits = 0;
-    int i, j;
-
-    if (    (s1 == floor(s1) && 0.0 <= s1 && s1 <= 4294967295.0)
-            && (s2 == floor(s2) && 0.0 <= s2 && s2 <= 4294967295.0) )
-    {
-        cong_state = ((unsigned int) s1) * 2 + 1;
-        shrgx = (unsigned int) s2;
-        for ( i = 0 ; i < N ; i++)
-        {
-            for ( j = 32 ; j > 0 ; j--)
-            {
-                cong_state = cong_state * 69069;
-                shrgx = shrgx ^ (shrgx >> 15);
-                shrgx = shrgx ^ (shrgx << 17);
-                tidbits = (tidbits >> 1) | (0x80000000 & (cong_state ^ shrgx));
-            }
-            swb_state[i] = tidbits;
-        }
-        swb_index = 0;
-        swb_flag = 0;
-        swb_index = advance_state_swb();  /* pour retrouver la même séquence que ds scilab V3.0 */
-        is_init = 1;
-        return swb_index;
-    }
-    else
-    {
-        sciprint(_("\nBad seed for fsultra, must be integers in [0, 2^32-1]\n"));
-        return -1;
-    }
-}
-
-int set_state_fsultra(double *s)
-{
-    double try;
-    int i;
-
-    try = s[0];
-    if ( floor(try) != try || try < 0.0  ||  try > (double) N)
-                {
-                    sciprint(_("\nThe first component of the fsultra state, must be an int in [0, %d]\n"), N);
-                    return 0;
-                }
-    swb_index = (int) try;
-
-    try = s[1];
-    if ( try != 0.0  &&  try != 1.0)
-            {
-                sciprint(_("\nThe second component of the fsultra state, must be 0 or 1\n"));
-                return 0;
-            }
-    swb_flag = (int) try;
-
-    try = s[2];
-    if ( floor(try) != try  ||  try <= 0 ||  try > 4294967295.0 )
-                {
-                    sciprint(_("\nThe third component of the fsultra state, must be an int in [1, 2^32-1]\n"));
-                    return 0;
-                }
-    cong_state = (unsigned int) try;
-
-    /* no verif here ... */
-    for (i = 0 ; i < N ; i++)
-    {
-        swb_state[i] = (int) (((unsigned int) s[i + 3]) & 0xffffffff);
-    }
-
-    is_init = 1;
-    return 1;
-}
-
-
-/*  to return the state at the scilab level  */
-void get_state_fsultra(double s[])
-{
-    int i;
-
-    if ( ! is_init )
-    {
-        set_state_fsultra_simple(DEFAULT_SEED1, DEFAULT_SEED2);
-    }
-
-    s[0] = (double)  swb_index;
-    s[1] = (double)  swb_flag;
-    s[2] = (double)  cong_state;
-    for (i = 0 ; i < N ; i++)
-    {
-        s[i + 3] = (double) (unsigned int) swb_state[i];
-    }
-}
-
-unsigned long int fsultra(void)
-{
-    if (swb_index >= N)  /* generate N words at one time */
-    {
-        if ( ! is_init )
-        {
-            swb_index = set_state_fsultra_simple(DEFAULT_SEED1, DEFAULT_SEED2);
-            if (swb_index < 0)
-            {
-                return 0;
-            }
-        }
-        else
-        {
-            swb_index = advance_state_swb();
-        }
-    }
-    return (swb_state[swb_index++] ^ (cong_state = cong_state * 69069));
-}
index 6694720..690ab20 100644 (file)
@@ -205,7 +205,6 @@ lib /DEF:"$(ProjectDir)elementary_functions_f_Import.def" /SUBSYSTEM:WINDOWS /MA
     <ClCompile Include="clcg2.c" />
     <ClCompile Include="clcg4.c" />
     <ClCompile Include="DllmainRandlib.c" />
-    <ClCompile Include="fsultra.c" />
     <ClCompile Include="genbet.c" />
     <ClCompile Include="ignbin.c" />
     <ClCompile Include="igngeom.c" />
@@ -261,4 +260,4 @@ lib /DEF:"$(ProjectDir)elementary_functions_f_Import.def" /SUBSYSTEM:WINDOWS /MA
   <Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
   <ImportGroup Label="ExtensionTargets">
   </ImportGroup>
-</Project>
\ No newline at end of file
+</Project>
index dd03944..23671a6 100644 (file)
@@ -32,9 +32,6 @@
     <ClCompile Include="DllmainRandlib.c">
       <Filter>Source Files</Filter>
     </ClCompile>
-    <ClCompile Include="fsultra.c">
-      <Filter>Source Files</Filter>
-    </ClCompile>
     <ClCompile Include="genbet.c">
       <Filter>Source Files</Filter>
     </ClCompile>
       <Filter>Resource Files</Filter>
     </ResourceCompile>
   </ItemGroup>
-</Project>
\ No newline at end of file
+</Project>
index 17c85ff..6b8c307 100644 (file)
@@ -34,7 +34,7 @@ unsigned long int clcg4_with_gen(void)
 double C2F(ranf)(void)
 {
     //  pointers onto the generators func
-    unsigned long int (*gen[NbGenInScilab])() = {randmt, kiss, clcg4_with_gen, clcg2, urandc, fsultra};
+    unsigned long int (*gen[NbGenInScilab])() = {randmt, kiss, clcg4_with_gen, clcg2, urandc};
 
     int current_gen = ConfigVariable::getCurrentBaseGen();
 
@@ -45,9 +45,8 @@ double C2F(ranf)(void)
         2.3283064365386963e-10,  // kiss
         4.6566128752457969e-10,  // clcg4
         4.6566130595601735e-10,  // clcg2
-        4.6566128730773926e-10,  // urand
-        2.3283064365386963e-10
-    }; // fsultra
+        4.6566128730773926e-10   // urand
+    };
 
     // random deviate from U[0,1)
     return ((double)gen[current_gen]() * factor[current_gen]);
@@ -56,7 +55,7 @@ double C2F(ranf)(void)
 double ignlgi(void)
 {
     //  pointers onto the generators func
-    unsigned long int (*gen[NbGenInScilab])() = {randmt, kiss, clcg4_with_gen, clcg2, urandc, fsultra};
+    unsigned long int (*gen[NbGenInScilab])() = {randmt, kiss, clcg4_with_gen, clcg2, urandc};
 
     int current_gen = ConfigVariable::getCurrentBaseGen();
 
@@ -93,9 +92,8 @@ double C2F(ignuin)(double *a, double *b)
         4294967295ul,  // kiss
         2147483646ul,  // clcg4
         2147483561ul,  // clcg2
-        2147483647ul,  // urand
-        4294967295ul
-    }; // fsultra
+        2147483647ul   // urand
+    };
 
     int current_gen = ConfigVariable::getCurrentBaseGen();
 
index 708f8af..8f53df6 100644 (file)
 // <-- CLI SHELL MODE -->
 // <-- ENGLISH IMPOSED -->
 function checkMeanVariance2arg ( m , n , name , A , B , mu , va , rtol )
-  // Check the mean and variance of random numbers.
-  //
-  // Parameters
-  // m : a 1-by-1 matrix of floating point integers, the number of rows
-  // n : a 1-by-1 matrix of floating point integers, the number of columns
-  // name: a 1-by-1 string, the name of the distribution function
-  // A : a 1-by-1 matrix of doubles, the value of the 1st parameter
-  // B : a 1-by-1 matrix of doubles, the value of the 2nd parameter
-  // mu : a 1-by-1 matrix of doubles, the expected mean
-  // va : a 1-by-1 matrix of doubles, the expected variance. If va==[], then the variance is not checked.
-  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
-  
-  R=grand(m,n,name,A,B);
-  assert_checkequal ( size(R) , [m,n] );
-  assert_checkequal ( typeof(R) , "constant" );
-  assert_checkalmostequal ( mean(R) , mu , rtol );
-  if ( va<>[] ) then
-    assert_checkalmostequal ( variance(R) , va , rtol );
-  end
+    // Check the mean and variance of random numbers.
+    //
+    // Parameters
+    // m : a 1-by-1 matrix of floating point integers, the number of rows
+    // n : a 1-by-1 matrix of floating point integers, the number of columns
+    // name: a 1-by-1 string, the name of the distribution function
+    // A : a 1-by-1 matrix of doubles, the value of the 1st parameter
+    // B : a 1-by-1 matrix of doubles, the value of the 2nd parameter
+    // mu : a 1-by-1 matrix of doubles, the expected mean
+    // va : a 1-by-1 matrix of doubles, the expected variance. If va==[], then the variance is not checked.
+    // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+    R=grand(m,n,name,A,B);
+    assert_checkequal ( size(R) , [m,n] );
+    assert_checkequal ( typeof(R) , "constant" );
+    assert_checkalmostequal ( mean(R) , mu , rtol );
+    if ( va<>[] ) then
+        assert_checkalmostequal ( variance(R) , va , rtol );
+    end
 endfunction
 function checkMeanVariance0arg ( m , n , name , mu , va , rtol )
-  // Check the mean and variance of random numbers.
-  //
-  // Parameters
-  // m : a 1-by-1 matrix of floating point integers, the number of rows
-  // n : a 1-by-1 matrix of floating point integers, the number of columns
-  // name: a 1-by-1 string, the name of the distribution function
-  // mu : a 1-by-1 matrix of doubles, the expected mean
-  // va : a 1-by-1 matrix of doubles, the expected variance. If va==[], then the variance is not checked.
-  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
-  
-  R=grand(m,n,name);
-  assert_checkequal ( size(R) , [m,n] );
-  assert_checkequal ( typeof(R) , "constant" );
-  assert_checkalmostequal ( mean(R) , mu , rtol );
-  if ( va<>[] ) then
-    assert_checkalmostequal ( variance(R) , va , rtol );
-  end
+    // Check the mean and variance of random numbers.
+    //
+    // Parameters
+    // m : a 1-by-1 matrix of floating point integers, the number of rows
+    // n : a 1-by-1 matrix of floating point integers, the number of columns
+    // name: a 1-by-1 string, the name of the distribution function
+    // mu : a 1-by-1 matrix of doubles, the expected mean
+    // va : a 1-by-1 matrix of doubles, the expected variance. If va==[], then the variance is not checked.
+    // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+    R=grand(m,n,name);
+    assert_checkequal ( size(R) , [m,n] );
+    assert_checkequal ( typeof(R) , "constant" );
+    assert_checkalmostequal ( mean(R) , mu , rtol );
+    if ( va<>[] ) then
+        assert_checkalmostequal ( variance(R) , va , rtol );
+    end
 endfunction
 function checkLaw0arg ( name , cdffun , N , NC , rtol )
-  //
-  // Check the random number generator for a continuous distribution function.
-  //
-  // Parameters
-  // name: a 1-by-1 string, the name of the distribution function
-  // cdffun : a function, the Cumulated Distribution Function
-  // NC : a 1-by-1 matrix of floating point integers, the number of classes
-  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
-  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
-  //
-  // Description
-  //  Compare the empirical histogram with the theoretical histogram.
-  R = grand(1,N,name);
-  [X,EmpiricalPDF] = histcompute(NC,R);
-  CDF = cdffun(X)
-  TheoricPDF = diff(CDF);
-  assert_checktrue ( abs(EmpiricalPDF-TheoricPDF) < rtol );
+    //
+    // Check the random number generator for a continuous distribution function.
+    //
+    // Parameters
+    // name: a 1-by-1 string, the name of the distribution function
+    // cdffun : a function, the Cumulated Distribution Function
+    // NC : a 1-by-1 matrix of floating point integers, the number of classes
+    // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+    // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+    //
+    // Description
+    //  Compare the empirical histogram with the theoretical histogram.
+    R = grand(1,N,name);
+    [X,EmpiricalPDF] = histcompute(NC,R);
+    CDF = cdffun(X)
+    TheoricPDF = diff(CDF);
+    assert_checktrue ( abs(EmpiricalPDF-TheoricPDF) < rtol );
     if ( %f ) then
-      plot(X(1:$-1),EmpiricalPDF,"bo-"); // Empirical Histogram
-      plot(X(1:$-1),TheoricPDF,"rox-"); // Theoretical Histogram
+        plot(X(1:$-1),EmpiricalPDF,"bo-"); // Empirical Histogram
+        plot(X(1:$-1),TheoricPDF,"rox-"); // Theoretical Histogram
     end
 endfunction
 function checkPieceLaw0arg ( name , cdffun , N , NC , rtol )
-  //
-  // Check the random number generator for a piecewise distribution function.
-  //
-  // Parameters
-  // name: a 1-by-1 string, the name of the distribution function
-  // cdffun : a function, the Cumulated Distribution Function
-  // NC : a 1-by-1 matrix of floating point integers, the number of classes
-  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
-  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
-  //
-  // Description
-  //  Compare the empirical histogram with the theoretical histogram.
-  //  The classes of the histogram are computed from the random samples, 
-  //  from which the unique entries are extracted.
-  R=grand(N,1,name);
-  X = unique(R);
-  for k = 1 : size(X,"*")
-    EmpiricalPDF(k) = length(find(R==X(k)));
-  end
-  EmpiricalPDF = EmpiricalPDF./N;
-  CDF = cdffun(X)
-  TheoricPDF=[CDF(1);diff(CDF)];
-  assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
-  if ( %f ) then
-    plot(X,EmpiricalPDF,"bo-"); // Empirical Histogram
-    plot(X,TheoricPDF,"rox-"); // Theoretical Histogram
-  end
+    //
+    // Check the random number generator for a piecewise distribution function.
+    //
+    // Parameters
+    // name: a 1-by-1 string, the name of the distribution function
+    // cdffun : a function, the Cumulated Distribution Function
+    // NC : a 1-by-1 matrix of floating point integers, the number of classes
+    // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+    // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+    //
+    // Description
+    //  Compare the empirical histogram with the theoretical histogram.
+    //  The classes of the histogram are computed from the random samples,
+    //  from which the unique entries are extracted.
+    R=grand(N,1,name);
+    X = unique(R);
+    for k = 1 : size(X,"*")
+        EmpiricalPDF(k) = length(find(R==X(k)));
+    end
+    EmpiricalPDF = EmpiricalPDF./N;
+    CDF = cdffun(X)
+    TheoricPDF=[CDF(1);diff(CDF)];
+    assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
+    if ( %f ) then
+        plot(X,EmpiricalPDF,"bo-"); // Empirical Histogram
+        plot(X,TheoricPDF,"rox-"); // Theoretical Histogram
+    end
 endfunction
 function checkPieceLaw2arg ( name , cdffun , N , NC , A , B , rtol )
-  //
-  // Check the random number generator for a piecewise distribution function.
-  //
-  // Parameters
-  // name: a 1-by-1 string, the name of the distribution function
-  // cdffun : a function, the Cumulated Distribution Function
-  // NC : a 1-by-1 matrix of floating point integers, the number of classes
-  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
-  // A : a 1-by-1 matrix of doubles, the value of the parameter
-  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
-  //
-  // Description
-  //  Compare the empirical histogram with the theoretical histogram.
-  //  The classes of the histogram are computed from the random samples, 
-  //  from which the unique entries are extracted.
-  R=grand(N,1,name,A,B);
-  X = unique(R);
-  for k = 1 : size(X,"*")
-    EmpiricalPDF(k) = length(find(R==X(k)));
-  end
-  EmpiricalPDF = EmpiricalPDF./N;
-  CDF = cdffun(X,A,B)
-  TheoricPDF=[CDF(1);diff(CDF)];
-  assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
-  if ( %f ) then
-    plot(X,EmpiricalPDF,"bo-"); // Empirical Histogram
-    plot(X,TheoricPDF,"rox-"); // Theoretical Histogram
-  end
+    //
+    // Check the random number generator for a piecewise distribution function.
+    //
+    // Parameters
+    // name: a 1-by-1 string, the name of the distribution function
+    // cdffun : a function, the Cumulated Distribution Function
+    // NC : a 1-by-1 matrix of floating point integers, the number of classes
+    // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+    // A : a 1-by-1 matrix of doubles, the value of the parameter
+    // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+    //
+    // Description
+    //  Compare the empirical histogram with the theoretical histogram.
+    //  The classes of the histogram are computed from the random samples,
+    //  from which the unique entries are extracted.
+    R=grand(N,1,name,A,B);
+    X = unique(R);
+    for k = 1 : size(X,"*")
+        EmpiricalPDF(k) = length(find(R==X(k)));
+    end
+    EmpiricalPDF = EmpiricalPDF./N;
+    CDF = cdffun(X,A,B)
+    TheoricPDF=[CDF(1);diff(CDF)];
+    assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
+    if ( %f ) then
+        plot(X,EmpiricalPDF,"bo-"); // Empirical Histogram
+        plot(X,TheoricPDF,"rox-"); // Theoretical Histogram
+    end
 endfunction
 function [x,y] = histcompute(n,data)
-   //
-   // Computes the histogram of a data.
-   //
-   // Parameters
-   // n : a 1-by-1 matrix of floating point integers, the number of classes.
-   // data : a matrix of doubles, the data
-   // x : 1-by-n+1 matrix of doubles, the classes of the histogram
-   // y : 1-by-n+1 matrix of doubles, the empirical probability that one value in data which fall in each class
-   x = linspace(min(data), max(data), n+1);
-   [ind , y] = dsearch(data, x)
-   y = y./length(data)
+    //
+    // Computes the histogram of a data.
+    //
+    // Parameters
+    // n : a 1-by-1 matrix of floating point integers, the number of classes.
+    // data : a matrix of doubles, the data
+    // x : 1-by-n+1 matrix of doubles, the classes of the histogram
+    // y : 1-by-n+1 matrix of doubles, the empirical probability that one value in data which fall in each class
+    x = linspace(min(data), max(data), n+1);
+    [ind , y] = dsearch(data, x)
+    y = y./length(data)
 endfunction
 function p = mycdfdef (X)
-  p = X
+    p = X
 endfunction
 function p = mycdfuin (X,A,B)
-  p = (floor(X)-A+1)/(B-A+1)
+    p = (floor(X)-A+1)/(B-A+1)
 endfunction
 function p = mycdflgi (X)
-  // Here, A and B are read from the environment
-  p = (floor(X)-A+1)/(B-A+1)
+    // Here, A and B are read from the environment
+    p = (floor(X)-A+1)/(B-A+1)
 endfunction
 //
 // These tests makes checks that all uniform generators are correct.
@@ -160,11 +158,11 @@ endfunction
 // MinInt : minimum of the uniform integer interval for random number generation
 // MaxInt : maximum of the uniform integer interval for random number generation
 //
-NGen = 6;
-generators = ["mt"   "kiss" "clcg2"    "clcg4" "fsultra" "urand"];
-seedsize =   [625    4      2          4       40        1];
-MinInt =     [0      0      0          0       0         0];
-MaxInt =     [2^32-1 2^32-1 2147483561 2^31-2  2^32-1    2^31-1];
+NGen = 5;
+generators = ["mt"   "kiss" "clcg2"    "clcg4" "urand"];
+seedsize =   [625    4      2          4       1];
+MinInt =     [0      0      0          0       0];
+MaxInt =     [2^32-1 2^32-1 2147483561 2^31-2  2^31-1];
 rtol = 1.e-2;
 // The number of classes in the histogram
 // NC must be even.
@@ -172,7 +170,7 @@ NC = 2*12;
 N=10000;
 //
 // The default generator must be Mersenne-Twister
-S=grand('getgen');
+S=grand("getgen");
 assert_checkequal ( S , "mt" );
 // The maximum integer generable with uin option
 UinMax = 2147483560;
@@ -183,39 +181,39 @@ UinMax = 2147483560;
 kgen = 1;
 gen = "mt";
 sdsize = seedsize(kgen);
-grand('setgen',gen);
-S=grand('getgen');
+grand("setgen",gen);
+S=grand("getgen");
 assert_checkequal ( S , gen );
 //
-grand('setsd',0);
-S=grand('getsd');
+grand("setsd",0);
+S=grand("getsd");
 assert_checkequal ( typeof(S) , "constant" );
 assert_checkequal ( size(S) , [sdsize 1] );
 //
-grand('setsd',123456);
-S=grand('getsd');
+grand("setsd",123456);
+S=grand("getsd");
 assert_checkequal ( typeof(S) , "constant" );
 assert_checkequal ( size(S) , [sdsize 1] );
 //
 // Check numbers
 expected = [
-0.5488135    0.6027634    0.4236548    0.4375872    0.9636628    0.7917250  
-0.5928446    0.8579456    0.6235637    0.2975346    0.2726563    0.8121687  
-0.7151894    0.5448832    0.6458941    0.891773     0.3834415    0.5288949  
-0.8442657    0.8472517    0.3843817    0.0567130    0.4776651    0.4799772  
+0.5488135    0.6027634    0.4236548    0.4375872    0.9636628    0.7917250
+0.5928446    0.8579456    0.6235637    0.2975346    0.2726563    0.8121687
+0.7151894    0.5448832    0.6458941    0.891773     0.3834415    0.5288949
+0.8442657    0.8472517    0.3843817    0.0567130    0.4776651    0.4799772
 ];
-grand('setsd',0);
+grand("setsd",0);
 computed = grand(4,6,"def");
 assert_checkalmostequal ( computed , expected, 1.e-6 );
 //
 // Check integers
 expected = [
-    2357136044.    2588848963.    1819583497.    1879422756.    4138900056.    3400433126.  
-    2546248239.    3684848379.    2678185683.    1277901399.    1171049868.    3488238119.  
-    3071714933.    2340255427.    2774094101.    3830135878.    1646868794.    2271586391.  
-    3626093760.    3638918503.    1650906866.    243580376.     2051556033.    2061486254.  
+2357136044.    2588848963.    1819583497.    1879422756.    4138900056.    3400433126.
+2546248239.    3684848379.    2678185683.    1277901399.    1171049868.    3488238119.
+3071714933.    2340255427.    2774094101.    3830135878.    1646868794.    2271586391.
+3626093760.    3638918503.    1650906866.    243580376.     2051556033.    2061486254.
 ];
-grand('setsd',0);
+grand("setsd",0);
 computed = grand(4,6,"lgi");
 assert_checkequal ( computed , expected );
 //
@@ -236,39 +234,39 @@ checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
 kgen = 2;
 gen = "kiss";
 sdsize = seedsize(kgen);
-grand('setgen',gen);
-S=grand('getgen');
+grand("setgen",gen);
+S=grand("getgen");
 assert_checkequal ( S , gen );
 //
-grand('setsd',0,0,0,0);
-S=grand('getsd');
+grand("setsd",0,0,0,0);
+S=grand("getsd");
 assert_checkequal ( typeof(S) , "constant" );
 assert_checkequal ( size(S) , [sdsize 1] );
 //
-grand('setsd',123456,123456,123456,123456);
-S=grand('getsd');
+grand("setsd",123456,123456,123456,123456);
+S=grand("getsd");
 assert_checkequal ( typeof(S) , "constant" );
 assert_checkequal ( size(S) , [sdsize 1] );
 //
 // Check numbers
 expected = [
-    2.874450D-04    9.423555D-01    7.707249D-02    2.078324D-02    4.746445D-01    1.895302D-01  
-    8.538282D-01    5.493145D-01    3.200836D-01    4.775516D-01    2.245108D-01    6.637360D-01  
-    5.815227D-02    6.006782D-01    8.569004D-01    1.235649D-02    7.357421D-01    5.837571D-01  
-    5.196679D-01    2.448867D-01    2.568304D-01    4.503826D-01    9.680347D-01    5.214808D-01  
+2.874450D-04    9.423555D-01    7.707249D-02    2.078324D-02    4.746445D-01    1.895302D-01
+8.538282D-01    5.493145D-01    3.200836D-01    4.775516D-01    2.245108D-01    6.637360D-01
+5.815227D-02    6.006782D-01    8.569004D-01    1.235649D-02    7.357421D-01    5.837571D-01
+5.196679D-01    2.448867D-01    2.568304D-01    4.503826D-01    9.680347D-01    5.214808D-01
 ];
-grand('setsd',0,0,0,0);
+grand("setsd",0,0,0,0);
 computed = grand(4,6,"def");
 assert_checkalmostequal ( computed , expected, 1.e-6 );
 //
 // Check integers
 expected = [
-    1234567.       4047385867.    331023823.     89263315.      2038582807.    814026139.   
-    3667164066.    2359287638.    1374748746.    2051068542.    964266482.     2850724518.  
-    249762113.     2579893349.    3680359369.    53070701.      3159988049.    2507217781.  
-    2231956628.    1051780200.    1103078268.    1934378448.    4157677412.    2239743032.  
+1234567.       4047385867.    331023823.     89263315.      2038582807.    814026139.
+3667164066.    2359287638.    1374748746.    2051068542.    964266482.     2850724518.
+249762113.     2579893349.    3680359369.    53070701.      3159988049.    2507217781.
+2231956628.    1051780200.    1103078268.    1934378448.    4157677412.    2239743032.
 ];
-grand('setsd',0,0,0,0);
+grand("setsd",0,0,0,0);
 computed = grand(4,6,"lgi");
 assert_checkequal ( computed , expected );
 //
@@ -289,39 +287,39 @@ checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
 kgen = 3;
 gen = "clcg2";
 sdsize = seedsize(kgen);
-grand('setgen',gen);
-S=grand('getgen');
+grand("setgen",gen);
+S=grand("getgen");
 assert_checkequal ( S , gen );
 //
-grand('setsd',1,1);
-S=grand('getsd');
+grand("setsd",1,1);
+S=grand("getsd");
 assert_checkequal ( typeof(S) , "constant" );
 assert_checkequal ( size(S) , [sdsize 1] );
 //
-grand('setsd',123456,123456);
-S=grand('getsd');
+grand("setsd",123456,123456);
+S=grand("getsd");
 assert_checkequal ( typeof(S) , "constant" );
 assert_checkequal ( size(S) , [sdsize 1] );
 //
 // Check numbers
 expected = [
-0.9999997    0.0369445    0.2041364    0.9100817    0.6998243    0.9596867  
-0.9745196    0.1617119    0.1673069    0.1117162    0.9502824    0.9149753  
-0.6474839    0.6646450    0.6549574    0.2990212    0.0918107    0.4411791  
-0.3330856    0.0846729    0.1288161    0.2654475    0.9023415    0.0735483  
+0.9999997    0.0369445    0.2041364    0.9100817    0.6998243    0.9596867
+0.9745196    0.1617119    0.1673069    0.1117162    0.9502824    0.9149753
+0.6474839    0.6646450    0.6549574    0.2990212    0.0918107    0.4411791
+0.3330856    0.0846729    0.1288161    0.2654475    0.9023415    0.0735483
 ];
-grand('setsd',1,1);
+grand("setsd",1,1);
 computed = grand(4,6,"def");
 assert_checkalmostequal ( computed , expected, 1.e-5 );
 //
 // Check integers
 expected = [
-    2147482884.    79337801.      438379562.     1954385533.    1502861091.    2060911403.  
-    2092764894.    347273588.     359288887.     239908781.     2040715732.    1964894377.  
-    1390461064.    1427314282.    1406510214.    642143055.     197161966.     947424871.   
-    715295839.     181833602.     276630551.     570044126.     1937763493.    157943826.   
+2147482884.    79337801.      438379562.     1954385533.    1502861091.    2060911403.
+2092764894.    347273588.     359288887.     239908781.     2040715732.    1964894377.
+1390461064.    1427314282.    1406510214.    642143055.     197161966.     947424871.
+715295839.     181833602.     276630551.     570044126.     1937763493.    157943826.
 ];
-grand('setsd',1,1);
+grand("setsd",1,1);
 computed = grand(4,6,"lgi");
 assert_checkequal ( computed , expected );
 //
@@ -342,46 +340,46 @@ checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
 kgen = 4;
 gen = "clcg4";
 sdsize = seedsize(kgen);
-grand('setgen',gen);
-S=grand('getgen');
+grand("setgen",gen);
+S=grand("getgen");
 assert_checkequal ( S , gen );
 //
 warning("off");
-grand('setsd',1,1,1,1);
+grand("setsd",1,1,1,1);
 warning("on");
-S=grand('getsd');
+S=grand("getsd");
 assert_checkequal ( typeof(S) , "constant" );
 assert_checkequal ( size(S) , [sdsize 1] );
 //
 warning("off");
-grand('setsd',123456,123456,123456,123456);
+grand("setsd",123456,123456,123456,123456);
 warning("on");
-S=grand('getsd');
+S=grand("getsd");
 assert_checkequal ( typeof(S) , "constant" );
 assert_checkequal ( size(S) , [sdsize 1] );
 //
 // Check numbers
 expected = [
-0.9999661    0.0552914    0.6345306    0.0640227    0.2885048    0.2781458  
-0.6852419    0.1806991    0.8665501    0.0981421    0.2660715    0.4279616  
-0.4370514    0.4956021    0.6870544    0.8501209    0.1271038    0.4554926  
-0.4202952    0.2903676    0.5712601    0.4764120    0.1818799    0.3121748  
+0.9999661    0.0552914    0.6345306    0.0640227    0.2885048    0.2781458
+0.6852419    0.1806991    0.8665501    0.0981421    0.2660715    0.4279616
+0.4370514    0.4956021    0.6870544    0.8501209    0.1271038    0.4554926
+0.4202952    0.2903676    0.5712601    0.4764120    0.1818799    0.3121748
 ];
 warning("off");
-grand('setsd',1,1,1,1);
+grand("setsd",1,1,1,1);
 warning("on");
 computed = grand(4,6,"def");
 assert_checkalmostequal ( computed , expected, 1.e-6 );
 //
 // Check integers
 expected = [
-    2147410798.    118737467.     1362644105.    137487708.     619559402.    597313629.  
-    1471545799.    388048305.     1860902184.    210758531.     571384155.    919040470.  
-    938560647.     1064297484.    1475437993.    1825620628.    272953383.    978162913.  
-    902576998.     623559632.     1226771622.    1023086907.    390584072.    670390361.  
+2147410798.    118737467.     1362644105.    137487708.     619559402.    597313629.
+1471545799.    388048305.     1860902184.    210758531.     571384155.    919040470.
+938560647.     1064297484.    1475437993.    1825620628.    272953383.    978162913.
+902576998.     623559632.     1226771622.    1023086907.    390584072.    670390361.
 ];
 warning("off");
-grand('setsd',1,1,1,1);
+grand("setsd",1,1,1,1);
 warning("on");
 computed = grand(4,6,"lgi");
 assert_checkequal ( computed , expected );
@@ -398,94 +396,41 @@ checkPieceLaw2arg ( "uin" , mycdfuin , N , NC , 0 , UinMax , rtol );
 checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
 ////////////////////////////////////////////////////////////////////
 //
-// "fsultra"
-//
-kgen = 5;
-gen = "fsultra";
-sdsize = seedsize(kgen);
-grand('setgen',gen);
-S=grand('getgen');
-assert_checkequal ( S , gen );
-//
-grand('setsd',1,1);
-S=grand('getsd');
-assert_checkequal ( typeof(S) , "constant" );
-assert_checkequal ( size(S) , [sdsize 1] );
-//
-grand('setsd',123456,123456);
-S=grand('getsd');
-assert_checkequal ( typeof(S) , "constant" );
-assert_checkequal ( size(S) , [sdsize 1] );
-//
-// Check numbers
-expected = [
-0.3314877    0.3699260    0.4383216    0.99706      0.0577929    0.4836669  
-0.5826624    0.9600475    0.2037475    0.6774254    0.4450278    0.3082941  
-0.1630857    0.2033307    0.4214824    0.6372521    0.0782678    0.4409892  
-0.7211611    0.1833922    0.8065496    0.6375251    0.2572713    0.8039582  
-];
-grand('setsd',1,1);
-computed = grand(4,6,"def");
-assert_checkalmostequal ( computed , expected, 1.e-6 );
-//
-// Check integers
-expected = [
-    1423728774.    1588820113.    1882577034.    4282340079.    248218608.     2077333598.  
-    2502516061.    4123372468.    875088704.     2909519830.    1911379739.    1324113135.  
-    700447838.     873298853.     1810253313.    2736976953.    336157762.     1894034123.  
-    3097363227.    787663378.     3464104206.    2738149622.    1104971606.    3452974260.  
-];
-grand('setsd',1,1);
-computed = grand(4,6,"lgi");
-assert_checkequal ( computed , expected );
-//
-// Check distribution of uniform numbers in [0,1[
-checkMeanVariance0arg ( 400 , 800 , "def" , 1/2 , 1/12 , rtol );
-checkLaw0arg ( "def" , mycdfdef , N , NC , rtol );
-//
-// Check distribution of uniform integers in [A,B]
-A = MinInt(kgen);
-B = MaxInt(kgen);
-checkMeanVariance0arg ( 400 , 800 , "lgi" , (A+B)/2 , ((B-A+1)^2-1)/12 , rtol );
-checkPieceLaw2arg ( "uin" , mycdfuin , N , NC , 0 , UinMax , rtol );
-checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
-////////////////////////////////////////////////////////////////////
-//
 // "urand"
 //
-kgen = 6;
+kgen = 5;
 gen = "urand";
-grand('setgen',gen);
-S=grand('getgen');
+grand("setgen",gen);
+S=grand("getgen");
 assert_checkequal ( S , gen );
 //
-grand('setsd',1);
-S=grand('getsd');
+grand("setsd",1);
+S=grand("getsd");
 assert_checkequal ( S , 1 );
 //
-grand('setsd',123456);
-S=grand('getsd');
+grand("setsd",123456);
+S=grand("getsd");
 assert_checkequal ( S , 123456 );
 //
 // Check numbers
 expected = [
-0.6040239    0.5321420    0.2276811    0.8979351    0.8925854    0.6928366  
-0.0079647    0.4138784    0.6656067    0.8274020    0.0848163    0.6776849  
-0.6643966    0.5036204    0.9694369    0.0664231    0.2566682    0.4284010  
-0.9832111    0.6850569    0.0775390    0.1099766    0.6507795    0.3725794  
+0.6040239    0.5321420    0.2276811    0.8979351    0.8925854    0.6928366
+0.0079647    0.4138784    0.6656067    0.8274020    0.0848163    0.6776849
+0.6643966    0.5036204    0.9694369    0.0664231    0.2566682    0.4284010
+0.9832111    0.6850569    0.0775390    0.1099766    0.6507795    0.3725794
 ];
-grand('setsd',1);
+grand("setsd",1);
 computed = grand(4,6,"def");
 assert_checkalmostequal ( computed , expected, 1.e-5 );
 //
 // Check integers
 expected = [
-    1297131554.    1142766270.    488941338.     1928300854.    1916812562.    1487855278.  
-    17103983.      888797147.     1429379591.    1776832243.    182141599.     1455317259.  
-    1426780792.    1081516660.    2081849904.    142642604.     551190760.     919984100.   
-    2111429773.    1471148505.    166513637.     236172977.     1397538365.    800108169.   
+1297131554.    1142766270.    488941338.     1928300854.    1916812562.    1487855278.
+17103983.      888797147.     1429379591.    1776832243.    182141599.     1455317259.
+1426780792.    1081516660.    2081849904.    142642604.     551190760.     919984100.
+2111429773.    1471148505.    166513637.     236172977.     1397538365.    800108169.
 ];
-grand('setsd',1);
+grand("setsd",1);
 computed = grand(4,6,"lgi");
 assert_checkequal ( computed , expected );
 //
index 0c73d3b..9537b74 100644 (file)
 // <-- ENGLISH IMPOSED -->
 
 function checkMeanVariance2arg ( m , n , name , A , B , mu , va , rtol )
-  // Check the mean and variance of random numbers.
-  //
-  // Parameters
-  // m : a 1-by-1 matrix of floating point integers, the number of rows
-  // n : a 1-by-1 matrix of floating point integers, the number of columns
-  // name: a 1-by-1 string, the name of the distribution function
-  // A : a 1-by-1 matrix of doubles, the value of the 1st parameter
-  // B : a 1-by-1 matrix of doubles, the value of the 2nd parameter
-  // mu : a 1-by-1 matrix of doubles, the expected mean
-  // va : a 1-by-1 matrix of doubles, the expected variance. If va==[], then the variance is not checked.
-  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
-  
-  R=grand(m,n,name,A,B);
-  assert_checkequal ( size(R) , [m,n] );
-  assert_checkequal ( typeof(R) , "constant" );
-  assert_checkalmostequal ( mean(R) , mu , rtol );
-  if ( va<>[] ) then
-    assert_checkalmostequal ( variance(R) , va , rtol );
-  end
+    // Check the mean and variance of random numbers.
+    //
+    // Parameters
+    // m : a 1-by-1 matrix of floating point integers, the number of rows
+    // n : a 1-by-1 matrix of floating point integers, the number of columns
+    // name: a 1-by-1 string, the name of the distribution function
+    // A : a 1-by-1 matrix of doubles, the value of the 1st parameter
+    // B : a 1-by-1 matrix of doubles, the value of the 2nd parameter
+    // mu : a 1-by-1 matrix of doubles, the expected mean
+    // va : a 1-by-1 matrix of doubles, the expected variance. If va==[], then the variance is not checked.
+    // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+
+    R=grand(m,n,name,A,B);
+    assert_checkequal ( size(R) , [m,n] );
+    assert_checkequal ( typeof(R) , "constant" );
+    assert_checkalmostequal ( mean(R) , mu , rtol );
+    if ( va<>[] ) then
+        assert_checkalmostequal ( variance(R) , va , rtol );
+    end
 endfunction
 
 function checkMeanVariance0arg ( m , n , name , mu , va , rtol )
-  // Check the mean and variance of random numbers.
-  //
-  // Parameters
-  // m : a 1-by-1 matrix of floating point integers, the number of rows
-  // n : a 1-by-1 matrix of floating point integers, the number of columns
-  // name: a 1-by-1 string, the name of the distribution function
-  // mu : a 1-by-1 matrix of doubles, the expected mean
-  // va : a 1-by-1 matrix of doubles, the expected variance. If va==[], then the variance is not checked.
-  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
-  
-  R=grand(m,n,name);
-  assert_checkequal ( size(R) , [m,n] );
-  assert_checkequal ( typeof(R) , "constant" );
-  assert_checkalmostequal ( mean(R) , mu , rtol );
-  if ( va<>[] ) then
-    assert_checkalmostequal ( variance(R) , va , rtol );
-  end
+    // Check the mean and variance of random numbers.
+    //
+    // Parameters
+    // m : a 1-by-1 matrix of floating point integers, the number of rows
+    // n : a 1-by-1 matrix of floating point integers, the number of columns
+    // name: a 1-by-1 string, the name of the distribution function
+    // mu : a 1-by-1 matrix of doubles, the expected mean
+    // va : a 1-by-1 matrix of doubles, the expected variance. If va==[], then the variance is not checked.
+    // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+
+    R=grand(m,n,name);
+    assert_checkequal ( size(R) , [m,n] );
+    assert_checkequal ( typeof(R) , "constant" );
+    assert_checkalmostequal ( mean(R) , mu , rtol );
+    if ( va<>[] ) then
+        assert_checkalmostequal ( variance(R) , va , rtol );
+    end
 endfunction
 
 function checkLaw0arg ( name , cdffun , N , NC , rtol )
-  //
-  // Check the random number generator for a continuous distribution function.
-  //
-  // Parameters
-  // name: a 1-by-1 string, the name of the distribution function
-  // cdffun : a function, the Cumulated Distribution Function
-  // NC : a 1-by-1 matrix of floating point integers, the number of classes
-  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
-  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
-  //
-  // Description
-  //  Compare the empirical histogram with the theoretical histogram.
+    //
+    // Check the random number generator for a continuous distribution function.
+    //
+    // Parameters
+    // name: a 1-by-1 string, the name of the distribution function
+    // cdffun : a function, the Cumulated Distribution Function
+    // NC : a 1-by-1 matrix of floating point integers, the number of classes
+    // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+    // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+    //
+    // Description
+    //  Compare the empirical histogram with the theoretical histogram.
 
-  R = grand(1,N,name);
-  [X,EmpiricalPDF] = histcompute(NC,R);
-  CDF = cdffun(X)
-  TheoricPDF = diff(CDF);
-  assert_checktrue ( abs(EmpiricalPDF-TheoricPDF) < rtol );
+    R = grand(1,N,name);
+    [X,EmpiricalPDF] = histcompute(NC,R);
+    CDF = cdffun(X)
+    TheoricPDF = diff(CDF);
+    assert_checktrue ( abs(EmpiricalPDF-TheoricPDF) < rtol );
     if ( %f ) then
-      plot(X(1:$-1),EmpiricalPDF,"bo-"); // Empirical Histogram
-      plot(X(1:$-1),TheoricPDF,"rox-"); // Theoretical Histogram
+        plot(X(1:$-1),EmpiricalPDF,"bo-"); // Empirical Histogram
+        plot(X(1:$-1),TheoricPDF,"rox-"); // Theoretical Histogram
     end
 endfunction
 
 function checkPieceLaw0arg ( name , cdffun , N , NC , rtol )
-  //
-  // Check the random number generator for a piecewise distribution function.
-  //
-  // Parameters
-  // name: a 1-by-1 string, the name of the distribution function
-  // cdffun : a function, the Cumulated Distribution Function
-  // NC : a 1-by-1 matrix of floating point integers, the number of classes
-  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
-  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
-  //
-  // Description
-  //  Compare the empirical histogram with the theoretical histogram.
-  //  The classes of the histogram are computed from the random samples, 
-  //  from which the unique entries are extracted.
+    //
+    // Check the random number generator for a piecewise distribution function.
+    //
+    // Parameters
+    // name: a 1-by-1 string, the name of the distribution function
+    // cdffun : a function, the Cumulated Distribution Function
+    // NC : a 1-by-1 matrix of floating point integers, the number of classes
+    // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+    // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+    //
+    // Description
+    //  Compare the empirical histogram with the theoretical histogram.
+    //  The classes of the histogram are computed from the random samples,
+    //  from which the unique entries are extracted.
 
-  R=grand(N,1,name);
-  X = unique(R);
-  for k = 1 : size(X,"*")
-    EmpiricalPDF(k) = length(find(R==X(k)));
-  end
-  EmpiricalPDF = EmpiricalPDF./N;
-  CDF = cdffun(X)
-  TheoricPDF=[CDF(1);diff(CDF)];
-  assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
-  if ( %f ) then
-    plot(X,EmpiricalPDF,"bo-"); // Empirical Histogram
-    plot(X,TheoricPDF,"rox-"); // Theoretical Histogram
-  end
+    R=grand(N,1,name);
+    X = unique(R);
+    for k = 1 : size(X,"*")
+        EmpiricalPDF(k) = length(find(R==X(k)));
+    end
+    EmpiricalPDF = EmpiricalPDF./N;
+    CDF = cdffun(X)
+    TheoricPDF=[CDF(1);diff(CDF)];
+    assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
+    if ( %f ) then
+        plot(X,EmpiricalPDF,"bo-"); // Empirical Histogram
+        plot(X,TheoricPDF,"rox-"); // Theoretical Histogram
+    end
 endfunction
 
 function checkPieceLaw2arg ( name , cdffun , N , NC , A , B , rtol )
-  //
-  // Check the random number generator for a piecewise distribution function.
-  //
-  // Parameters
-  // name: a 1-by-1 string, the name of the distribution function
-  // cdffun : a function, the Cumulated Distribution Function
-  // NC : a 1-by-1 matrix of floating point integers, the number of classes
-  // N : a 1-by-1 matrix of floating point integers, the number of random values to test
-  // A : a 1-by-1 matrix of doubles, the value of the parameter
-  // rtol : a 1-by-1 matrix of doubles, the relative tolerance
-  //
-  // Description
-  //  Compare the empirical histogram with the theoretical histogram.
-  //  The classes of the histogram are computed from the random samples, 
-  //  from which the unique entries are extracted.
+    //
+    // Check the random number generator for a piecewise distribution function.
+    //
+    // Parameters
+    // name: a 1-by-1 string, the name of the distribution function
+    // cdffun : a function, the Cumulated Distribution Function
+    // NC : a 1-by-1 matrix of floating point integers, the number of classes
+    // N : a 1-by-1 matrix of floating point integers, the number of random values to test
+    // A : a 1-by-1 matrix of doubles, the value of the parameter
+    // rtol : a 1-by-1 matrix of doubles, the relative tolerance
+    //
+    // Description
+    //  Compare the empirical histogram with the theoretical histogram.
+    //  The classes of the histogram are computed from the random samples,
+    //  from which the unique entries are extracted.
 
-  R=grand(N,1,name,A,B);
-  X = unique(R);
-  for k = 1 : size(X,"*")
-    EmpiricalPDF(k) = length(find(R==X(k)));
-  end
-  EmpiricalPDF = EmpiricalPDF./N;
-  CDF = cdffun(X,A,B)
-  TheoricPDF=[CDF(1);diff(CDF)];
-  assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
-  if ( %f ) then
-    plot(X,EmpiricalPDF,"bo-"); // Empirical Histogram
-    plot(X,TheoricPDF,"rox-"); // Theoretical Histogram
-  end
+    R=grand(N,1,name,A,B);
+    X = unique(R);
+    for k = 1 : size(X,"*")
+        EmpiricalPDF(k) = length(find(R==X(k)));
+    end
+    EmpiricalPDF = EmpiricalPDF./N;
+    CDF = cdffun(X,A,B)
+    TheoricPDF=[CDF(1);diff(CDF)];
+    assert_checktrue( abs(EmpiricalPDF-TheoricPDF) < rtol );
+    if ( %f ) then
+        plot(X,EmpiricalPDF,"bo-"); // Empirical Histogram
+        plot(X,TheoricPDF,"rox-"); // Theoretical Histogram
+    end
 endfunction
 
 function [x,y] = histcompute(n,data)
-   //
-   // Computes the histogram of a data.
-   //
-   // Parameters
-   // n : a 1-by-1 matrix of floating point integers, the number of classes.
-   // data : a matrix of doubles, the data
-   // x : 1-by-n+1 matrix of doubles, the classes of the histogram
-   // y : 1-by-n+1 matrix of doubles, the empirical probability that one value in data which fall in each class
-   x = linspace(min(data), max(data), n+1);
-   [ind , y] = dsearch(data, x)
-   y = y./length(data)
+    //
+    // Computes the histogram of a data.
+    //
+    // Parameters
+    // n : a 1-by-1 matrix of floating point integers, the number of classes.
+    // data : a matrix of doubles, the data
+    // x : 1-by-n+1 matrix of doubles, the classes of the histogram
+    // y : 1-by-n+1 matrix of doubles, the empirical probability that one value in data which fall in each class
+    x = linspace(min(data), max(data), n+1);
+    [ind , y] = dsearch(data, x)
+    y = y./length(data)
 endfunction
 
 function p = mycdfdef (X)
-  p = X
+    p = X
 endfunction
 
 function p = mycdfuin (X,A,B)
-  p = (floor(X)-A+1)/(B-A+1)
+    p = (floor(X)-A+1)/(B-A+1)
 endfunction
 
 function p = mycdflgi (X)
-  // Here, A and B are read from the environment
-  p = (floor(X)-A+1)/(B-A+1)
+    // Here, A and B are read from the environment
+    p = (floor(X)-A+1)/(B-A+1)
 endfunction
 
 
@@ -174,11 +174,11 @@ endfunction
 // MinInt : minimum of the uniform integer interval for random number generation
 // MaxInt : maximum of the uniform integer interval for random number generation
 //
-NGen = 6;
-generators = ["mt"   "kiss" "clcg2"    "clcg4" "fsultra" "urand"];
-seedsize =   [625    4      2          4       40        1];
-MinInt =     [0      0      0          0       0         0];
-MaxInt =     [2^32-1 2^32-1 2147483561 2^31-2  2^32-1    2^31-1];
+NGen = 5;
+generators = ["mt"   "kiss" "clcg2"    "clcg4" "urand"];
+seedsize =   [625    4      2          4       1];
+MinInt =     [0      0      0          0       0];
+MaxInt =     [2^32-1 2^32-1 2147483561 2^31-2  2^31-1];
 
 rtol = 1.e-2;
 
@@ -189,7 +189,7 @@ N=10000;
 
 //
 // The default generator must be Mersenne-Twister
-S=grand('getgen');
+S=grand("getgen");
 assert_checkequal ( S , "mt" );
 
 // The maximum integer generable with uin option
@@ -202,39 +202,39 @@ UinMax = 2147483560;
 kgen = 1;
 gen = "mt";
 sdsize = seedsize(kgen);
-grand('setgen',gen);
-S=grand('getgen');
+grand("setgen",gen);
+S=grand("getgen");
 assert_checkequal ( S , gen );
 //
-grand('setsd',0);
-S=grand('getsd');
+grand("setsd",0);
+S=grand("getsd");
 assert_checkequal ( typeof(S) , "constant" );
 assert_checkequal ( size(S) , [sdsize 1] );
 //
-grand('setsd',123456);
-S=grand('getsd');
+grand("setsd",123456);
+S=grand("getsd");
 assert_checkequal ( typeof(S) , "constant" );
 assert_checkequal ( size(S) , [sdsize 1] );
 //
 // Check numbers
 expected = [
-0.5488135    0.6027634    0.4236548    0.4375872    0.9636628    0.7917250  
-0.5928446    0.8579456    0.6235637    0.2975346    0.2726563    0.8121687  
-0.7151894    0.5448832    0.6458941    0.891773     0.3834415    0.5288949  
-0.8442657    0.8472517    0.3843817    0.0567130    0.4776651    0.4799772  
+0.5488135    0.6027634    0.4236548    0.4375872    0.9636628    0.7917250
+0.5928446    0.8579456    0.6235637    0.2975346    0.2726563    0.8121687
+0.7151894    0.5448832    0.6458941    0.891773     0.3834415    0.5288949
+0.8442657    0.8472517    0.3843817    0.0567130    0.4776651    0.4799772
 ];
-grand('setsd',0);
+grand("setsd",0);
 computed = grand(4,6,"def");
 assert_checkalmostequal ( computed , expected, 1.e-6 );
 //
 // Check integers
 expected = [
-    2357136044.    2588848963.    1819583497.    1879422756.    4138900056.    3400433126.  
-    2546248239.    3684848379.    2678185683.    1277901399.    1171049868.    3488238119.  
-    3071714933.    2340255427.    2774094101.    3830135878.    1646868794.    2271586391.  
-    3626093760.    3638918503.    1650906866.    243580376.     2051556033.    2061486254.  
+2357136044.    2588848963.    1819583497.    1879422756.    4138900056.    3400433126.
+2546248239.    3684848379.    2678185683.    1277901399.    1171049868.    3488238119.
+3071714933.    2340255427.    2774094101.    3830135878.    1646868794.    2271586391.
+3626093760.    3638918503.    1650906866.    243580376.     2051556033.    2061486254.
 ];
-grand('setsd',0);
+grand("setsd",0);
 computed = grand(4,6,"lgi");
 assert_checkequal ( computed , expected );
 //
@@ -256,39 +256,39 @@ checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
 kgen = 2;
 gen = "kiss";
 sdsize = seedsize(kgen);
-grand('setgen',gen);
-S=grand('getgen');
+grand("setgen",gen);
+S=grand("getgen");
 assert_checkequal ( S , gen );
 //
-grand('setsd',0,0,0,0);
-S=grand('getsd');
+grand("setsd",0,0,0,0);
+S=grand("getsd");
 assert_checkequal ( typeof(S) , "constant" );
 assert_checkequal ( size(S) , [sdsize 1] );
 //
-grand('setsd',123456,123456,123456,123456);
-S=grand('getsd');
+grand("setsd",123456,123456,123456,123456);
+S=grand("getsd");
 assert_checkequal ( typeof(S) , "constant" );
 assert_checkequal ( size(S) , [sdsize 1] );
 //
 // Check numbers
 expected = [
-    2.874450D-04    9.423555D-01    7.707249D-02    2.078324D-02    4.746445D-01    1.895302D-01  
-    8.538282D-01    5.493145D-01    3.200836D-01    4.775516D-01    2.245108D-01    6.637360D-01  
-    5.815227D-02    6.006782D-01    8.569004D-01    1.235649D-02    7.357421D-01    5.837571D-01  
-    5.196679D-01    2.448867D-01    2.568304D-01    4.503826D-01    9.680347D-01    5.214808D-01  
+2.874450D-04    9.423555D-01    7.707249D-02    2.078324D-02    4.746445D-01    1.895302D-01
+8.538282D-01    5.493145D-01    3.200836D-01    4.775516D-01    2.245108D-01    6.637360D-01
+5.815227D-02    6.006782D-01    8.569004D-01    1.235649D-02    7.357421D-01    5.837571D-01
+5.196679D-01    2.448867D-01    2.568304D-01    4.503826D-01    9.680347D-01    5.214808D-01
 ];
-grand('setsd',0,0,0,0);
+grand("setsd",0,0,0,0);
 computed = grand(4,6,"def");
 assert_checkalmostequal ( computed , expected, 1.e-6 );
 //
 // Check integers
 expected = [
-    1234567.       4047385867.    331023823.     89263315.      2038582807.    814026139.   
-    3667164066.    2359287638.    1374748746.    2051068542.    964266482.     2850724518.  
-    249762113.     2579893349.    3680359369.    53070701.      3159988049.    2507217781.  
-    2231956628.    1051780200.    1103078268.    1934378448.    4157677412.    2239743032.  
+1234567.       4047385867.    331023823.     89263315.      2038582807.    814026139.
+3667164066.    2359287638.    1374748746.    2051068542.    964266482.     2850724518.
+249762113.     2579893349.    3680359369.    53070701.      3159988049.    2507217781.
+2231956628.    1051780200.    1103078268.    1934378448.    4157677412.    2239743032.
 ];
-grand('setsd',0,0,0,0);
+grand("setsd",0,0,0,0);
 computed = grand(4,6,"lgi");
 assert_checkequal ( computed , expected );
 //
@@ -309,39 +309,39 @@ checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
 kgen = 3;
 gen = "clcg2";
 sdsize = seedsize(kgen);
-grand('setgen',gen);
-S=grand('getgen');
+grand("setgen",gen);
+S=grand("getgen");
 assert_checkequal ( S , gen );
 //
-grand('setsd',1,1);
-S=grand('getsd');
+grand("setsd",1,1);
+S=grand("getsd");
 assert_checkequal ( typeof(S) , "constant" );
 assert_checkequal ( size(S) , [sdsize 1] );
 //
-grand('setsd',123456,123456);
-S=grand('getsd');
+grand("setsd",123456,123456);
+S=grand("getsd");
 assert_checkequal ( typeof(S) , "constant" );
 assert_checkequal ( size(S) , [sdsize 1] );
 //
 // Check numbers
 expected = [
-0.9999997    0.0369445    0.2041364    0.9100817    0.6998243    0.9596867  
-0.9745196    0.1617119    0.1673069    0.1117162    0.9502824    0.9149753  
-0.6474839    0.6646450    0.6549574    0.2990212    0.0918107    0.4411791  
-0.3330856    0.0846729    0.1288161    0.2654475    0.9023415    0.0735483  
+0.9999997    0.0369445    0.2041364    0.9100817    0.6998243    0.9596867
+0.9745196    0.1617119    0.1673069    0.1117162    0.9502824    0.9149753
+0.6474839    0.6646450    0.6549574    0.2990212    0.0918107    0.4411791
+0.3330856    0.0846729    0.1288161    0.2654475    0.9023415    0.0735483
 ];
-grand('setsd',1,1);
+grand("setsd",1,1);
 computed = grand(4,6,"def");
 assert_checkalmostequal ( computed , expected, 1.e-5 );
 //
 // Check integers
 expected = [
-    2147482884.    79337801.      438379562.     1954385533.    1502861091.    2060911403.  
-    2092764894.    347273588.     359288887.     239908781.     2040715732.    1964894377.  
-    1390461064.    1427314282.    1406510214.    642143055.     197161966.     947424871.   
-    715295839.     181833602.     276630551.     570044126.     1937763493.    157943826.   
+2147482884.    79337801.      438379562.     1954385533.    1502861091.    2060911403.
+2092764894.    347273588.     359288887.     239908781.     2040715732.    1964894377.
+1390461064.    1427314282.    1406510214.    642143055.     197161966.     947424871.
+715295839.     181833602.     276630551.     570044126.     1937763493.    157943826.
 ];
-grand('setsd',1,1);
+grand("setsd",1,1);
 computed = grand(4,6,"lgi");
 assert_checkequal ( computed , expected );
 //
@@ -362,46 +362,46 @@ checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
 kgen = 4;
 gen = "clcg4";
 sdsize = seedsize(kgen);
-grand('setgen',gen);
-S=grand('getgen');
+grand("setgen",gen);
+S=grand("getgen");
 assert_checkequal ( S , gen );
 //
 warning("off");
-grand('setsd',1,1,1,1);
+grand("setsd",1,1,1,1);
 warning("on");
-S=grand('getsd');
+S=grand("getsd");
 assert_checkequal ( typeof(S) , "constant" );
 assert_checkequal ( size(S) , [sdsize 1] );
 //
 warning("off");
-grand('setsd',123456,123456,123456,123456);
+grand("setsd",123456,123456,123456,123456);
 warning("on");
-S=grand('getsd');
+S=grand("getsd");
 assert_checkequal ( typeof(S) , "constant" );
 assert_checkequal ( size(S) , [sdsize 1] );
 //
 // Check numbers
 expected = [
-0.9999661    0.0552914    0.6345306    0.0640227    0.2885048    0.2781458  
-0.6852419    0.1806991    0.8665501    0.0981421    0.2660715    0.4279616  
-0.4370514    0.4956021    0.6870544    0.8501209    0.1271038    0.4554926  
-0.4202952    0.2903676    0.5712601    0.4764120    0.1818799    0.3121748  
+0.9999661    0.0552914    0.6345306    0.0640227    0.2885048    0.2781458
+0.6852419    0.1806991    0.8665501    0.0981421    0.2660715    0.4279616
+0.4370514    0.4956021    0.6870544    0.8501209    0.1271038    0.4554926
+0.4202952    0.2903676    0.5712601    0.4764120    0.1818799    0.3121748
 ];
 warning("off");
-grand('setsd',1,1,1,1);
+grand("setsd",1,1,1,1);
 warning("on");
 computed = grand(4,6,"def");
 assert_checkalmostequal ( computed , expected, 1.e-6 );
 //
 // Check integers
 expected = [
-    2147410798.    118737467.     1362644105.    137487708.     619559402.    597313629.  
-    1471545799.    388048305.     1860902184.    210758531.     571384155.    919040470.  
-    938560647.     1064297484.    1475437993.    1825620628.    272953383.    978162913.  
-    902576998.     623559632.     1226771622.    1023086907.    390584072.    670390361.  
+2147410798.    118737467.     1362644105.    137487708.     619559402.    597313629.
+1471545799.    388048305.     1860902184.    210758531.     571384155.    919040470.
+938560647.     1064297484.    1475437993.    1825620628.    272953383.    978162913.
+902576998.     623559632.     1226771622.    1023086907.    390584072.    670390361.
 ];
 warning("off");
-grand('setsd',1,1,1,1);
+grand("setsd",1,1,1,1);
 warning("on");
 computed = grand(4,6,"lgi");
 assert_checkequal ( computed , expected );
@@ -418,94 +418,41 @@ checkPieceLaw2arg ( "uin" , mycdfuin , N , NC , 0 , UinMax , rtol );
 checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
 ////////////////////////////////////////////////////////////////////
 //
-// "fsultra"
-//
-kgen = 5;
-gen = "fsultra";
-sdsize = seedsize(kgen);
-grand('setgen',gen);
-S=grand('getgen');
-assert_checkequal ( S , gen );
-//
-grand('setsd',1,1);
-S=grand('getsd');
-assert_checkequal ( typeof(S) , "constant" );
-assert_checkequal ( size(S) , [sdsize 1] );
-//
-grand('setsd',123456,123456);
-S=grand('getsd');
-assert_checkequal ( typeof(S) , "constant" );
-assert_checkequal ( size(S) , [sdsize 1] );
-//
-// Check numbers
-expected = [
-0.3314877    0.3699260    0.4383216    0.99706      0.0577929    0.4836669  
-0.5826624    0.9600475    0.2037475    0.6774254    0.4450278    0.3082941  
-0.1630857    0.2033307    0.4214824    0.6372521    0.0782678    0.4409892  
-0.7211611    0.1833922    0.8065496    0.6375251    0.2572713    0.8039582  
-];
-grand('setsd',1,1);
-computed = grand(4,6,"def");
-assert_checkalmostequal ( computed , expected, 1.e-6 );
-//
-// Check integers
-expected = [
-    1423728774.    1588820113.    1882577034.    4282340079.    248218608.     2077333598.  
-    2502516061.    4123372468.    875088704.     2909519830.    1911379739.    1324113135.  
-    700447838.     873298853.     1810253313.    2736976953.    336157762.     1894034123.  
-    3097363227.    787663378.     3464104206.    2738149622.    1104971606.    3452974260.  
-];
-grand('setsd',1,1);
-computed = grand(4,6,"lgi");
-assert_checkequal ( computed , expected );
-//
-// Check distribution of uniform numbers in [0,1[
-checkMeanVariance0arg ( 400 , 800 , "def" , 1/2 , 1/12 , rtol );
-checkLaw0arg ( "def" , mycdfdef , N , NC , rtol );
-//
-// Check distribution of uniform integers in [A,B]
-A = MinInt(kgen);
-B = MaxInt(kgen);
-checkMeanVariance0arg ( 400 , 800 , "lgi" , (A+B)/2 , ((B-A+1)^2-1)/12 , rtol );
-checkPieceLaw2arg ( "uin" , mycdfuin , N , NC , 0 , UinMax , rtol );
-checkPieceLaw0arg ( "lgi" , mycdflgi , N , NC , rtol );
-////////////////////////////////////////////////////////////////////
-//
 // "urand"
 //
-kgen = 6;
+kgen = 5;
 gen = "urand";
-grand('setgen',gen);
-S=grand('getgen');
+grand("setgen",gen);
+S=grand("getgen");
 assert_checkequal ( S , gen );
 //
-grand('setsd',1);
-S=grand('getsd');
+grand("setsd",1);
+S=grand("getsd");
 assert_checkequal ( S , 1 );
 //
-grand('setsd',123456);
-S=grand('getsd');
+grand("setsd",123456);
+S=grand("getsd");
 assert_checkequal ( S , 123456 );
 //
 // Check numbers
 expected = [
-0.6040239    0.5321420    0.2276811    0.8979351    0.8925854    0.6928366  
-0.0079647    0.4138784    0.6656067    0.8274020    0.0848163    0.6776849  
-0.6643966    0.5036204    0.9694369    0.0664231    0.2566682    0.4284010  
-0.9832111    0.6850569    0.0775390    0.1099766    0.6507795    0.3725794  
+0.6040239    0.5321420    0.2276811    0.8979351    0.8925854    0.6928366
+0.0079647    0.4138784    0.6656067    0.8274020    0.0848163    0.6776849
+0.6643966    0.5036204    0.9694369    0.0664231    0.2566682    0.4284010
+0.9832111    0.6850569    0.0775390    0.1099766    0.6507795    0.3725794
 ];
-grand('setsd',1);
+grand("setsd",1);
 computed = grand(4,6,"def");
 assert_checkalmostequal ( computed , expected, 1.e-5 );
 //
 // Check integers
 expected = [
-    1297131554.    1142766270.    488941338.     1928300854.    1916812562.    1487855278.  
-    17103983.      888797147.     1429379591.    1776832243.    182141599.     1455317259.  
-    1426780792.    1081516660.    2081849904.    142642604.     551190760.     919984100.   
-    2111429773.    1471148505.    166513637.     236172977.     1397538365.    800108169.   
+1297131554.    1142766270.    488941338.     1928300854.    1916812562.    1487855278.
+17103983.      888797147.     1429379591.    1776832243.    182141599.     1455317259.
+1426780792.    1081516660.    2081849904.    142642604.     551190760.     919984100.
+2111429773.    1471148505.    166513637.     236172977.     1397538365.    800108169.
 ];
-grand('setsd',1);
+grand("setsd",1);
 computed = grand(4,6,"lgi");
 assert_checkequal ( computed , expected );
 //
index 2995e10..5612480 100644 (file)
 ///////////////////////////////////////////////////////////////////////////////
 //
 // Dimensions
-mat = grand(100, 101, 102, 'unf', 0, 1);
+mat = grand(100, 101, 102, "unf", 0, 1);
 assert_checktrue(size(mat) == [100 101 102]);
 ///////////////////////////////////////////////////////////////////////////////
 //
 // Generators
 // The grand(i, j, ...) should be equal to the first element of grand(i, j, k, ...).
 // mt generator
-grand('setgen', 'mt');
-grand('setsd', 0);
-expected = grand(4, 6, 'def');
-grand('setsd', 0); // Resetting the seed to obtain the same results
-computed = grand(4, 6, 5, 'def');
+grand("setgen", "mt");
+grand("setsd", 0);
+expected = grand(4, 6, "def");
+grand("setsd", 0); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, "def");
 assert_checkequal(expected, computed(:, :, 1));
-grand('setsd', 0);
-expected = grand(4, 6, 'lgi');
-grand('setsd', 0); // Resetting the seed to obtain the same results
-computed = grand(4, 6, 5, 'lgi');
+grand("setsd", 0);
+expected = grand(4, 6, "lgi");
+grand("setsd", 0); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, "lgi");
 assert_checkequal(expected, computed(:, :, 1));
 // kiss generator
-grand('setgen', 'kiss');
-grand('setsd', 0, 0, 0, 0);
-expected = grand(4, 6, 'def');
-grand('setsd', 0, 0, 0, 0); // Resetting the seed to obtain the same results
-computed = grand(4, 6, 5, 'def');
+grand("setgen", "kiss");
+grand("setsd", 0, 0, 0, 0);
+expected = grand(4, 6, "def");
+grand("setsd", 0, 0, 0, 0); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, "def");
 assert_checkequal(expected, computed(:, :, 1));
-grand('setsd', 0, 0, 0, 0);
-expected = grand(4, 6, 'lgi');
-grand('setsd', 0, 0, 0, 0); // Resetting the seed to obtain the same results
-computed = grand(4, 6, 5, 'lgi');
+grand("setsd", 0, 0, 0, 0);
+expected = grand(4, 6, "lgi");
+grand("setsd", 0, 0, 0, 0); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, "lgi");
 assert_checkequal(expected, computed(:, :, 1));
 // clcg2 generator
-grand('setgen', 'clcg2');
-grand('setsd', 1, 1);
-expected = grand(4, 6, 'def');
-grand('setsd', 1, 1); // Resetting the seed to obtain the same results
-computed = grand(4, 6, 5, 'def');
+grand("setgen", "clcg2");
+grand("setsd", 1, 1);
+expected = grand(4, 6, "def");
+grand("setsd", 1, 1); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, "def");
 assert_checkequal(expected, computed(:, :, 1));
-grand('setsd', 1, 1);
-expected = grand(4, 6, 'lgi');
-grand('setsd', 1, 1); // Resetting the seed to obtain the same results
-computed = grand(4, 6, 5, 'lgi');
+grand("setsd", 1, 1);
+expected = grand(4, 6, "lgi");
+grand("setsd", 1, 1); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, "lgi");
 assert_checkequal(expected, computed(:, :, 1));
 // clcg4 generator
-grand('setgen', 'clcg4');
-warning('off');
-grand('setsd',1,1,1,1);
-warning('on');
-expected = grand(4, 6, 'def');
-warning('off');
-grand('setsd',1,1,1,1); // Resetting the seed to obtain the same results
-warning('on');
-computed = grand(4, 6, 5, 'def');
+grand("setgen", "clcg4");
+warning("off");
+grand("setsd",1,1,1,1);
+warning("on");
+expected = grand(4, 6, "def");
+warning("off");
+grand("setsd",1,1,1,1); // Resetting the seed to obtain the same results
+warning("on");
+computed = grand(4, 6, 5, "def");
 assert_checkequal(expected, computed(:, :, 1));
-warning('off');
-grand('setsd',1,1,1,1); // Resetting the seed to obtain the same results
-warning('on');
-expected = grand(4, 6, 'lgi');
-warning('off');
-grand('setsd',1,1,1,1); // Resetting the seed to obtain the same results
-warning('on');
-computed = grand(4, 6, 5, 'lgi');
-assert_checkequal(expected, computed(:, :, 1));
-// fsultra generator
-grand('setgen', 'fsultra');
-grand('setsd', 1, 1);
-expected = grand(4, 6, 'def');
-grand('setsd', 1, 1); // Resetting the seed to obtain the same results
-computed = grand(4, 6, 5, 'def');
-assert_checkequal(expected, computed(:, :, 1));
-grand('setsd', 1, 1);
-expected = grand(4, 6, 'lgi');
-grand('setsd', 1, 1); // Resetting the seed to obtain the same results
-computed = grand(4, 6, 5, 'lgi');
+warning("off");
+grand("setsd",1,1,1,1); // Resetting the seed to obtain the same results
+warning("on");
+expected = grand(4, 6, "lgi");
+warning("off");
+grand("setsd",1,1,1,1); // Resetting the seed to obtain the same results
+warning("on");
+computed = grand(4, 6, 5, "lgi");
 assert_checkequal(expected, computed(:, :, 1));
 // urand generator
-grand('setgen', 'urand');
-grand('setsd', 1);
-expected = grand(4, 6, 'def');
-grand('setsd', 1); // Resetting the seed to obtain the same results
-computed = grand(4, 6, 5, 'def');
+grand("setgen", "urand");
+grand("setsd", 1);
+expected = grand(4, 6, "def");
+grand("setsd", 1); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, "def");
 assert_checkequal(expected, computed(:, :, 1));
-grand('setsd', 1);
-expected = grand(4, 6, 'lgi');
-grand('setsd', 1); // Resetting the seed to obtain the same results
-computed = grand(4, 6, 5, 'lgi');
+grand("setsd", 1);
+expected = grand(4, 6, "lgi");
+grand("setsd", 1); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, "lgi");
 assert_checkequal(expected, computed(:, :, 1));
index 08702c6..3d760da 100644 (file)
@@ -14,7 +14,7 @@
 ///////////////////////////////////////////////////////////////////////////////
 //
 // Dimensions
-mat = grand(100, 101, 102, 'unf', 0, 1);
+mat = grand(100, 101, 102, "unf", 0, 1);
 assert_checktrue(size(mat) == [100 101 102]);
 
 ///////////////////////////////////////////////////////////////////////////////
@@ -23,87 +23,74 @@ assert_checktrue(size(mat) == [100 101 102]);
 // The grand(i, j, ...) should be equal to the first element of grand(i, j, k, ...).
 
 // mt generator
-grand('setgen', 'mt');
-grand('setsd', 0);
-expected = grand(4, 6, 'def');
-grand('setsd', 0); // Resetting the seed to obtain the same results
-computed = grand(4, 6, 5, 'def');
+grand("setgen", "mt");
+grand("setsd", 0);
+expected = grand(4, 6, "def");
+grand("setsd", 0); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, "def");
 assert_checkequal(expected, computed(:, :, 1));
-grand('setsd', 0);
-expected = grand(4, 6, 'lgi');
-grand('setsd', 0); // Resetting the seed to obtain the same results
-computed = grand(4, 6, 5, 'lgi');
+grand("setsd", 0);
+expected = grand(4, 6, "lgi");
+grand("setsd", 0); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, "lgi");
 assert_checkequal(expected, computed(:, :, 1));
 
 // kiss generator
-grand('setgen', 'kiss');
-grand('setsd', 0, 0, 0, 0);
-expected = grand(4, 6, 'def');
-grand('setsd', 0, 0, 0, 0); // Resetting the seed to obtain the same results
-computed = grand(4, 6, 5, 'def');
+grand("setgen", "kiss");
+grand("setsd", 0, 0, 0, 0);
+expected = grand(4, 6, "def");
+grand("setsd", 0, 0, 0, 0); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, "def");
 assert_checkequal(expected, computed(:, :, 1));
-grand('setsd', 0, 0, 0, 0);
-expected = grand(4, 6, 'lgi');
-grand('setsd', 0, 0, 0, 0); // Resetting the seed to obtain the same results
-computed = grand(4, 6, 5, 'lgi');
+grand("setsd", 0, 0, 0, 0);
+expected = grand(4, 6, "lgi");
+grand("setsd", 0, 0, 0, 0); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, "lgi");
 assert_checkequal(expected, computed(:, :, 1));
 
 // clcg2 generator
-grand('setgen', 'clcg2');
-grand('setsd', 1, 1);
-expected = grand(4, 6, 'def');
-grand('setsd', 1, 1); // Resetting the seed to obtain the same results
-computed = grand(4, 6, 5, 'def');
+grand("setgen", "clcg2");
+grand("setsd", 1, 1);
+expected = grand(4, 6, "def");
+grand("setsd", 1, 1); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, "def");
 assert_checkequal(expected, computed(:, :, 1));
-grand('setsd', 1, 1);
-expected = grand(4, 6, 'lgi');
-grand('setsd', 1, 1); // Resetting the seed to obtain the same results
-computed = grand(4, 6, 5, 'lgi');
+grand("setsd", 1, 1);
+expected = grand(4, 6, "lgi");
+grand("setsd", 1, 1); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, "lgi");
 assert_checkequal(expected, computed(:, :, 1));
 
 // clcg4 generator
-grand('setgen', 'clcg4');
-warning('off');
-grand('setsd',1,1,1,1);
-warning('on');
-expected = grand(4, 6, 'def');
-warning('off');
-grand('setsd',1,1,1,1); // Resetting the seed to obtain the same results
-warning('on');
-computed = grand(4, 6, 5, 'def');
+grand("setgen", "clcg4");
+warning("off");
+grand("setsd",1,1,1,1);
+warning("on");
+expected = grand(4, 6, "def");
+warning("off");
+grand("setsd",1,1,1,1); // Resetting the seed to obtain the same results
+warning("on");
+computed = grand(4, 6, 5, "def");
 assert_checkequal(expected, computed(:, :, 1));
-warning('off');
-grand('setsd',1,1,1,1); // Resetting the seed to obtain the same results
-warning('on');
-expected = grand(4, 6, 'lgi');
-warning('off');
-grand('setsd',1,1,1,1); // Resetting the seed to obtain the same results
-warning('on');
-computed = grand(4, 6, 5, 'lgi');
-assert_checkequal(expected, computed(:, :, 1));
-
-// fsultra generator
-grand('setgen', 'fsultra');
-grand('setsd', 1, 1);
-expected = grand(4, 6, 'def');
-grand('setsd', 1, 1); // Resetting the seed to obtain the same results
-computed = grand(4, 6, 5, 'def');
-assert_checkequal(expected, computed(:, :, 1));
-grand('setsd', 1, 1);
-expected = grand(4, 6, 'lgi');
-grand('setsd', 1, 1); // Resetting the seed to obtain the same results
-computed = grand(4, 6, 5, 'lgi');
+warning("off");
+grand("setsd",1,1,1,1); // Resetting the seed to obtain the same results
+warning("on");
+expected = grand(4, 6, "lgi");
+warning("off");
+grand("setsd",1,1,1,1); // Resetting the seed to obtain the same results
+warning("on");
+computed = grand(4, 6, 5, "lgi");
 assert_checkequal(expected, computed(:, :, 1));
 
 // urand generator
-grand('setgen', 'urand');
-grand('setsd', 1);
-expected = grand(4, 6, 'def');
-grand('setsd', 1); // Resetting the seed to obtain the same results
-computed = grand(4, 6, 5, 'def');
+grand("setgen", "urand");
+grand("setsd", 1);
+expected = grand(4, 6, "def");
+grand("setsd", 1); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, "def");
 assert_checkequal(expected, computed(:, :, 1));
-grand('setsd', 1);
-expected = grand(4, 6, 'lgi');
-grand('setsd', 1); // Resetting the seed to obtain the same results
-computed = grand(4, 6, 5, 'lgi');
+grand("setsd", 1);
+expected = grand(4, 6, "lgi");
+grand("setsd", 1); // Resetting the seed to obtain the same results
+computed = grand(4, 6, 5, "lgi");
 assert_checkequal(expected, computed(:, :, 1));