GetAllGIs(minA, maxA, b, p, onlyDirect=0) = { if(onlyDirect,allGIs=Vec(null,0);numAllGIs=0,numAllGIs = 1;allGIs = Vec(null,1);allGIs[1]=0); for(a=minA,maxA, if(a<4,maxY=floor(b^(a/p)-0.75), maxY=floor(b^(a/p)-1)); for(y=0,maxY, if(Mod(y,b)!=0 || y==0, g=FindGIs(a,y,b,p,onlyDirect); if(g!=null, numAllGIs=numAllGIs+length(g); allGIs = Vec(allGIs,numAllGIs); for(i=1,length(g),allGIs[numAllGIs - length(g) + i]=g[i]); ) ) ) ); if(numAllGIs==0,return(null),return(vecsort(allGIs))) } FindGIs(a,y,b,p,onlyDirect)={ x=GetGN(a,y,b,p); if(onlyDirect&&(y==0||x>=1-10^-precision(0.1)), return(null)); if(x>1&&a>=p, return(null)); searchX=x; if(onlyDirect, numTests = 1; while(searchX < 1/b, searchX=searchX*b; numTests=numTests+1); if (numTests > floor((precision(x) - a) / p) - 1, numTests = floor((precision(x) - a) / p) - 1); if (a >= precision(x) / 2, numTests = 1) , if(y==0&&a==0,numTests=1,numTests = floor((precision(0.1)-a) / p) - 1); ); GIs=null; numGIsFound=0; for(n=1,numTests, gis = SearchGNForGIs(x,a,y,b,p,n); if(gis != null, for(ind=1, length(gis), if(Mod(gis[ind],b^a) != 0 || y==0, numGIsFound=numGIsFound+1; GIs = Vec(GIs,numGIsFound); GIs[numGIsFound]=gis[ind]; ) ) ) ); return(GIs); } GetGN(a,y,b,p,solNum=1) = { if(a==0&&y==0,return(1), if(y==0,return(b^(a/(p-1))), if(a==0,return(0), if(solNum==1, if(a+y<20, return(sumpos(n=1,ypC(n,y,p)/b^(a*n))), return(suminf(n=1,ypC(n,y,p)/b^(a*n)))) , if(solNum==2, return(b^(a/(p-1))-p*y/(p-1)-sumpos(n=1,GetMultiBinomial(p*n,n,p-1)*y^(n+1)/((n+1)*(p-1)*b^(a*n/(p-1))))) ) ) ) ) ) } ypC(n,y,p) = binomial(p*n,n)*y^((p-1)*n+1)/((p-1)*n+1) GetMultiBinomial(n,k,m) = return(GetMultiFact(n,m) / (GetMultiFact(k,m)*GetMultiFact(n-k,m))) GetMultiFact(n,m) = return(if(n<m,1,n*GetMultiFact(n-m,m))) SearchGNForGIs(x,a,y,b,p,n) = { z=x*b^(a+p*(n-1)); ceilZ = ceil(z); floorZ = floor(z); r = GetRValue(x,y,p); i = GetIValue(r); retGIs = Vec(null,1); if(r<0, retGIs[1] = floorZ; numSGNFGI = 1; d = floorZ-z - 1; while(i<=d && floorZ/numSGNFGI > 2, retGIs = Vec(retGIs,numSGNFGI+1); retGIs[numSGNFGI+1] = floorZ-numSGNFGI; numSGNFGI=numSGNFGI+1; d=d-1; ); return(retGIs); , d=z-floorZ; numSGNFGI = 0; while(d>=i, retGIs = Vec(retGIs,numSGNFGI+1); retGIs[numSGNFGI+1] = ceilZ+numSGNFGI; numSGNFGI = numSGNFGI+1; d=d-1; if(floorZ/numSGNFGI <= 1, d=i-1); ); if(numSGNFGI==0,return(null),return(retGIs)) ); } GetRValue(x,y,p) = (y/x - p + 1)/p GetIValue(r) = if(r<0,1/r,1-1/r) /*Sample Run Using a = [0,14], b=10, p=4*/ GetAllGIs(0,14,10,4,1)