\\addhelp(upto, "Computes the sequence upto term max(n, 3).") upto(n) = {my(res = List([1, 1, 2])); slopes = List([List([0,1]), List([1])]); for(i=4,n, listput(res, nxt(res)) );res } addhelp(nxt, "Given the sequence upto some term m, computes the next term.") nxt(l) = {my(n = #l+1, taboo = List(), m, c); for(i=#slopes+1,n,listput(slopes, List())); \\find the 'forbidden numbers'; the numbers that create a contradiction. for(i=1,#l, fordiv(n-i, X, for(j=1,#slopes[X], c = l[i]+(n-i)*slopes[X][j]/X; if(c>0, listput(taboo, c) ) ) ) ); listsort(taboo,1); \\find the least positive integer that doesn't give a contradiction (= the next term). m = missing(taboo); \\put in the new slopes. for(i=1,#l, f = (m-l[i])/(n-i); listput(slopes[denominator(f)], numerator(f)); );m } addhelp(missing, "Finds the least positive integer not present in a sorted list l of positive integers.") missing(l, {start=1}, {end=#l}) = {my(mid); if(start > end, return(end+1)); if(start!=l[start], return(start)); mid = (start + end)>>1; if(l[mid]==mid, return(missing(l, mid+1, end)) , return(missing(l, start, mid)) ) }