def FindFixed(RF, seed=None, which = 1, size=64):
    prec = RF.precision()+size
    CF = ComplexField(prec + max(64, prec/4))
    if seed is None:
        a0 = CF(7/22, which*4/3)
    else:
        a0 = CF(seed)
    delta = CF(2^(-(prec+2)))
    perturb = delta/(2^8)
    mag = 1
    maxout = 0
    while (mag > delta) and (maxout < 50):
        a0 = a0 + perturb
        a1 = CF(log(a0))
        diff = a1-a0
        axial = diff/(a0-1)
        mag = abs(diff)
        maxout = maxout + 1
        a0 = a1 + axial
    return a0
    
def calc_sing_deriv(RF, afix, size):
    CF = ComplexField(RF.precision()+2*size)
    # The upper and lower logarithms are complex conjugates, so adding
    # them cancels the imaginary parts. Thus, we only need twice the
    # real part of either logarithm.
    return vector(RealField(RF.precision()+2*size), size,
        [2*CF((afix - pi*I)/afix).real()] +
        [2*CF(factorial(k-1)/(-afix*(afix)^k)).real() for k in xrange(1, size)])

def calc_sing_taylor(RF, afix, size):
    return vector(RF, size,
        [(2/(-afix*k*(afix)^k)).real() for k in xrange(1, size+1)])
    
def calc_sing_exp_deriv(RF, afix, size):
    CF = ComplexField(RF.precision()+2*size)
    reta = []
    reta.append(2*(ln(1-afix)/afix).real())
    ram1 = 1/(afix-1)  # reciprocal of afix minus 1
    ram1k = 1/afix
    rks = []
    SN = [1]
    for k in xrange(1, size):
        ram1k *= ram1
        rks.append(ram1k.real())
        if k > 1:
            SN.append((k-1)*SN[k-2])
            for j in xrange(k-2, 0, -1):
                SN[j] = (j+1)*SN[j] + j*SN[j-1]
        snv = vector(RF, k, SN)
        rkv = vector(RF, k, rks)
        reta.append(-2*(snv*rkv))
    return vector(RealField(RF.precision()+2*size), size, reta)
    
def Vec1(RF, size):
    return vector(RealField(RF.precision()+2*size), [(k==0) for k in xrange(size)])
        
def VecAccel(RF, afix, size):
    s0 = calc_sing_deriv(RF, afix, size)
    s1 = calc_sing_exp_deriv(RF, afix, size)
    v1 = Vec1(RF, size)
    return (v1 - (s1-s0)).change_ring(RF)
    
def Solve_System(RF, vec, vec2 = None):
    size = vec.degree()
    if vec2 is not None:
        vv = matrix(2, size, [vec, vec2]).transpose()
    else:
        vv = vec
    m = matrix(RF, size, size, lambda r, c: ((c+1)^r) - (factorial(c+1) if (c==(r-1)) else 0))
    ret = m.solve_right(vv)
    return ret

def Solve_System_Unsafe(RF, vec, vec2 = None):
    size = vec.degree()
    if vec2 is not None:
        m = matrix(RF, size, size+2, [[((c+1)^r) -
            factorial(c+1)*(c==(r-1))
            for c in xrange(size)] + [vec[r], vec2[r]] for r in xrange(size)])
    else:
        m = matrix(RF, size, size+2, [[((c+1)^r) -
            factorial(c+1)*(c==(r-1))
            for c in xrange(size)] + [vec[r]] for r in xrange(size)])
    m.echelonize()
    return m.matrix_from_columns(range(size,m.ncols()))

def Solve_System_Gaussian(RF, vec, vec2 = None):
    size = vec.degree()
    m = []
    RF2 = RealField(RF.precision()+size//4)
    ret = vec.change_ring(RF2)
    if vec2 is not None:
        ret2 = vec2.change_ring(RF2)
    for rr in xrange(size):
        m.append(vector(RF, size, [((c+1)^rr) - (factorial(rr) if ((c+1)==rr) else 0) for c in xrange(size)]))
    for j in xrange(size):
        jrow = m[j].change_ring(RF2)
        ret[j] /= jrow[j]
        if ret2 is not None:
            ret2[j] /= jrow[j]
        jrow /= jrow[j]
        m[j] = jrow.change_ring(RF)
        for k in xrange(j+1, size):
            krow = m[k].change_ring(RF2)
            ret[k] -= krow[j]*ret[j]
            if ret2 is not None:
                ret2[k] -= krow[j]*ret2[j]
            krow -= krow[j]*jrow
            m[k] = krow.change_ring(RF)
    for j in xrange(size-1, 0, -1):
        for k in xrange(j-1, -1, -1):
            ret[k] -= RF2(m[k][j])*ret[j]
            if ret2 is not None:
                ret2[k] -= RF2(m[k][j])*ret2[j]
    if ret2 is not None:
        return ret, ret2
    else:
        return ret

def Solve_System_GaussElim(RF, vecs):
    if ((type(vecs) is list) or (type(vecs) is tuple)):
        rets = [vecs[k].change_ring(RF).__copy__() for k in xrange(len(vecs))]
    else:   # assume it's a vector at this point
        rets = [vecs.change_ring(RF).__copy__()]
    size = rets[0].degree()
    m = []
    for rr in xrange(size):
        m.append(vector(RF, size, [((c+1)^rr) - (factorial(rr) if ((c+1)==rr) else 0) for c in xrange(size)]))
    for j in xrange(size):
        for ret in rets:
            ret[j] /= m[j][j]
        m[j] /= m[j][j]
        for k in xrange(j+1, size):
            for ret in rets:
                ret[k] -= m[k][j]*ret[j]
            m[k] -= m[k][j]*m[j]
    #print m
    for j in xrange(size-1, 0, -1):
        for k in xrange(j-1, -1, -1):
            for ret in rets:
                ret[k] -= m[k][j]*ret[j]
    return rets

def Solve_System_Gauss_SpaceSaver(RF, vecs):
    if ((type(vecs) is list) or (type(vecs) is tuple)):
        rets = [vecs[k].change_ring(RF).__copy__() for k in xrange(len(vecs))]
    else:   # assume it's a vector at this point
        rets = [vecs.change_ring(RF).__copy__()]
    size = rets[0].degree()
    m = []
    for rr in xrange(size//2+1):
        m.append(vector(RF, size, [((c+1)^rr) - (factorial(rr) if ((c+1)==rr) else 0) for c in xrange(size)]))
    for rr in xrange(size):
        rw = rr
        if rr >= size//2:
            rw = size//2
            for c in xrange(size):
                m[rw][c] = ((c+1)^rr) - (factorial(rr) if ((c+1)==rr) else 0)
        for k in xrange(rr):
            for ret in rets:
                ret[rr] -= m[rw][k]*ret[k]
            for p in xrange(k+1, size):
                if k < size//2:
                    m[rw][p] -= m[rw][k]*m[k][p]
                else:
                    m[rw][p] -= m[rw][k]*m[size-1-k][p-k-1]
        recip = RF(1/m[rw][rr])
        for ret in rets:
            ret[rr] *= recip
        for k in xrange(rr+1, size):
            m[rw][k] *= recip
        if rr >= size//2:
            for p in xrange(rr+1, size):
                m[size-1-rr][p-rr-1] = m[rw][p]
    for rr in xrange(size-1, 0, -1):
        for k in xrange(rr-1, -1, -1):
            for ret in rets:
                if k < size//2:
                    ret[k] -= m[k][rr]*ret[rr]
                else:
                    ret[k] -= m[size-1-k][rr-k-1]*ret[rr]
    return rets

def Solve_System_Gauss_SpaceSaver2(RF, vecs):
    if ((type(vecs) is list) or (type(vecs) is tuple)):
        rets = [vecs[k].change_ring(RF).__copy__() for k in xrange(len(vecs))]
    else:   # assume it's a vector at this point
        rets = [vecs.change_ring(RF).__copy__()]
    size = rets[0].degree()
    m = []
    for rr in xrange(size//2+1):
        m.append(vector(RF, size, [((c+1)^rr) - (factorial(rr) if ((c+1)==rr) else 0) for c in xrange(size)]))
    for rr in xrange(size):
        rw = rr
        if rr >= size//2:
            rw = size//2
            for c in xrange(size):
                m[rw][c] = ((c+1)^rr) - (factorial(rr) if ((c+1)==rr) else 0)
        for k in xrange(rr):
            for ret in rets:
                ret[rr] -= m[rw][k]*ret[k]
            if k < size//2:
                m[rw][k+1:size] -= m[rw][k]*m[k][k+1:size]
            else:
                m[rw][k+1:size] -= m[rw][k]*m[size-1-k][0:size-k-1]
        recip = RF(1/m[rw][rr])
        for ret in rets:
            ret[rr] *= recip
        m[rw][rr+1:size] *= recip
        if rr >= size//2:
            m[size-1-rr][0:size-rr-1] = m[rw][rr+1:size]
    for rr in xrange(size-1, 0, -1):
        for k in xrange(rr-1, -1, -1):
            for ret in rets:
                if k < size//2:
                    ret[k] -= m[k][rr]*ret[rr]
                else:
                    ret[k] -= m[size-1-k][rr-k-1]*ret[rr]
    return rets

# Out of Core (OOC) Gaussian solver
def Solve_System_Gauss_OOC(RF, vecs, chunkSize = 32, startChunk = 0):
    if ((type(vecs) is list) or (type(vecs) is tuple)):
        rets = [vecs[k].change_ring(RF).__copy__() for k in xrange(len(vecs))]
    else:   # assume it's a vector at this point
        rets = [vecs.change_ring(RF).__copy__()]
    size = rets[0].degree()
    numChunks = size // chunkSize
    #chunks = []
    RF2 = RealField(24)
    m = []
    path = "slog_" + str(size) + "_" + str(chunkSize) + "_" + str(RF.precision())
    if not (os.access(path, os.F_OK)):
        os.mkdir(path)
    for rr in xrange(chunkSize):
        m.append(vector(RF, size))
    t0 = cputime()
    w0 = walltime()
    t1 = t0
    w1 = w0
    for ch in xrange(startChunk, numChunks):
        chRow = ch * chunkSize
        # Initialize current chunk
        for rr in xrange(chRow, chRow + chunkSize):
            rw = rr - chRow
            for c in xrange(size):
                m[rw][size-1-c] = ((c+1)^rr) - (factorial(rr) if ((c+1)==rr) else 0)
        # Perform Gaussian elimination from previous chunks
        for cc in xrange(ch):
            pc = load(path + "/chunk_" + str(cc))
            pcr = cc * chunkSize
            for rr in xrange(chRow, chRow + chunkSize):
                rw = rr - chRow
                for k in xrange(pcr, pcr + chunkSize):
                    rk = k - pcr
                    for ret in rets:
                        ret[rr] -= m[rw][size-1-k]*ret[k]
                    m[rw][0:size-1-k] -= m[rw][size-1-k]*pc[rk][0:size-1-k]
        # Perform Gaussian elimination within current chunk
        for rr in xrange(chRow, chRow + chunkSize):
            rw = rr - chRow
            for k in xrange(chRow, rr):
                rk = k - chRow
                for ret in rets:
                    ret[rr] -= m[rw][size-1-k]*ret[k]
                m[rw][0:size-k-1] -= m[rw][size-1-k]*m[rk][0:size-1-k]
            recip = RF(1/m[rw][size-1-rr])
            for ret in rets:
                ret[rr] *= recip
            m[rw][0:size-1-rr] *= recip
        m2 = []
        for rr in xrange(chRow, chRow + chunkSize):
            m2.append(m[rr-chRow][0:size-1-rr].__copy__())
        print "Finished chunk " + str(ch+1) + " of " + str(numChunks) + "; Saving to disk, do not abort until finished!"
        save(m2, path + "/chunk_" + str(ch) + ".sobj")
        save(rets, path + "/rets_" + str(ch % 4) + ".sobj")
        with open(path+"/chunk.txt", "w") as chFile:
            chFile.write(str(ch)+"\n")
        with open(path+"/chunk_backup.txt", "w") as chFile:
            chFile.write(str(ch)+"\n")
        t2 = cputime()
        w2 = walltime()
        print "Save completed.  Elapsed seconds for chunk. CPU: " + str(t2-t1) + "; Wall: " + str(int(1/2+w2-w1)) + "; Total Elapsed Wall: " + str(int(1/2+w2-w0))
        t1 = t2
        w1 = w2
    save(rets, path + "/rets_back_sub.sobj")
    #rets = Solve_BackSub_OOC(RF, rets, chunkSize, size)
    #for ch in xrange(numChunks-1, -1, -1):
    #    chRow = ch * chunkSize
    #    m = chunks[ch]
    #    for rr in xrange(chRow+chunkSize-1, chRow-1, -1):
    #        rw = rr - chRow
    #        for k in xrange(size-1, rr, -1):
    #            for ret in rets:
    #                ret[rr] -= m[rw][size-1-k]*ret[k]
    #for ch in xrange(len(chunks)):
        #chRow = ch * chunkSize
        #print "Chunk: " + str(ch)
        #for rr in xrange(chRow, chRow+chunkSize):
            #print "Row: " + str(rr)
            #print chunks[ch][rr-chRow].change_ring(RF2)
    return rets

def Solve_BackSub_OOC(RF, vecs, chunkSize, size = None):
    if ((type(vecs) is list) or (type(vecs) is tuple)):
        vsize = vecs[0].degree()
        if size is None:
            size = vsize
        rets = [vecs[k][0:size].change_ring(RF).__copy__() for k in xrange(len(vecs))]
    else:   # assume it's a vector at this point
        vsize = vecs.degree()
        if size is None:
            size = vsize
        rets = [vecs[0:size].change_ring(RF).__copy__()]
    numChunks = 1 + ((size-1) // chunkSize)
    path = "slog_" + str(vsize) + "_" + str(chunkSize) + "_" + str(RF.precision())
    for ch in xrange(numChunks-1, -1, -1):
        chRow = ch * chunkSize
        m = load(path + "/chunk_" + str(ch))
        for rr in xrange(min(size-1, chRow+chunkSize-1), chRow-1, -1):
            rw = rr - chRow
            for k in xrange(size-1, rr, -1):
                for ret in rets:
                    ret[rr] -= m[rw][vsize-1-k]*ret[k]
    return rets

def Init_OOC_Solver(RF, size, chunkSize):
    ch = 0
    path = "slog_" + str(size) + "_" + str(chunkSize) + "_" + str(RF.precision())
    if os.access(path, os.F_OK) and os.access(path+"/chunk.txt", os.R_OK):
        with open(path+"/chunk.txt", "r") as chFile:
            ch = int(chFile.readline().strip())
        vecs = load(path + "/rets_" + str(ch % 4) + ".sobj")
        ch += 1
    else:
        afix = FindFixed(RF, size=size)
        vecs = [Vec1(RF, size), VecAccel(RF, afix, size)]
    return ch, vecs