DataMuseum.dk

top - metrics - download

    Length: 2304 (0x900)
    Types: TextFile
    Names: »extstudent«

└─⟦621cfb9a2⟧ Bits:30002817 RC8000 Dump tape fra HCØ.  Detaljer om "HC8000" projekt.
    └─⟦0364f57e3⟧ 
        └─⟦58ca399f1⟧ »extbib« 
            └─⟦this⟧ 

TextFile

student=algol list.yes index.no
external
real procedure student(t1,t2,m,eps);
value t1,t2,m,eps; integer m; real t1,t2,eps;
begin
real t,pt,dpt,norm;
real procedure fract(n);
  value n; integer n;
  begin integer i; real f;
  comment calculates gamma(n/2)/gamma(n/2-1/2);
  f:= 1; i:= n+2;
  for i:=i-2 while i>1 do f:= f*i/(i-1);
  fract:= if i=0 then f/sqrt(pi) else f*sqrt(pi)/2
end fract;
real procedure f(t,m);
  value t,m; integer m; real t;
  f:= (1+t*t/m)**(-(m+1)/2);
real procedure integr(f,t1,t2,m,eps);
  value t1,t2,m,eps; integer m; real t1,t2,eps; real procedure f;
  begin integer i,n; real p,q,s,x,d;
  q:= (f(t1,m)+f(t2,m))/2;
  s:= q/2; n:= 1; d:= t2-t1;
  for p:=0 while abs(q-s) > eps do begin
    s:= q; x:= t1+d/2;
    for i:=1 step 1 until n do begin
      p:=p+f(x,m); x:= x+d
    end;
    q:= (s+p/n)/2; n:= n+n; d:= d/2;
  end;
  integr:= s*(t2-t1)
end integration;
norm:= fract(m-1)/sqrt(pi*m);
if t2=0 then begin
  t:=t1; pt:= 0; dpt:= 1;
  for t:= t+2 while dpt*norm > eps do begin
    dpt:= integr(f,t-2,t,m,eps/norm); pt:= pt+dpt end;
end else pt:= integr(f,t1,t2,m,eps/norm);
student:= pt*norm
end student-distribution;
end

testst=algol list.yes
begin
integer m; real t1,t2,eps,s;
rep:
read(in,t1,t2,m,eps);
if m>0 then begin
  s:= student(t1,t2,m,eps);
  write(out,<:<10>:>,<<ddd.dd>,t1,t2,<<dddd>,m,
        <<   d'-d>,eps,<<-ddd.dddddd>,s);
  goto rep
end end;
testst
0 1 20 '-4
0 2 20 '-4
0 5 20 '-4
5 0 20 '-4
0 1 10 '-4
0 2 10 '-4
0 5 10 '-4
5 0 10 '-4
0 1 9 '-4
0 2 9 '-4
0 5 9 '-4
5 0 9 '-4
0 0 0 0
▶EOF◀