home *** CD-ROM | disk | FTP | other *** search
/ ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng) / ProfitPress-MegaCDROM2.B6I / PROG / PASCAL / NRPAS13.ZIP / BNLDEV.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1991-04-29  |  1.4 KB  |  57 lines
  1. FUNCTION bnldev(pp: real; n: integer; VAR idum: integer): real;
  2. LABEL 1;
  3. CONST
  4.    pi=3.141592654;
  5. VAR
  6.    am,em,en,g,angle: real;
  7.    oldg,p,pc,bnl: real;
  8.    pclog,plog,pold,sq,t,y: real;
  9.    j,nold: integer;
  10. BEGIN
  11.    nold := -1;
  12.    pold := -1.0;
  13.    IF (pp <= 0.5) THEN p := pp ELSE p := 1.0-pp;
  14.    am := n*p;
  15.    IF (n < 25) THEN BEGIN
  16.       bnl := 0.0;
  17.       FOR j := 1 TO n DO BEGIN
  18.          IF (ran3(idum) < p) THEN bnl := bnl+1.0
  19.       END
  20.    END ELSE IF (am < 1.0) THEN BEGIN
  21.       g := exp(-am);
  22.       t := 1.0;
  23.       FOR j := 0 TO n DO BEGIN
  24.          t := t*ran3(idum);
  25.          IF (t < g) THEN GOTO 1
  26.       END;
  27.       j := n;
  28. 1:      bnl := j
  29.    END ELSE BEGIN
  30.       IF (n <> nold) THEN BEGIN
  31.          en := n;
  32.          oldg := gammln(en+1.0);
  33.          nold := n
  34.       END;
  35.       IF (p <> pold) THEN BEGIN
  36.          pc := 1.0-p;
  37.          plog := ln(p);
  38.          pclog := ln(pc);
  39.          pold := p
  40.       END;
  41.       sq := sqrt(2.0*am*pc);
  42.       REPEAT
  43.          REPEAT
  44.             angle := pi*ran3(idum);
  45.             y := sin(angle)/cos(angle);
  46.             em := sq*y+am;
  47.          UNTIL ((em >= 0.0) AND (em < en+1.0));
  48.          em := trunc(em);
  49.          t := 1.2*sq*(1.0+sqr(y))*exp(oldg-gammln(em+1.0)
  50.             -gammln(en-em+1.0)+em*plog+(en-em)*pclog);
  51.       UNTIL (ran3(idum) <= t);
  52.       bnl := em
  53.    END;
  54.    IF (p <> pp) THEN bnl := n-bnl;
  55.    bnldev := bnl
  56. END;
  57.