# -*- coding: utf-8 -*-
"""
Spyder Editor

This temporary script file is located here:
/home/alkis/.spyder2/.temp.py
"""
from sympy import (var, degree, diff, Poly, expand, sturm,
                   Matrix, pprint, rem, det)

from sympy.polys.subresultants_qq_zz import *

x = var('x')

def variations(num_list):
    """
    num_list is a numerical list;
    return the number of sign variations in the list.
    """
    nlp = []   # copy num_list into nlp
    [nlp.append(num_list[i]) for i in range(len(num_list))] 
    while 0 in nlp: nlp.remove(0)
    counter = 0
    for j in range(len(nlp) - 1):
        if nlp[j] * nlp[j+1] < 0:
            counter = counter + 1
    return counter


#### FOURIER's theorem for max number of real roots in (0,2)
print('FOURIER theorem ', '\n')
f = x**3 - 7*x + 7

Fseq = []
[Fseq.append(diff(f,x,i)) for i in range(degree(f)+1)]

print('Fseq = ', Fseq, '\n')
    

# evaluate at left end
fL = []
[fL.append(Fseq[i].subs(x,0)) for i in range(degree(f)+1)]
print('fL = ', fL, '\n')

vfL=variations(fL)

print('vfL = ', vfL, '\n')

# evaluate at right end
fR = []
[fR.append(Fseq[i].subs(x,2)) for i in range(0, degree(f)+1)]
print('fR = ', fR, '\n')
 
    
vfR=variations(fR)

print('vfR = ', vfR, '\n')

# max number of real roots in (0,2) is vL-vR

print('possible # of real roots inside (0, 2) is = vfL - vfR = ', vfL - vfR, '\n')


#### BUDAN's theorem theorem for max number of real roots in (0, 2)
print('BUDAN theorem ', '\n')
f = x**3 - 7*x + 7

# at left end poly has 2 sign variations

bL = Poly(f).all_coeffs()
print('coeffs at left point = ', bL, '\n')
vbL = variations(bL)


# at right end poly has 0 sign variations
g = expand(f.subs(x,x+2))
print(g)

bR=Poly(g).all_coeffs()
print('coeffs at right point = ', bR, '\n')
vbR = variations(bR)

# max number of real roots in (0,2) is vLL-vRR
print('max number of real roots inside (0,2) is vbL-vbR = ')
print(vbL - vbR, '\n')


## STURM's theorem for exact number of real roots in (0,2)
print('STURM theorem ', '\n')
Sseq = sturm(f)
print(Sseq, '\n')

# evaluate at left end
sL = []
[sL.append(Sseq[i].subs(x,0)) for i in range(degree(f)+1)]
print('sL = ', sL, '\n')

vsL=variations(sL)

print('vsL = ', vsL, '\n')

# evaluate at right end
sR = []
[sR.append(Sseq[i].subs(x,2)) for i in range(0, degree(f)+1)]
print('sR = ', sR, '\n')
 
    
vsR=variations(sR)

print('vsR = ', vsR, '\n')

# EXACT number of real roots in (0,2) is vsL-vsR

print('EXACT # of real roots inside (0, 2) is = vL - vR = ', vsL - vsR, '\n')



# we can construct this sequence using rem_z()

def sturmSeqZ(f, x):
    ## sequence in ZZ[x];  use rem_z
    Sseq = []
    a0 = f
    a1 = diff(f,x,1)
    a2 = -rem_z(a0, a1, x)
    Sseq.append(a0)
    Sseq.append(a1)
    Sseq.append(a2)
    while degree(a2,x) >= 1:
        a0 = a1
        a1 = a2
        a2 = -rem_z(a0, a1, x)
        Sseq.append(a2)
    return Sseq
    
def sturmSeqQ(f, x):
    ## sequence in Q[x]; use rem
    Sseq = []
    a0 = f
    a1 = diff(f,x,1)
    a2 = -rem(a0, a1, x)
    Sseq.append(a0)
    Sseq.append(a1)
    Sseq.append(a2)
    while degree(a2,x) >= 1:
        a0 = a1
        a1 = a2
        a2 = -rem(a0, a1, x)
        Sseq.append(a2)
    return Sseq

      
print('POLYNOMIAL DIVISION by matrix triangularization ', '\n')
# polynomial remainders with matrix triangularization
M = [[3,0,-7,0],[0,3,0,-7],[1,0,-7,7]]


MS = Matrix(M) 
pprint(MS)
print('\n')
 

pprint(MS.LUdecompositionFF()[3])
print('\n')




f = x**5 + 5*x**2 - 7*x + 7
g = diff(f, x)
print('bezout','\n')
m = bezout(f, g, x)
pprint(m)
print('sylvester 2','\n')
m = sylvester(f, g, x, 2)
print('LUdecompositionFF()[3]','\n')
pprint(m.LUdecompositionFF()[3])
print('\n')
print('sturm = ', sturm_amv(f, g, x),'\n')
#print('subresultants_vv(f, g, x, 2)','\n')
#pprint(subresultants_vv(f, g, x, 2))
print('\n')
#print('euclid = ', euclid_amv(f, g, x),'\n')
print('\n')
#print('subresultants_vv(f, g, x)','\n')
#pprint(subresultants_vv(f, g, x))
print('\n')
#pprint(subresultants_vv_2(f, g, x))
print('\n')
print('sylvester 1','\n')
m = sylvester(f, g, x, 1)
pprint(m)
print('\n')
print('sylvester 2','\n')
m = sylvester(f, g, x, 2)
pprint(m)
print('\n')
#print('LUdecompositionFF()[3]','\n')
#pprint(m.LUdecompositionFF()[3])
#print('\n')
m2 = m[0:4,:]
pprint(m2)
print('\n')
print('coeff = ', det(m2[:,0:4]))
print('\n')
m2.col_swap(3,4)
print('coeff = ', det(m2[:,0:4]))
print('\n')
m2.col_swap(3,5)
print('coeff = ', det(m2[:,0:4]))
print('\n')
m2.col_swap(3,6)
print('coeff = ', det(m2[:,0:4]))
print('\n')
print('sturm = ', sturm_amv(f, g, x),'\n')

m2 = m[0:6,:]
pprint(m2)
print('\n')
print('coeff = ', det(m2[:,0:6]))
print('\n')
m2.col_swap(5,6)
print('coeff = ', det(m2[:,0:6]))
print('\n')
m2.col_swap(5,7)
print('coeff = ', det(m2[:,0:6]))
print('\n')


m2 = m[0:8,:]
pprint(m2)
print('\n')
print('coeff = ', det(m2[:,0:8]))
print('\n')
m2.col_swap(7,8)
print('coeff = ', det(m2[:,0:8]))
print('\n')

print('sturm = ', sturm_amv(f, g, x),'\n')

print('det = ',det(m))

f = x**10 - 10*x**8 + 12*x**6 - 5*x**3 + 1
g = diff(f, x)

print('euclid = ', sign_seq(euclid_amv(f, g, x),x),'\n')
print('sturm = ', sign_seq(sturm_amv(f, g, x),x),'\n')