Correspondingly, here is the improved version of a1,

a1:=proc(b)
local A,i,j,k,arem,rem,dig;
arem:=proc() rem:=irem(b*j,i); dig:=convert(j,base,b); 
     seq(`if`(member(k*i-rem,dig),NULL,b*j+k*i-rem),
          k = 1 .. iquo(b-1+rem,i)) end;
A[1] := [$1..b-1];
for i from 2 do
     A[i] := [seq](arem(), j = A[i-1]);
     if A[i] = [] then break fi od;
i-1 = A[i-1]
end:

For decimal system, it works very fast,

tt:=time(): a1(10); time()-tt;

                           9 = [381654729]

                                0.016

In addition to the values in the worksheet, I calculated

a1(21);
                   18 = [509504810627259318940020]

a1(22);

                   19 = [5213608783930183587520297]

a1(23);

  21 = [3017921136920002608782486643, 3202513921881602294966195322,

        10912277070078955588719411381, 11506928512609494037599575154,

        19136629049078922486640055262, 19563154864313291881366262910,

        20927442077211608842020972981, 22095893110934309880723905997,

        24479055281091281744429147151, 26034588683899462487295851544,

        27944901128718982841182304058, 28236074834439966722409861288,

        30327538677373615144106864301, 30327538677373615144106864322,

        33286709815522787435710859625, 33637644479370238638607793874,

        39311846134353288073052119134]

s:=proc(L) local n,a,b,i;
n:=nops(L)-1;
a:=add(1/i,i=L);
b:=add(1/i^2,i=L);
((n-1)*(n/2-a)+(a^2-b)/2)*mul(i,i=L)
end;

s([0.5,0.5,0.5,0.5,0.5]);
                             0.5000000000
s([.500,.621,.656,.528,.595]);
                             0.6484707154
s([.500,.621,.656,.528,.595,.111]);
                             0.3398817684

For small b, say from 2 to 20, or 30, K(b) ~ (b-1)*e gives a good approximation. More precise asymptotics that is better for large b, such as b close to 1,000,000 and larger, seems to be K(b) ~ b*e - ln(b)*γ

Looking at my user number, 135, Joe Riel's 84, and Thomas Richard's 50, I found out that we have something in common. gcd(50,84)=2 gcd(50,135)=5 gcd(84,135)=3

The same calculation for even bases shows that the maximal number of digits cannot be b -2 for even bases - either it is b -1, or ≤b -3, because for a (b -2)-digit number with all digits different and non-zero, adding the remaining non-zero digit at the end automatically produces a number divisible by b -1. P.S. I called the previous comment Even Bases so I didn't have other choice as to call this Odd Bases. They are odd as in the well-known phrase: "All prime numbers are odd and 2 is the oddest."

Thank you, Joe. Both things are good. The next improvement can be achieved by separating the loop in 2 parts,

a:=proc(b)
local A,i,j,k,arem,brem,rem;
arem:=proc() rem:= irem(b*j,i); 
     seq(b*j+k*i-rem,k = `if`(rem=0,0,1) .. iquo(b-1+rem,i)) end;
brem:=proc() rem:=irem(b*j,i);
     if rem=0 then b*j elif b-1+rem<i then NULL else b*j+i-rem fi end;
A[1] := [$1..b-1];
for i from 2 to b-1 do
     A[i] := [seq](arem(), j = A[i-1]) od;
for i from b do 
     A[i] := [seq](brem(), j = A[i-1]);
     if A[i] = [] then break fi od;
i-1 = A[i-1]
end:

Meanwhile, I calculated a(16): 39 = [18872900738885736149574055538327802527212537551, 60753927368683934227793588395570842550542338031, 89515749136034833729775437005460258167590093634] Converted to hex, the numbers look like 34E4A468166CD8604EC0F8106AB4326098286CF AA44CE207C78FC30003C3CC0D8382E2078D07EF FAE06678C2E884607EB8B4E0B0A0F0603420342

Yes, there is a simple proof. A number written in base b is divisible by b -1 if and only if the sum of digits is divisible by b -1 because a[0]*b^k + a[1]*b^(k-1) + a[k] mod (b-1) = a[0] + a[1]+ a[k] mod (b-1). If all b -1 digits are different and non-zero, then their sum is 1+2+...+(b -1) = b*(b -1)/2. This number, divided by b -1 is b/2 that is not integer for odd b. Thus, there are no (b -1)-digit numbers in base b with all digits different and non-zero, divisible by b -1, q.e.d.

The probability that a K-digit number is divisible by K is approximately 1/K. The probability that first (K-1) digit of it is divisible by (K-1) is approximately 1/(K-1), etc... That gives a probability that a K-digit number satisfies the problem conditions, approximately 1/K! Now, there are (b-1)*b^(K-1) K-digit numbers in base b. Thus, the expected number of them, satisfying the problem conditions, is (b-1)*b^(K-1)/K! Approximating factorial by Stirling formula, we get the following asymptotics of K with expected number 1, K(b) ~ (b-1)*e that is pretty close to the calculated maximal possible values of K.