---------- ---------- ---------- ---------- ---------- ---------- ---------- ----------
# With divp you can divide a polynomial by a polynomial. Some examples:
#
# (x^4-1)/(x+1)
A=c(1,0,0,0,-1); B=c(1,1); divp(A,B)
# 1 0 0 0         -1 0 0 -1      
#   -1 0 0        1 0 -1         
#     1 0         -1 -1  
#       -1        0      STOP
#  x^3 - x^2 + x -1  and remainder  0 
#
#  Interpretation of the outputs:
#  1 0 0 0  ->  x^3 + 0*x^2 + 0*x + 0
#   -1 0 0  ->       -1*x^2 + 0*x + 0
#      1 0  ->                1*x + 0
#       -1  ->                     -1
#
# The parts on the right describe the algorithm:
# x^4     - 1  :  x + 1
# ----------------------
# x^4 + x^3       x^3
# -------------
# - x^3   - 1    -x^2
# - x^3 - x^2
# -------------
#  x^2    - 1     x
#  x^2 + x
# -------------
# - x - 1        -1
# - x - 1
# -------------
# 0
#
# (x^3-3*x^2+x-1)/(2*x^2+1)
divp( c(1,-3,1,-1), c(2,0,1) )
# 0.5 0     -3 0.5 -1      
#   -1.5     0.5 0.5        STOP
#  0.5*x - 1.5  and remainder  0.5*x + 0.5
#
# (x^5 + x^3 - 1) / (x^2 + x - 5)
divp( c(1,0,1,0,0,-1), c(1,1,-5) )
#  1 0 0 0         -1 6 0 0 -1    
#    -1 0 0        7 -5 0 -1      
#      7 0         -12 35 -1      
#        -12       47 -61         STOP
# x^3 - x^2 + 7*x - 12   remainder:  47*x - 61
#
# (2*x^8+7/4*x^7-33/4*x^6+3*x^5-41/4*x^4-3/4*x^3-7/4*x^2+ 41/4*x-5/4) / ( (x^2+x-5)*(2*x-1/4) )
divp(c(2,7/4,-33/4,3,-41/4,-3/4,-7/4,41/4,- 5/4), prodp( c(2,-1/4),c(1,1,-5) ) )
# 1 0 0 0 0 0     2 1.75 -10.25 -0.75 -1.75 10.25 -1.25  
#   1 0 0 0       -2 -1.75 10.25 -1.25   
#     -1          0      STOP
#  x^5 + x^3 - 1  and remainder  0