Thanks. The large S-unit generator is indeed produced by pari (but I'm not certain whether these are incorrect generators, or whether the initial bug comes from testing them for being S-units, and whether this also explains other inconsistencies between coecions between L and OLSstar and back in other examples):
L = bnfinit(x^6 - 68463*x^4 - 5120808*x^3 + 1250774892*x^2 + 192368273328*x + 7520491439712,1)
S2 = idealprimedec(L,2)
S2 = [S2[3],S2[2],S2[1]]
S7 = idealprimedec(L,7)
S13 = idealprimedec(L,13)
S13 = [S13[1],S13[3],S13[2]]
S5 = idealprimedec(L,5)
S = concat([S2,S7,S13,S5])
#The order matters, this one is compatible with sage's S = L.primes_above(2*5*7*13)
sfu = bnfsunit(L,S)
In GP 2.11.2 (via Sage 9.0), the last two generators in sfu are indeed very large.
In GP 2.9.4, it runs for a minute before raising a stack overflow.
Running the same example with a different order of input primes raises another bug (perhaps they are related):
S2 = idealprimedec(L,2)
S2 = [S2[3],S2[2],S2[1]]
S5 = idealprimedec(L,5)
S7 = idealprimedec(L,7)
S13 = idealprimedec(L,13)
S13 = [S13[1],S13[3],S13[2]]
S = concat([S2,S5,S7,S13]) #
same order as what L.primes_above(2),...,L.primes_above(13) yields in sage.
sfu = bnfsunit(L,S)
In GP 2.11.2 (via Sage 9.0) this immediately raises an error:
*** bnfsunit: impossible inverse in ZM_inv: [60, 40, 2, 0, 0, 0, 0, 54; 0, 20, 2, 0, 0, 0, 0, 0; 0, 0, 2, 0, 0, 0, 0, 0; 0, 0, 0, 6, 1, 2, 0, 0; 0, 0, 0, 0, 0, 6, 3, 0; 0, 0, 0, 0, 0, 0, 3, 0; 0, 0, 0, 0, 0, 0, 0, 42; 0, 0, 0, 0, 0, 0, 0, 0].
In GP 2.9.4 it again runs for a minute before raising a stack overflow.
Benjamin