Xor parity constraint with more than two variables

Question

I am using the following Python function to model Xor constraints with an arbitrary number of variables for Microsoft's Z3 solver with Python API:

#  parity function
#  break up long input lists in two smaller lists
def odd(solver, lits):
    length = len(lits)
    check(length > 0, "Odd needs argument(s)")
    if length == 1:
        solver.add(lits[0])
    elif length == 2:
        solver.add(Xor(lits[0], lits[1]))
    elif length == 3:
        solver.add(Xor(lits[0], Xor(lits[1], lits[2])))
    elif length == 4:
        #  this symmetric form is much faster than the chained forms
        solver.add(Xor(Xor(lits[0], lits[1]), Xor(lits[2], lits[3])))
    else:
        aux = get_aux_variable()

        #  cf. http://www.gregorybard.com/papers/bard_thesis.pdf
        cut_len = 3
        odd(solver, lits[:cut_len] + [aux])
        odd(solver, [aux] + lits[cut_len:])

auxVariableIdx = 0

def get_aux_variable():
    global auxVariableIdx

    auxVariableIdx += 1
    aux = Bool('t%s' % auxVariableIdx)
    return aux

I've noticed that the overall solver performance depends quite a lot on the way, the smaller Xor sub-cases are modelled. Example: I am using "Xor(Xor(a,b), Xor(c,d))" to declare an Xor with four inputs. The alternative "Xor(a,Xor(b,Xor(c,d)))" turned out to be slower by a factor of ten for my examples.

Another example: It makes a big difference if the auxiliary switching variable created to link sub-expressions is added to the inputs lists at the front or at the back of the list.

What is the recommended way to model Xor constraints with many variables in z3py?