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,12], b=10, p=2*/
GetAllGIs(0,12,10,2,1)