Mercurial > hg > octave-nkf
annotate liboctave/numeric/randmtzig.c @ 15271:648dabbb4c6b
build: Refactor liboctave into multiple subdirectories. Move libcruft into liboctave.
* array/Array-C.cc, array/Array-b.cc, array/Array-ch.cc, array/Array-d.cc,
array/Array-f.cc, array/Array-fC.cc, array/Array-i.cc, array/Array-idx-vec.cc,
array/Array-s.cc, array/Array-str.cc, array/Array-util.cc, array/Array-util.h,
array/Array-voidp.cc, array/Array.cc, array/Array.h, array/Array2.h,
array/Array3.h, array/ArrayN.h, array/CColVector.cc, array/CColVector.h,
array/CDiagMatrix.cc, array/CDiagMatrix.h, array/CMatrix.cc, array/CMatrix.h,
array/CNDArray.cc, array/CNDArray.h, array/CRowVector.cc, array/CRowVector.h,
array/CSparse.cc, array/CSparse.h, array/DiagArray2.cc, array/DiagArray2.h,
array/MArray-C.cc, array/MArray-d.cc, array/MArray-decl.h, array/MArray-defs.h,
array/MArray-f.cc, array/MArray-fC.cc, array/MArray-i.cc, array/MArray-s.cc,
array/MArray.cc, array/MArray.h, array/MArray2.h, array/MArrayN.h,
array/MDiagArray2.cc, array/MDiagArray2.h, array/MSparse-C.cc,
array/MSparse-d.cc, array/MSparse-defs.h, array/MSparse.cc, array/MSparse.h,
array/Matrix.h, array/MatrixType.cc, array/MatrixType.h, array/PermMatrix.cc,
array/PermMatrix.h, array/Range.cc, array/Range.h, array/Sparse-C.cc,
array/Sparse-b.cc, array/Sparse-d.cc, array/Sparse.cc, array/Sparse.h,
array/boolMatrix.cc, array/boolMatrix.h, array/boolNDArray.cc,
array/boolNDArray.h, array/boolSparse.cc, array/boolSparse.h,
array/chMatrix.cc, array/chMatrix.h, array/chNDArray.cc, array/chNDArray.h,
array/dColVector.cc, array/dColVector.h, array/dDiagMatrix.cc,
array/dDiagMatrix.h, array/dMatrix.cc, array/dMatrix.h, array/dNDArray.cc,
array/dNDArray.h, array/dRowVector.cc, array/dRowVector.h, array/dSparse.cc,
array/dSparse.h, array/dim-vector.cc, array/dim-vector.h, array/fCColVector.cc,
array/fCColVector.h, array/fCDiagMatrix.cc, array/fCDiagMatrix.h,
array/fCMatrix.cc, array/fCMatrix.h, array/fCNDArray.cc, array/fCNDArray.h,
array/fCRowVector.cc, array/fCRowVector.h, array/fColVector.cc,
array/fColVector.h, array/fDiagMatrix.cc, array/fDiagMatrix.h,
array/fMatrix.cc, array/fMatrix.h, array/fNDArray.cc, array/fNDArray.h,
array/fRowVector.cc, array/fRowVector.h, array/idx-vector.cc,
array/idx-vector.h, array/int16NDArray.cc, array/int16NDArray.h,
array/int32NDArray.cc, array/int32NDArray.h, array/int64NDArray.cc,
array/int64NDArray.h, array/int8NDArray.cc, array/int8NDArray.h,
array/intNDArray.cc, array/intNDArray.h, array/module.mk,
array/uint16NDArray.cc, array/uint16NDArray.h, array/uint32NDArray.cc,
array/uint32NDArray.h, array/uint64NDArray.cc, array/uint64NDArray.h,
array/uint8NDArray.cc, array/uint8NDArray.h:
Moved from liboctave dir to array subdirectory.
* cruft/Makefile.am, cruft/amos/README, cruft/amos/cacai.f, cruft/amos/cacon.f,
cruft/amos/cairy.f, cruft/amos/casyi.f, cruft/amos/cbesh.f, cruft/amos/cbesi.f,
cruft/amos/cbesj.f, cruft/amos/cbesk.f, cruft/amos/cbesy.f, cruft/amos/cbinu.f,
cruft/amos/cbiry.f, cruft/amos/cbknu.f, cruft/amos/cbuni.f, cruft/amos/cbunk.f,
cruft/amos/ckscl.f, cruft/amos/cmlri.f, cruft/amos/crati.f, cruft/amos/cs1s2.f,
cruft/amos/cseri.f, cruft/amos/cshch.f, cruft/amos/cuchk.f, cruft/amos/cunhj.f,
cruft/amos/cuni1.f, cruft/amos/cuni2.f, cruft/amos/cunik.f, cruft/amos/cunk1.f,
cruft/amos/cunk2.f, cruft/amos/cuoik.f, cruft/amos/cwrsk.f,
cruft/amos/dgamln.f, cruft/amos/gamln.f, cruft/amos/module.mk,
cruft/amos/xzabs.f, cruft/amos/xzexp.f, cruft/amos/xzlog.f,
cruft/amos/xzsqrt.f, cruft/amos/zacai.f, cruft/amos/zacon.f,
cruft/amos/zairy.f, cruft/amos/zasyi.f, cruft/amos/zbesh.f, cruft/amos/zbesi.f,
cruft/amos/zbesj.f, cruft/amos/zbesk.f, cruft/amos/zbesy.f, cruft/amos/zbinu.f,
cruft/amos/zbiry.f, cruft/amos/zbknu.f, cruft/amos/zbuni.f, cruft/amos/zbunk.f,
cruft/amos/zdiv.f, cruft/amos/zkscl.f, cruft/amos/zmlri.f, cruft/amos/zmlt.f,
cruft/amos/zrati.f, cruft/amos/zs1s2.f, cruft/amos/zseri.f, cruft/amos/zshch.f,
cruft/amos/zuchk.f, cruft/amos/zunhj.f, cruft/amos/zuni1.f, cruft/amos/zuni2.f,
cruft/amos/zunik.f, cruft/amos/zunk1.f, cruft/amos/zunk2.f, cruft/amos/zuoik.f,
cruft/amos/zwrsk.f, cruft/blas-xtra/cconv2.f, cruft/blas-xtra/cdotc3.f,
cruft/blas-xtra/cmatm3.f, cruft/blas-xtra/csconv2.f, cruft/blas-xtra/dconv2.f,
cruft/blas-xtra/ddot3.f, cruft/blas-xtra/dmatm3.f, cruft/blas-xtra/module.mk,
cruft/blas-xtra/sconv2.f, cruft/blas-xtra/sdot3.f, cruft/blas-xtra/smatm3.f,
cruft/blas-xtra/xcdotc.f, cruft/blas-xtra/xcdotu.f, cruft/blas-xtra/xddot.f,
cruft/blas-xtra/xdnrm2.f, cruft/blas-xtra/xdznrm2.f, cruft/blas-xtra/xerbla.f,
cruft/blas-xtra/xscnrm2.f, cruft/blas-xtra/xsdot.f, cruft/blas-xtra/xsnrm2.f,
cruft/blas-xtra/xzdotc.f, cruft/blas-xtra/xzdotu.f, cruft/blas-xtra/zconv2.f,
cruft/blas-xtra/zdconv2.f, cruft/blas-xtra/zdotc3.f, cruft/blas-xtra/zmatm3.f,
cruft/daspk/datv.f, cruft/daspk/dcnst0.f, cruft/daspk/dcnstr.f,
cruft/daspk/ddasic.f, cruft/daspk/ddasid.f, cruft/daspk/ddasik.f,
cruft/daspk/ddaspk.f, cruft/daspk/ddstp.f, cruft/daspk/ddwnrm.f,
cruft/daspk/dfnrmd.f, cruft/daspk/dfnrmk.f, cruft/daspk/dhels.f,
cruft/daspk/dheqr.f, cruft/daspk/dinvwt.f, cruft/daspk/dlinsd.f,
cruft/daspk/dlinsk.f, cruft/daspk/dmatd.f, cruft/daspk/dnedd.f,
cruft/daspk/dnedk.f, cruft/daspk/dnsd.f, cruft/daspk/dnsid.f,
cruft/daspk/dnsik.f, cruft/daspk/dnsk.f, cruft/daspk/dorth.f,
cruft/daspk/dslvd.f, cruft/daspk/dslvk.f, cruft/daspk/dspigm.f,
cruft/daspk/dyypnw.f, cruft/daspk/module.mk, cruft/dasrt/ddasrt.f,
cruft/dasrt/drchek.f, cruft/dasrt/droots.f, cruft/dasrt/module.mk,
cruft/dassl/ddaini.f, cruft/dassl/ddajac.f, cruft/dassl/ddanrm.f,
cruft/dassl/ddaslv.f, cruft/dassl/ddassl.f, cruft/dassl/ddastp.f,
cruft/dassl/ddatrp.f, cruft/dassl/ddawts.f, cruft/dassl/module.mk,
cruft/fftpack/cfftb.f, cruft/fftpack/cfftb1.f, cruft/fftpack/cfftf.f,
cruft/fftpack/cfftf1.f, cruft/fftpack/cffti.f, cruft/fftpack/cffti1.f,
cruft/fftpack/fftpack.doc, cruft/fftpack/module.mk, cruft/fftpack/passb.f,
cruft/fftpack/passb2.f, cruft/fftpack/passb3.f, cruft/fftpack/passb4.f,
cruft/fftpack/passb5.f, cruft/fftpack/passf.f, cruft/fftpack/passf2.f,
cruft/fftpack/passf3.f, cruft/fftpack/passf4.f, cruft/fftpack/passf5.f,
cruft/fftpack/zfftb.f, cruft/fftpack/zfftb1.f, cruft/fftpack/zfftf.f,
cruft/fftpack/zfftf1.f, cruft/fftpack/zffti.f, cruft/fftpack/zffti1.f,
cruft/fftpack/zpassb.f, cruft/fftpack/zpassb2.f, cruft/fftpack/zpassb3.f,
cruft/fftpack/zpassb4.f, cruft/fftpack/zpassb5.f, cruft/fftpack/zpassf.f,
cruft/fftpack/zpassf2.f, cruft/fftpack/zpassf3.f, cruft/fftpack/zpassf4.f,
cruft/fftpack/zpassf5.f, cruft/lapack-xtra/crsf2csf.f,
cruft/lapack-xtra/module.mk, cruft/lapack-xtra/xclange.f,
cruft/lapack-xtra/xdlamch.f, cruft/lapack-xtra/xdlange.f,
cruft/lapack-xtra/xilaenv.f, cruft/lapack-xtra/xslamch.f,
cruft/lapack-xtra/xslange.f, cruft/lapack-xtra/xzlange.f,
cruft/lapack-xtra/zrsf2csf.f, cruft/link-deps.mk, cruft/misc/blaswrap.c,
cruft/misc/cquit.c, cruft/misc/d1mach-tst.for, cruft/misc/d1mach.f,
cruft/misc/f77-extern.cc, cruft/misc/f77-fcn.c, cruft/misc/f77-fcn.h,
cruft/misc/i1mach.f, cruft/misc/lo-error.c, cruft/misc/lo-error.h,
cruft/misc/module.mk, cruft/misc/quit.cc, cruft/misc/quit.h,
cruft/misc/r1mach.f, cruft/mkf77def.in, cruft/odepack/cfode.f,
cruft/odepack/dlsode.f, cruft/odepack/ewset.f, cruft/odepack/intdy.f,
cruft/odepack/module.mk, cruft/odepack/prepj.f, cruft/odepack/scfode.f,
cruft/odepack/sewset.f, cruft/odepack/sintdy.f, cruft/odepack/slsode.f,
cruft/odepack/solsy.f, cruft/odepack/sprepj.f, cruft/odepack/ssolsy.f,
cruft/odepack/sstode.f, cruft/odepack/stode.f, cruft/odepack/svnorm.f,
cruft/odepack/vnorm.f, cruft/ordered-qz/README, cruft/ordered-qz/dsubsp.f,
cruft/ordered-qz/exchqz.f, cruft/ordered-qz/module.mk,
cruft/ordered-qz/sexchqz.f, cruft/ordered-qz/ssubsp.f, cruft/quadpack/dqagi.f,
cruft/quadpack/dqagie.f, cruft/quadpack/dqagp.f, cruft/quadpack/dqagpe.f,
cruft/quadpack/dqelg.f, cruft/quadpack/dqk15i.f, cruft/quadpack/dqk21.f,
cruft/quadpack/dqpsrt.f, cruft/quadpack/module.mk, cruft/quadpack/qagi.f,
cruft/quadpack/qagie.f, cruft/quadpack/qagp.f, cruft/quadpack/qagpe.f,
cruft/quadpack/qelg.f, cruft/quadpack/qk15i.f, cruft/quadpack/qk21.f,
cruft/quadpack/qpsrt.f, cruft/quadpack/xerror.f, cruft/ranlib/Basegen.doc,
cruft/ranlib/HOWTOGET, cruft/ranlib/README, cruft/ranlib/advnst.f,
cruft/ranlib/genbet.f, cruft/ranlib/genchi.f, cruft/ranlib/genexp.f,
cruft/ranlib/genf.f, cruft/ranlib/gengam.f, cruft/ranlib/genmn.f,
cruft/ranlib/genmul.f, cruft/ranlib/gennch.f, cruft/ranlib/gennf.f,
cruft/ranlib/gennor.f, cruft/ranlib/genprm.f, cruft/ranlib/genunf.f,
cruft/ranlib/getcgn.f, cruft/ranlib/getsd.f, cruft/ranlib/ignbin.f,
cruft/ranlib/ignlgi.f, cruft/ranlib/ignnbn.f, cruft/ranlib/ignpoi.f,
cruft/ranlib/ignuin.f, cruft/ranlib/initgn.f, cruft/ranlib/inrgcm.f,
cruft/ranlib/lennob.f, cruft/ranlib/mltmod.f, cruft/ranlib/module.mk,
cruft/ranlib/phrtsd.f, cruft/ranlib/qrgnin.f, cruft/ranlib/randlib.chs,
cruft/ranlib/randlib.fdoc, cruft/ranlib/ranf.f, cruft/ranlib/setall.f,
cruft/ranlib/setant.f, cruft/ranlib/setgmn.f, cruft/ranlib/setsd.f,
cruft/ranlib/sexpo.f, cruft/ranlib/sgamma.f, cruft/ranlib/snorm.f,
cruft/ranlib/tstbot.for, cruft/ranlib/tstgmn.for, cruft/ranlib/tstmid.for,
cruft/ranlib/wrap.f, cruft/slatec-err/fdump.f, cruft/slatec-err/ixsav.f,
cruft/slatec-err/j4save.f, cruft/slatec-err/module.mk,
cruft/slatec-err/xerclr.f, cruft/slatec-err/xercnt.f,
cruft/slatec-err/xerhlt.f, cruft/slatec-err/xermsg.f,
cruft/slatec-err/xerprn.f, cruft/slatec-err/xerrwd.f,
cruft/slatec-err/xersve.f, cruft/slatec-err/xgetf.f, cruft/slatec-err/xgetua.f,
cruft/slatec-err/xsetf.f, cruft/slatec-err/xsetua.f, cruft/slatec-fn/acosh.f,
cruft/slatec-fn/albeta.f, cruft/slatec-fn/algams.f, cruft/slatec-fn/alngam.f,
cruft/slatec-fn/alnrel.f, cruft/slatec-fn/asinh.f, cruft/slatec-fn/atanh.f,
cruft/slatec-fn/betai.f, cruft/slatec-fn/csevl.f, cruft/slatec-fn/d9gmit.f,
cruft/slatec-fn/d9lgic.f, cruft/slatec-fn/d9lgit.f, cruft/slatec-fn/d9lgmc.f,
cruft/slatec-fn/dacosh.f, cruft/slatec-fn/dasinh.f, cruft/slatec-fn/datanh.f,
cruft/slatec-fn/dbetai.f, cruft/slatec-fn/dcsevl.f, cruft/slatec-fn/derf.f,
cruft/slatec-fn/derfc.in.f, cruft/slatec-fn/dgami.f, cruft/slatec-fn/dgamit.f,
cruft/slatec-fn/dgamlm.f, cruft/slatec-fn/dgamma.f, cruft/slatec-fn/dgamr.f,
cruft/slatec-fn/dlbeta.f, cruft/slatec-fn/dlgams.f, cruft/slatec-fn/dlngam.f,
cruft/slatec-fn/dlnrel.f, cruft/slatec-fn/dpchim.f, cruft/slatec-fn/dpchst.f,
cruft/slatec-fn/erf.f, cruft/slatec-fn/erfc.in.f, cruft/slatec-fn/gami.f,
cruft/slatec-fn/gamit.f, cruft/slatec-fn/gamlim.f, cruft/slatec-fn/gamma.f,
cruft/slatec-fn/gamr.f, cruft/slatec-fn/initds.f, cruft/slatec-fn/inits.f,
cruft/slatec-fn/module.mk, cruft/slatec-fn/pchim.f, cruft/slatec-fn/pchst.f,
cruft/slatec-fn/r9gmit.f, cruft/slatec-fn/r9lgic.f, cruft/slatec-fn/r9lgit.f,
cruft/slatec-fn/r9lgmc.f, cruft/slatec-fn/xacosh.f, cruft/slatec-fn/xasinh.f,
cruft/slatec-fn/xatanh.f, cruft/slatec-fn/xbetai.f, cruft/slatec-fn/xdacosh.f,
cruft/slatec-fn/xdasinh.f, cruft/slatec-fn/xdatanh.f,
cruft/slatec-fn/xdbetai.f, cruft/slatec-fn/xderf.f, cruft/slatec-fn/xderfc.f,
cruft/slatec-fn/xdgami.f, cruft/slatec-fn/xdgamit.f, cruft/slatec-fn/xdgamma.f,
cruft/slatec-fn/xerf.f, cruft/slatec-fn/xerfc.f, cruft/slatec-fn/xgamma.f,
cruft/slatec-fn/xgmainc.f, cruft/slatec-fn/xsgmainc.f:
Moved from top-level libcruft to cruft directory below liboctave.
* numeric/CmplxAEPBAL.cc, numeric/CmplxAEPBAL.h, numeric/CmplxCHOL.cc,
numeric/CmplxCHOL.h, numeric/CmplxGEPBAL.cc, numeric/CmplxGEPBAL.h,
numeric/CmplxHESS.cc, numeric/CmplxHESS.h, numeric/CmplxLU.cc,
numeric/CmplxLU.h, numeric/CmplxQR.cc, numeric/CmplxQR.h, numeric/CmplxQRP.cc,
numeric/CmplxQRP.h, numeric/CmplxSCHUR.cc, numeric/CmplxSCHUR.h,
numeric/CmplxSVD.cc, numeric/CmplxSVD.h, numeric/CollocWt.cc,
numeric/CollocWt.h, numeric/DAE.h, numeric/DAEFunc.h, numeric/DAERT.h,
numeric/DAERTFunc.h, numeric/DASPK-opts.in, numeric/DASPK.cc, numeric/DASPK.h,
numeric/DASRT-opts.in, numeric/DASRT.cc, numeric/DASRT.h,
numeric/DASSL-opts.in, numeric/DASSL.cc, numeric/DASSL.h, numeric/DET.h,
numeric/EIG.cc, numeric/EIG.h, numeric/LSODE-opts.in, numeric/LSODE.cc,
numeric/LSODE.h, numeric/ODE.h, numeric/ODEFunc.h, numeric/ODES.cc,
numeric/ODES.h, numeric/ODESFunc.h, numeric/Quad-opts.in, numeric/Quad.cc,
numeric/Quad.h, numeric/SparseCmplxCHOL.cc, numeric/SparseCmplxCHOL.h,
numeric/SparseCmplxLU.cc, numeric/SparseCmplxLU.h, numeric/SparseCmplxQR.cc,
numeric/SparseCmplxQR.h, numeric/SparseQR.cc, numeric/SparseQR.h,
numeric/SparsedbleCHOL.cc, numeric/SparsedbleCHOL.h, numeric/SparsedbleLU.cc,
numeric/SparsedbleLU.h, numeric/base-aepbal.h, numeric/base-dae.h,
numeric/base-de.h, numeric/base-lu.cc, numeric/base-lu.h, numeric/base-min.h,
numeric/base-qr.cc, numeric/base-qr.h, numeric/bsxfun-decl.h,
numeric/bsxfun-defs.cc, numeric/bsxfun.h, numeric/dbleAEPBAL.cc,
numeric/dbleAEPBAL.h, numeric/dbleCHOL.cc, numeric/dbleCHOL.h,
numeric/dbleGEPBAL.cc, numeric/dbleGEPBAL.h, numeric/dbleHESS.cc,
numeric/dbleHESS.h, numeric/dbleLU.cc, numeric/dbleLU.h, numeric/dbleQR.cc,
numeric/dbleQR.h, numeric/dbleQRP.cc, numeric/dbleQRP.h, numeric/dbleSCHUR.cc,
numeric/dbleSCHUR.h, numeric/dbleSVD.cc, numeric/dbleSVD.h,
numeric/eigs-base.cc, numeric/fCmplxAEPBAL.cc, numeric/fCmplxAEPBAL.h,
numeric/fCmplxCHOL.cc, numeric/fCmplxCHOL.h, numeric/fCmplxGEPBAL.cc,
numeric/fCmplxGEPBAL.h, numeric/fCmplxHESS.cc, numeric/fCmplxHESS.h,
numeric/fCmplxLU.cc, numeric/fCmplxLU.h, numeric/fCmplxQR.cc,
numeric/fCmplxQR.h, numeric/fCmplxQRP.cc, numeric/fCmplxQRP.h,
numeric/fCmplxSCHUR.cc, numeric/fCmplxSCHUR.h, numeric/fCmplxSVD.cc,
numeric/fCmplxSVD.h, numeric/fEIG.cc, numeric/fEIG.h, numeric/floatAEPBAL.cc,
numeric/floatAEPBAL.h, numeric/floatCHOL.cc, numeric/floatCHOL.h,
numeric/floatGEPBAL.cc, numeric/floatGEPBAL.h, numeric/floatHESS.cc,
numeric/floatHESS.h, numeric/floatLU.cc, numeric/floatLU.h, numeric/floatQR.cc,
numeric/floatQR.h, numeric/floatQRP.cc, numeric/floatQRP.h,
numeric/floatSCHUR.cc, numeric/floatSCHUR.h, numeric/floatSVD.cc,
numeric/floatSVD.h, numeric/lo-mappers.cc, numeric/lo-mappers.h,
numeric/lo-specfun.cc, numeric/lo-specfun.h, numeric/module.mk,
numeric/oct-convn.cc, numeric/oct-convn.h, numeric/oct-fftw.cc,
numeric/oct-fftw.h, numeric/oct-norm.cc, numeric/oct-norm.h,
numeric/oct-rand.cc, numeric/oct-rand.h, numeric/oct-spparms.cc,
numeric/oct-spparms.h, numeric/randgamma.c, numeric/randgamma.h,
numeric/randmtzig.c, numeric/randmtzig.h, numeric/randpoisson.c,
numeric/randpoisson.h, numeric/sparse-base-chol.cc, numeric/sparse-base-chol.h,
numeric/sparse-base-lu.cc, numeric/sparse-base-lu.h, numeric/sparse-dmsolve.cc:
Moved from liboctave dir to numeric subdirectory.
* operators/Sparse-diag-op-defs.h, operators/Sparse-op-defs.h,
operators/Sparse-perm-op-defs.h, operators/config-ops.sh, operators/mk-ops.awk,
operators/module.mk, operators/mx-base.h, operators/mx-defs.h,
operators/mx-ext.h, operators/mx-inlines.cc, operators/mx-op-decl.h,
operators/mx-op-defs.h, operators/mx-ops, operators/sparse-mk-ops.awk,
operators/sparse-mx-ops, operators/vx-ops:
Moved from liboctave dir to operators subdirectory.
* system/dir-ops.cc, system/dir-ops.h, system/file-ops.cc, system/file-ops.h,
system/file-stat.cc, system/file-stat.h, system/lo-sysdep.cc,
system/lo-sysdep.h, system/mach-info.cc, system/mach-info.h, system/module.mk,
system/oct-env.cc, system/oct-env.h, system/oct-group.cc, system/oct-group.h,
system/oct-openmp.h, system/oct-passwd.cc, system/oct-passwd.h,
system/oct-syscalls.cc, system/oct-syscalls.h, system/oct-time.cc,
system/oct-time.h, system/oct-uname.cc, system/oct-uname.h, system/pathlen.h,
system/sysdir.h, system/syswait.h, system/tempnam.c, system/tempname.c:
Moved from liboctave dir to system subdirectory.
* util/base-list.h, util/byte-swap.h, util/caseless-str.h, util/cmd-edit.cc,
util/cmd-edit.h, util/cmd-hist.cc, util/cmd-hist.h, util/data-conv.cc,
util/data-conv.h, util/f2c-main.c, util/functor.h, util/glob-match.cc,
util/glob-match.h, util/kpse.cc, util/lo-array-gripes.cc,
util/lo-array-gripes.h, util/lo-cieee.c, util/lo-cutils.c, util/lo-cutils.h,
util/lo-ieee.cc, util/lo-ieee.h, util/lo-macros.h, util/lo-math.h,
util/lo-traits.h, util/lo-utils.cc, util/lo-utils.h, util/module.mk,
util/oct-alloc.cc, util/oct-alloc.h, util/oct-base64.cc, util/oct-base64.h,
util/oct-binmap.h, util/oct-cmplx.h, util/oct-glob.cc, util/oct-glob.h,
util/oct-inttypes.cc, util/oct-inttypes.h, util/oct-locbuf.cc,
util/oct-locbuf.h, util/oct-md5.cc, util/oct-md5.h, util/oct-mem.h,
util/oct-mutex.cc, util/oct-mutex.h, util/oct-refcount.h, util/oct-rl-edit.c,
util/oct-rl-edit.h, util/oct-rl-hist.c, util/oct-rl-hist.h, util/oct-shlib.cc,
util/oct-shlib.h, util/oct-sort.cc, util/oct-sort.h, util/oct-sparse.h,
util/pathsearch.cc, util/pathsearch.h, util/regexp.cc, util/regexp.h,
util/singleton-cleanup.cc, util/singleton-cleanup.h, util/sparse-sort.cc,
util/sparse-sort.h, util/sparse-util.cc, util/sparse-util.h, util/statdefs.h,
util/str-vec.cc, util/str-vec.h, util/sun-utils.h:
Moved from liboctave dir to util subdirectory.
* Makefile.am: Eliminate reference to top-level liboctave directory.
* autogen.sh: cd to new liboctave/operators directory to run config-ops.sh.
* build-aux/common.mk: Eliminate LIBCRUFT references.
* configure.ac: Eliminate libcruft top-level references. Switch test
programs to find files in liboctave/cruft subdirectory.
* OctaveFAQ.texi, install.txi, mkoctfile.1: Eliminate references to libcruft in
docs.
* libgui/src/Makefile.am, libinterp/Makefile.am, src/Makefile.am: Update
include file locations. Stop linking against libcruft.
* libinterp/corefcn/module.mk: Update location of OPT_INC files which are
now in numeric/ subdirectory.
* libinterp/dldfcn/config-module.awk: Stop linking against libcruft.
* libinterp/interpfcn/toplev.cc: Remove reference to LIBCRUFT.
* libinterp/link-deps.mk, liboctave/link-deps.mk:
Add GNULIB_LINK_DEPS to link dependencies.
* libinterp/oct-conf.in.h: Remove reference to OCTAVE_CONF_LIBCRUFT.
* liboctave/Makefile.am: Overhaul to use convenience libraries in
subdirectories.
* scripts/miscellaneous/mkoctfile.m: Eliminate reference to LIBCRUFT.
* src/mkoctfile.in.cc, src/mkoctfile.in.sh: Stop linking againt libcruft.
Eliminate references to LIBCRUFT.
| author | Rik <rik@octave.org> |
|---|---|
| date | Fri, 31 Aug 2012 20:00:20 -0700 |
| parents | liboctave/randmtzig.c@3d8ace26c5b4 |
| children | c0dd3ce6d1e2 |
| rev | line source |
|---|---|
| 7019 | 1 /* |
| 2 | |
|
14138
72c96de7a403
maint: update copyright notices for 2012
John W. Eaton <jwe@octave.org>
parents:
11586
diff
changeset
|
3 Copyright (C) 2006-2012 John W. Eaton |
| 7019 | 4 |
| 5 This file is part of Octave. | |
| 6 | |
| 7 Octave is free software; you can redistribute it and/or modify it | |
| 8 under the terms of the GNU General Public License as published by the | |
| 9 Free Software Foundation; either version 3 of the License, or (at your | |
| 10 option) any later version. | |
| 11 | |
| 12 Octave is distributed in the hope that it will be useful, but WITHOUT | |
| 13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
| 14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
| 15 for more details. | |
| 16 | |
| 17 You should have received a copy of the GNU General Public License | |
| 18 along with Octave; see the file COPYING. If not, see | |
| 19 <http://www.gnu.org/licenses/>. | |
| 20 | |
| 21 */ | |
| 22 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
23 /* |
| 5742 | 24 A C-program for MT19937, with initialization improved 2002/2/10. |
| 25 Coded by Takuji Nishimura and Makoto Matsumoto. | |
| 26 This is a faster version by taking Shawn Cokus's optimization, | |
| 27 Matthe Bellew's simplification, Isaku Wada's real version. | |
| 28 David Bateman added normal and exponential distributions following | |
| 29 Marsaglia and Tang's Ziggurat algorithm. | |
| 30 | |
| 31 Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, | |
| 32 Copyright (C) 2004, David Bateman | |
| 33 All rights reserved. | |
| 34 | |
| 35 Redistribution and use in source and binary forms, with or without | |
| 36 modification, are permitted provided that the following conditions | |
| 37 are met: | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
38 |
| 5742 | 39 1. Redistributions of source code must retain the above copyright |
| 40 notice, this list of conditions and the following disclaimer. | |
| 41 | |
| 42 2. Redistributions in binary form must reproduce the above copyright | |
| 43 notice, this list of conditions and the following disclaimer in the | |
| 44 documentation and/or other materials provided with the distribution. | |
| 45 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
46 3. The names of its contributors may not be used to endorse or promote |
|
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
47 products derived from this software without specific prior written |
| 5742 | 48 permission. |
| 49 | |
| 50 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS | |
| 51 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT | |
| 52 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
53 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER |
| 5742 | 54 OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, |
| 55 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, | |
| 56 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR | |
| 57 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF | |
| 58 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING | |
| 59 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS | |
| 60 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. | |
| 61 | |
| 62 | |
| 63 Any feedback is very welcome. | |
| 64 http://www.math.keio.ac.jp/matumoto/emt.html | |
| 65 email: matumoto@math.keio.ac.jp | |
| 66 | |
| 67 * 2004-01-19 Paul Kienzle | |
| 68 * * comment out main | |
| 69 * add init_by_entropy, get_state, set_state | |
| 70 * * converted to allow compiling by C++ compiler | |
| 71 * | |
| 72 * 2004-01-25 David Bateman | |
| 73 * * Add Marsaglia and Tsang Ziggurat code | |
| 74 * | |
| 75 * 2004-07-13 Paul Kienzle | |
| 76 * * make into an independent library with some docs. | |
| 77 * * introduce new main and test code. | |
| 78 * | |
| 79 * 2004-07-28 Paul Kienzle & David Bateman | |
| 80 * * add -DALLBITS flag for 32 vs. 53 bits of randomness in mantissa | |
| 81 * * make the naming scheme more uniform | |
| 82 * * add -DHAVE_X86 for faster support of 53 bit mantissa on x86 arch. | |
| 83 * | |
| 84 * 2005-02-23 Paul Kienzle | |
| 85 * * fix -DHAVE_X86_32 flag and add -DUSE_X86_32=0|1 for explicit control | |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
86 * |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
87 * 2006-04-01 David Bateman |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
88 * * convert for use in octave, declaring static functions only used |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
89 * here and adding oct_ to functions visible externally |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
90 * * inverse sense of ALLBITS |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
91 * |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
92 * 2012-05-18 David Bateman |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
93 * * Remove randu64 and ALLBIT option |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
94 * * Add the single precision generators |
| 5742 | 95 */ |
| 96 | |
| 97 /* | |
| 98 === Build instructions === | |
| 99 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
100 Compile with -DHAVE_GETTIMEOFDAY if the gettimeofday function is |
| 5742 | 101 available. This is not necessary if your architecture has |
| 102 /dev/urandom defined. | |
| 103 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
104 Uses implicit -Di386 or explicit -DHAVE_X86_32 to determine if CPU=x86. |
|
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
105 You can force X86 behaviour with -DUSE_X86_32=1, or suppress it with |
|
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
106 -DUSE_X86_32=0. You should also consider -march=i686 or similar for |
| 5742 | 107 extra performance. Check whether -DUSE_X86_32=0 is faster on 64-bit |
| 108 x86 architectures. | |
| 109 | |
| 110 If you want to replace the Mersenne Twister with another | |
| 111 generator then redefine randi32 appropriately. | |
| 112 | |
| 113 === Usage instructions === | |
| 114 Before using any of the generators, initialize the state with one of | |
| 115 oct_init_by_int, oct_init_by_array or oct_init_by_entropy. | |
| 116 | |
| 117 All generators share the same state vector. | |
| 118 | |
| 119 === Mersenne Twister === | |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
120 void oct_init_by_int (uint32_t s) 32-bit initial state |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
121 void oct_init_by_array (uint32_t k[],int m) m*32-bit initial state |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
122 void oct_init_by_entropy (void) random initial state |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
123 void oct_get_state (uint32_t save[MT_N+1]) saves state in array |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
124 void oct_set_state (uint32_t save[MT_N+1]) restores state from array |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
125 static uint32_t randmt (void) returns 32-bit unsigned int |
| 5742 | 126 |
| 127 === inline generators === | |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
128 static uint32_t randi32 (void) returns 32-bit unsigned int |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
129 static uint64_t randi53 (void) returns 53-bit unsigned int |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
130 static uint64_t randi54 (void) returns 54-bit unsigned int |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
131 static float randu32 (void) returns 32-bit uniform in (0,1) |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
132 static double randu53 (void) returns 53-bit uniform in (0,1) |
| 5742 | 133 |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
134 double oct_randu (void) returns M-bit uniform in (0,1) |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
135 double oct_randn (void) returns M-bit standard normal |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
136 double oct_rande (void) returns N-bit standard exponential |
| 5742 | 137 |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
138 float oct_float_randu (void) returns M-bit uniform in (0,1) |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
139 float oct_float_randn (void) returns M-bit standard normal |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
140 float oct_float_rande (void) returns N-bit standard exponential |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
141 |
| 5742 | 142 === Array generators === |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
143 void oct_fill_randi32 (octave_idx_type, uint32_t []) |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
144 void oct_fill_randi64 (octave_idx_type, uint64_t []) |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
145 |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
146 void oct_fill_randu (octave_idx_type, double []) |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
147 void oct_fill_randn (octave_idx_type, double []) |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
148 void oct_fill_rande (octave_idx_type, double []) |
| 5742 | 149 |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
150 void oct_fill_float_randu (octave_idx_type, float []) |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
151 void oct_fill_float_randn (octave_idx_type, float []) |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
152 void oct_fill_float_rande (octave_idx_type, float []) |
| 5742 | 153 */ |
| 154 | |
| 155 #if defined (HAVE_CONFIG_H) | |
| 156 #include <config.h> | |
| 157 #endif | |
| 158 | |
| 159 #include <stdio.h> | |
| 160 #include <time.h> | |
| 161 | |
| 162 #ifdef HAVE_GETTIMEOFDAY | |
| 163 #include <sys/time.h> | |
| 164 #endif | |
| 165 | |
| 7231 | 166 #include "lo-math.h" |
| 5742 | 167 #include "randmtzig.h" |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
168 |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
169 /* FIXME may want to suppress X86 if sizeof(long) > 4 */ |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
170 #if !defined (USE_X86_32) |
|
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
171 # if defined (i386) || defined (HAVE_X86_32) |
| 5742 | 172 # define USE_X86_32 1 |
| 173 # else | |
| 174 # define USE_X86_32 0 | |
| 175 # endif | |
| 176 #endif | |
| 177 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
178 /* ===== Mersenne Twister 32-bit generator ===== */ |
| 5742 | 179 |
| 180 #define MT_M 397 | |
| 181 #define MATRIX_A 0x9908b0dfUL /* constant vector a */ | |
| 182 #define UMASK 0x80000000UL /* most significant w-r bits */ | |
| 183 #define LMASK 0x7fffffffUL /* least significant r bits */ | |
| 184 #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) ) | |
| 185 #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL)) | |
| 186 | |
| 187 static uint32_t *next; | |
| 188 static uint32_t state[MT_N]; /* the array for the state vector */ | |
| 189 static int left = 1; | |
| 190 static int initf = 0; | |
| 191 static int initt = 1; | |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
192 static int inittf = 1; |
| 5742 | 193 |
| 194 /* initializes state[MT_N] with a seed */ | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
195 void |
| 5742 | 196 oct_init_by_int (uint32_t s) |
| 197 { | |
| 198 int j; | |
| 199 state[0] = s & 0xffffffffUL; | |
| 200 for (j = 1; j < MT_N; j++) { | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
201 state[j] = (1812433253UL * (state[j-1] ^ (state[j-1] >> 30)) + j); |
| 5742 | 202 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ |
| 203 /* In the previous versions, MSBs of the seed affect */ | |
| 204 /* only MSBs of the array state[]. */ | |
| 205 /* 2002/01/09 modified by Makoto Matsumoto */ | |
| 206 state[j] &= 0xffffffffUL; /* for >32 bit machines */ | |
| 207 } | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
208 left = 1; |
| 5742 | 209 initf = 1; |
| 210 } | |
| 211 | |
| 212 /* initialize by an array with array-length */ | |
| 213 /* init_key is the array for initializing keys */ | |
| 214 /* key_length is its length */ | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
215 void |
|
7533
ff52243af934
save state separately for each MT random number generator
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
216 oct_init_by_array (uint32_t *init_key, int key_length) |
| 5742 | 217 { |
| 218 int i, j, k; | |
| 219 oct_init_by_int (19650218UL); | |
| 220 i = 1; | |
| 221 j = 0; | |
| 222 k = (MT_N > key_length ? MT_N : key_length); | |
| 223 for (; k; k--) | |
| 224 { | |
| 225 state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1664525UL)) | |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
226 + init_key[j] + j; /* non linear */ |
| 5742 | 227 state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ |
| 228 i++; | |
| 229 j++; | |
| 230 if (i >= MT_N) | |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
231 { |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
232 state[0] = state[MT_N-1]; |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
233 i = 1; |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
234 } |
| 5742 | 235 if (j >= key_length) |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
236 j = 0; |
| 5742 | 237 } |
| 238 for (k = MT_N - 1; k; k--) | |
| 239 { | |
| 240 state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL)) | |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
241 - i; /* non linear */ |
| 5742 | 242 state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ |
| 243 i++; | |
| 244 if (i >= MT_N) | |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
245 { |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
246 state[0] = state[MT_N-1]; |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
247 i = 1; |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
248 } |
| 5742 | 249 } |
| 250 | |
| 251 state[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ | |
| 252 left = 1; | |
| 253 initf = 1; | |
| 254 } | |
| 255 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
256 void |
| 5742 | 257 oct_init_by_entropy (void) |
| 258 { | |
| 259 uint32_t entropy[MT_N]; | |
| 260 int n = 0; | |
| 261 | |
| 262 /* Look for entropy in /dev/urandom */ | |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
263 FILE* urandom =fopen ("/dev/urandom", "rb"); |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
264 if (urandom) |
| 5742 | 265 { |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
266 while (n < MT_N) |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
267 { |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
268 unsigned char word[4]; |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
269 if (fread (word, 4, 1, urandom) != 1) |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
270 break; |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
271 entropy[n++] = word[0]+(word[1]<<8)+(word[2]<<16)+(word[3]<<24); |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
272 } |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
273 fclose (urandom); |
| 5742 | 274 } |
| 275 | |
| 276 /* If there isn't enough entropy, gather some from various sources */ | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
277 if (n < MT_N) |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
278 entropy[n++] = time (NULL); /* Current time in seconds */ |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
279 if (n < MT_N) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
280 entropy[n++] = clock (); /* CPU time used (usec) */ |
| 5742 | 281 #ifdef HAVE_GETTIMEOFDAY |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
282 if (n < MT_N) |
| 5742 | 283 { |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
284 struct timeval tv; |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
285 if (gettimeofday (&tv, NULL) != -1) |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
286 entropy[n++] = tv.tv_usec; /* Fractional part of current time */ |
| 5742 | 287 } |
| 288 #endif | |
| 289 /* Send all the entropy into the initial state vector */ | |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
290 oct_init_by_array (entropy,n); |
| 5742 | 291 } |
| 292 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
293 void |
|
7533
ff52243af934
save state separately for each MT random number generator
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
294 oct_set_state (uint32_t *save) |
| 5742 | 295 { |
| 296 int i; | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
297 for (i = 0; i < MT_N; i++) |
| 5742 | 298 state[i] = save[i]; |
| 299 left = save[MT_N]; | |
| 300 next = state + (MT_N - left + 1); | |
| 301 } | |
| 302 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
303 void |
|
7533
ff52243af934
save state separately for each MT random number generator
John W. Eaton <jwe@octave.org>
parents:
7231
diff
changeset
|
304 oct_get_state (uint32_t *save) |
| 5742 | 305 { |
| 306 int i; | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
307 for (i = 0; i < MT_N; i++) |
| 5742 | 308 save[i] = state[i]; |
| 309 save[MT_N] = left; | |
| 310 } | |
| 311 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
312 static void |
| 5742 | 313 next_state (void) |
| 314 { | |
| 315 uint32_t *p = state; | |
| 316 int j; | |
| 317 | |
| 318 /* if init_by_int() has not been called, */ | |
| 319 /* a default initial seed is used */ | |
| 320 /* if (initf==0) init_by_int(5489UL); */ | |
| 321 /* Or better yet, a random seed! */ | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
322 if (initf == 0) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
323 oct_init_by_entropy (); |
| 5742 | 324 |
| 325 left = MT_N; | |
| 326 next = state; | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
327 |
|
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
328 for (j = MT_N - MT_M + 1; --j; p++) |
| 5742 | 329 *p = p[MT_M] ^ TWIST(p[0], p[1]); |
| 330 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
331 for (j = MT_M; --j; p++) |
| 5742 | 332 *p = p[MT_M-MT_N] ^ TWIST(p[0], p[1]); |
| 333 | |
| 334 *p = p[MT_M-MT_N] ^ TWIST(p[0], state[0]); | |
| 335 } | |
| 336 | |
| 337 /* generates a random number on [0,0xffffffff]-interval */ | |
| 5766 | 338 static uint32_t |
| 5742 | 339 randmt (void) |
| 340 { | |
| 341 register uint32_t y; | |
| 342 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
343 if (--left == 0) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
344 next_state (); |
| 5742 | 345 y = *next++; |
| 346 | |
| 347 /* Tempering */ | |
| 348 y ^= (y >> 11); | |
| 349 y ^= (y << 7) & 0x9d2c5680UL; | |
| 350 y ^= (y << 15) & 0xefc60000UL; | |
| 351 return (y ^ (y >> 18)); | |
| 352 } | |
| 353 | |
| 354 /* ===== Uniform generators ===== */ | |
| 355 | |
| 356 /* Select which 32 bit generator to use */ | |
| 357 #define randi32 randmt | |
| 358 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
359 static uint64_t |
| 5742 | 360 randi53 (void) |
| 361 { | |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
362 const uint32_t lo = randi32 (); |
|
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
363 const uint32_t hi = randi32 ()&0x1FFFFF; |
| 5742 | 364 #if HAVE_X86_32 |
| 365 uint64_t u; | |
| 366 uint32_t *p = (uint32_t *)&u; | |
| 367 p[0] = lo; | |
| 368 p[1] = hi; | |
| 369 return u; | |
| 370 #else | |
| 371 return (((uint64_t)hi<<32)|lo); | |
| 372 #endif | |
| 373 } | |
| 374 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
375 static uint64_t |
| 5742 | 376 randi54 (void) |
| 377 { | |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
378 const uint32_t lo = randi32 (); |
|
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
379 const uint32_t hi = randi32 ()&0x3FFFFF; |
| 5742 | 380 #if HAVE_X86_32 |
| 381 uint64_t u; | |
| 382 uint32_t *p = (uint32_t *)&u; | |
| 383 p[0] = lo; | |
| 384 p[1] = hi; | |
| 385 return u; | |
| 386 #else | |
| 387 return (((uint64_t)hi<<32)|lo); | |
| 388 #endif | |
| 389 } | |
| 390 | |
| 391 /* generates a random number on (0,1)-real-interval */ | |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
392 static float |
| 5742 | 393 randu32 (void) |
| 394 { | |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
395 return ((float)randi32 () + 0.5) * (1.0/4294967296.0); |
| 5742 | 396 /* divided by 2^32 */ |
| 397 } | |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
398 |
| 5742 | 399 /* generates a random number on (0,1) with 53-bit resolution */ |
| 5766 | 400 static double |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
401 randu53 (void) |
|
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
402 { |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
403 const uint32_t a=randi32 ()>>5; |
|
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
404 const uint32_t b=randi32 ()>>6; |
| 6959 | 405 return (a*67108864.0+b+0.4) * (1.0/9007199254740992.0); |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
406 } |
| 5742 | 407 |
| 408 /* Determine mantissa for uniform doubles */ | |
| 409 double | |
| 410 oct_randu (void) | |
| 411 { | |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
412 return randu53 (); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
413 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
414 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
415 /* Determine mantissa for uniform floats */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
416 float |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
417 oct_float_randu (void) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
418 { |
| 6959 | 419 return randu32 (); |
| 5742 | 420 } |
| 421 | |
| 422 /* ===== Ziggurat normal and exponential generators ===== */ | |
| 423 | |
| 424 #define ZIGGURAT_TABLE_SIZE 256 | |
| 425 | |
| 426 #define ZIGGURAT_NOR_R 3.6541528853610088 | |
| 427 #define ZIGGURAT_NOR_INV_R 0.27366123732975828 | |
| 428 #define NOR_SECTION_AREA 0.00492867323399 | |
| 429 | |
| 430 #define ZIGGURAT_EXP_R 7.69711747013104972 | |
| 431 #define ZIGGURAT_EXP_INV_R 0.129918765548341586 | |
| 432 #define EXP_SECTION_AREA 0.0039496598225815571993 | |
| 433 | |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
434 #define ZIGINT uint64_t |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
435 #define EMANTISSA 9007199254740992.0 /* 53 bit mantissa */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
436 #define ERANDI randi53() /* 53 bits for mantissa */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
437 #define NMANTISSA EMANTISSA |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
438 #define NRANDI randi54() /* 53 bits for mantissa + 1 bit sign */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
439 #define RANDU randu53() |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
440 |
| 5742 | 441 static ZIGINT ki[ZIGGURAT_TABLE_SIZE]; |
| 442 static double wi[ZIGGURAT_TABLE_SIZE], fi[ZIGGURAT_TABLE_SIZE]; | |
| 443 static ZIGINT ke[ZIGGURAT_TABLE_SIZE]; | |
| 444 static double we[ZIGGURAT_TABLE_SIZE], fe[ZIGGURAT_TABLE_SIZE]; | |
| 445 | |
| 446 /* | |
| 447 This code is based on the paper Marsaglia and Tsang, "The ziggurat method | |
| 448 for generating random variables", Journ. Statistical Software. Code was | |
| 449 presented in this paper for a Ziggurat of 127 levels and using a 32 bit | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
450 integer random number generator. This version of the code, uses the |
|
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
451 Mersenne Twister as the integer generator and uses 256 levels in the |
| 5742 | 452 Ziggurat. This has several advantages. |
| 453 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
454 1) As Marsaglia and Tsang themselves states, the more levels the few |
| 5742 | 455 times the expensive tail algorithm must be called |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
456 2) The cycle time of the generator is determined by the integer |
|
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
457 generator, thus the use of a Mersenne Twister for the core random |
| 5742 | 458 generator makes this cycle extremely long. |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
459 3) The license on the original code was unclear, thus rewriting the code |
| 5742 | 460 from the article means we are free of copyright issues. |
| 461 4) Compile flag for full 53-bit random mantissa. | |
| 462 | |
| 463 It should be stated that the authors made my life easier, by the fact that | |
| 464 the algorithm developed in the text of the article is for a 256 level | |
| 465 ziggurat, even if the code itself isn't... | |
| 466 | |
| 467 One modification to the algorithm developed in the article, is that it is | |
| 468 assumed that 0 <= x < Inf, and "unsigned long"s are used, thus resulting in | |
| 469 terms like 2^32 in the code. As the normal distribution is defined between | |
| 470 -Inf < x < Inf, we effectively only have 31 bit integers plus a sign. Thus | |
| 471 in Marsaglia and Tsang, terms like 2^32 become 2^31. We use NMANTISSA for | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
472 this term. The exponential distribution is one sided so we use the |
| 5742 | 473 full 32 bits. We use EMANTISSA for this term. |
| 474 | |
| 475 It appears that I'm slightly slower than the code in the article, this | |
| 476 is partially due to a better generator of random integers than they | |
| 477 use. But might also be that the case of rapid return was optimized by | |
| 478 inlining the relevant code with a #define. As the basic Mersenne | |
| 479 Twister is only 25% faster than this code I suspect that the main | |
| 480 reason is just the use of the Mersenne Twister and not the inlining, | |
| 481 so I'm not going to try and optimize further. | |
| 482 */ | |
| 483 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
484 static void |
| 5742 | 485 create_ziggurat_tables (void) |
| 486 { | |
| 487 int i; | |
| 488 double x, x1; | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
489 |
| 5742 | 490 /* Ziggurat tables for the normal distribution */ |
| 491 x1 = ZIGGURAT_NOR_R; | |
| 492 wi[255] = x1 / NMANTISSA; | |
| 493 fi[255] = exp (-0.5 * x1 * x1); | |
| 494 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
495 /* Index zero is special for tail strip, where Marsaglia and Tsang |
|
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
496 * defines this as |
| 5742 | 497 * k_0 = 2^31 * r * f(r) / v, w_0 = 0.5^31 * v / f(r), f_0 = 1, |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
498 * where v is the area of each strip of the ziggurat. |
| 5742 | 499 */ |
| 500 ki[0] = (ZIGINT) (x1 * fi[255] / NOR_SECTION_AREA * NMANTISSA); | |
| 501 wi[0] = NOR_SECTION_AREA / fi[255] / NMANTISSA; | |
| 502 fi[0] = 1.; | |
| 503 | |
| 504 for (i = 254; i > 0; i--) | |
| 505 { | |
| 506 /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus | |
| 507 * need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y)) | |
| 508 */ | |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
509 x = sqrt (-2. * log (NOR_SECTION_AREA / x1 + fi[i+1])); |
| 5742 | 510 ki[i+1] = (ZIGINT)(x / x1 * NMANTISSA); |
| 511 wi[i] = x / NMANTISSA; | |
| 512 fi[i] = exp (-0.5 * x * x); | |
| 513 x1 = x; | |
| 514 } | |
| 515 | |
| 516 ki[1] = 0; | |
| 517 | |
| 518 /* Zigurrat tables for the exponential distribution */ | |
| 519 x1 = ZIGGURAT_EXP_R; | |
| 520 we[255] = x1 / EMANTISSA; | |
| 521 fe[255] = exp (-x1); | |
| 522 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
523 /* Index zero is special for tail strip, where Marsaglia and Tsang |
|
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
524 * defines this as |
| 5742 | 525 * k_0 = 2^32 * r * f(r) / v, w_0 = 0.5^32 * v / f(r), f_0 = 1, |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
526 * where v is the area of each strip of the ziggurat. |
| 5742 | 527 */ |
| 528 ke[0] = (ZIGINT) (x1 * fe[255] / EXP_SECTION_AREA * EMANTISSA); | |
| 529 we[0] = EXP_SECTION_AREA / fe[255] / EMANTISSA; | |
| 530 fe[0] = 1.; | |
| 531 | |
| 532 for (i = 254; i > 0; i--) | |
| 533 { | |
| 534 /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus | |
| 535 * need inverse operator of y = exp(-x) -> x = -ln(y) | |
| 536 */ | |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
537 x = - log (EXP_SECTION_AREA / x1 + fe[i+1]); |
| 5742 | 538 ke[i+1] = (ZIGINT)(x / x1 * EMANTISSA); |
| 539 we[i] = x / EMANTISSA; | |
| 540 fe[i] = exp (-x); | |
| 541 x1 = x; | |
| 542 } | |
| 543 ke[1] = 0; | |
| 544 | |
| 545 initt = 0; | |
| 546 } | |
| 547 | |
| 548 /* | |
| 549 * Here is the guts of the algorithm. As Marsaglia and Tsang state the | |
| 550 * algorithm in their paper | |
| 551 * | |
| 552 * 1) Calculate a random signed integer j and let i be the index | |
| 553 * provided by the rightmost 8-bits of j | |
| 554 * 2) Set x = j * w_i. If j < k_i return x | |
| 555 * 3) If i = 0, then return x from the tail | |
| 556 * 4) If [f(x_{i-1}) - f(x_i)] * U < f(x) - f(x_i), return x | |
| 557 * 5) goto step 1 | |
| 558 * | |
| 559 * Where f is the functional form of the distribution, which for a normal | |
| 560 * distribution is exp(-0.5*x*x) | |
| 561 */ | |
| 562 | |
| 563 double | |
| 564 oct_randn (void) | |
| 565 { | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
566 if (initt) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
567 create_ziggurat_tables (); |
| 5742 | 568 |
| 569 while (1) | |
| 570 { | |
| 571 /* The following code is specialized for 32-bit mantissa. | |
| 572 * Compared to the arbitrary mantissa code, there is a performance | |
| 573 * gain for 32-bits: PPC: 2%, MIPS: 8%, x86: 40% | |
| 574 * There is a bigger performance gain compared to using a full | |
| 575 * 53-bit mantissa: PPC: 60%, MIPS: 65%, x86: 240% | |
| 576 * Of course, different compilers and operating systems may | |
| 577 * have something to do with this. | |
| 578 */ | |
| 579 # if HAVE_X86_32 | |
| 580 /* 53-bit mantissa, 1-bit sign, x86 32-bit architecture */ | |
| 581 double x; | |
| 582 int si,idx; | |
| 583 register uint32_t lo, hi; | |
| 584 int64_t rabs; | |
| 585 uint32_t *p = (uint32_t *)&rabs; | |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
586 lo = randi32 (); |
| 5742 | 587 idx = lo&0xFF; |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
588 hi = randi32 (); |
| 5742 | 589 si = hi&UMASK; |
| 590 p[0] = lo; | |
| 591 p[1] = hi&0x1FFFFF; | |
| 592 x = ( si ? -rabs : rabs ) * wi[idx]; | |
| 593 # else /* !HAVE_X86_32 */ | |
| 594 /* arbitrary mantissa (selected by NRANDI, with 1 bit for sign) */ | |
| 595 const uint64_t r = NRANDI; | |
| 596 const int64_t rabs=r>>1; | |
| 597 const int idx = (int)(rabs&0xFF); | |
| 598 const double x = ( r&1 ? -rabs : rabs) * wi[idx]; | |
| 599 # endif /* !HAVE_X86_32 */ | |
| 600 if (rabs < (int64_t)ki[idx]) | |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
601 return x; /* 99.3% of the time we return here 1st try */ |
| 5742 | 602 else if (idx == 0) |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
603 { |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
604 /* As stated in Marsaglia and Tsang |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
605 * |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
606 * For the normal tail, the method of Marsaglia[5] provides: |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
607 * generate x = -ln(U_1)/r, y = -ln(U_2), until y+y > x*x, |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
608 * then return r+x. Except that r+x is always in the positive |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
609 * tail!!!! Any thing random might be used to determine the |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
610 * sign, but as we already have r we might as well use it |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
611 * |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
612 * [PAK] but not the bottom 8 bits, since they are all 0 here! |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
613 */ |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
614 double xx, yy; |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
615 do |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
616 { |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
617 xx = - ZIGGURAT_NOR_INV_R * log (RANDU); |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
618 yy = - log (RANDU); |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
619 } |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
620 while ( yy+yy <= xx*xx); |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
621 return (rabs&0x100 ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx); |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
622 } |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
623 else if ((fi[idx-1] - fi[idx]) * RANDU + fi[idx] < exp (-0.5*x*x)) |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
624 return x; |
| 5742 | 625 } |
| 626 } | |
| 627 | |
| 628 double | |
| 629 oct_rande (void) | |
| 630 { | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
631 if (initt) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
632 create_ziggurat_tables (); |
| 5742 | 633 |
| 634 while (1) | |
| 635 { | |
| 636 ZIGINT ri = ERANDI; | |
| 637 const int idx = (int)(ri & 0xFF); | |
| 638 const double x = ri * we[idx]; | |
| 639 if (ri < ke[idx]) | |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
640 return x; // 98.9% of the time we return here 1st try |
| 5742 | 641 else if (idx == 0) |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
642 { |
|
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
643 /* As stated in Marsaglia and Tsang |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
644 * |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
645 * For the exponential tail, the method of Marsaglia[5] provides: |
| 5742 | 646 * x = r - ln(U); |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
647 */ |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
648 return ZIGGURAT_EXP_R - log (RANDU); |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
649 } |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
650 else if ((fe[idx-1] - fe[idx]) * RANDU + fe[idx] < exp (-x)) |
|
10317
42d098307c30
untabify additional source files
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
651 return x; |
| 5742 | 652 } |
| 653 } | |
| 654 | |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
655 #undef ZIGINT |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
656 #undef EMANTISSA |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
657 #undef ERANDI |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
658 #undef NMANTISSA |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
659 #undef NRANDI |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
660 #undef RANDU |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
661 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
662 #define ZIGINT uint32_t |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
663 #define EMANTISSA 4294967296.0 /* 32 bit mantissa */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
664 #define ERANDI randi32() /* 32 bits for mantissa */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
665 #define NMANTISSA 2147483648.0 /* 31 bit mantissa */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
666 #define NRANDI randi32() /* 31 bits for mantissa + 1 bit sign */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
667 #define RANDU randu32() |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
668 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
669 static ZIGINT fki[ZIGGURAT_TABLE_SIZE]; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
670 static float fwi[ZIGGURAT_TABLE_SIZE], ffi[ZIGGURAT_TABLE_SIZE]; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
671 static ZIGINT fke[ZIGGURAT_TABLE_SIZE]; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
672 static float fwe[ZIGGURAT_TABLE_SIZE], ffe[ZIGGURAT_TABLE_SIZE]; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
673 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
674 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
675 static void |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
676 create_ziggurat_float_tables (void) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
677 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
678 int i; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
679 float x, x1; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
680 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
681 /* Ziggurat tables for the normal distribution */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
682 x1 = ZIGGURAT_NOR_R; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
683 fwi[255] = x1 / NMANTISSA; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
684 ffi[255] = exp (-0.5 * x1 * x1); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
685 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
686 /* Index zero is special for tail strip, where Marsaglia and Tsang |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
687 * defines this as |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
688 * k_0 = 2^31 * r * f(r) / v, w_0 = 0.5^31 * v / f(r), f_0 = 1, |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
689 * where v is the area of each strip of the ziggurat. |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
690 */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
691 fki[0] = (ZIGINT) (x1 * ffi[255] / NOR_SECTION_AREA * NMANTISSA); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
692 fwi[0] = NOR_SECTION_AREA / ffi[255] / NMANTISSA; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
693 ffi[0] = 1.; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
694 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
695 for (i = 254; i > 0; i--) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
696 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
697 /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
698 * need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y)) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
699 */ |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
700 x = sqrt (-2. * log (NOR_SECTION_AREA / x1 + ffi[i+1])); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
701 fki[i+1] = (ZIGINT)(x / x1 * NMANTISSA); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
702 fwi[i] = x / NMANTISSA; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
703 ffi[i] = exp (-0.5 * x * x); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
704 x1 = x; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
705 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
706 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
707 fki[1] = 0; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
708 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
709 /* Zigurrat tables for the exponential distribution */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
710 x1 = ZIGGURAT_EXP_R; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
711 fwe[255] = x1 / EMANTISSA; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
712 ffe[255] = exp (-x1); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
713 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
714 /* Index zero is special for tail strip, where Marsaglia and Tsang |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
715 * defines this as |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
716 * k_0 = 2^32 * r * f(r) / v, w_0 = 0.5^32 * v / f(r), f_0 = 1, |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
717 * where v is the area of each strip of the ziggurat. |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
718 */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
719 fke[0] = (ZIGINT) (x1 * ffe[255] / EXP_SECTION_AREA * EMANTISSA); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
720 fwe[0] = EXP_SECTION_AREA / ffe[255] / EMANTISSA; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
721 ffe[0] = 1.; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
722 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
723 for (i = 254; i > 0; i--) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
724 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
725 /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
726 * need inverse operator of y = exp(-x) -> x = -ln(y) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
727 */ |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
728 x = - log (EXP_SECTION_AREA / x1 + ffe[i+1]); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
729 fke[i+1] = (ZIGINT)(x / x1 * EMANTISSA); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
730 fwe[i] = x / EMANTISSA; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
731 ffe[i] = exp (-x); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
732 x1 = x; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
733 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
734 fke[1] = 0; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
735 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
736 inittf = 0; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
737 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
738 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
739 /* |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
740 * Here is the guts of the algorithm. As Marsaglia and Tsang state the |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
741 * algorithm in their paper |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
742 * |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
743 * 1) Calculate a random signed integer j and let i be the index |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
744 * provided by the rightmost 8-bits of j |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
745 * 2) Set x = j * w_i. If j < k_i return x |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
746 * 3) If i = 0, then return x from the tail |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
747 * 4) If [f(x_{i-1}) - f(x_i)] * U < f(x) - f(x_i), return x |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
748 * 5) goto step 1 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
749 * |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
750 * Where f is the functional form of the distribution, which for a normal |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
751 * distribution is exp(-0.5*x*x) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
752 */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
753 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
754 float |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
755 oct_float_randn (void) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
756 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
757 if (inittf) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
758 create_ziggurat_float_tables (); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
759 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
760 while (1) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
761 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
762 /* 32-bit mantissa */ |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
763 const uint32_t r = randi32 (); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
764 const uint32_t rabs = r&LMASK; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
765 const int idx = (int)(r&0xFF); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
766 const float x = ((int32_t)r) * fwi[idx]; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
767 if (rabs < fki[idx]) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
768 return x; /* 99.3% of the time we return here 1st try */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
769 else if (idx == 0) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
770 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
771 /* As stated in Marsaglia and Tsang |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
772 * |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
773 * For the normal tail, the method of Marsaglia[5] provides: |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
774 * generate x = -ln(U_1)/r, y = -ln(U_2), until y+y > x*x, |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
775 * then return r+x. Except that r+x is always in the positive |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
776 * tail!!!! Any thing random might be used to determine the |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
777 * sign, but as we already have r we might as well use it |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
778 * |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
779 * [PAK] but not the bottom 8 bits, since they are all 0 here! |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
780 */ |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
781 float xx, yy; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
782 do |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
783 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
784 xx = - ZIGGURAT_NOR_INV_R * log (RANDU); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
785 yy = - log (RANDU); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
786 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
787 while ( yy+yy <= xx*xx); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
788 return (rabs&0x100 ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
789 } |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
790 else if ((ffi[idx-1] - ffi[idx]) * RANDU + ffi[idx] < exp (-0.5*x*x)) |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
791 return x; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
792 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
793 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
794 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
795 float |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
796 oct_float_rande (void) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
797 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
798 if (inittf) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
799 create_ziggurat_float_tables (); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
800 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
801 while (1) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
802 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
803 ZIGINT ri = ERANDI; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
804 const int idx = (int)(ri & 0xFF); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
805 const float x = ri * fwe[idx]; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
806 if (ri < fke[idx]) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
807 return x; // 98.9% of the time we return here 1st try |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
808 else if (idx == 0) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
809 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
810 /* As stated in Marsaglia and Tsang |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
811 * |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
812 * For the exponential tail, the method of Marsaglia[5] provides: |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
813 * x = r - ln(U); |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
814 */ |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
815 return ZIGGURAT_EXP_R - log (RANDU); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
816 } |
|
15018
3d8ace26c5b4
maint: Use Octave coding conventions for cuddled parentheses in liboctave/.
Rik <rik@octave.org>
parents:
14846
diff
changeset
|
817 else if ((ffe[idx-1] - ffe[idx]) * RANDU + ffe[idx] < exp (-x)) |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
818 return x; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
819 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
820 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
821 |
| 5742 | 822 /* Array generators */ |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
823 void |
| 5742 | 824 oct_fill_randu (octave_idx_type n, double *p) |
| 825 { | |
| 826 octave_idx_type i; | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
827 for (i = 0; i < n; i++) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
828 p[i] = oct_randu (); |
| 5742 | 829 } |
| 830 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
831 void |
| 5742 | 832 oct_fill_randn (octave_idx_type n, double *p) |
| 833 { | |
| 834 octave_idx_type i; | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
835 for (i = 0; i < n; i++) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
836 p[i] = oct_randn (); |
| 5742 | 837 } |
| 838 | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
839 void |
| 5742 | 840 oct_fill_rande (octave_idx_type n, double *p) |
| 841 { | |
| 842 octave_idx_type i; | |
|
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
843 for (i = 0; i < n; i++) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
844 p[i] = oct_rande (); |
| 5742 | 845 } |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
846 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
847 void |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
848 oct_fill_float_randu (octave_idx_type n, float *p) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
849 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
850 octave_idx_type i; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
851 for (i = 0; i < n; i++) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
852 p[i] = oct_float_randu (); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
853 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
854 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
855 void |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
856 oct_fill_float_randn (octave_idx_type n, float *p) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
857 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
858 octave_idx_type i; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
859 for (i = 0; i < n; i++) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
860 p[i] = oct_float_randn (); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
861 } |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
862 |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
863 void |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
864 oct_fill_float_rande (octave_idx_type n, float *p) |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
865 { |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
866 octave_idx_type i; |
|
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
867 for (i = 0; i < n; i++) |
|
14846
460a3c6d8bf1
maint: Use Octave coding convention for cuddled parenthis in function calls with empty argument lists.
Rik <octave@nomad.inbox5.com>
parents:
14655
diff
changeset
|
868 p[i] = oct_float_rande (); |
|
14655
43db83eff9db
Implement single precision rand, randn, rande, randg and randp generators (bug #34351, #36293)
David Bateman <dbateman@free.fr>
parents:
14138
diff
changeset
|
869 } |