code 35193;
  comment COMPUTATION OF NONEXPONENTIAL MODIFIED BESSEL
  FUNCTIONS OF FRACTIONAL ORDER;
  procedure NONEXP BESS IAPLUSN(A, X, N, IA); value A, X, N;
  real X, A; integer N; array IA;
  if X= 0 then 
  begin IA[0]:= if A= 0 then 1 else 0;
    for N:= N step -1 until 1 do IA[N]:= 0 end 
  else if A= 0 then 
  begin 
    NONEXP BESSI(X, N, IA)
  end else if A= .5 then 
  begin real C;
    C:= .797 884 560 802 865 * SQRT(X);
    NONEXP SPHER BESSI(X, N, IA);
    for N:= N step -1 until 0 do IA[N]:= C * IA[N]
  end else 
  begin integer M, NU; real R, S, LABDA, L, A2, X2;
    A2:= A+A; X2:= 2/X; L:=1;
    NU:= START(X,N,1) ; R:= S:= 0;
    for M:= 1 step 1 until NU do L:= L * (M+A2)/(M+1);
    for M:= NU step -1 until 1 do 
    begin R:= 1/(X2 *(A+M)+R); L:= L*(M+1)/(M+A2);
      LABDA:= L*(M+A) * 2; S:= R * (LABDA + S);
      if M  <=  N then IA[M]:= R
    end;
    IA[0]:= R:= 1/(1+S)/GAMMA(1+A)/X2 **A;
    for M:= 1 step 1 until N do IA[M]:= R:= IA[M] * R;
  end;
        eop