Subversion Repositories wimsdev

Rev

Rev 13586 | Rev 18195 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
20 reyssat 1
!if $wims_read_parm!=slib_header
13586 bpr 2
  !goto proc
20 reyssat 3
!endif
4161 bpr 4
 
20 reyssat 5
slib_parms=3\
6
 ,m\
7
 ,n\
8
 ,p
4330 bpr 9
slib_author=Adeline, Grelot
4161 bpr 10
 
18193 bpr 11
slib_example=5,10,0.5\
12
5,10,1\
13
5,11,1.5
14
 
3265 bpr 15
slib_require=pari
3307 bpr 16
!exit
17
 
20 reyssat 18
:proc
19
slib_out=
20
 
21
!distribute item $wims_read_parm into slib_m, slib_n, slib_p
18193 bpr 22
!if $slib_p>1 or $slib_p<0
23
  slib_out=
24
  !exit
25
!endif
26
!if $slib_p=0
27
  slib_out=!makelist 0 for x=1 to $slib_m
28
  !exit
29
!endif
30
!if $slib_p=1
31
  slib_out=!makelist $slib_n for x=1 to $slib_m
32
  !exit
33
!endif
20 reyssat 34
slib_u=!random 0,1 repeat $slib_m
35
slib_C=!values x for x=0 to $[$slib_n]
36
slib_out=!exec pari {slib_g(p,x,m)=local(L,u,q,i,k) ;L=listcreate(m);\
18193 bpr 37
  for(k=1,m,u=[$slib_u][k];q=p[1];i=1;\
38
  if(u<q,listput(L,x[1]),while(u>=q,i=i+1;q+=p[i]);listput(L,x[i])));L};\
39
  {slib_bino(n,p,m)=local(q,b,prob);b=listcreate(n+1);\
40
  q=p/(1-p);prob=(1-p)^n;listput(b,prob);\
41
  for(i=1,n,prob=prob*(n-i+1)/i*q;listput(b,prob));\
42
  Vec(slib_g(b,[$slib_C],m))};\
43
  slib_bino($slib_n,$slib_p,$slib_m)
20 reyssat 44
 
45
slib_out=!declosing $slib_out
46
slib_out=!trim $slib_out