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.


 
AKS: procedure options (main, reorder); /* 16 September 2015, derived from Fortran */
 
  /* Coefficients of polynomial expansion */
  declare coeffs(*) fixed (31) controlled;
  declare n fixed(3);
 
 
  /* Point #2 */
  do n = 0 to 7;
    call polynomial_expansion(n, coeffs);
    put edit ( '(x - 1)^', trim(n), ' =' ) (a);
    call print_polynomial (coeffs);
  end;
 
  /* Point #4 */
   put skip;
  do n = 2 to 35;
    if is_prime(n) then put edit ( trim (n) ) (x(1), a);
  end;
 
  /* Point #5 */
  put skip;
  do n = 2 to 97;
    if is_prime(n) then put edit ( trim (n) ) (x(1), a);
  end;
  put skip;
 
 
 
  /* Calculate coefficients of (x - 1)^n using binomial theorem */
polynomial_expansion: procedure (n, coeffs);
    declare n fixed binary;
    declare coeffs (*) fixed (31) controlled;
    declare i fixed binary;
 
    if allocation(coeffs) > 0 then free coeffs;
    allocate coeffs (n+1);
 
    do i = 1 to n + 1;
      coeffs(i) = binomial(n, i - 1);
      if iand(n - i - 1, 1) = 1 then coeffs(i) = -coeffs(i);
    end;
  end polynomial_expansion;
 
  /* Calculate binomial coefficient using recurrent relation, as calculation	*/
  /* using factorial overflows too quickly.					*/
binomial: procedure (n, k) returns (fixed(31));
    declare (n, k) fixed;
    declare i fixed;
    declare result fixed (31) initial (n);
 
    if k = 0 then return (1);
 
    do i = 1 to k - 1;
      result = (result*(n - i))/(i + 1);
    end;
    return (result);
  end binomial;
 
  /* Outputs polynomial with given coefficients */
print_polynomial: procedure (coeffs);
    declare coeffs (*) fixed (31) controlled;
    declare ( i, p ) fixed binary;
    declare non_zero bit (1) aligned;
    declare (true initial ('1'b), false initial ('0'b)) bit (1);
 
    if allocation(coeffs) = 0 then return;
 
    non_zero = false;
 
    do i = 1 to hbound(coeffs);
      if coeffs(i) = 0 then iterate;
 
      p = i - 1;
 
      if non_zero then
         do;
            if coeffs(i) > 0 then
               put edit ( ' + ' ) (a);
            else
               put edit ( ' - ' ) (a);
         end;
      else
        do;
           if coeffs(i) > 0 then
              put edit ( '   ' ) (a);
           else
              put edit ( ' - ' ) (a);
        end;
 
      if p = 0 then
        put edit ( trim(abs(coeffs(i))) ) (a);
      else if p = 1 then
         do;
            if coeffs(i) = 1 then
               put edit ( 'x' ) (a);
            else
               put edit ( trim(abs(coeffs(i))), 'x' ) (a);
        end;
      else
        do;
            if coeffs(i) = 1 then
               put edit ( 'x^', trim(p) ) (a);
            else
               put edit (  trim(abs(coeffs(i)) ), 'x^', trim(p)) (a);
        end;
 
      non_zero = true;
    end;
 
    put skip;
  end print_polynomial;
 
  /* Test if n is prime using AKS test. Point #3. */
is_prime: procedure (n) returns (bit (1));
    declare n fixed (15);
    declare result bit (1) aligned;
    declare coeffs (*) fixed (31) controlled;
    declare i fixed binary;
 
    call polynomial_expansion(n, coeffs);
    coeffs(1) = coeffs(1) + 1;
    coeffs(n + 1) = coeffs(n + 1) - 1;
 
    result = '1'b;
 
    do i = 1 to n + 1;
      result = result & (mod(coeffs(i), n) = 0);
    end;
 
    if allocation(coeffs) > 0 then free coeffs;
    return (result);
  end is_prime;
 
end AKS;
 