Old code under the ACM license (the last one) replaced by some code found in freemat...
Delphine GASC [Wed, 9 Jul 2008 13:11:35 +0000 (13:11 +0000)]
14 files changed:
scilab/modules/randlib/Makefile.am
scilab/modules/randlib/Makefile.in
scilab/modules/randlib/src/c/genbet.c [new file with mode: 0644]
scilab/modules/randlib/src/c/ignbin.c [new file with mode: 0644]
scilab/modules/randlib/src/c/ignpoi.c [new file with mode: 0644]
scilab/modules/randlib/src/c/sexpo.c [new file with mode: 0644]
scilab/modules/randlib/src/c/sgamma.c [new file with mode: 0644]
scilab/modules/randlib/src/fortran/genbet.f [deleted file]
scilab/modules/randlib/src/fortran/ignbin.f [deleted file]
scilab/modules/randlib/src/fortran/ignpoi.f [deleted file]
scilab/modules/randlib/src/fortran/sexpo.f [deleted file]
scilab/modules/randlib/src/fortran/sgamma.f [deleted file]
scilab/modules/randlib/tests/unit_tests/grand.dia.ref
scilab/modules/randlib/tests/unit_tests/grand.tst

index 0391d3b..bd2b214 100644 (file)
@@ -9,32 +9,37 @@ src/c/igngeom.c \
 src/c/kiss.c \
 src/c/urand.c \
 src/c/clcg2.c \
-src/c/clcg4.c
+src/c/clcg4.c \
+src/c/sexpo.c \
+src/c/ignbin.c \
+src/c/ignpoi.c \
+src/c/sgamma.c \
+src/c/genbet.c 
 
 RANDLIB_FORTRAN_SOURCES = \
 src/fortran/gennf.f \
 src/fortran/genchi.f \
 src/fortran/setgmn.f \
 src/fortran/gengam.f \
-src/fortran/sexpo.f \
 src/fortran/snorm.f \
 src/fortran/gennch.f \
 src/fortran/genmn.f \
-src/fortran/genbet.f \
 src/fortran/phrtsd.f \
-src/fortran/sgamma.f \
-src/fortran/ignbin.f \
 src/fortran/spofa.f \
 src/fortran/sdot.f \
 src/fortran/ignnbn.f \
 src/fortran/lennob.f \
 src/fortran/genf.f \
 src/fortran/genunf.f \
-src/fortran/ignpoi.f \
 src/fortran/genexp.f \
 src/fortran/genmul.f \
 src/fortran/genprm.f \
 src/fortran/gennor.f
+#src/fortran/genbet.f
+#src/fortran/sexpo.f 
+#src/fortran/ignbin.f
+#src/fortran/ignpoi.f
+#src/fortran/sgamma.f
 
 GATEWAY_C_SOURCES = sci_gateway/c/gw_randlib.c \
 sci_gateway/c/sci_grand.c
index 6762b4a..5aec260 100644 (file)
@@ -100,13 +100,15 @@ libscirandlib_la_DEPENDENCIES =  \
 am__objects_1 = libscirandlib_la-fsultra.lo libscirandlib_la-mt.lo \
        libscirandlib_la-igngeom.lo libscirandlib_la-kiss.lo \
        libscirandlib_la-urand.lo libscirandlib_la-clcg2.lo \
-       libscirandlib_la-clcg4.lo
+       libscirandlib_la-clcg4.lo libscirandlib_la-sexpo.lo \
+       libscirandlib_la-ignbin.lo libscirandlib_la-ignpoi.lo \
+       libscirandlib_la-sgamma.lo libscirandlib_la-genbet.lo
 am__objects_2 = libscirandlib_la-gw_randlib.lo \
        libscirandlib_la-sci_grand.lo
-am__objects_3 = gennf.lo genchi.lo setgmn.lo gengam.lo sexpo.lo \
-       snorm.lo gennch.lo genmn.lo genbet.lo phrtsd.lo sgamma.lo \
-       ignbin.lo spofa.lo sdot.lo ignnbn.lo lennob.lo genf.lo \
-       genunf.lo ignpoi.lo genexp.lo genmul.lo genprm.lo gennor.lo
+am__objects_3 = gennf.lo genchi.lo setgmn.lo gengam.lo snorm.lo \
+       gennch.lo genmn.lo phrtsd.lo spofa.lo sdot.lo ignnbn.lo \
+       lennob.lo genf.lo genunf.lo genexp.lo genmul.lo genprm.lo \
+       gennor.lo
 am_libscirandlib_la_OBJECTS = $(am__objects_1) $(am__objects_2) \
        $(am__objects_3)
 libscirandlib_la_OBJECTS = $(am_libscirandlib_la_OBJECTS)
@@ -339,33 +341,38 @@ src/c/igngeom.c \
 src/c/kiss.c \
 src/c/urand.c \
 src/c/clcg2.c \
-src/c/clcg4.c
+src/c/clcg4.c \
+src/c/sexpo.c \
+src/c/ignbin.c \
+src/c/ignpoi.c \
+src/c/sgamma.c \
+src/c/genbet.c 
 
 RANDLIB_FORTRAN_SOURCES = \
 src/fortran/gennf.f \
 src/fortran/genchi.f \
 src/fortran/setgmn.f \
 src/fortran/gengam.f \
-src/fortran/sexpo.f \
 src/fortran/snorm.f \
 src/fortran/gennch.f \
 src/fortran/genmn.f \
-src/fortran/genbet.f \
 src/fortran/phrtsd.f \
-src/fortran/sgamma.f \
-src/fortran/ignbin.f \
 src/fortran/spofa.f \
 src/fortran/sdot.f \
 src/fortran/ignnbn.f \
 src/fortran/lennob.f \
 src/fortran/genf.f \
 src/fortran/genunf.f \
-src/fortran/ignpoi.f \
 src/fortran/genexp.f \
 src/fortran/genmul.f \
 src/fortran/genprm.f \
 src/fortran/gennor.f
 
+#src/fortran/genbet.f
+#src/fortran/sexpo.f 
+#src/fortran/ignbin.f
+#src/fortran/ignpoi.f
+#src/fortran/sgamma.f
 GATEWAY_C_SOURCES = sci_gateway/c/gw_randlib.c \
 sci_gateway/c/sci_grand.c
 
@@ -546,11 +553,16 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libscirandlib_la-clcg2.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libscirandlib_la-clcg4.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libscirandlib_la-fsultra.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libscirandlib_la-genbet.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libscirandlib_la-gw_randlib.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libscirandlib_la-ignbin.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libscirandlib_la-igngeom.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libscirandlib_la-ignpoi.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libscirandlib_la-kiss.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libscirandlib_la-mt.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libscirandlib_la-sci_grand.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libscirandlib_la-sexpo.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libscirandlib_la-sgamma.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libscirandlib_la-urand.Plo@am__quote@
 
 .c.o:
@@ -623,6 +635,41 @@ libscirandlib_la-clcg4.lo: src/c/clcg4.c
 @AMDEP_TRUE@@am__fastdepCC_FALSE@      DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCC_FALSE@  $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libscirandlib_la_CFLAGS) $(CFLAGS) -c -o libscirandlib_la-clcg4.lo `test -f 'src/c/clcg4.c' || echo '$(srcdir)/'`src/c/clcg4.c
 
+libscirandlib_la-sexpo.lo: src/c/sexpo.c
+@am__fastdepCC_TRUE@   $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libscirandlib_la_CFLAGS) $(CFLAGS) -MT libscirandlib_la-sexpo.lo -MD -MP -MF $(DEPDIR)/libscirandlib_la-sexpo.Tpo -c -o libscirandlib_la-sexpo.lo `test -f 'src/c/sexpo.c' || echo '$(srcdir)/'`src/c/sexpo.c
+@am__fastdepCC_TRUE@   mv -f $(DEPDIR)/libscirandlib_la-sexpo.Tpo $(DEPDIR)/libscirandlib_la-sexpo.Plo
+@AMDEP_TRUE@@am__fastdepCC_FALSE@      source='src/c/sexpo.c' object='libscirandlib_la-sexpo.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCC_FALSE@      DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCC_FALSE@  $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libscirandlib_la_CFLAGS) $(CFLAGS) -c -o libscirandlib_la-sexpo.lo `test -f 'src/c/sexpo.c' || echo '$(srcdir)/'`src/c/sexpo.c
+
+libscirandlib_la-ignbin.lo: src/c/ignbin.c
+@am__fastdepCC_TRUE@   $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libscirandlib_la_CFLAGS) $(CFLAGS) -MT libscirandlib_la-ignbin.lo -MD -MP -MF $(DEPDIR)/libscirandlib_la-ignbin.Tpo -c -o libscirandlib_la-ignbin.lo `test -f 'src/c/ignbin.c' || echo '$(srcdir)/'`src/c/ignbin.c
+@am__fastdepCC_TRUE@   mv -f $(DEPDIR)/libscirandlib_la-ignbin.Tpo $(DEPDIR)/libscirandlib_la-ignbin.Plo
+@AMDEP_TRUE@@am__fastdepCC_FALSE@      source='src/c/ignbin.c' object='libscirandlib_la-ignbin.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCC_FALSE@      DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCC_FALSE@  $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libscirandlib_la_CFLAGS) $(CFLAGS) -c -o libscirandlib_la-ignbin.lo `test -f 'src/c/ignbin.c' || echo '$(srcdir)/'`src/c/ignbin.c
+
+libscirandlib_la-ignpoi.lo: src/c/ignpoi.c
+@am__fastdepCC_TRUE@   $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libscirandlib_la_CFLAGS) $(CFLAGS) -MT libscirandlib_la-ignpoi.lo -MD -MP -MF $(DEPDIR)/libscirandlib_la-ignpoi.Tpo -c -o libscirandlib_la-ignpoi.lo `test -f 'src/c/ignpoi.c' || echo '$(srcdir)/'`src/c/ignpoi.c
+@am__fastdepCC_TRUE@   mv -f $(DEPDIR)/libscirandlib_la-ignpoi.Tpo $(DEPDIR)/libscirandlib_la-ignpoi.Plo
+@AMDEP_TRUE@@am__fastdepCC_FALSE@      source='src/c/ignpoi.c' object='libscirandlib_la-ignpoi.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCC_FALSE@      DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCC_FALSE@  $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libscirandlib_la_CFLAGS) $(CFLAGS) -c -o libscirandlib_la-ignpoi.lo `test -f 'src/c/ignpoi.c' || echo '$(srcdir)/'`src/c/ignpoi.c
+
+libscirandlib_la-sgamma.lo: src/c/sgamma.c
+@am__fastdepCC_TRUE@   $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libscirandlib_la_CFLAGS) $(CFLAGS) -MT libscirandlib_la-sgamma.lo -MD -MP -MF $(DEPDIR)/libscirandlib_la-sgamma.Tpo -c -o libscirandlib_la-sgamma.lo `test -f 'src/c/sgamma.c' || echo '$(srcdir)/'`src/c/sgamma.c
+@am__fastdepCC_TRUE@   mv -f $(DEPDIR)/libscirandlib_la-sgamma.Tpo $(DEPDIR)/libscirandlib_la-sgamma.Plo
+@AMDEP_TRUE@@am__fastdepCC_FALSE@      source='src/c/sgamma.c' object='libscirandlib_la-sgamma.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCC_FALSE@      DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCC_FALSE@  $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libscirandlib_la_CFLAGS) $(CFLAGS) -c -o libscirandlib_la-sgamma.lo `test -f 'src/c/sgamma.c' || echo '$(srcdir)/'`src/c/sgamma.c
+
+libscirandlib_la-genbet.lo: src/c/genbet.c
+@am__fastdepCC_TRUE@   $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libscirandlib_la_CFLAGS) $(CFLAGS) -MT libscirandlib_la-genbet.lo -MD -MP -MF $(DEPDIR)/libscirandlib_la-genbet.Tpo -c -o libscirandlib_la-genbet.lo `test -f 'src/c/genbet.c' || echo '$(srcdir)/'`src/c/genbet.c
+@am__fastdepCC_TRUE@   mv -f $(DEPDIR)/libscirandlib_la-genbet.Tpo $(DEPDIR)/libscirandlib_la-genbet.Plo
+@AMDEP_TRUE@@am__fastdepCC_FALSE@      source='src/c/genbet.c' object='libscirandlib_la-genbet.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCC_FALSE@      DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCC_FALSE@  $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libscirandlib_la_CFLAGS) $(CFLAGS) -c -o libscirandlib_la-genbet.lo `test -f 'src/c/genbet.c' || echo '$(srcdir)/'`src/c/genbet.c
+
 libscirandlib_la-gw_randlib.lo: sci_gateway/c/gw_randlib.c
 @am__fastdepCC_TRUE@   $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(libscirandlib_la_CFLAGS) $(CFLAGS) -MT libscirandlib_la-gw_randlib.lo -MD -MP -MF $(DEPDIR)/libscirandlib_la-gw_randlib.Tpo -c -o libscirandlib_la-gw_randlib.lo `test -f 'sci_gateway/c/gw_randlib.c' || echo '$(srcdir)/'`sci_gateway/c/gw_randlib.c
 @am__fastdepCC_TRUE@   mv -f $(DEPDIR)/libscirandlib_la-gw_randlib.Tpo $(DEPDIR)/libscirandlib_la-gw_randlib.Plo
@@ -658,9 +705,6 @@ setgmn.lo: src/fortran/setgmn.f
 gengam.lo: src/fortran/gengam.f
        $(LIBTOOL) --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o gengam.lo `test -f 'src/fortran/gengam.f' || echo '$(srcdir)/'`src/fortran/gengam.f
 
-sexpo.lo: src/fortran/sexpo.f
-       $(LIBTOOL) --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o sexpo.lo `test -f 'src/fortran/sexpo.f' || echo '$(srcdir)/'`src/fortran/sexpo.f
-
 snorm.lo: src/fortran/snorm.f
        $(LIBTOOL) --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o snorm.lo `test -f 'src/fortran/snorm.f' || echo '$(srcdir)/'`src/fortran/snorm.f
 
@@ -670,18 +714,9 @@ gennch.lo: src/fortran/gennch.f
 genmn.lo: src/fortran/genmn.f
        $(LIBTOOL) --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o genmn.lo `test -f 'src/fortran/genmn.f' || echo '$(srcdir)/'`src/fortran/genmn.f
 
-genbet.lo: src/fortran/genbet.f
-       $(LIBTOOL) --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o genbet.lo `test -f 'src/fortran/genbet.f' || echo '$(srcdir)/'`src/fortran/genbet.f
-
 phrtsd.lo: src/fortran/phrtsd.f
        $(LIBTOOL) --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o phrtsd.lo `test -f 'src/fortran/phrtsd.f' || echo '$(srcdir)/'`src/fortran/phrtsd.f
 
-sgamma.lo: src/fortran/sgamma.f
-       $(LIBTOOL) --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o sgamma.lo `test -f 'src/fortran/sgamma.f' || echo '$(srcdir)/'`src/fortran/sgamma.f
-
-ignbin.lo: src/fortran/ignbin.f
-       $(LIBTOOL) --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o ignbin.lo `test -f 'src/fortran/ignbin.f' || echo '$(srcdir)/'`src/fortran/ignbin.f
-
 spofa.lo: src/fortran/spofa.f
        $(LIBTOOL) --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o spofa.lo `test -f 'src/fortran/spofa.f' || echo '$(srcdir)/'`src/fortran/spofa.f
 
@@ -700,9 +735,6 @@ genf.lo: src/fortran/genf.f
 genunf.lo: src/fortran/genunf.f
        $(LIBTOOL) --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o genunf.lo `test -f 'src/fortran/genunf.f' || echo '$(srcdir)/'`src/fortran/genunf.f
 
-ignpoi.lo: src/fortran/ignpoi.f
-       $(LIBTOOL) --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o ignpoi.lo `test -f 'src/fortran/ignpoi.f' || echo '$(srcdir)/'`src/fortran/ignpoi.f
-
 genexp.lo: src/fortran/genexp.f
        $(LIBTOOL) --tag=F77 $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(F77) $(AM_FFLAGS) $(FFLAGS) -c -o genexp.lo `test -f 'src/fortran/genexp.f' || echo '$(srcdir)/'`src/fortran/genexp.f
 
diff --git a/scilab/modules/randlib/src/c/genbet.c b/scilab/modules/randlib/src/c/genbet.c
new file mode 100644 (file)
index 0000000..c20e4bc
--- /dev/null
@@ -0,0 +1,237 @@
+#include "stack-c.h"
+#include <math.h>
+#include "grand.h"
+#include "machine.h" 
+#include "core_math.h"
+
+double C2F(genbet)(double *aa,double *bb)
+/*
+**********************************************************************
+This source code was taken in the project "freemat"(BSD license)
+This source code was modified by Gaüzère Sabine according to the 
+modifications done by JJV
+
+     float genbet(float aa,float bb)
+               GeNerate BETa random deviate
+                              Function
+     Returns a single random deviate from the beta distribution with
+     parameters A and B.  The density of the beta is
+               x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 < x < 1
+
+                              Arguments
+
+     A --> First parameter of the beta distribution
+                         DOUBLE PRECISION A
+     JJV                 (A > 1.0E-37)
+
+     B --> Second parameter of the beta distribution
+                         DOUBLE PRECISION B
+     JJV                 (B > 1.0E-37)
+
+                              Method
+     R. C. H. Cheng
+     Generating Beta Variatew with Nonintegral Shape Parameters
+     Communications of the ACM, 21:317-322  (1978)
+     (Algorithms BB and BC)
+**********************************************************************
+*/
+{
+#define expmax 87.498222998046875
+#define infnty 1.0E38
+#define minlog 1.0E-37
+
+static double olda = -1.0E37;
+static double oldb = -1.0E37;
+static double genbet,a,alpha,b,beta,delta,gamma,k1,k2,r,s,t,u1,u2,v,w,y,z;
+static int qsame;
+
+    qsame = olda == *aa && oldb == *bb;
+    if(qsame) goto S20;
+    //if(!(aa <= 0.0 || bb <= 0.0)) goto S10;
+    //    Scierror(999, "RANLIB Error: AA or BB <= 0 in GENBET - Abort!");
+//  JJV added small minimum for small log problem in calc of W
+//  in Rand.c
+S10:
+    olda = *aa;
+    oldb = *bb;
+S20:
+    if(!(Min(*aa,*bb) > 1.0)) goto S100;
+/*
+     Alborithm BB
+     Initialize
+*/
+    if(qsame) goto S30;
+    a = Min(*aa,*bb);
+    b = Max(*aa,*bb);
+    alpha = a+b;
+    beta = sqrt((alpha-2.0)/(2.0*a*b-alpha));
+    gamma = a+1.0/beta;
+S30:
+S40:
+    u1 = C2F(ranf)();
+/*
+     Step 1
+*/
+    u2 = C2F(ranf)();
+    v = beta*log(u1/(1.0-u1));
+//  JJV altered this
+    if(v > expmax) goto S55;
+//  JJV added checker to see if a*exp(v) will overflow
+//  JJV 50 _was_ w = a*exp(v); also note here a > 1.0
+//    if(!(v > expmax)) goto S50;
+//    w = infnty;
+//    goto S60;
+S50:
+    w= exp(v);
+    if(w > (infnty/a)) goto S55;
+    w *= a;
+    goto S60;
+S55:
+    w = infnty;
+S60:
+    z = pow(u1,2.0f)*u2;
+    r = gamma*v-1.3862943649291992188;
+    s = a+r-w;
+/*
+     Step 2
+*/
+    if(s+2.6094379425048828125 >= 5.0*z) goto S70;
+/*
+     Step 3
+*/
+    t = log(z);
+    if(s > t) goto S70;
+/*
+     Step 4
+*/
+//    JJV added checker to see if log(alpha/(b+w)) will 
+//    JJV overflow.  If so, we count the log as -INF, and
+//    JJV consequently evaluate conditional as true, i.e.
+//    JJV the algorithm rejects the trial and starts over
+//    JJV May not need this here since ALPHA > 2.0
+    if(alpha/(b+w) < minlog) goto S40;
+    if(r+alpha*log(alpha/(b+w)) < t) goto S40;
+S70:
+/*
+     Step 5
+*/
+    if(!(*aa == a)) goto S80;
+    genbet = w/(b+w);
+    goto S90;
+S80:
+    genbet = b/(b+w);
+S90:
+    goto S230;
+S100:
+/*
+     Algorithm BC
+     Initialize
+*/
+    if(qsame) goto S110;
+    a = Max(*aa,*bb);
+    b = Min(*aa,*bb);
+    alpha = a+b;
+    beta = 1.0/b;
+    delta = 1.0+a-b;
+    k1 = delta*(1.3888900168240070343E-2+4.1666701436042785645E-2*b)/(a*beta-0.77777802944183349609);
+    k2 = 0.25+(0.5+0.25/delta)*b;
+S110:
+S120:
+    u1 = C2F(ranf)();
+/*
+     Step 1
+*/
+    u2 = C2F(ranf)();
+    if(u1 >= 0.5) goto S130;
+/*
+     Step 2
+*/
+    y = u1*u2;
+    z = u1*y;
+    if(0.25*u2+z-y >= k1) goto S120;
+    goto S170;
+S130:
+/*
+     Step 3
+*/
+    z = pow(u1,2.0)*u2;
+    if(!(z <= 0.25)) goto S160;
+    v = beta*log(u1/(1.0-u1));
+
+//  JJV instead of checking v > expmax at top, I will check
+//  JJV if a < 1, then check the appropriate values
+
+    if(a > 1.0) goto S135;
+//  JJV A < 1 so it can help out if EXP(V) would overflow
+    if(v > expmax) goto S132;
+    w = a*exp(v);
+    goto S200;
+S132:
+    w = v + log(a);
+    if(w > expmax) goto S140;
+    w = exp(w);
+    goto S200;
+
+//  JJV in this case A > 1
+S135:
+    if(v > expmax) goto S140;
+    w = exp(v);
+    if(w > infnty/a) goto S140;
+    w *= a;
+    goto S200;
+S140:
+    w = infnty;
+    goto S200;
+//S150:
+//    goto S200;
+S160:
+    if(z >= k2) goto S120;
+S170:
+/*
+     Step 4
+     Step 5
+*/
+    v = beta*log(u1/(1.0-u1));
+//  JJV same kind of checking as above
+    if(a > 1.0) goto S175;
+//  JJV A < 1 so it can help out if EXP(V) would overflow
+    if(v > expmax) goto S172;
+    w = a*exp(v);
+    goto S190;
+S172:
+    w = v + log(a);
+    if(w > expmax) goto S180;
+    w = exp(w);
+    goto S190;
+//  JJV in this case A > 1
+S175:
+    if(v > expmax) goto S180;
+    w = exp(v);
+    if(w > infnty/a) goto S180;
+    w *= a;
+    goto S190;
+
+S180:
+    w = infnty;
+
+//  JJV here we also check to see if log overlows; if so, we treat it
+//  JJV as -INF, which means condition is true, i.e. restart
+S190:
+    if(alpha/(b+w) < minlog) goto S120;
+    if(alpha*(log(alpha/(b+w))+v)-1.3862943649291992188 < log(z)) goto S120;
+S200:
+/*
+     Step 6
+*/
+    if(!(a == *aa)) goto S210;
+    genbet = w/(b+w);
+    goto S220;
+S210:
+    genbet = b/(b+w);
+S220:
+S230:
+    return genbet;
+#undef expmax
+#undef infnty
+#undef minlog
+}
diff --git a/scilab/modules/randlib/src/c/ignbin.c b/scilab/modules/randlib/src/c/ignbin.c
new file mode 100644 (file)
index 0000000..70c10b8
--- /dev/null
@@ -0,0 +1,273 @@
+#include "stack-c.h"
+#include <math.h>
+#include "grand.h"
+#include "machine.h" 
+#include "core_math.h"
+
+int C2F(ignbin)(int *n,double *pp)
+/*
+**********************************************************************
+This source code was taken in the project "freemat"(BSD license)
+This source code was modified by Gaüzère Sabine according to the 
+modifications done by JJV
+
+     long ignbin(long n,float pp)
+                    GENerate BINomial random deviate
+                              Function
+     Generates a single random deviate from a binomial
+     distribution whose number of trials is N and whose
+     probability of an event in each trial is P.
+                              Arguments
+     n  --> The number of trials in the binomial distribution
+            from which a random deviate is to be generated.
+     p  --> The probability of an event in each trial of the
+            binomial distribution from which a random deviate
+            is to be generated.
+     ignbin <-- A random deviate yielding the number of events
+                from N independent trials, each of which has
+                a probability of event P.
+                              Method
+     This is algorithm BTPE from:
+         Kachitvichyanukul, V. and Schmeiser, B. W.
+         Binomial Random Variate Generation.
+         Communications of the ACM, 31, 2
+         (February, 1988) 216.
+**********************************************************************
+     SUBROUTINE BTPEC(N,PP,ISEED,JX)
+     BINOMIAL RANDOM VARIATE GENERATOR
+     MEAN .LT. 30 -- INVERSE CDF
+       MEAN .GE. 30 -- ALGORITHM BTPE:  ACCEPTANCE-REJECTION VIA
+       FOUR REGION COMPOSITION.  THE FOUR REGIONS ARE A TRIANGLE
+       (SYMMETRIC IN THE CENTER), A PAIR OF PARALLELOGRAMS (ABOVE
+       THE TRIANGLE), AND EXPONENTIAL LEFT AND RIGHT TAILS.
+     BTPE REFERS TO BINOMIAL-TRIANGLE-PARALLELOGRAM-EXPONENTIAL.
+     BTPEC REFERS TO BTPE AND "COMBINED."  THUS BTPE IS THE
+       RESEARCH AND BTPEC IS THE IMPLEMENTATION OF A COMPLETE
+       USABLE ALGORITHM.
+     REFERENCE:  VORATAS KACHITVICHYANUKUL AND BRUCE SCHMEISER,
+       "BINOMIAL RANDOM VARIATE GENERATION,"
+       COMMUNICATIONS OF THE ACM, FORTHCOMING
+     WRITTEN:  SEPTEMBER 1980.
+       LAST REVISED:  MAY 1985, JULY 1987
+     REQUIRED SUBPROGRAM:  RAND() -- A UNIFORM (0,1) RANDOM NUMBER
+                           GENERATOR
+     ARGUMENTS
+       N : NUMBER OF BERNOULLI TRIALS            (INPUT)
+       PP : PROBABILITY OF SUCCESS IN EACH TRIAL (INPUT)
+       ISEED:  RANDOM NUMBER SEED                (INPUT AND OUTPUT)
+       JX:  RANDOMLY GENERATED OBSERVATION       (OUTPUT)
+     VARIABLES
+       PSAVE: VALUE OF PP FROM THE LAST CALL TO BTPEC
+       NSAVE: VALUE OF N FROM THE LAST CALL TO BTPEC
+       XNP:  VALUE OF THE MEAN FROM THE LAST CALL TO BTPEC
+       P: PROBABILITY USED IN THE GENERATION PHASE OF BTPEC
+       FFM: TEMPORARY VARIABLE EQUAL TO XNP + P
+       M:  INTEGER VALUE OF THE CURRENT MODE
+       FM:  FLOATING POINT VALUE OF THE CURRENT MODE
+       XNPQ: TEMPORARY VARIABLE USED IN SETUP AND SQUEEZING STEPS
+       P1:  AREA OF THE TRIANGLE
+       C:  HEIGHT OF THE PARALLELOGRAMS
+       XM:  CENTER OF THE TRIANGLE
+       XL:  LEFT END OF THE TRIANGLE
+       XR:  RIGHT END OF THE TRIANGLE
+       AL:  TEMPORARY VARIABLE
+       XLL:  RATE FOR THE LEFT EXPONENTIAL TAIL
+       XLR:  RATE FOR THE RIGHT EXPONENTIAL TAIL
+       P2:  AREA OF THE PARALLELOGRAMS
+       P3:  AREA OF THE LEFT EXPONENTIAL TAIL
+       P4:  AREA OF THE RIGHT EXPONENTIAL TAIL
+       U:  A U(0,P4) RANDOM VARIATE USED FIRST TO SELECT ONE OF THE
+           FOUR REGIONS AND THEN CONDITIONALLY TO GENERATE A VALUE
+           FROM THE REGION
+       V:  A U(0,1) RANDOM NUMBER USED TO GENERATE THE RANDOM VALUE
+           (REGION 1) OR TRANSFORMED INTO THE VARIATE TO ACCEPT OR
+           REJECT THE CANDIDATE VALUE
+       IX:  INTEGER CANDIDATE VALUE
+       X:  PRELIMINARY CONTINUOUS CANDIDATE VALUE IN REGION 2 LOGIC
+           AND A FLOATING POINT IX IN THE ACCEPT/REJECT LOGIC
+       K:  ABSOLUTE VALUE OF (IX-M)
+       F:  THE HEIGHT OF THE SCALED DENSITY FUNCTION USED IN THE
+           ACCEPT/REJECT DECISION WHEN BOTH M AND IX ARE SMALL
+           ALSO USED IN THE INVERSE TRANSFORMATION
+       R: THE RATIO P/Q
+       G: CONSTANT USED IN CALCULATION OF PROBABILITY
+       MP:  MODE PLUS ONE, THE LOWER INDEX FOR EXPLICIT CALCULATION
+            OF F WHEN IX IS GREATER THAN M
+       IX1:  CANDIDATE VALUE PLUS ONE, THE LOWER INDEX FOR EXPLICIT
+             CALCULATION OF F WHEN IX IS LESS THAN M
+       I:  INDEX FOR EXPLICIT CALCULATION OF F FOR BTPE
+       AMAXP: MAXIMUM ERROR OF THE LOGARITHM OF NORMAL BOUND
+       YNORM: LOGARITHM OF NORMAL BOUND
+       ALV:  NATURAL LOGARITHM OF THE ACCEPT/REJECT VARIATE V
+       X1,F1,Z,W,Z2,X2,F2, AND W2 ARE TEMPORARY VARIABLES TO BE
+       USED IN THE FINAL ACCEPT/REJECT TEST
+       QN: PROBABILITY OF NO SUCCESS IN N TRIALS
+     REMARK
+       IX AND JX COULD LOGICALLY BE THE SAME VARIABLE, WHICH WOULD
+       SAVE A MEMORY POSITION AND A LINE OF CODE.  HOWEVER, SOME
+       COMPILERS (E.G.,CDC MNF) OPTIMIZE BETTER WHEN THE ARGUMENTS
+       ARE NOT INVOLVED.
+     ISEED NEEDS TO BE DOUBLE PRECISION IF THE IMSL ROUTINE
+     GGUBFS IS USED TO GENERATE UNIFORM RANDOM NUMBER, OTHERWISE
+     TYPE OF ISEED SHOULD BE DICTATED BY THE UNIFORM GENERATOR
+**********************************************************************
+*****DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY
+*/
+{
+static double psave = -1.0E37;
+static int nsave = -214748365;
+static int ignbin,i,ix,ix1,k,m,mp,T1;
+static double al,alv,amaxp,c,f,f1,f2,ffm,fm,g,p,p1,p2,p3,p4,q,qn,r,u,v,w,w2,x,x1,
+    x2,xl,xll,xlr,xm,xnp,xnpq,xr,ynorm,z,z2;
+
+    if(*pp != psave) goto S10;
+    if(*n != nsave) goto S20;
+    if(xnp < 30.0) goto S150;
+    goto S30;
+S10:
+/*
+*****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE
+*/
+/*   JJV added the argument checker - involved only renaming 10
+     JJV and 20 to the checkers and adding checkers
+     JJV Only remaining problem - if called initially with the
+     JJV initial values of psave and nsave, it will hang
+*/
+    psave = *pp;
+    p = Min(psave,1.0-psave);
+    q = 1.0-p;
+S20:
+    xnp = *n * p;
+    nsave = *n;
+    if(xnp < 30.0) goto S140;
+    ffm = xnp+p;
+    m = ffm;
+    fm = m;
+    xnpq = xnp*q;
+    p1 = (int) (2.195*sqrt(xnpq)-4.6*q)+0.5;
+    xm = fm+0.5;
+    xl = xm-p1;
+    xr = xm+p1;
+    c = 0.134+20.5/(15.3+fm);
+    al = (ffm-xl)/(ffm-xl*p);
+    xll = al*(1.0+0.5*al);
+    al = (xr-ffm)/(xr*q);
+    xlr = al*(1.0+0.5*al);
+    p2 = p1*(1.0+c+c);
+    p3 = p2+c/xll;
+    p4 = p3+c/xlr;
+S30:
+/*
+*****GENERATE VARIATE
+*/
+    u = C2F(ranf)()*p4;
+    v = C2F(ranf)();
+/*
+     TRIANGULAR REGION
+*/
+    if(u > p1) goto S40;
+    ix = xm-p1*v+u;
+    goto S170;
+S40:
+/*
+     PARALLELOGRAM REGION
+*/
+    if(u > p2) goto S50;
+    x = xl+(u-p1)/c;
+    v = v*c+1.0-Abs(xm-x)/p1;
+    if(v > 1.0 || v <= 0.0) goto S30;
+    ix = x;
+    goto S70;
+S50:
+/*
+     LEFT TAIL
+*/
+    if(u > p3) goto S60;
+    ix = xl+log(v)/xll;
+    if(ix < 0) goto S30;
+    v *= ((u-p2)*xll);
+    goto S70;
+S60:
+/*
+     RIGHT TAIL
+*/
+    ix = xr-log(v)/xlr;
+    if(ix > *n) goto S30;
+    v *= ((u-p3)*xlr);
+S70:
+/*
+*****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST
+*/
+    k = Abs(ix-m);
+    if(k > 20 && k < xnpq/2-1) goto S130;
+/*
+     EXPLICIT EVALUATION
+*/
+    f = 1.0;
+    r = p/q;
+    g = (*n+1)*r;
+    T1 = m-ix;
+    if(T1 < 0) goto S80;
+    else if(T1 == 0) goto S120;
+    else  goto S100;
+S80:
+    mp = m+1;
+    for(i=mp; i<=ix; i++) f *= (g/i-r);
+    goto S120;
+S100:
+    ix1 = ix+1;
+    for(i=ix1; i<=m; i++) f /= (g/i-r);
+S120:
+    if(v <= f) goto S170;
+    goto S30;
+S130:
+/*
+     SQUEEZING USING UPPER AND LOWER BOUNDS ON ALOG(F(X))
+*/
+    amaxp = k/xnpq*((k*(k/3.0+0.625)+0.1666666666666)/xnpq+0.5);
+    ynorm = -(k*k/(2.0*xnpq));
+    alv = log(v);
+    if(alv < ynorm-amaxp) goto S170;
+    if(alv > ynorm+amaxp) goto S30;
+/*
+     STIRLING'S FORMULA TO MACHINE ACCURACY FOR
+     THE FINAL ACCEPTANCE/REJECTION TEST
+*/
+    x1 = ix+1.0;
+    f1 = fm+1.0;
+    z = *n+1.0-fm;
+    w = *n-ix+1.0;
+    z2 = z*z;
+    x2 = x1*x1;
+    f2 = f1*f1;
+    w2 = w*w;
+    if(alv <= xm*log(f1/x1)+(*n-m+0.5)*log(z/w)+(ix-m)*log(w*p/(x1*q))+(13860.0-
+      (462.0-(132.0-(99.0-140.0/f2)/f2)/f2)/f2)/f1/166320.0+(13860.0-(462.0-
+      (132.0-(99.0-140.0/z2)/z2)/z2)/z2)/z/166320.0+(13860.0-(462.0-(132.0-
+      (99.0-140.0/x2)/x2)/x2)/x2)/x1/166320.0+(13860.0-(462.0-(132.0-(99.0
+      -140.0/w2)/w2)/w2)/w2)/w/166320.0) goto S170;
+    goto S30;
+S140:
+/*
+     INVERSE CDF LOGIC FOR MEAN LESS THAN 30
+*/
+    qn = pow((double)q,(double)*n);
+    r = p/q;
+    g = r*(*n+1);
+S150:
+    ix = 0;
+    f = qn;
+    u = C2F(ranf)();
+S160:
+    if(u < f) goto S170;
+    if(ix > 110) goto S150;
+    u -= f;
+    ix += 1;
+    f *= (g/ix-r);
+    goto S160;
+S170:
+    if(psave > 0.5) ix = *n-ix;
+    ignbin = ix;
+    return ignbin;
+}
+
diff --git a/scilab/modules/randlib/src/c/ignpoi.c b/scilab/modules/randlib/src/c/ignpoi.c
new file mode 100644 (file)
index 0000000..4573d98
--- /dev/null
@@ -0,0 +1,267 @@
+#include "stack-c.h"
+#include <math.h>
+#include "grand.h"
+#include "machine.h" 
+#include "core_math.h"
+
+double fsign (double x, double y)
+{
+if (y >= 0.0)
+return Abs(x);
+else 
+return -Abs(x);
+}
+
+int C2F(ignpoi)(double *mu)
+/*
+**********************************************************************
+This source code was taken in the project "freemat"(BSD license)
+This source code was modified by Gaüzère Sabine according to the 
+modifications done by JJV
+
+     long ignpoi(float mu)
+                    GENerate POIsson random deviate
+                              Function
+     Generates a single random deviate from a Poisson
+     distribution with mean AV.
+                              Arguments
+     av --> The mean of the Poisson distribution from which
+            a random deviate is to be generated.
+     genexp <-- The random deviate.
+                              Method
+     Renames KPOIS from TOMS as slightly modified by BWB to use RANF
+     instead of SUNIF.
+     For details see:
+               Ahrens, J.H. and Dieter, U.
+               Computer Generation of Poisson Deviates
+               From Modified Normal Distributions.
+               ACM Trans. Math. Software, 8, 2
+               (June 1982),163-179
+**********************************************************************
+**********************************************************************
+
+
+     P O I S S O N  DISTRIBUTION
+
+
+**********************************************************************
+**********************************************************************
+
+     FOR DETAILS SEE:
+
+               AHRENS, J.H. AND DIETER, U.
+               COMPUTER GENERATION OF POISSON DEVIATES
+               FROM MODIFIED NORMAL DISTRIBUTIONS.
+               ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179.
+
+     (SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE)
+
+**********************************************************************
+      INTEGER FUNCTION IGNPOI(IR,MU)
+     INPUT:  IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR
+             MU=MEAN MU OF THE POISSON DISTRIBUTION
+     OUTPUT: IGNPOI=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION
+     MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR B.
+     TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
+     COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
+     SEPARATION OF CASES A AND B
+*/
+{
+//extern float sign( float num, float sign );
+static double a0 = -0.5;
+static double a1 = 0.3333333;
+static double a2 = -0.2500068;
+static double a3 = 0.2000118;
+static double a4 = -0.1661269;
+static double a5 = 0.1421878;
+static double a6 = -0.1384794;
+static double a7 = 0.125006;
+static double muold = -1.0E37;
+static double muprev = -1.0E37;
+static double fact[10] = {
+    1.0,1.0,2.0,6.0,24.0,120.0,720.0,5040.0,40320.0,362880.0
+};
+static int ignpoi,j,k,kflag,l,m,ll;     
+// JJV I added a variable 'll' here - it is the 'l' for CASE A
+static double b1,b2,c,c0,c1,c2,c3,d,del,difmuk,e,fk,fx,fy,g,omega,p,p0,px,py,q,s,
+    t,u,v,x,xx,pp[35];
+
+    if(*mu == muprev) goto S10;
+    if(*mu < 10.0) goto S120;
+/*
+     C A S E  A. (RECALCULATION OF S,D,L IF MU HAS CHANGED)
+*/
+//   JJV This is the case where I changed 'l' to 'll'
+//   JJV Here 'll' is set once and used in a comparison once
+    muprev = *mu;
+    s = sqrt(*mu);
+    d = 6.0* *mu * *mu;
+/*
+             THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
+             PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484)
+             IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
+*/
+    ll = (int) (*mu-1.1484);
+S10:
+/*
+     STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE
+*/
+    g = *mu+s*C2F(snorm)();
+    if(g < 0.0) goto S20;
+    ignpoi = (int) (g);
+/*
+     STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH
+*/
+    if(ignpoi >= ll) return ignpoi;
+/*
+     STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U
+*/
+    fk = (double)ignpoi;
+    difmuk = *mu-fk;
+    u = C2F(ranf)();
+    if(d*u >= difmuk*difmuk*difmuk) return ignpoi;
+S20:
+/*
+     STEP P. PREPARATIONS FOR STEPS Q AND H.
+             (RECALCULATIONS OF PARAMETERS IF NECESSARY)
+             .3989423=(2*PI)**(-.5)  .416667E-1=1./24.  .1428571=1./7.
+             THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
+             APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
+             C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
+*/
+    if(*mu == muold) goto S30;
+    muold = *mu;
+    omega = 0.3989423/s;
+    b1 = 4.166667E-2/ *mu;
+    b2 = 0.3*b1*b1;
+    c3 = 0.1428571*b1*b2;
+    c2 = b2-15.0*c3;
+    c1 = b1-6.0*b2+45.0*c3;
+    c0 = 1.0-b1+3.0*b2-15.0*c3;
+    c = 0.1069/ *mu;
+S30:
+    if(g < 0.0) goto S50;
+/*
+             'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
+*/
+    kflag = 0;
+    goto S70;
+S40:
+/*
+     STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
+*/
+    if(fy-u*fy <= py*exp(px-fx)) return ignpoi;
+S50:
+/*
+     STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL
+             DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
+             (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
+*/
+    e = C2F(sexpo)();
+    u = C2F(ranf)();
+    u += (u-1.0);
+    t = 1.8+fsign(e,u);
+    if(t <= -0.6744) goto S50;
+    ignpoi = (int) (*mu+s*t);
+    fk = (double)ignpoi;
+    difmuk = *mu-fk;
+/*
+             'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
+*/
+    kflag = 1;
+    goto S70;
+S60:
+/*
+     STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)
+*/
+    if(c*Abs(u) > py*exp(px+e)-fy*exp(fx+e)) goto S50;
+    return ignpoi;
+S70:
+/*
+     STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY.
+             CASE IGNPOI .LT. 10 USES FACTORIALS FROM TABLE FACT
+*/
+    if(ignpoi >= 10) goto S80;
+    px = - *mu;
+    py = pow((double) *mu,(double)ignpoi)/ *(fact+ignpoi);
+    goto S110;
+S80:
+/*
+             CASE IGNPOI .GE. 10 USES POLYNOMIAL APPROXIMATION
+             A0-A7 FOR ACCURACY WHEN ADVISABLE
+             .8333333E-1=1./12.  .3989423=(2*PI)**(-.5)
+*/
+    del = 8.333333E-2/fk;
+    del -= (4.8*del*del*del);
+    v = difmuk/fk;
+    if(Abs(v) <= 0.25) goto S90;
+    px = fk*log(1.0+v)-difmuk-del;
+    goto S100;
+S90:
+    px = fk*v*v*(((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0)-del;
+S100:
+    py = 0.3989423/sqrt(fk);
+S110:
+    x = (0.5-difmuk)/s;
+    xx = x*x;
+    fx = -0.5*xx;
+    fy = omega*(((c3*xx+c2)*xx+c1)*xx+c0);
+    if(kflag <= 0) goto S40;
+    goto S60;
+S120:
+/*
+     C A S E  B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY)
+*/
+//  JJV changed MUPREV assignment from 0.0 to initial value
+    muprev = -1.0E37;
+//     Jpc 1999: the next lines seams to produce a bug 
+//     and I finaly commented them out 
+//     IF (mu.EQ.muold) GO TO 130
+//     JJV added argument checker here
+//     JJV added line label here
+// 125  muold = mu
+//    if(mu == muold) goto S130;
+//    muold = mu;
+    m = Max(1L,(int) (*mu));
+    l = 0;
+    p = exp(- *mu);
+    q = p0 = p;
+S130:
+/*
+     STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
+*/
+    u = C2F(ranf)();
+    ignpoi = 0;
+    if(u <= p0) return ignpoi;
+/*
+     STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
+             PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
+             (0.458=PP(9) FOR MU=10)
+*/
+    if(l == 0) goto S150;
+    j = 1;
+    if(u > 0.458) j = Min(l,m);
+    for(k=j; k<=l; k++) {
+        if(u <= *(pp+k-1)) goto S180;
+    }
+    if(l == 35) goto S130;
+S150:
+/*
+     STEP C. CREATION OF NEW POISSON PROBABILITIES P
+             AND THEIR CUMULATIVES Q=PP(K)
+*/
+    l += 1;
+    for(k=l; k<=35; k++) {
+        p = p* *mu/(double)k;
+        q += p;
+        *(pp+k-1) = q;
+        if(u <= q) goto S170;
+    }
+    l = 35;
+    goto S130;
+S170:
+    l = k;
+S180:
+    ignpoi = k;
+    return ignpoi;
+}
diff --git a/scilab/modules/randlib/src/c/sexpo.c b/scilab/modules/randlib/src/c/sexpo.c
new file mode 100644 (file)
index 0000000..f556722
--- /dev/null
@@ -0,0 +1,86 @@
+#include "stack-c.h"
+#include <math.h>
+#include "grand.h"
+#include "machine.h" 
+#include "core_math.h"
+
+double C2F(sexpo)(void)
+/*
+**********************************************************************
+
+
+     (STANDARD-)  E X P O N E N T I A L   DISTRIBUTION
+
+
+**********************************************************************
+**********************************************************************
+
+This source code was taken in the project "freemat"(BSD license)
+This source code was modified by Gaüzère Sabine according to the 
+modifications done by JJV
+
+     FOR DETAILS SEE:
+
+               AHRENS, J.H. AND DIETER, U.
+               COMPUTER METHODS FOR SAMPLING FROM THE
+               EXPONENTIAL AND NORMAL DISTRIBUTIONS.
+               COMM. ACM, 15,10 (OCT. 1972), 873 - 882.
+
+     ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM
+     'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)
+
+     Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
+     SUNIF.  The argument IR thus goes away.
+
+**********************************************************************
+     Q(N) = SUM(ALOG(2.0)**K/K!)    K=1,..,N ,      THE HIGHEST N
+     (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION
+*/
+{
+/* ** OLD VALUES
+static double q[8] = {
+    0.69314718055995, 0.93337368751905, 0.98887779618387, 0.99849592529150,
+    0.99982928110614, 0.99998331641007, 0.99999856914388, 0.99999989069256
+//0.6931472,0.9333737,0.9888778,0.9984959,0.9998293,0.9999833,0.9999986,0.9999999
+};
+*/
+static double q[8] = {
+ 0.69314718246459960938,
+ 0.93337368965148925781,
+ 0.98887777328491210938,
+ 0.99849587678909301758,
+ 0.99982929229736328125,
+ 0.99998331069946289062,
+ 0.99999862909317016602,
+ 0.99999988079071044922
+};
+static int i;
+static double sexpo,a,u,ustar,umin;
+
+    a = 0.0;
+    u = C2F(ranf)();
+    goto S30;
+S20:
+    a += q[0];
+S30:
+    u += u;
+//  JJV changed the following to reflect the true algorithm and
+//  JJV prevent unpredictable behavior if U is initially 0.5.
+//  IF (u.LE.1.0) GO TO 20
+    if(u < 1.0) goto S20;
+    u -= 1.0;
+    if(u > q[0]) goto S60;
+    sexpo = a+u;
+    return sexpo;
+S60:
+    i = 1;
+    ustar = C2F(ranf)();
+    umin = ustar;
+S70:
+    ustar = C2F(ranf)();
+    if(ustar < umin) umin = ustar;
+    i += 1;
+    if(u > q[i-1]) goto S70;
+    sexpo = a+umin*q[0];
+    return sexpo;
+}
diff --git a/scilab/modules/randlib/src/c/sgamma.c b/scilab/modules/randlib/src/c/sgamma.c
new file mode 100644 (file)
index 0000000..25e1015
--- /dev/null
@@ -0,0 +1,250 @@
+#include "stack-c.h"
+#include <math.h>
+#include "grand.h"
+#include "machine.h" 
+#include "core_math.h"
+
+
+double fsign1 (double x, double y)
+{
+if (y >= 0.0)
+return Abs(x);
+else 
+return -Abs(x);
+}
+
+double C2F(sgamma)(double *a)
+/*
+**********************************************************************
+
+
+     (STANDARD-)  G A M M A  DISTRIBUTION
+
+
+**********************************************************************
+**********************************************************************
+
+               PARAMETER  A >= 1.0  !
+
+**********************************************************************
+This source code was taken in the project "freemat"(BSD license)
+This source code was modified by Gaüzère Sabine according to the 
+modifications done by JJV
+
+     FOR DETAILS SEE:
+
+               AHRENS, J.H. AND DIETER, U.
+               GENERATING GAMMA VARIATES BY A
+               MODIFIED REJECTION TECHNIQUE.
+               COMM. ACM, 25,1 (JAN. 1982), 47 - 54.
+
+     STEP NUMBERS CORRESPOND TO ALGORITHM 'GD' IN THE ABOVE PAPER
+                                 (STRAIGHTFORWARD IMPLEMENTATION)
+
+     Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of
+     SUNIF.  The argument IR thus goes away.
+
+**********************************************************************
+
+               PARAMETER  0.0 < A < 1.0  !
+
+**********************************************************************
+
+     FOR DETAILS SEE:
+
+               AHRENS, J.H. AND DIETER, U.
+               COMPUTER METHODS FOR SAMPLING FROM GAMMA,
+               BETA, POISSON AND BINOMIAL DISTRIBUTIONS.
+               COMPUTING, 12 (1974), 223 - 246.
+
+     (ADAPTED IMPLEMENTATION OF ALGORITHM 'GS' IN THE ABOVE PAPER)
+
+**********************************************************************
+     INPUT: A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION
+     OUTPUT: SGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION
+     COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K))
+     COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K)
+     COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K)
+     PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A"
+     SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380
+*/
+{
+//extern float sign( float num, float sign );
+static double q1 = 4.166669026017189026E-2;
+static double q2 = 2.083148062229156494E-2;
+static double q3 = 8.01191013306379318E-3;
+static double q4 = 1.44121004268527031E-3;
+static double q5 = -7.388000085484236E-5;
+static double q6 = 2.4510998628102243E-4;
+static double q7 = 2.4239999765995890E-4;
+static double a1 = 0.33333331346511840820;
+static double a2 = -0.25000301003456115723;
+static double a3 = 0.20000620186328887939;
+static double a4 = -0.16629210114479064941;
+static double a5 = 0.14236569404602050781;
+static double a6 = -0.13671770691871643066;
+static double a7 = 0.12337949872016906738;
+static double e1 = 1.0;
+static double e2 = 0.49998968839645385742;
+static double e3 = 0.16682900488376617432;
+static double e4 = 4.077529907226562500E-2;
+static double e5 = 1.029300037771463394E-2;
+static double aa = 0.0;
+static double aaa = 0.0;
+static double sqrt32 = 5.65685415267944335938;
+static double sgamma,s2,s,d,t,x,u,r,q0,b,si,c,v,q,e,w,p,b0;
+    if(*a == aa) goto S10;
+    if(*a < 1.0) goto S130;
+/*
+     STEP  1:  RECALCULATIONS OF S2,S,D IF A HAS CHANGED
+*/
+    aa = *a;
+    s2 = *a -0.5;
+    s = sqrt(s2);
+    d = sqrt32-12.0*s;
+S10:
+/*
+     STEP  2:  T=STANDARD NORMAL DEVIATE,
+               X=(S,1/2)-NORMAL DEVIATE.
+               IMMEDIATE ACCEPTANCE (I)
+*/
+    t = C2F(snorm)();
+    x = s+0.5*t;
+    sgamma = x*x;
+    if(t >= 0.0) return sgamma;
+/*
+     STEP  3:  U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S)
+*/
+    u = C2F(ranf)();
+    if(d*u <= t*t*t) return sgamma;
+/*
+     STEP  4:  RECALCULATIONS OF Q0,B,SI,C IF NECESSARY
+*/
+    if(*a == aaa) goto S40;
+    aaa = *a;
+    r = 1.0/ *a;
+    q0 = ((((((q7*r+q6)*r+q5)*r+q4)*r+q3)*r+q2)*r+q1)*r;
+/*
+               APPROXIMATION DEPENDING ON SIZE OF PARAMETER A
+               THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND
+               C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS
+*/
+    if(*a <= 3.686) goto S30; //3.68600010871887207031
+    if(*a <= 13.022) goto S20; //13.02200031280517578125
+/*
+               CASE 3:  A .GT. 13.022
+*/
+    b = 1.76999998092651367188;
+    si = 0.75;
+    c = 0.15150000154972076416/s;
+    goto S40;
+S20:
+/*
+               CASE 2:  3.686 .LT. A .LE. 13.022
+*/
+    b = 1.65400004386901855469+7.60000012814998627E-3*s2;
+    si = 1.67999994754791259766/s+0.27500000596046447754;
+    c = 6.199999898672103882E-2/s+2.400000020861625671E-2;
+    goto S40;
+S30:
+/*
+               CASE 1:  A .LE. 3.686
+*/
+    b = 0.46299999952316284180+s+0.17800000309944152832*s2;
+    si = 1.23500001430511474609;
+    c = 0.19499999284744262695/s-7.900000363588333130E-2+1.5999999642372131348E-1*s;
+S40:
+/*
+     STEP  5:  NO QUOTIENT TEST IF X NOT POSITIVE
+*/
+    if(x <= 0.0) goto S70;
+/*
+     STEP  6:  CALCULATION OF V AND QUOTIENT Q
+*/
+    v = t/(s+s);
+    if(Abs(v) <= 0.25) goto S50;
+    q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
+    goto S60;
+S50:
+    q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
+S60:
+/*
+     STEP  7:  QUOTIENT ACCEPTANCE (Q)
+*/
+    if(log(1.0-u) <= q) return sgamma;
+S70:
+/*
+     STEP  8:  E=STANDARD EXPONENTIAL DEVIATE
+               U= 0,1 -UNIFORM DEVIATE
+               T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE
+*/
+    e = C2F(sexpo)();
+    u = C2F(ranf)();
+    u += (u-1.0);
+    t = b+fsign1(si*e,u);
+/*
+     STEP  9:  REJECTION IF T .LT. TAU(1) = -.71874483771719
+*/
+    if(t < -0.7187449) goto S70; //-0.71874487400054931641
+/*
+     STEP 10:  CALCULATION OF V AND QUOTIENT Q
+*/
+    v = t/(s+s);
+    if(Abs(v) <= 0.25) goto S90;
+    q = q0-s*t+0.25*t*t+(s2+s2)*log(1.0+v);
+    goto S100;
+S90:
+    q = q0+0.5*t*t*((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v;
+S100:
+/*
+     STEP 11:  HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8)
+*/
+    if(q <= 0.0) goto S70;
+    if(q <= 0.5) goto S110;
+//  JJV modified the code through line 125 to handle large Q case
+    if(q < 15.0) goto S105;
+//  JJV Here Q is large enough that Q = log(exp(Q) - 1.0) (for DOUBLE PRECISION Q)
+//  JJV so reformulate test at 120 in terms of one EXP, if not too big
+//  JJV 87.49823 is close to the largest DOUBLE PRECISION which can be
+//  JJV exponentiated (87.49823 = log(1.0E38))
+    if((q+e-0.5*t*t) > 87.49823) goto S125;//87.498222998046875
+    if(c*Abs(u) > exp(q+e-0.5*t*t)) goto S70;
+    goto S125;
+S105:
+    w = exp(q)-1.0;
+    goto S120;
+S110:
+    w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*q;
+S120:
+/*
+               IF T IS REJECTED, SAMPLE AGAIN AT STEP 8
+*/
+    if(c*Abs(u) > w*exp(e-0.5*t*t)) goto S70;
+S125:
+    x = s+0.5*t;
+    sgamma = x*x;
+    return sgamma;
+/*
+     ALTERNATE METHOD FOR PARAMETERS A BELOW 1  (.3678794=EXP(-1.))
+*/
+//     JJV changed B to B0 (which was added to declarations for this)
+//     JJV in 130 to END to fix rare and subtle bug.
+//     JJV Line: '130 aa = 0.0' was removed (unnecessary, wasteful).
+//     JJV Reasons: the state of AA only serves to tell the A .GE. 1.0
+//     JJV case if certain A-dependant constants need to be recalculated.
+//     JJV The A .LT. 1.0 case (here) no longer changes any of these, and
+//     JJV the recalculation of B (which used to change with an
+//     JJV A .LT. 1.0 call) is governed by the state of AAA anyway.
+S130:
+    b0 = 1.0+0.3678794 * (*a); //0.36787939071655273438
+S140:
+    p = b0*C2F(ranf)();
+    if(p >= 1.0) goto S150;
+    sgamma = exp(log(p)/ *a);
+    if(C2F(sexpo)() < sgamma) goto S140;
+    return sgamma;
+S150:
+    sgamma = -log((b0-p)/ *a);
+    if(C2F(sexpo)() < (1.0- *a)*log(sgamma)) goto S140;
+    return sgamma;
+}
diff --git a/scilab/modules/randlib/src/fortran/genbet.f b/scilab/modules/randlib/src/fortran/genbet.f
deleted file mode 100644 (file)
index 646a3e8..0000000
+++ /dev/null
@@ -1,245 +0,0 @@
-      DOUBLE PRECISION FUNCTION genbet(aa,bb)
-C**********************************************************************
-C
-C     DOUBLE PRECISION FUNCTION GENBET( A, B )
-C               GeNerate BETa random deviate
-C
-C
-C                              Function
-C
-C
-C     Returns a single random deviate from the beta distribution with
-C     parameters A and B.  The density of the beta is
-C               x^(a-1) * (1-x)^(b-1) / B(a,b) for 0 < x < 1
-C
-C
-C                              Arguments
-C
-C
-C     A --> First parameter of the beta distribution
-C                         DOUBLE PRECISION A
-C     JJV                 (A > 1.0E-37)
-C
-C     B --> Second parameter of the beta distribution
-C                         DOUBLE PRECISION B
-C     JJV                 (B > 1.0E-37)
-C
-C
-C                              Method
-C
-C
-C     R. C. H. Cheng
-C     Generating Beta Variates with Nonintegral Shape Parameters
-C     Communications of the ACM, 21:317-322  (1978)
-C     (Algorithms BB and BC)
-C
-C**********************************************************************
-C     .. Parameters ..
-C     Close to the largest number that can be exponentiated
-      DOUBLE PRECISION expmax
-C     JJV changed this - 89 was too high, and LOG(1.0E38) = 87.49823
-      PARAMETER (expmax=87.49823)
-C     Close to the largest representable single precision number
-      DOUBLE PRECISION infnty
-      PARAMETER (infnty=1.0E38)
-C     JJV added the parameter minlog
-C     Close to the smallest number of which a LOG can be taken.
-      DOUBLE PRECISION minlog
-      PARAMETER (minlog=1.0E-37)
-C     ..
-C     .. Scalar Arguments ..
-      DOUBLE PRECISION aa,bb
-C     ..
-C     .. Local Scalars ..
-      DOUBLE PRECISION a,alpha,b,beta,delta,gamma,k1,k2,olda,oldb,r,s,t,
-     +     u1,u2,v,w,y,z
-      LOGICAL qsame
-C     ..
-C     .. External Functions ..
-      DOUBLE PRECISION ranf
-      EXTERNAL ranf
-C     ..
-C     .. Intrinsic Functions ..
-      INTRINSIC exp,log,max,min,sqrt
-C     ..
-C     .. Save statement ..
-C     JJV added a,b
-      SAVE olda,oldb,alpha,beta,gamma,k1,k2,a,b
-C     ..
-C     .. Data statements ..
-C     JJV changed these to ridiculous values
-      DATA olda,oldb/-1.0E37,-1.0E37/
-C     ..
-C     .. Executable Statements ..
-      qsame = (olda.EQ.aa) .AND. (oldb.EQ.bb)
-      IF (qsame) GO TO 20
-C     JJV added small minimum for small log problem in calc of W
-C     in Rand.c
-   10 olda = aa
-      oldb = bb
-   20 IF (.NOT. (min(aa,bb).GT.1.0)) GO TO 100
-
-
-C     Alborithm BB
-
-C
-C     Initialize
-C
-      IF (qsame) GO TO 30
-      a = min(aa,bb)
-      b = max(aa,bb)
-      alpha = a + b
-      beta = sqrt((alpha-2.0)/ (2.0*a*b-alpha))
-      gamma = a + 1.0/beta
-   30 CONTINUE
-   40 u1 = ranf()
-C
-C     Step 1
-C
-      u2 = ranf()
-      v = beta*log(u1/ (1.0-u1))
-C     JJV altered this
-      IF (v.GT.expmax) GO TO 55
-C     JJV added checker to see if a*exp(v) will overflow
-C     JJV 50 _was_ w = a*exp(v); also note here a > 1.0
-   50 w = exp(v)
-      IF (w.GT.infnty/a) GO TO 55
-      w = a*w
-      GO TO 60
- 55   w = infnty
-
-   60 z = u1**2*u2
-      r = gamma*v - 1.3862944
-      s = a + r - w
-C
-C     Step 2
-C
-      IF ((s+2.609438).GE. (5.0*z)) GO TO 70
-C
-C     Step 3
-C
-      t = log(z)
-      IF (s.GT.t) GO TO 70
-C
-C     Step 4
-C
-C     JJV added checker to see if log(alpha/(b+w)) will 
-C     JJV overflow.  If so, we count the log as -INF, and
-C     JJV consequently evaluate conditional as true, i.e.
-C     JJV the algorithm rejects the trial and starts over
-C     JJV May not need this here since ALPHA > 2.0
-      IF (alpha/(b+w).LT.minlog) GO TO 40
-
-      IF ((r+alpha*log(alpha/ (b+w))).LT.t) GO TO 40
-C
-C     Step 5
-C
-   70 IF (.NOT. (aa.EQ.a)) GO TO 80
-      genbet = w/ (b+w)
-      GO TO 90
-
-   80 genbet = b/ (b+w)
-   90 GO TO 230
-
-
-C     Algorithm BC
-
-C
-C     Initialize
-C
-  100 IF (qsame) GO TO 110
-      a = max(aa,bb)
-      b = min(aa,bb)
-      alpha = a + b
-      beta = 1.0/b
-      delta = 1.0 + a - b
-      k1 = delta* (0.0138889+0.0416667*b)/ (a*beta-0.777778)
-      k2 = 0.25 + (0.5+0.25/delta)*b
-  110 CONTINUE
-  120 u1 = ranf()
-C
-C     Step 1
-C
-      u2 = ranf()
-      IF (u1.GE.0.5) GO TO 130
-C
-C     Step 2
-C
-      y = u1*u2
-      z = u1*y
-      IF ((0.25*u2+z-y).GE.k1) GO TO 120
-      GO TO 170
-C
-C     Step 3
-C
-  130 z = u1**2*u2
-      IF (.NOT. (z.LE.0.25)) GO TO 160
-      v = beta*log(u1/ (1.0-u1))
-
-C     JJV instead of checking v > expmax at top, I will check
-C     JJV if a < 1, then check the appropriate values
-
-      IF (a.GT.1.0) GO TO 135
-C     JJV A < 1 so it can help out if EXP(V) would overflow
-      IF (v.GT.expmax) GO TO 132
-      w = a*exp(v)
-      GO TO 200
- 132  w = v + log(a)
-      IF (w.GT.expmax) GO TO 140
-      w = exp(w)
-      GO TO 200
-
-C     JJV in this case A > 1
- 135  IF (v.GT.expmax) GO TO 140
-      w = exp(v)
-      IF (w.GT.infnty/a) GO TO 140
-      w = a*w
-      GO TO 200
- 140  w = infnty
-      GO TO 200
-
-  160 IF (z.GE.k2) GO TO 120
-C
-C     Step 4
-C
-C
-C     Step 5
-C
-  170 v = beta*log(u1/ (1.0-u1))
-
-C     JJV same kind of checking as above
-      IF (a.GT.1.0) GO TO 175
-C     JJV A < 1 so it can help out if EXP(V) would overflow
-      IF (v.GT.expmax) GO TO 172
-      w = a*exp(v)
-      GO TO 190
- 172  w = v + log(a)
-      IF (w.GT.expmax) GO TO 180
-      w = exp(w)
-      GO TO 190
-
-C     JJV in this case A > 1
- 175  IF (v.GT.expmax) GO TO 180
-      w = exp(v)
-      IF (w.GT.infnty/a) GO TO 180
-      w = a*w
-      GO TO 190
-
-  180 w = infnty
-
-C     JJV here we also check to see if log overlows; if so, we treat it
-C     JJV as -INF, which means condition is true, i.e. restart
-  190 IF (alpha/(b+w).LT.minlog) GO TO 120
-      IF ((alpha* (log(alpha/ (b+w))+v)-1.3862944).LT.log(z)) GO TO 120
-C
-C     Step 6
-C
-  200 IF (.NOT. (a.EQ.aa)) GO TO 210
-      genbet = w/ (b+w)
-      GO TO 220
-
-  210 genbet = b/ (b+w)
-  220 CONTINUE
-  230 RETURN
-
-      END
diff --git a/scilab/modules/randlib/src/fortran/ignbin.f b/scilab/modules/randlib/src/fortran/ignbin.f
deleted file mode 100644 (file)
index 659de4f..0000000
+++ /dev/null
@@ -1,339 +0,0 @@
-      INTEGER FUNCTION ignbin(n,pp)
-C**********************************************************************
-C
-C     INTEGER FUNCTION IGNBIN( N, PP )
-C
-C                    GENerate BINomial random deviate
-C
-C
-C                              Function
-C
-C
-C     Generates a single random deviate from a binomial
-C     distribution whose number of trials is N and whose
-C     probability of an event in each trial is P.
-C
-C
-C                              Arguments
-C
-C
-C     N  --> The number of trials in the binomial distribution
-C            from which a random deviate is to be generated.
-C                              INTEGER N
-C     JJV                      (N >= 0)
-C
-C     PP --> The probability of an event in each trial of the
-C            binomial distribution from which a random deviate
-C            is to be generated.
-C                              DOUBLE PRECISION PP
-C     JJV                      (0.0 <= pp <= 1.0)
-C
-C     IGNBIN <-- A random deviate yielding the number of events
-C                from N independent trials, each of which has
-C                a probability of event P.
-C                              INTEGER IGNBIN
-C
-C
-C                              Note
-C
-C
-C     Uses RANF so the value of the seeds, ISEED1 and ISEED2 must be set
-C     by a call similar to the following
-C          DUM = RANSET( ISEED1, ISEED2 )
-C
-C
-C                              Method
-C
-C
-C     This is algorithm BTPE from:
-C
-C         Kachitvichyanukul, V. and Schmeiser, B. W.
-C
-C         Binomial Random Variate Generation.
-C         Communications of the ACM, 31, 2
-C         (February, 1988) 216.
-C
-C**********************************************************************
-C     SUBROUTINE BTPEC(N,PP,ISEED,JX)
-C
-C     BINOMIAL RANDOM VARIATE GENERATOR
-C     MEAN .LT. 30 -- INVERSE CDF
-C       MEAN .GE. 30 -- ALGORITHM BTPE:  ACCEPTANCE-REJECTION VIA
-C       FOUR REGION COMPOSITION.  THE FOUR REGIONS ARE A TRIANGLE
-C       (SYMMETRIC IN THE CENTER), A PAIR OF PARALLELOGRAMS (ABOVE
-C       THE TRIANGLE), AND EXPONENTIAL LEFT AND RIGHT TAILS.
-C
-C     BTPE REFERS TO BINOMIAL-TRIANGLE-PARALLELOGRAM-EXPONENTIAL.
-C     BTPEC REFERS TO BTPE AND "COMBINED."  THUS BTPE IS THE
-C       RESEARCH AND BTPEC IS THE IMPLEMENTATION OF A COMPLETE
-C       USABLE ALGORITHM.
-C     REFERENCE:  VORATAS KACHITVICHYANUKUL AND BRUCE SCHMEISER,
-C       "BINOMIAL RANDOM VARIATE GENERATION,"
-C       COMMUNICATIONS OF THE ACM, FORTHCOMING
-C     WRITTEN:  SEPTEMBER 1980.
-C       LAST REVISED:  MAY 1985, JULY 1987
-C     REQUIRED SUBPROGRAM:  RAND() -- A UNIFORM (0,1) RANDOM NUMBER
-C                           GENERATOR
-C     ARGUMENTS
-C
-C       N : NUMBER OF BERNOULLI TRIALS            (INPUT)
-C       PP : PROBABILITY OF SUCCESS IN EACH TRIAL (INPUT)
-C       ISEED:  RANDOM NUMBER SEED                (INPUT AND OUTPUT)
-C       JX:  RANDOMLY GENERATED OBSERVATION       (OUTPUT)
-C
-C     VARIABLES
-C       PSAVE: VALUE OF PP FROM THE LAST CALL TO BTPEC
-C       NSAVE: VALUE OF N FROM THE LAST CALL TO BTPEC
-C       XNP:  VALUE OF THE MEAN FROM THE LAST CALL TO BTPEC
-C
-C       P: PROBABILITY USED IN THE GENERATION PHASE OF BTPEC
-C       FFM: TEMPORARY VARIABLE EQUAL TO XNP + P
-C       M:  INTEGER VALUE OF THE CURRENT MODE
-C       FM:  FLOATING POINT VALUE OF THE CURRENT MODE
-C       XNPQ: TEMPORARY VARIABLE USED IN SETUP AND SQUEEZING STEPS
-C       P1:  AREA OF THE TRIANGLE
-C       C:  HEIGHT OF THE PARALLELOGRAMS
-C       XM:  CENTER OF THE TRIANGLE
-C       XL:  LEFT END OF THE TRIANGLE
-C       XR:  RIGHT END OF THE TRIANGLE
-C       AL:  TEMPORARY VARIABLE
-C       XLL:  RATE FOR THE LEFT EXPONENTIAL TAIL
-C       XLR:  RATE FOR THE RIGHT EXPONENTIAL TAIL
-C       P2:  AREA OF THE PARALLELOGRAMS
-C       P3:  AREA OF THE LEFT EXPONENTIAL TAIL
-C       P4:  AREA OF THE RIGHT EXPONENTIAL TAIL
-C       U:  A U(0,P4) RANDOM VARIATE USED FIRST TO SELECT ONE OF THE
-C           FOUR REGIONS AND THEN CONDITIONALLY TO GENERATE A VALUE
-C           FROM THE REGION
-C       V:  A U(0,1) RANDOM NUMBER USED TO GENERATE THE RANDOM VALUE
-C           (REGION 1) OR TRANSFORMED INTO THE VARIATE TO ACCEPT OR
-C           REJECT THE CANDIDATE VALUE
-C       IX:  INTEGER CANDIDATE VALUE
-C       X:  PRELIMINARY CONTINUOUS CANDIDATE VALUE IN REGION 2 LOGIC
-C           AND A FLOATING POINT IX IN THE ACCEPT/REJECT LOGIC
-C       K:  ABSOLUTE VALUE OF (IX-M)
-C       F:  THE HEIGHT OF THE SCALED DENSITY FUNCTION USED IN THE
-C           ACCEPT/REJECT DECISION WHEN BOTH M AND IX ARE SMALL
-C           ALSO USED IN THE INVERSE TRANSFORMATION
-C       R: THE RATIO P/Q
-C       G: CONSTANT USED IN CALCULATION OF PROBABILITY
-C       MP:  MODE PLUS ONE, THE LOWER INDEX FOR EXPLICIT CALCULATION
-C            OF F WHEN IX IS GREATER THAN M
-C       IX1:  CANDIDATE VALUE PLUS ONE, THE LOWER INDEX FOR EXPLICIT
-C             CALCULATION OF F WHEN IX IS LESS THAN M
-C       I:  INDEX FOR EXPLICIT CALCULATION OF F FOR BTPE
-C       AMAXP: MAXIMUM ERROR OF THE LOGARITHM OF NORMAL BOUND
-C       YNORM: LOGARITHM OF NORMAL BOUND
-C       ALV:  NATURAL LOGARITHM OF THE ACCEPT/REJECT VARIATE V
-C
-C       X1,F1,Z,W,Z2,X2,F2, AND W2 ARE TEMPORARY VARIABLES TO BE
-C       USED IN THE FINAL ACCEPT/REJECT TEST
-C
-C       QN: PROBABILITY OF NO SUCCESS IN N TRIALS
-C
-C     REMARK
-C       IX AND JX COULD LOGICALLY BE THE SAME VARIABLE, WHICH WOULD
-C       SAVE A MEMORY POSITION AND A LINE OF CODE.  HOWEVER, SOME
-C       COMPILERS (E.G.,CDC MNF) OPTIMIZE BETTER WHEN THE ARGUMENTS
-C       ARE NOT INVOLVED.
-C
-C     ISEED NEEDS TO BE DOUBLE PRECISION IF THE IMSL ROUTINE
-C     GGUBFS IS USED TO GENERATE UNIFORM RANDOM NUMBER, OTHERWISE
-C     TYPE OF ISEED SHOULD BE DICTATED BY THE UNIFORM GENERATOR
-C
-C**********************************************************************
-
-C
-C
-C
-C*****DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY
-C
-C     ..
-C     .. Scalar Arguments ..
-      DOUBLE PRECISION pp
-      INTEGER n
-C     ..
-C     .. Local Scalars ..
-      DOUBLE PRECISION al,alv,amaxp,c,f,f1,f2,ffm,fm,g,p,p1,p2,p3,p4,
-     +     psave,q,qn,r,u,
-     +     v,w,w2,x,x1,x2,xl,xll,xlr,xm,xnp,xnpq,xr,ynorm,z,z2
-      INTEGER i,ix,ix1,k,m,mp,nsave
-C     ..
-C     .. External Functions ..
-      DOUBLE PRECISION ranf
-      EXTERNAL ranf
-C     ..
-C     .. Intrinsic Functions ..
-      INTRINSIC abs,log,dmin1,iabs,int,sqrt
-C     JJV ..
-C     JJV .. Save statement ..
-      SAVE p,q,m,fm,xnp,xnpq,p1,xm,xl,xr,c,xll,xlr,p2,p3,p4,qn,r,g,
-     +     psave,nsave
-C     JJV I am including the variables in data statements
-C     ..
-C     .. Data statements ..
-C     JJV made these ridiculous starting values - the hope is that
-C     JJV no one will call this the first time with them as args
-      DATA psave,nsave/-1.0E37,-214748365/
-C     ..
-C     .. Executable Statements ..
-      IF (pp.NE.psave) GO TO 10
-      IF (n.NE.nsave) GO TO 20
-      IF ((xnp-30.) .lt. 0) then 
-         goto 150
-      else
-         goto 30
-      endif
-C
-C*****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE
-C
-
-C     JJV added the argument checker - involved only renaming 10
-C     JJV and 20 to the checkers and adding checkers
-C     JJV Only remaining problem - if called initially with the
-C     JJV initial values of psave and nsave, it will hang
- 10   continue
-      psave = pp
-      p = dmin1(psave,1.-psave)
-      q = 1. - p
- 20   continue
-      xnp = n*p
-      nsave = n
-      IF (xnp.LT.30.) GO TO 140
-      ffm = xnp + p
-      m = ffm
-      fm = m
-      xnpq = xnp*q
-      p1 = int(2.195*sqrt(xnpq)-4.6*q) + 0.5
-      xm = fm + 0.5
-      xl = xm - p1
-      xr = xm + p1
-      c = 0.134 + 20.5/ (15.3+fm)
-      al = (ffm-xl)/ (ffm-xl*p)
-      xll = al* (1.+.5*al)
-      al = (xr-ffm)/ (xr*q)
-      xlr = al* (1.+.5*al)
-      p2 = p1* (1.+c+c)
-      p3 = p2 + c/xll
-      p4 = p3 + c/xlr
-C      WRITE(6,100) N,P,P1,P2,P3,P4,XL,XR,XM,FM
-C  100 FORMAT(I15,4F18.7/5F18.7)
-C
-C*****GENERATE VARIATE
-C
-   30 u = ranf()*p4
-      v = ranf()
-C
-C     TRIANGULAR REGION
-C
-      IF (u.GT.p1) GO TO 40
-      ix = xm - p1*v + u
-      GO TO 170
-C
-C     PARALLELOGRAM REGION
-C
-   40 IF (u.GT.p2) GO TO 50
-      x = xl + (u-p1)/c
-      v = v*c + 1. - abs(xm-x)/p1
-      IF (v.GT.1. .OR. v.LE.0.) GO TO 30
-      ix = x
-      GO TO 70
-C
-C     LEFT TAIL
-C
-   50 IF (u.GT.p3) GO TO 60
-      ix = xl + log(v)/xll
-      IF (ix.LT.0) GO TO 30
-      v = v* (u-p2)*xll
-      GO TO 70
-C
-C     RIGHT TAIL
-C
-   60 ix = xr - log(v)/xlr
-      IF (ix.GT.n) GO TO 30
-      v = v* (u-p3)*xlr
-C
-C*****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST
-C
-   70 k = iabs(ix-m)
-      IF (k.GT.20 .AND. k.LT.xnpq/2-1) GO TO 130
-C
-C     EXPLICIT EVALUATION
-C
-      f = 1.0
-      r = p/q
-      g = (n+1)*r
-      CRES=m-ix
-      IF (CRES .lt. 0) then 
-         goto 80
-      elseif (CRES .eq. 0) then 
-         goto 120
-      else
-         goto 100
-      endif
-   80 mp = m + 1
-      DO 90 i = mp,ix
-          f = f* (g/i-r)
-   90 CONTINUE
-      GO TO 120
-
-  100 ix1 = ix + 1
-      DO 110 i = ix1,m
-          f = f/ (g/i-r)
-  110 CONTINUE
-  120 IF ((v-f) .le. 0) then 
-         goto 170
-      else
-         goto 30
-      endif
-C
-C     SQUEEZING USING UPPER AND LOWER BOUNDS ON LOG(F(X))
-C
-  130 amaxp = (k/xnpq)* ((k* (k/3.+.625)+.1666666666666)/xnpq+.5)
-      ynorm = -k*k/ (2.*xnpq)
-      alv = log(v)
-      IF (alv.LT.ynorm-amaxp) GO TO 170
-      IF (alv.GT.ynorm+amaxp) GO TO 30
-C
-C     STIRLING'S FORMULA TO MACHINE ACCURACY FOR
-C     THE FINAL ACCEPTANCE/REJECTION TEST
-C
-      x1 = ix + 1
-      f1 = fm + 1.
-      z = n + 1 - fm
-      w = n - ix + 1.
-      z2 = z*z
-      x2 = x1*x1
-      f2 = f1*f1
-      w2 = w*w
-      IF ((alv- (xm*log(f1/x1)+ (n-m+.5)*log(z/w)+ (ix-
-     +    m)*log(w*p/ (x1*q))+ (13860.- (462.- (132.- (99.-
-     +    140./f2)/f2)/f2)/f2)/f1/166320.+ (13860.- (462.- (132.- (99.-
-     +    140./z2)/z2)/z2)/z2)/z/166320.+ (13860.- (462.- (132.- (99.-
-     +    140./x2)/x2)/x2)/x2)/x1/166320.+ (13860.- (462.- (132.- (99.-
-     +    140./w2)/w2)/w2)/w2)/w/166320.)) .le. 0) then
-         goto 170
-      else
-         goto 30
-      endif
-C
-C     INVERSE CDF LOGIC FOR MEAN LESS THAN 30
-C
-  140 qn = q**n
-      r = p/q
-      g = r* (n+1)
-  150 ix = 0
-      f = qn
-      u = ranf()
-  160 IF (u.LT.f) GO TO 170
-      IF (ix.GT.110) GO TO 150
-      u = u - f
-      ix = ix + 1
-      f = f* (g/ix-r)
-      GO TO 160
-
-  170 IF (psave.GT.0.5) ix = n - ix
-      ignbin = ix
-      RETURN
-
-      END
diff --git a/scilab/modules/randlib/src/fortran/ignpoi.f b/scilab/modules/randlib/src/fortran/ignpoi.f
deleted file mode 100644 (file)
index 9ab2f08..0000000
+++ /dev/null
@@ -1,292 +0,0 @@
-      INTEGER FUNCTION ignpoi(mu)
-C**********************************************************************
-C
-C     INTEGER FUNCTION IGNPOI( MU )
-C
-C                    GENerate POIsson random deviate
-C
-C
-C                              Function
-C
-C
-C     Generates a single random deviate from a Poisson
-C     distribution with mean MU.
-C
-C
-C                              Arguments
-C
-C
-C     MU --> The mean of the Poisson distribution from which
-C            a random deviate is to be generated.
-C                              DOUBLE PRECISION MU
-C     JJV                    (MU >= 0.0)
-C
-C     IGNPOI <-- The random deviate.
-C                              INTEGER IGNPOI (non-negative)
-C
-C
-C                              Method
-C
-C
-C     Renames KPOIS from TOMS as slightly modified by BWB to use RANF
-C     instead of SUNIF.
-C
-C     For details see:
-C
-C               Ahrens, J.H. and Dieter, U.
-C               Computer Generation of Poisson Deviates
-C               From Modified Normal Distributions.
-C               ACM Trans. Math. Software, 8, 2
-C               (June 1982),163-179
-C
-C**********************************************************************
-C**********************************************************************C
-C**********************************************************************C
-C                                                                      C
-C                                                                      C
-C     P O I S S O N  DISTRIBUTION                                      C
-C                                                                      C
-C                                                                      C
-C**********************************************************************C
-C**********************************************************************C
-C                                                                      C
-C     FOR DETAILS SEE:                                                 C
-C                                                                      C
-C               AHRENS, J.H. AND DIETER, U.                            C
-C               COMPUTER GENERATION OF POISSON DEVIATES                C
-C               FROM MODIFIED NORMAL DISTRIBUTIONS.                    C
-C               ACM TRANS. MATH. SOFTWARE, 8,2 (JUNE 1982), 163 - 179. C
-C                                                                      C
-C     (SLIGHTLY MODIFIED VERSION OF THE PROGRAM IN THE ABOVE ARTICLE)  C
-C                                                                      C
-C**********************************************************************C
-C
-C      INTEGER FUNCTION IGNPOI(IR,MU)
-C
-C     INPUT:  IR=CURRENT STATE OF BASIC RANDOM NUMBER GENERATOR
-C             MU=MEAN MU OF THE POISSON DISTRIBUTION
-C     OUTPUT: IGNPOI=SAMPLE FROM THE POISSON-(MU)-DISTRIBUTION
-C
-C
-C
-C     MUPREV=PREVIOUS MU, MUOLD=MU AT LAST EXECUTION OF STEP P OR CASE B
-C     TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
-C     COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
-C
-C
-C
-C     SEPARATION OF CASES A AND B
-C
-C     .. Scalar Arguments ..
-      DOUBLE PRECISION mu
-C     ..
-C     .. Local Scalars ..
-      DOUBLE PRECISION a0,a1,a2,a3,a4,a5,a6,a7,b1,b2,c,c0,c1,c2,c3,
-     $     d,del,difmuk,e,
-     $     fk,fx,fy,g,muold,muprev,omega,p,p0,px,py,q,s,t,u,v,x,xx
-C     JJV I added a variable 'll' here - it is the 'l' for CASE A
-      INTEGER j,k,kflag,l,ll,m
-C     ..
-C     .. Local Arrays ..
-      DOUBLE PRECISION fact(10),pp(35)
-C     ..
-C     .. External Functions ..
-      DOUBLE PRECISION ranf,sexpo,snorm
-      EXTERNAL ranf,sexpo,snorm
-C     ..
-C     .. Intrinsic Functions ..
-      INTRINSIC abs,log,exp,float,idint,max0,min0,dsign,sqrt
-C     ..
-C     JJV added this for case: mu unchanged
-C     .. Save statement ..
-      SAVE s, d, l, ll, omega, c3, c2, c1, c0, c, m, p, q, p0,
-     +     a0, a1, a2, a3, a4, a5, a6, a7, fact, muprev, muold
-C     ..
-C     JJV end addition - I am including vars in Data statements
-C     .. Data statements ..
-C     JJV changed initial values of MUPREV and MUOLD to -1.0E37
-C     JJV if no one calls IGNPOI with MU = -1.0E37 the first time,
-C     JJV the code shouldn't break
-      DATA muprev,muold/-1.0E37,-1.0E37/
-      DATA a0,a1,a2,a3,a4,a5,a6,a7/-.5,.3333333,-.2500068,.2000118,
-     +     -.1661269,.1421878,-.1384794,.1250060/
-      DATA fact/1.,1.,2.,6.,24.,120.,720.,5040.,40320.,362880./
-C     ..
-C     .. Executable Statements ..
-
-      IF (mu.EQ.muprev) GO TO 10
-      IF (mu.LT.10.0) GO TO 120
-C
-C     C A S E  A. (RECALCULATION OF S,D,LL IF MU HAS CHANGED)
-C
-C     JJV This is the case where I changed 'l' to 'll'
-C     JJV Here 'll' is set once and used in a comparison once
-
-      muprev = mu
-      s = sqrt(mu)
-      d = 6.0*mu*mu
-C
-C             THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
-C             PROBABILITIES FK WHENEVER K >= M(MU). LL=IDINT(MU-1.1484)
-C             IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
-C
-      ll = idint(mu-1.1484)
-C
-C     STEP N. NORMAL SAMPLE - SNORM(IR) FOR STANDARD NORMAL DEVIATE
-C
-   10 g = mu + s*snorm()
-      IF (g.LT.0.0) GO TO 20
-      ignpoi = idint(g)
-C
-C     STEP I. IMMEDIATE ACCEPTANCE IF IGNPOI IS LARGE ENOUGH
-C
-      IF (ignpoi.GE.ll) RETURN
-C
-C     STEP S. SQUEEZE ACCEPTANCE - SUNIF(IR) FOR (0,1)-SAMPLE U
-C
-      fk = dble(ignpoi)
-      difmuk = mu - fk
-      u = ranf()
-      IF (d*u.GE.difmuk*difmuk*difmuk) RETURN
-C
-C     STEP P. PREPARATIONS FOR STEPS Q AND H.
-C             (RECALCULATIONS OF PARAMETERS IF NECESSARY)
-C             .3989423=(2*PI)**(-.5)  .416667E-1=1./24.  .1428571=1./7.
-C             THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
-C             APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
-C             C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
-C
-   20 IF (mu.EQ.muold) GO TO 30
-      muold = mu
-      omega = .3989423/s
-      b1 = .4166667E-1/mu
-      b2 = .3*b1*b1
-      c3 = .1428571*b1*b2
-      c2 = b2 - 15.*c3
-      c1 = b1 - 6.*b2 + 45.*c3
-      c0 = 1. - b1 + 3.*b2 - 15.*c3
-      c = .1069/mu
-   30 IF (g.LT.0.0) GO TO 50
-C
-C             'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
-C
-      kflag = 0
-      GO TO 70
-C
-C     STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
-C
-   40 IF (fy-u*fy.LE.py*exp(px-fx)) RETURN
-C
-C     STEP E. EXPONENTIAL SAMPLE - SEXPO(IR) FOR STANDARD EXPONENTIAL
-C             DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
-C             (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
-C
-   50 e = sexpo()
-      u = ranf()
-      u = u + u - 1.0
-      t = 1.8d0 + dsign(e,u)
-      IF (t.LE. (-.6744)) GO TO 50
-      ignpoi = idint(mu+s*t)
-      fk = dble(ignpoi)
-      difmuk = mu - fk
-C
-C             'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
-C
-      kflag = 1
-      GO TO 70
-C
-C     STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)
-C
-   60 IF (c*abs(u).GT.py*exp(px+e)-fy*exp(fx+e)) GO TO 50
-      RETURN
-C
-C     STEP F. 'SUBROUTINE' F. CALCULATION OF PX,PY,FX,FY.
-C             CASE IGNPOI .LT. 10 USES FACTORIALS FROM TABLE FACT
-C
-   70 IF (ignpoi.GE.10) GO TO 80
-      px = -mu
-      py = mu**ignpoi/fact(ignpoi+1)
-      GO TO 110
-C
-C             CASE IGNPOI .GE. 10 USES POLYNOMIAL APPROXIMATION
-C             A0-A7 FOR ACCURACY WHEN ADVISABLE
-C             .8333333E-1=1./12.  .3989423=(2*PI)**(-.5)
-C
-   80 del = .8333333E-1/fk
-      del = del - 4.8*del*del*del
-      v = difmuk/fk
-      IF (abs(v).LE.0.25) GO TO 90
-      px = fk*log(1.0+v) - difmuk - del
-      GO TO 100
-
-   90 px = fk*v*v* (((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0) -
-     +     del
-  100 py = .3989423/sqrt(fk)
-  110 x = (0.5-difmuk)/s
-      xx = x*x
-      fx = -0.5*xx
-      fy = omega* (((c3*xx+c2)*xx+c1)*xx+c0)
-      IF (kflag .le. 0) then
-         goto 40
-      else
-         goto 60
-      endif
-      
-C
-C     C A S E  B. (START NEW TABLE AND CALCULATE P0 IF NECESSARY)
-C
-C     JJV changed MUPREV assignment from 0.0 to initial value
-  120 muprev = -1.0E37
-C     Jpc 1999: the next lines seams to produce a bug 
-C     and I finaly commented them out 
-C     IF (mu.EQ.muold) GO TO 130
-C     JJV added argument checker here
-C     JJV added line label here
-C 125  muold = mu
-      m = max0(1,idint(mu))
-      l = 0
-      p = exp(-mu)
-      q = p
-      p0 = p
-C
-C     STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
-C
-  130 u = ranf()
-      ignpoi = 0
-      IF (u.LE.p0) RETURN
-C
-C     STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
-C             PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
-C             (0.458=PP(9) FOR MU=10)
-C
-      IF (l.EQ.0) GO TO 150
-      j = 1
-      IF (u.GT.0.458) j = min0(l,m)
-      DO 140 k = j,l
-          IF (u.LE.pp(k)) GO TO 180
-  140 CONTINUE
-      IF (l.EQ.35) GO TO 130
-C
-C     STEP C. CREATION OF NEW POISSON PROBABILITIES P
-C             AND THEIR CUMULATIVES Q=PP(K)
-C
-  150 l = l + 1
-      DO 160 k = l,35
-          p = p*mu/float(k)
-          q = q + p
-          pp(k) = q
-          IF (u.LE.q) GO TO 170
-  160 CONTINUE
-      l = 35
-      GO TO 130
-
-  170 l = k
-  180 ignpoi = k
-      RETURN
-
-      END
-
-
-
-
-
diff --git a/scilab/modules/randlib/src/fortran/sexpo.f b/scilab/modules/randlib/src/fortran/sexpo.f
deleted file mode 100644 (file)
index 976ae33..0000000
+++ /dev/null
@@ -1,78 +0,0 @@
-      DOUBLE PRECISION FUNCTION sexpo()
-C**********************************************************************C
-C                                                                      C
-C                                                                      C
-C     (STANDARD-)  E X P O N E N T I A L   DISTRIBUTION                C
-C                                                                      C
-C                                                                      C
-C**********************************************************************C
-C**********************************************************************C
-C                                                                      C
-C     FOR DETAILS SEE:                                                 C
-C                                                                      C
-C               AHRENS, J.H. AND DIETER, U.                            C
-C               COMPUTER METHODS FOR SAMPLING FROM THE                 C
-C               EXPONENTIAL AND NORMAL DISTRIBUTIONS.                  C
-C               COMM. ACM, 15,10 (OCT. 1972), 873 - 882.               C
-C                                                                      C
-C     ALL STATEMENT NUMBERS CORRESPOND TO THE STEPS OF ALGORITHM       C
-C     'SA' IN THE ABOVE PAPER (SLIGHTLY MODIFIED IMPLEMENTATION)       C
-C                                                                      C
-C     Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of   C
-C     SUNIF.  The argument IR thus goes away.                          C
-C                                                                      C
-C**********************************************************************C
-C
-C
-C     Q(N) = SUM(ALOG(2.0)**K/K!)    K=1,..,N ,      THE HIGHEST N
-C     (HERE 8) IS DETERMINED BY Q(N)=1.0 WITHIN STANDARD PRECISION
-C
-C     JJV added a Save statement for q (in Data statement)
-C     .. Local Scalars ..
-      DOUBLE PRECISION a,q1,u,umin,ustar
-      INTEGER i
-C     ..
-C     .. Local Arrays ..
-      DOUBLE PRECISION q(8)
-C     ..
-C     .. External Functions ..
-      DOUBLE PRECISION ranf
-      EXTERNAL ranf
-C     ..
-C     .. Equivalences ..
-      EQUIVALENCE (q(1),q1)
-C     ..
-C     .. Save statement ..
-      SAVE q
-C     ..
-C     .. Data statements ..
-      DATA q/.6931472,.9333737,.9888778,.9984959,.9998293,.9999833,
-     +     .9999986,.9999999/
-C     ..
-C
-   10 a = 0.0
-      u = ranf()
-      GO TO 30
-
-   20 a = a + q1
-   30 u = u + u
-C     JJV changed the following to reflect the true algorithm and
-C     JJV prevent unpredictable behavior if U is initially 0.5.
-C      IF (u.LE.1.0) GO TO 20
-      IF (u.LT.1.0) GO TO 20
-   40 u = u - 1.0
-      IF (u.GT.q1) GO TO 60
-   50 sexpo = a + u
-      RETURN
-
-   60 i = 1
-      ustar = ranf()
-      umin = ustar
-   70 ustar = ranf()
-      IF (ustar.LT.umin) umin = ustar
-   80 i = i + 1
-      IF (u.GT.q(i)) GO TO 70
-   90 sexpo = a + umin*q1
-      RETURN
-
-      END
diff --git a/scilab/modules/randlib/src/fortran/sgamma.f b/scilab/modules/randlib/src/fortran/sgamma.f
deleted file mode 100644 (file)
index c2528a0..0000000
+++ /dev/null
@@ -1,236 +0,0 @@
-      DOUBLE PRECISION FUNCTION sgamma(a)
-C**********************************************************************C
-C                                                                      C
-C                                                                      C
-C     (STANDARD-)  G A M M A  DISTRIBUTION                             C
-C                                                                      C
-C                                                                      C
-C**********************************************************************C
-C**********************************************************************C
-C                                                                      C
-C               PARAMETER  A >= 1.0  !                                 C
-C                                                                      C
-C**********************************************************************C
-C                                                                      C
-C     FOR DETAILS SEE:                                                 C
-C                                                                      C
-C               AHRENS, J.H. AND DIETER, U.                            C
-C               GENERATING GAMMA VARIATES BY A                         C
-C               MODIFIED REJECTION TECHNIQUE.                          C
-C               COMM. ACM, 25,1 (JAN. 1982), 47 - 54.                  C
-C                                                                      C
-C     STEP NUMBERS CORRESPOND TO ALGORITHM 'GD' IN THE ABOVE PAPER     C
-C                                 (STRAIGHTFORWARD IMPLEMENTATION)     C
-C                                                                      C
-C     Modified by Barry W. Brown, Feb 3, 1988 to use RANF instead of   C
-C     SUNIF.  The argument IR thus goes away.                          C
-C                                                                      C
-C**********************************************************************C
-C                                                                      C
-C               PARAMETER  0.0 < A < 1.0  !                            C
-C                                                                      C
-C**********************************************************************C
-C                                                                      C
-C     FOR DETAILS SEE:                                                 C
-C                                                                      C
-C               AHRENS, J.H. AND DIETER, U.                            C
-C               COMPUTER METHODS FOR SAMPLING FROM GAMMA,              C
-C               BETA, POISSON AND BINOMIAL DISTRIBUTIONS.              C
-C               COMPUTING, 12 (1974), 223 - 246.                       C
-C                                                                      C
-C     (ADAPTED IMPLEMENTATION OF ALGORITHM 'GS' IN THE ABOVE PAPER)    C
-C                                                                      C
-C**********************************************************************C
-C
-C
-C     INPUT: A =PARAMETER (MEAN) OF THE STANDARD GAMMA DISTRIBUTION
-C     OUTPUT: SGAMMA = SAMPLE FROM THE GAMMA-(A)-DISTRIBUTION
-C
-C     COEFFICIENTS Q(K) - FOR Q0 = SUM(Q(K)*A**(-K))
-C     COEFFICIENTS A(K) - FOR Q = Q0+(T*T/2)*SUM(A(K)*V**K)
-C     COEFFICIENTS E(K) - FOR EXP(Q)-1 = SUM(E(K)*Q**K)
-C
-C     .. Scalar Arguments ..
-      DOUBLE PRECISION a
-C     ..
-C     .. Local Scalars .. (JJV added B0 to fix rare and subtle bug)
-      DOUBLE PRECISION a1,a2,a3,a4,a5,a6,a7,aa,aaa,b,b0,c,d,e,e1,e2,e3,
-     +     e4,e5,p,q,q0,
-     +     q1,q2,q3,q4,q5,q6,q7,r,s,s2,si,sqrt32,t,u,v,w,x
-C     ..
-C     .. External Functions ..
-      DOUBLE PRECISION ranf,sexpo,snorm
-      EXTERNAL ranf,sexpo,snorm
-C     ..
-C     .. Intrinsic Functions ..
-      INTRINSIC abs,log,exp,sign,sqrt
-C     ..
-C     .. Save statement ..
-C     JJV added Save statement for vars in Data satatements
-      SAVE aa,aaa,s2,s,d,q0,b,si,c,q1,q2,q3,q4,q5,q6,q7,a1,a2,a3,a4,a5,
-     +     a6,a7,e1,e2,e3,e4,e5,sqrt32
-C     ..
-C     .. Data statements ..
-C
-C     PREVIOUS A PRE-SET TO ZERO - AA IS A', AAA IS A"
-C     SQRT32 IS THE SQUAREROOT OF 32 = 5.656854249492380
-C
-      DATA q1,q2,q3,q4,q5,q6,q7/.04166669,.02083148,.00801191,.00144121,
-     +     -.00007388,.00024511,.00024240/
-      DATA a1,a2,a3,a4,a5,a6,a7/.3333333,-.2500030,.2000062,-.1662921,
-     +     .1423657,-.1367177,.1233795/
-      DATA e1,e2,e3,e4,e5/1.,.4999897,.1668290,.0407753,.0102930/
-      DATA aa/0.0/,aaa/0.0/,sqrt32/5.656854/
-C     ..
-C     .. Executable Statements ..
-C
-      IF (a.EQ.aa) GO TO 10
-      IF (a.LT.1.0) GO TO 130
-C
-C     STEP  1:  RECALCULATIONS OF S2,S,D IF A HAS CHANGED
-C
-      aa = a
-      s2 = a - 0.5
-      s = sqrt(s2)
-      d = sqrt32 - 12.0*s
-C
-C     STEP  2:  T=STANDARD NORMAL DEVIATE,
-C               X=(S,1/2)-NORMAL DEVIATE.
-C               IMMEDIATE ACCEPTANCE (I)
-C
-   10 t = snorm()
-      x = s + 0.5*t
-      sgamma = x*x
-      IF (t.GE.0.0) RETURN
-C
-C     STEP  3:  U= 0,1 -UNIFORM SAMPLE. SQUEEZE ACCEPTANCE (S)
-C
-      u = ranf()
-      IF (d*u.LE.t*t*t) RETURN
-C
-C     STEP  4:  RECALCULATIONS OF Q0,B,SI,C IF NECESSARY
-C
-      IF (a.EQ.aaa) GO TO 40
-      aaa = a
-      r = 1.0/a
-      q0 = ((((((q7*r+q6)*r+q5)*r+q4)*r+q3)*r+q2)*r+q1)*r
-C
-C               APPROXIMATION DEPENDING ON SIZE OF PARAMETER A
-C               THE CONSTANTS IN THE EXPRESSIONS FOR B, SI AND
-C               C WERE ESTABLISHED BY NUMERICAL EXPERIMENTS
-C
-      IF (a.LE.3.686) GO TO 30
-      IF (a.LE.13.022) GO TO 20
-C
-C               CASE 3:  A .GT. 13.022
-C
-      b = 1.77
-      si = .75
-      c = .1515/s
-      GO TO 40
-C
-C               CASE 2:  3.686 .LT. A .LE. 13.022
-C
-   20 b = 1.654 + .0076*s2
-      si = 1.68/s + .275
-      c = .062/s + .024
-      GO TO 40
-C
-C               CASE 1:  A .LE. 3.686
-C
-   30 b = .463 + s + .178*s2
-      si = 1.235
-      c = .195/s - .079 + .16*s
-C
-C     STEP  5:  NO QUOTIENT TEST IF X NOT POSITIVE
-C
-   40 IF (x.LE.0.0) GO TO 70
-C
-C     STEP  6:  CALCULATION OF V AND QUOTIENT Q
-C
-      v = t/ (s+s)
-      IF (abs(v).LE.0.25) GO TO 50
-      q = q0 - s*t + 0.25*t*t + (s2+s2)*log(1.0+v)
-      GO TO 60
-
-   50 q = q0 + 0.5*t*t* ((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v
-C
-C     STEP  7:  QUOTIENT ACCEPTANCE (Q)
-C
-   60 IF (log(1.0-u).LE.q) RETURN
-C
-C     STEP  8:  E=STANDARD EXPONENTIAL DEVIATE
-C               U= 0,1 -UNIFORM DEVIATE
-C               T=(B,SI)-DOUBLE EXPONENTIAL (LAPLACE) SAMPLE
-C
-   70 e = sexpo()
-      u = ranf()
-      u = u + u - 1.0
-      t = b + sign(si*e,u)
-C
-C     STEP  9:  REJECTION IF T .LT. TAU(1) = -.71874483771719
-C
-   80 IF (t.LT. (-.7187449)) GO TO 70
-C
-C     STEP 10:  CALCULATION OF V AND QUOTIENT Q
-C
-      v = t/ (s+s)
-      IF (abs(v).LE.0.25) GO TO 90
-      q = q0 - s*t + 0.25*t*t + (s2+s2)*log(1.0+v)
-      GO TO 100
-
-   90 q = q0 + 0.5*t*t* ((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v
-C
-C     STEP 11:  HAT ACCEPTANCE (H) (IF Q NOT POSITIVE GO TO STEP 8)
-C
-  100 IF (q.LE.0.0) GO TO 70
-      IF (q.LE.0.5) GO TO 110
-C
-C     JJV modified the code through line 125 to handle large Q case
-C
-      IF (q.LT.15.0) GO TO 105
-C
-C     JJV Here Q is large enough that Q = log(exp(Q) - 1.0) (for DOUBLE PRECISION Q)
-C     JJV so reformulate test at 120 in terms of one EXP, if not too big
-C     JJV 87.49823 is close to the largest DOUBLE PRECISION which can be
-C     JJV exponentiated (87.49823 = log(1.0E38))
-C
-      IF ((q+e-0.5*t*t).GT.87.49823) GO TO 125
-      IF (c*abs(u).GT.exp(q+e-0.5*t*t)) GO TO 70
-      GO TO 125
-
- 105  w = exp(q) - 1.0
-      GO TO 120
-
-  110 w = ((((e5*q+e4)*q+e3)*q+e2)*q+e1)*q
-C
-C               IF T IS REJECTED, SAMPLE AGAIN AT STEP 8
-C
-  120 IF (c*abs(u).GT.w*exp(e-0.5*t*t)) GO TO 70
- 125  x = s + 0.5*t
-      sgamma = x*x
-      RETURN
-C
-C     ALTERNATE METHOD FOR PARAMETERS A BELOW 1  (.3678794=EXP(-1.))
-C
-C     JJV changed B to B0 (which was added to declarations for this)
-C     JJV in 130 to END to fix rare and subtle bug.
-C     JJV Line: '130 aa = 0.0' was removed (unnecessary, wasteful).
-C     JJV Reasons: the state of AA only serves to tell the A .GE. 1.0
-C     JJV case if certain A-dependant constants need to be recalculated.
-C     JJV The A .LT. 1.0 case (here) no longer changes any of these, and
-C     JJV the recalculation of B (which used to change with an
-C     JJV A .LT. 1.0 call) is governed by the state of AAA anyway.
-C
- 130  b0 = 1.0 + .3678794*a
-  140 p = b0*ranf()
-      IF (p.GE.1.0) GO TO 150
-      sgamma = exp(log(p)/a)
-      IF (sexpo().LT.sgamma) GO TO 140
-      RETURN
-
-  150 sgamma = -log((b0-p)/a)
-      IF (sexpo().LT. (1.0-a)*log(sgamma)) GO TO 140
-      RETURN
-
-      END
index 21d22d7..86386ab 100644 (file)
@@ -1,9 +1,204 @@
-// Copyright INRIA
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - INRIA
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- TEST WITH GRAPHIC -->
+//---------------------Tests added by Sabine Gaüzère----------------------------
+//Same tests as below but with loops in order to have more examples.
 // test for grand
 prec = 1
  prec  =
  
     1.  
+// test for exp
+for i=0.1:0.1:50,
+  N=10000;A=i;
+  Rdev=grand(1,N,'exp',A);
+  RdevS=sort(Rdev);RdevS=RdevS($:-1:1)';
+  PS=(1:N)'/N;
+//plot2d(RdevS,PS);
+// theorical result
+  P=1-exp(-RdevS/A);
+//plot2d(RdevS,P,2,"000")
+  if norm(P-PS) > prec then bugmes();quit;end
+end
+// test for gamma
+for i=1:1:15,
+  for j=1:1:15,
+    N=10000;A=i;B=j;
+    Rdev=grand(1,N,'gam',A,B);
+    RdevS=sort(Rdev);RdevS=RdevS($:-1:1)';
+    PS=(1:N)'/N;
+//plot2d(RdevS,PS);
+// theorical result
+    [P]=cdfgam("PQ",RdevS,A*ones(RdevS),B*ones(RdevS));
+//plot2d(RdevS,P,2,"000")
+    if norm(P-PS)> prec then bugmes();quit;end
+  end
+end
+// test for beta random deviate
+for i=1:1:20,
+  for j=1:1:20,
+    N=10000;A=i;B=j;
+    Rdev=grand(1,N,'bet',A,B);
+    RdevS=sort(Rdev);RdevS=RdevS($:-1:1)';
+    PS=(1:N)'/N;
+//plot2d(RdevS,PS);
+// theorical result
+    [P]=cdfbet("PQ",RdevS,1-RdevS,A*ones(RdevS),B*ones(RdevS));
+//plot2d(RdevS,P,2,"000")
+    if norm(P-PS)> prec then bugmes();quit;end
+  end
+end
+// test for poi
+for i=50:1:70,
+  N=10000;A=i;
+  Rdev=grand(1,N,'poi',A);
+  RdevS=sort(Rdev);RdevS=RdevS($:-1:1)';
+  PS=(1:N)'/N;
+//plot2d(RdevS,PS);
+// theorical result
+  [P]=cdfpoi("PQ",RdevS,A*ones(RdevS));
+//plot2d(RdevS,P,2,"000")
+// idem need an other test P is piecewize linear and PS
+// linear
+//if norm(P-PS) > prec then bugmes();quit;end
+end
+//Tests to verify the compatibility with a matrix
+N=100;A=1;B=3;
+Rdev=grand(N,N,'bet',A,B);
+Rdev=matrix(Rdev,1,N^2);
+RdevS=sort(Rdev);RdevS=RdevS($:-1:1)';
+PS=(1:(N^2))'/(N^2);
+[P]=cdfbet("PQ",RdevS,1-RdevS,A*ones(RdevS),B*ones(RdevS));
+if norm(P-PS) > prec then bugmes();quit;end
+N=100;A=1;B=3;
+Rdev=grand(N,N,'gam',A,B);
+Rdev=matrix(Rdev,1,N^2);
+RdevS=sort(Rdev);RdevS=RdevS($:-1:1)';
+PS=(1:(N^2))'/(N^2);
+[P]=cdfgam("PQ",RdevS,A*ones(RdevS),B*ones(RdevS));
+if norm(P-PS) > prec then bugmes();quit;end
+N=100;A=5;B=0.7;
+Rdev=grand(N,N,'bin',A,B);
+Rdev=matrix(Rdev,1,N^2);
+RdevS=sort(Rdev);RdevS=RdevS($:-1:1)';
+PS=(1:(N^2))'/(N^2);
+[P]=cdfbin("PQ",RdevS,A*ones(RdevS),B*ones(RdevS),(1-B)*ones(RdevS));
+//if norm(P-PS) > prec then bugmes();quit;end
+N=100;A=50;
+Rdev=grand(N,N,'poi',A);
+Rdev=matrix(Rdev,1,N^2);
+RdevS=sort(Rdev);RdevS=RdevS($:-1:1)';
+PS=(1:(N^2))'/(N^2);
+[P]=cdfpoi("PQ",RdevS,A*ones(RdevS));
+//if norm(P-PS) > prec then bugmes();quit;end
+N=100;A=2;
+Rdev=grand(N,N,'exp',A);
+Rdev=matrix(Rdev,1,N^2);
+RdevS=sort(Rdev);RdevS=RdevS($:-1:1)';
+PS=(1:(N^2))'/(N^2);
+P=1-exp(-RdevS/A);
+if norm(P-PS) > prec then bugmes();quit;end
+//Tests with graphic
+//Comparison of pseudo-random numbers following an exponential distribution
+//and the density of this distribution
+//Parameter of the distribution which can be modified
+lambda=1.6;
+// Number of random variable generated
+N=100000;
+//Generation of a vector of numbers following an exponential distribution
+X = grand(1,N,"exp",lambda);
+xbasc();
+//Discretisation of the abscisses in classes
+classes = linspace(0,12,25);
+//Draw in histogram
+histplot(classes,X)
+//Draw the density
+x=linspace(0,12,25);
+y = (1/lambda)*exp(-(1/lambda)*x);
+plot2d(x,y,3);
+f=gcf();
+delete(f);
+//Comparison of pseudo-random numbers following a beta distribution
+//and the density of this distribution
+//Parameters of the distribution which can be modified
+A=1;B=3;
+// Number of random variable generated
+N=100000;
+//Generation of a vector of numbers following a beta distribution
+X = grand(1,N,"bet",A,B);
+xbasc();
+//Discretisation of the abscisses in classes
+classes = linspace(0,1,50);
+//Draw in histogram
+histplot(classes,X)
+//Draw the density
+x=linspace(0,1,50);
+y = (1/(beta(A,B))).*(x^(A-1)).*((1-x)^(B-1)) ;
+plot2d(x,y,2);
+f=gcf();
+delete(f);
+//Comparison of pseudo-random numbers following a gamma distribution
+//and the density of this distribution
+//Parameters of the distribution which can be modified
+A=2;B=1;
+// Number of random variable generated
+N=100000;
+//Generation of a vector of numbers following a gamma distribution
+X = grand(1,N,"gam",A,B);
+xbasc();
+//Discretisation of the abscisses in classes
+classes = linspace(0,2,50);
+//Draw in histogram
+histplot(classes,X)
+//Draw the density
+x=linspace(0,2,50);
+y = (B/(gamma(A))).*exp(-B*x).*(B*x)^(A-1);
+plot2d(x,y,2);
+f=gcf();
+delete(f);
+//Comparison of pseudo-random numbers following a binomial distribution
+//and the density of this distribution
+//Parameters of the distribution which can be modified
+n=50;p=0.3;
+// Number of random variable generated
+N=100000;
+//Generation of a vector of numbers following a binomial distribution
+X = grand(1,N,"bin",n,p);
+xbasc();
+//Discretisation of the abscisses in classes
+classes = linspace(0,n,n+1);
+//Draw in histogram
+histplot(classes,X)
+//Draw the density
+x=linspace(0,n,n+1);
+y = binomial(p,n);
+plot2d(x,y,2);
+f=gcf();
+delete(f);
+//Comparison of pseudo-random numbers following a poisson distribution
+//and the density of this distribution
+//Parameters of the distribution which can be modified
+mu=50;
+// Number of random variable generated
+N=100000;
+//Generation of a vector of numbers following a poisson distribution
+X = grand(1,N,"poi",mu);
+xbasc();
+//Discretisation of the abscisses in classes
+classes = linspace(0,2*mu,101);
+//Draw in histogram
+histplot(classes,X)
+//Draw the density
+[x]=linspace(0,2*mu,101);
+y = exp(-mu).*(mu^x)./factorial(x);
+plot2d(x,y,2);
+f=gcf();
+delete(f);
+//---------------------------------------------------------------------------
 // test for beta random deviate
 N=10000;A=1;B=3;
 Rdev=grand(1,N,'bet',A,B);
index 6898df3..f6247a2 100644 (file)
 //
 //  This file is distributed under the same license as the Scilab package.
 // =============================================================================
+// <-- TEST WITH GRAPHIC -->
 
+//---------------------Tests added by Sabine Gaüzère----------------------------
+
+
+//Same tests as below but with loops in order to have more examples.
 // test for grand
 prec = 1
 
+// test for exp 
+
+for i=0.1:0.1:50,
+  N=10000;A=i;
+  Rdev=grand(1,N,'exp',A);
+  RdevS=sort(Rdev);RdevS=RdevS($:-1:1)';
+  PS=(1:N)'/N;
+//plot2d(RdevS,PS);
+// theorical result 
+  P=1-exp(-RdevS/A);
+//plot2d(RdevS,P,2,"000")
+  if norm(P-PS) > prec then pause,end
+end
+
+
+// test for gamma 
+for i=1:1:15,
+  for j=1:1:15,
+    N=10000;A=i;B=j;
+    Rdev=grand(1,N,'gam',A,B); 
+    RdevS=sort(Rdev);RdevS=RdevS($:-1:1)';
+    PS=(1:N)'/N;
+//plot2d(RdevS,PS);
+// theorical result 
+    [P]=cdfgam("PQ",RdevS,A*ones(RdevS),B*ones(RdevS));
+//plot2d(RdevS,P,2,"000")
+    if norm(P-PS)> prec then pause,end
+  end
+end
+
+
+// test for beta random deviate 
+for i=1:1:20,
+  for j=1:1:20,
+    N=10000;A=i;B=j;
+    Rdev=grand(1,N,'bet',A,B); 
+    RdevS=sort(Rdev);RdevS=RdevS($:-1:1)';
+    PS=(1:N)'/N;
+//plot2d(RdevS,PS);
+// theorical result 
+    [P]=cdfbet("PQ",RdevS,1-RdevS,A*ones(RdevS),B*ones(RdevS));
+//plot2d(RdevS,P,2,"000")
+    if norm(P-PS)> prec then pause,end
+  end
+end
+
+// test for poi 
+for i=50:1:70,
+  N=10000;A=i;
+  Rdev=grand(1,N,'poi',A);
+  RdevS=sort(Rdev);RdevS=RdevS($:-1:1)';
+  PS=(1:N)'/N;
+//plot2d(RdevS,PS);
+// theorical result 
+  [P]=cdfpoi("PQ",RdevS,A*ones(RdevS));
+//plot2d(RdevS,P,2,"000")
+// idem need an other test P is piecewize linear and PS 
+// linear 
+//if norm(P-PS) > prec then pause,end
+end
+
+//Tests to verify the compatibility with a matrix
+  
+N=100;A=1;B=3;
+Rdev=grand(N,N,'bet',A,B); 
+Rdev=matrix(Rdev,1,N^2);
+RdevS=sort(Rdev);RdevS=RdevS($:-1:1)';
+PS=(1:(N^2))'/(N^2);
+[P]=cdfbet("PQ",RdevS,1-RdevS,A*ones(RdevS),B*ones(RdevS));
+if norm(P-PS) > prec then pause,end
+  
+N=100;A=1;B=3;
+Rdev=grand(N,N,'gam',A,B); 
+Rdev=matrix(Rdev,1,N^2);
+RdevS=sort(Rdev);RdevS=RdevS($:-1:1)';
+PS=(1:(N^2))'/(N^2);
+[P]=cdfgam("PQ",RdevS,A*ones(RdevS),B*ones(RdevS));
+if norm(P-PS) > prec then pause,end
+  
+N=100;A=5;B=0.7;
+Rdev=grand(N,N,'bin',A,B); 
+Rdev=matrix(Rdev,1,N^2);
+RdevS=sort(Rdev);RdevS=RdevS($:-1:1)';
+PS=(1:(N^2))'/(N^2);
+[P]=cdfbin("PQ",RdevS,A*ones(RdevS),B*ones(RdevS),(1-B)*ones(RdevS));
+//if norm(P-PS) > prec then pause,end
+
+N=100;A=50;
+Rdev=grand(N,N,'poi',A); 
+Rdev=matrix(Rdev,1,N^2);
+RdevS=sort(Rdev);RdevS=RdevS($:-1:1)';
+PS=(1:(N^2))'/(N^2);
+[P]=cdfpoi("PQ",RdevS,A*ones(RdevS));
+//if norm(P-PS) > prec then pause,end
+
+N=100;A=2;
+Rdev=grand(N,N,'exp',A); 
+Rdev=matrix(Rdev,1,N^2);
+RdevS=sort(Rdev);RdevS=RdevS($:-1:1)';
+PS=(1:(N^2))'/(N^2);
+P=1-exp(-RdevS/A);
+if norm(P-PS) > prec then pause,end
+
+//Tests with graphic
+
+//Comparison of pseudo-random numbers following an exponential distribution 
+//and the density of this distribution
+//Parameter of the distribution which can be modified
+lambda=1.6;
+// Number of random variable generated
+N=100000;
+//Generation of a vector of numbers following an exponential distribution
+X = grand(1,N,"exp",lambda);
+xbasc();
+//Discretisation of the abscisses in classes
+classes = linspace(0,12,25);
+//Draw in histogram
+histplot(classes,X)
+//Draw the density 
+x=linspace(0,12,25);
+y = (1/lambda)*exp(-(1/lambda)*x);
+plot2d(x,y,3);
+f=gcf();
+delete(f);
+
+//Comparison of pseudo-random numbers following a beta distribution 
+//and the density of this distribution
+//Parameters of the distribution which can be modified
+A=1;B=3;
+// Number of random variable generated
+N=100000;
+//Generation of a vector of numbers following a beta distribution
+X = grand(1,N,"bet",A,B);
+xbasc();
+//Discretisation of the abscisses in classes
+classes = linspace(0,1,50);
+//Draw in histogram
+histplot(classes,X)
+//Draw the density 
+x=linspace(0,1,50);
+y = (1/(beta(A,B))).*(x^(A-1)).*((1-x)^(B-1)) ;
+plot2d(x,y,2);
+f=gcf();
+delete(f);
+
+//Comparison of pseudo-random numbers following a gamma distribution 
+//and the density of this distribution
+//Parameters of the distribution which can be modified
+A=2;B=1;
+// Number of random variable generated
+N=100000;
+//Generation of a vector of numbers following a gamma distribution
+X = grand(1,N,"gam",A,B);
+xbasc();
+//Discretisation of the abscisses in classes
+classes = linspace(0,2,50);
+//Draw in histogram
+histplot(classes,X)
+//Draw the density 
+x=linspace(0,2,50);
+y = (B/(gamma(A))).*exp(-B*x).*(B*x)^(A-1);
+plot2d(x,y,2);
+f=gcf();
+delete(f);
+
+
+//Comparison of pseudo-random numbers following a binomial distribution 
+//and the density of this distribution
+//Parameters of the distribution which can be modified
+n=50;p=0.3;
+// Number of random variable generated
+N=100000;
+//Generation of a vector of numbers following a binomial distribution
+X = grand(1,N,"bin",n,p);
+xbasc();
+//Discretisation of the abscisses in classes
+classes = linspace(0,n,n+1);
+//Draw in histogram
+histplot(classes,X)
+//Draw the density
+x=linspace(0,n,n+1);
+y = binomial(p,n);
+plot2d(x,y,2);
+f=gcf();
+delete(f);
+
+//Comparison of pseudo-random numbers following a poisson distribution 
+//and the density of this distribution
+//Parameters of the distribution which can be modified
+mu=50;
+// Number of random variable generated
+N=100000;
+//Generation of a vector of numbers following a poisson distribution
+X = grand(1,N,"poi",mu);
+xbasc();
+//Discretisation of the abscisses in classes
+classes = linspace(0,2*mu,101);
+//Draw in histogram
+histplot(classes,X)
+//Draw the density
+[x]=linspace(0,2*mu,101);
+y = exp(-mu).*(mu^x)./factorial(x);
+plot2d(x,y,2);
+f=gcf();
+delete(f);
+
+//---------------------------------------------------------------------------
 
 // test for beta random deviate