abc1bN:=proc(N) ## list of a(n)=A216323(n) for b values in the range 1..N local list1pos,sqfkn,isabctriplex01,bN: list1pos:=proc(L::list) ## picks 1's in a list L and gives position local k,Lm: Lm:=proc(m) if L[m]=1 then m fi: end proc: [seq(Lm(k),k=1..nops(L))] end proc: sqfkn:=proc(n) ##(adapted from OEIS A007947) local j,t1,nrt1,t2: if n=1 then 1 else t1:=ifactors(n)[2]: nrt1:=nops(t1): t2 := product('op(1,t1[j])','j'=1..nrt1): fi: end proc: ## sqfkn isabctriplex01:=proc(a,b) local c,rad: c:=a+b: rad:=proc(a,b) sqfkn(a)*sqfkn(b)*sqfkn(c) end proc: ## rad if not type(a,posint) then print(`a not natural number`) elif not type(b,posint) then print(`b not natural number`) elif not igcd(a,b)=1 then x ## if GCD(1,b) not 1 elif rad(a,b)<c then 1 else 0 fi: end proc: ## isabctriplex01 bN:=[seq(isabctriplex01(1,b),b=1..N)]: list1pos(bN(N)): end proc: ## abc1bN ## A216323:=abc1bN(30000);