A = 4  (x^6) + 3  (x^5) + (x^4) + 7  (x^3) + 6  (x^2) + 2  x + 4

p = 13

n = Exponent[A, x]

lcoeff = Coefficient[A, x, n]

a = PolynomialMod[A PowerMod[lcoeff, -1, p], p]

reste = 
      Table[PolynomialMod[PolynomialRemainder[(x^(j p)), a, x], p], {j, 0, 
          n - 1}]

CoefficientList[(∑_(k = 0)^(n - 1))(Subscript[β, k]) (reste[[k + 1]] - (x^k)), x]

sol = 
      Solve[Union[
          Map[# == 0 &, 
            CoefficientList[
              Sum[(Subscript[β, k]) (reste[[k + 1]] - (x^k)), {k, 0, 
                  n - 1}], x]], {Modulus == p}], 
        Table[(Subscript[β, k]), {k, 0, n - 1}]]

b = 
      PolynomialMod[(∑_(k = 0)^(n - 1))(Subscript[β, k]) (x^k) /. 
            sol[[1]] /. (Subscript[β, l])_ -> 1, p]

tab = 
      Table[PolynomialGCD[a, b - k, Modulus -> p], {k, 0, 
          p - 1}]

Apply[Times, tab]

Factor[A, Modulus -> 13]

Clear[Berlekamp]
\(Berlekamp[A_, x_, p_] := 
      Module[{n, lcoeff, a, reste, j, k, sol, β, b, 
          tab}, 
n = 
          Exponent[A, x]; 
lcoeff = 
          Coefficient[A, x, n]; 
a = 
          PolynomialMod[A\ PowerMod[lcoeff, \(-1\), p], 
            p]; 
reste = 
          Table[PolynomialMod[PolynomialRemainder[x\^\(j\ p\), a, x], p], {j, 
              0, n - 1}]; 
sol = 
          Solve[Union[
              Map[# == 0 &, 
                CoefficientList[
                  Sum[\(β\_k\) \((reste[\([k + 1]\)] - x\^k)\), {k, 0, 
                      n - 1}], x]], {Modulus == p}], 
            Table[β\_k, {k, 0, n - 1}]]; 

        b = PolynomialMod[
            Sum[\(\(β\_k\) x\^k /. 
                  sol[\([1]\)]\) /. β\_l_ -> 1, {k, 0, n - 1}], 
            p]; 
tab = 
          Table[PolynomialGCD[a, b - k, Modulus -> p], {k, 0, 
              p - 1}]; 
lcoeff\ Apply[Times, 
            tab]
]\)

Berlekamp[A, x, p]

Berlekamp[(∑_(k = 0)^7) (x^k), x, 13]

Berlekamp[(x^4) + 1, x, 13]

Factor[(∑_(k = 0)^7) (x^k), Modulus -> 13]

Timing[Berlekamp[(∑_(k = 0)^7) (x^k), x, 10007]]

Timing[
      Factor[(∑_(k = 0)^7) (x^k), Modulus -> 10007]]