Given a number field $K$, I would like to compute (in Sage) generators for the group of totally positive units of $K$.
Update: I've tried some code (details below), which I've received some help on in the comments. To clarify, my ideal would be to write up my own algorithm (rather than use code I found, as I discuss below).
I see reference to a strategy in the comments of the related question magma generators for unit group/ sage totally positive. The comment says that "it takes a bit of linear algebra mod 2 to construct a basis of totally positive units" from an existing set of fundamental units. I see some matrix algebra modulo 2 in the Pari code I found, but I don't know Pari well.
Question: Can anyone clarify the linear algebra to be done modulo 2, because then I could just go ahead and implement it in Sage myself?
What I had previously tried:
Option 1: (existing Sage functions) The Sage function
K.S_units(S = [])
computes a set of (not necessarily positive) generators. For $K$ of degree 3 or 4, I have successfully used this to get to a system of totally positive units. But I want an algorithm I can write for general degrees. I searched for other Sage functions that might help, without success. (Maybe I missed something?)
Option 2: (existing GP/Pari functions) I'm aware that I can use Pari in Sage. I've tried to adapt code on page 103 of the "User's Guide to PARI", link at https://math.mit.edu/~brubaker/PARI/PARIusers.pdf to arrive at the following (where I use %%gp to tell Sage to use GP for the cell):
%%gp
bnf = bnfinit(x^3 + x^2 - 1)
S = bnfsignunit(bnf)
d = matsize(S)
S = matrix(d[1],d[2], i,j, if (S[i,j] < 0, 1,0))
S = concat(S, vectorv(d[1],i,1))
K = lift(matker(S * Mod(1,2)))
if(K, ex = mathnfmodid(K, 2), ex = 2*matid(d[1]))
units = concat(bnf.tu[2], bnf.fu)
result = vector(length(ex)-1, i, factorback([units, ex[,i+1]]))
print(result)
The code runs without errors. But I'm not getting expected results. For example, the shown polynomial of $x^3 + x^2 - 1$ gives me a result of $2x$ -- I note that this is just 2 times the unit returned by bnf.fu (which happens to already be totally positive). Meanwhile, I tried it on $x^4 + x - 1$ with result $(-2x, -x^3 - 1)$, neither of which is totally positive.