%%HP: T(3)A(D)F(.);
DIR

@ POLY, a polynomial utilities collection.
@ by Wayne H Scott

  ZTRIM  @ remove leading zeros
    \<< LIST\-> \-> n
      \<< n
        WHILE ROLL DUP 0 ==
        REPEAT DROP n 1 - DUP 'n' STO
        END n ROLLD
        IF n 0 ==
        THEN { 0 }
        ELSE n \->LIST
        END
      \>>
    \>>

  RDER  @ derivative of a rational function
    \<< \-> F G
      \<< G F PDER PMUL G PDER { -1 } PMUL F PMUL PADD G G PMUL
      \>>
    \>>

  IRT  @ invert root program (roots -> polynomial)
    \<< LIST\-> \-> n
      \<<
        IF n 0 >
        THEN 1 n
          START n ROLL { 1 } SWAP NEG +
          NEXT
        ELSE { 1 }
        END
        IF n 1 >
        THEN 2 n
          START PMUL
          NEXT
        END
      \>>
    \>>

  PDER  @ derivative of a polynomial
    \<< LIST\-> \-> n
      \<< 1 n
        FOR i n ROLL n i - *
        NEXT DROP
        IF n 1 ==
        THEN { 0 }
        ELSE n 1 - \->LIST
        END
      \>>
    \>>

  PF  @ partial fractions; handles multiple roots
    \<< MAXR \->NUM { } \-> Z P OLD LAST
      \<< 1 P SIZE
        FOR I P I GET \-> p1
          \<<
            IF p1 OLD \=/
            THEN Z p1 EVPLY 1 P SIZE
              FOR J
                IF P J GET P I GET \=/
                THEN p1 P J GET - /
                END
              NEXT p1 'OLD' STO { } 'LAST' STO
            ELSE
              IF { } LAST SAME
              THEN 1 { } 1 P SIZE
                FOR K P K GET
                  IF DUP p1 ==
                  THEN DROP
                  ELSE +
                  END
                NEXT IRT Z SWAP
              ELSE LAST LIST\-> DROP
              END RDER DUP2 5 PICK 1 + 3 ROLLD 3 \->LIST 'LAST' STO p1
              EVPLY SWAP p1 EVPLY SWAP / SWAP FACT /
            END
          \>>
        NEXT P SIZE \->LIST
      \>>
    \>>

  FCTP  @ factor polynomial
    \<<
      IF DUP SIZE 3 >
      THEN DUP BAIRS SWAP OVER PDIV DROP FCTP
      END
    \>>

  RT  @ find roots of any polynomial
    \<< ZTRIM DUP SIZE \-> n
      \<<
        IF n 3 >
        THEN DUP BAIRS SWAP OVER PDIV DROP \-> A B
          \<< A RT B RT 
          \>>
        ELSE
          IF n 2 >
          THEN QUD 
          ELSE LIST\-> DROP NEG SWAP /
          END
        END 
      \>>
    \>>

  L2A  @ LIST->ARRY and ARRY->LIST
    \<<
      IF DUP TYPE 5 ==
      THEN LIST\-> \->ARRY
      ELSE ARRY\-> 1 GET \->LIST
      END
    \>>

  PADD  @ add two polynomials
    \<< DUP2 SIZE SWAP SIZE \-> A B nB nA
      \<< A L2A B L2A
        IF nA nB <
        THEN SWAP
        END
        IF nA nB \=/
        THEN 1 nA nB - ABS
          START 0
          NEXT
        END nA nB - ABS 1 + ROLL ARRY\-> 1 GET nA nB - ABS + \->ARRY + L2A
      \>>
    \>>

  PMUL  @ multiply two polynomials by doing a convolution
    \<< DUP2 SIZE SWAP SIZE \-> A B nB nA
      \<< { }
        IF nB 1 >
        THEN 2 nB
          START 0 +
          NEXT
        END DUP A + SWAP + 'A' STO  @ pad A with zeros
        A LIST\-> \->ARRY B LIST\-> DROP
        IF nB 1 >
        THEN 2 nB
          FOR J J ROLL  @ reverse B
          NEXT
        END
        IF 3 nA nB + \<=
        THEN 3 nA nB +
          START 0
          NEXT
        END nA nB 1 - 2 * + \->ARRY  @ and pad with zeros
        2 nA nB +
        START DUP2 DOT  @ do dot product and then rotate B
          3 ROLLD ARRY\-> SWAP nA nB 1 - 2 * + 1 + ROLLD \->ARRY
        NEXT DROP2 nA nB + 1 - \->LIST  @ put dot products into list
      \>>
    \>>

  PDIV  @ divide two polynomials
    \<< DUP SIZE 3 ROLLD LIST\-> \->ARRY SWAP LIST\-> \->ARRY \-> c b a
      \<< a b
        WHILE OVER SIZE 1 GET c \>=
        REPEAT DIVV
        END DROP \-> d
        \<< a SIZE 1 GET c 1 - - \->LIST d ARRY\-> LIST\-> DROP \->LIST
        \>>
      \>>
    \>>

  EVPLY  @ evaluate a polynomial at a point
    \<< OVER
      IF DUP TYPE 5 ==
      THEN SIZE
      ELSE SIZE 1 GET
      END \-> a r n
      \<< a 1 GET
        IF n 1 >
        THEN 2 n
          FOR i r * a i GET +
          NEXT
        END
      \>>
    \>>

  COEF  @ equation -> polynomial list  (coefficient extractor)
    \<< \-> E n
      \<< 0 n
        FOR I 0 'X' STO E EVAL 'X' PURGE E 'X' \.d 'E' STO I FACT /
        NEXT 2 n 1 +
        FOR I I ROLL
        NEXT n 1 + \->LIST
      \>>
    \>>

  DIVV  @ subroutine used by PDIV
    \<< DUP 1 GET \-> a b c
      \<< a 1 GET c / DUP b * a SIZE RDM a SWAP - ARRY\-> 1 GETI 1 - PUT
         \->ARRY SWAP DROP b
      \>>
    \>>

  QUD  @ quadratic equation
    \<< LIST\-> \->ARRY DUP 1 GET / ARRY\-> DROP ROT DROP SWAP 2 / NEG DUP SQ
      ROT - \v/ DUP2 + 3 ROLLD -
    \>>

  BAIRS  @ Bairstow's method for quadratic factors
    \<< LIST\-> 1 1 \-> n R S
      \<<
        DO 0 n 1 + PICK 0 0 0 4 PICK 5 n + 7
          FOR J J PICK R 7 PICK * + S 8 PICK * + 7 ROLL DROP DUP 6 ROLLD R
            3 PICK * + S 4 PICK * + 5 ROLL DROP -1 
          STEP 3 PICK SQ 3 PICK 6 PICK * -
          IF DUP 0 ==
          THEN DROP 1 1
          ELSE 6 PICK 6 PICK * 5 PICK 9 PICK * - OVER / 4 PICK 9 PICK * 8 PICK
            7 PICK * - ROT /
          END DUP S + 'S' STO SWAP DUP R + 'R' STO
        UNTIL R\->C ABS .000000001 < 7 ROLLD 6 DROPN
        END n DROPN 1 R NEG S NEG 3 \->LIST
      \>>
    \>>
END