The AKS algorithm for testing whether a number is prime is a polynomial-time algorithm based on an elementary theorem about Pascal triangles.

The theorem on which the test is based can be stated as follows:

   a number   



p


{\displaystyle p}

   is prime   if and only if   all the coefficients of the polynomial expansion of
 



(
x
−
1

)

p


−
(

x

p


−
1
)


{\displaystyle (x-1)^{p}-(x^{p}-1)}


are divisible by   



p


{\displaystyle p}

.



Example

Using   



p
=
3


{\displaystyle p=3}

:

         (x-1)^3 - (x^3 - 1)
            = (x^3 - 3x^2 + 3x - 1) - (x^3 - 1)
            = -3x^2 + 3x



And all the coefficients are divisible by 3,   so 3 is prime.





 Note:
This task is not the AKS primality test.   It is an inefficient exponential time algorithm discovered in the late 1600s and used as an introductory lemma in the AKS derivation.



Task



 Create a function/subroutine/method that given   



p


{\displaystyle p}

   generates the coefficients of the expanded polynomial representation of   



(
x
−
1

)

p




{\displaystyle (x-1)^{p}}

.
 Use the function to show here the polynomial expansions of   



(
x
−
1

)

p




{\displaystyle (x-1)^{p}}

   for   



p


{\displaystyle p}

   in the range   0   to at least   7,   inclusive.
 Use the previous function in creating another function that when given   



p


{\displaystyle p}

   returns whether   



p


{\displaystyle p}

   is prime using the theorem.
 Use your test to generate a list of all primes under   35.
 As a stretch goal,   generate all primes under   50   (needs integers larger than 31-bit).


References

 Agrawal-Kayal-Saxena (AKS) primality test (Wikipedia) 
 Fool-Proof Test for Primes - Numberphile (Video). The accuracy of this video is disputed -- at best it is an oversimplification.


 
program aks
  implicit none
 
  ! Coefficients of polynomial expansion
  integer(kind=16), dimension(:), allocatable :: coeffs
  integer(kind=16) :: n
  ! Character variable for I/O
  character(len=40) :: tmp
 
  ! Point #2
  do n = 0, 7
    write(tmp, *) n
    call polynomial_expansion(n, coeffs)
    write(*, fmt='(A)', advance='no') '(x - 1)^'//trim(adjustl(tmp))//' ='
    call print_polynom(coeffs)
  end do
 
  ! Point #4
  do n = 2, 35
    if (is_prime(n)) write(*, '(I4)', advance='no') n
  end do
  write(*, *)
 
  ! Point #5
  do n = 2, 124
    if (is_prime(n)) write(*, '(I4)', advance='no') n
  end do
  write(*, *)
 
  if (allocated(coeffs)) deallocate(coeffs)
contains
  ! Calculate coefficients of (x - 1)^n using binomial theorem
  subroutine polynomial_expansion(n, coeffs)
    integer(kind=16), intent(in) :: n
    integer(kind=16), dimension(:), allocatable, intent(out) :: coeffs
    integer(kind=16) :: i, j
 
    if (allocated(coeffs)) deallocate(coeffs)
 
    allocate(coeffs(n + 1))
 
    do i = 1, n + 1
      coeffs(i) = binomial(n, i - 1)*(-1)**(n - i - 1)
    end do
  end subroutine
 
  ! Calculate binomial coefficient using recurrent relation, as calculation
  ! using factorial overflows too quickly.
  function binomial(n, k) result (res)
    integer(kind=16), intent(in) :: n, k
    integer(kind=16) :: res
    integer(kind=16) :: i
 
    if (k == 0) then
      res = 1
      return
    end if
 
    res = 1
    do i = 0, k - 1
      res = res*(n - i)/(i + 1)
    end do
  end function
 
  ! Outputs polynomial with given coefficients
  subroutine print_polynom(coeffs)
    integer(kind=16), dimension(:), allocatable, intent(in) :: coeffs
    integer(kind=4) :: i, p
    character(len=40) :: cbuf, pbuf
    logical(kind=1) :: non_zero
 
    if (.not. allocated(coeffs)) return
 
    non_zero = .false.
 
    do i = 1, size(coeffs)
      if (coeffs(i) .eq. 0) cycle
 
      p = i - 1
      write(cbuf, '(I40)') abs(coeffs(i))
      write(pbuf, '(I40)') p 
 
      if (non_zero) then
        if (coeffs(i) .gt. 0) then
          write(*, fmt='(A)', advance='no') ' + '
        else
          write(*, fmt='(A)', advance='no') ' - '
        endif
      else
        if (coeffs(i) .gt. 0) then
          write(*, fmt='(A)', advance='no') '   '
        else
          write(*, fmt='(A)', advance='no') ' - '
        endif
      endif
 
      if (p .eq. 0) then
        write(*, fmt='(A)', advance='no') trim(adjustl(cbuf))
      elseif (p .eq. 1) then
        if (coeffs(i) .eq. 1) then
          write(*, fmt='(A)', advance='no') 'x'
        else
          write(*, fmt='(A)', advance='no') trim(adjustl(cbuf))//'x'
        end if
      else
        if (coeffs(i) .eq. 1) then
          write(*, fmt='(A)', advance='no') 'x^'//trim(adjustl(pbuf))
        else
          write(*, fmt='(A)', advance='no') &
            trim(adjustl(cbuf))//'x^'//trim(adjustl(pbuf))
        end if
      end if
      non_zero = .true.
    end do
 
    write(*, *)
  end subroutine
 
  ! Test if n is prime using AKS test. Point #3.
  function is_prime(n) result (res)
    integer(kind=16), intent (in) :: n
    logical(kind=1) :: res
    integer(kind=16), dimension(:), allocatable :: coeffs
    integer(kind=16) :: i
 
    call polynomial_expansion(n, coeffs)
    coeffs(1) = coeffs(1) + 1
    coeffs(n + 1) = coeffs(n + 1) - 1
 
    res = .true.
 
    do i = 1, n + 1
      res = res .and. (mod(coeffs(i), n) == 0)
    end do
 
    if (allocated(coeffs)) deallocate(coeffs)
  end function
end program aks
 