code 41533;
real procedure NCSTUDENT(X, DF, DELTA);
value X, DF, DELTA; real X, DF, DELTA;
begin real A, B, A2, WB, D2, TOL, TOLI, H, HELP, RESULT;
    Boolean DFEVEN;
    
      real procedure INTEGRATE(Y0, Y4, F0, F2, F4);
      value Y0, Y4, F0, F2, F4; real Y0, Y4, F0, F2, F4;
      begin real F1, F3, Y2, TEE, Y;
          Y2:=(Y0 + Y4)/2;
          Y :=(Y0 + Y2)/2; F1:= EXP(H×(1 + Y×Y))/(1 + Y×Y);
          Y :=(Y2 + Y4)/2; F3:= EXP(H×(1 + Y×Y))/(1 + Y×Y);
          TEE:=6×F2 - 4×(F1 + F3) + F0 + F4;
          INTEGRATE:=if ABS(TEE) < TOLI
               then (Y4 - Y0)×(4×(F1 + F3) + 2×F2 +
                      F0 + F4 - TEE/15)
               else INTEGRATE(Y0, Y2, F0, F1, F2) +
                      INTEGRATE(Y2, Y4, F2, F3, F4);
      end INTEGRATE;
   
      real procedure SUMMATION OF FACTORS M;
      begin integer I;
          real MSUM, COEF, MIMIN2, MIMIN1, MI;
          Boolean ADD;
          MSUM:= 0;
          if DF > 1 then
          begin
              MIMIN2:= A×WB × EXP(H) × PHI(HELP×WB) ×
                       .3989422804;
               if DFEVEN then MSUM:= MSUM + MIMIN2;
               if DF > 2 then
               begin COEF:= 1;
                   MIMIN1:= B×(HELP×MIMIN2 +
                                   A×.1591549431×EXP(-.5×D2));
                     if ¬ DFEVEN then MSUM:= MSUM + MIMIN1;
                     ADD:= DFEVEN;
                     for I:= 2 step 1 until DF - 2 do
                     begin MI:= (I - 1)/I×B×
                                   (COEF×HELP×MIMIN1 + MIMIN2);
                         if ADD then MSUM:= MSUM + MI;
                         ADD:=¬ ADD; COEF:= 1/(I - 1)/COEF;
                         MIMIN2:= MIMIN1; MIMIN1:= MI;
                     end I;
              end DF>2;
          end DF>1;
          SUMMATION OF FACTORS M:= MSUM;
      end SUMMATION OF FACTORS M;
      
      procedure INITIALISATION;
      begin TOL:= 10-8;
          if DF < 1 ∨ ENTIER(DF) ≠ DF then
          STATAL3 ERROR(“NCSTUDENT”, 2, DF);
          DFEVEN:= ENTIER(DF/2) = DF/2;
          A:= X/SQRT(DF); A2:= A×A; D2:= DELTA×DELTA;
          HELP:= DELTA×A;
          B:=DF/(DF + X×X); WB:= SQRT(B); H:=-D2×B×.5;
          if ABS(A) > TOL then TOLI:= 180 × TOL / ABS(A);
      end INITIALISATION;
      
      INITIALISATION;
      
      RESULT :=
      if DFEVEN then
          PHI(-DELTA) + SUMMATION OF FACTORS M × 2.5066282746
      else
          PHI(-DELTA×WB) + SUMMATION OF FACTORS M × 2 +
          (if ABS(A) ≤ TOL then 0 else
          .31830 98862 × INTEGRATE(0, A, EXP(H),
          EXP(H×(1 + A2/4))/(1 + A2/4),
          EXP(H×(1 + A2))/(1 + A2)) / 12);
          
      NCSTUDENT:=
    if TOL ≤ RESULT ∧ RESULT ≤ 1 - TOL then RESULT
    else
    if ABS(RESULT) < TOL then 0 else
    if ABS(RESULT - 1) < TOL then 1 else
    STATAL3 ERROR(“NCSTUDENT”, 0, RESULT);
end NCSTUDENT;
eop