// AKS Test for Primes task
// https://rosettacode.org/wiki/AKS_test_for_primes


#build ShowMoreWarnings NO

begin globals
  sInt64 c(100)
  //double n
end globals

local fn coef(nx as short)
  // out-by-1, ie coef(1)==^0, coef(2)==^1, coef(3)==^2 etc.
  c(nx) = 1
  short i
  for i = nx-1 to 2 step -1
    c(i) = c(i) + c(i-1)
  next
end fn

local fn is_prime(nx as short) as boolean
  short i
  bool result = _false
  fn coef(nx+1) // (I said it was out-by-1)
  for i = 2 to nx-1   // (technically "to n" is more correct)
    if int(c(i)/nx) <> c(i)/nx
      return _false
    end if
  next

  result = _true
end fn = result

local fn show(nx as short)
  // (As per coef, this is (working) out-by-1)

  double ci
  str255 cix
  short i

  for i = nx to 1 step -1
    ci = c(i)
    if ci = 1
      if (nx-i) mod 2 = 0
        if i = 1
          if nx = 1
            cix = "    1"
          else
            cix = "+1"
          end if
        else
          cix = ""
        end if
      else
        cix = "-1"
      end if
    else
      if (nx-i) mod 2 = 0
        cix = "+" + str$(ci)
      else
        cix = "-" + str$(ci)
      end if
    end if

    if i = 1 // ie ^0
      print cix;
    else
      if i = 2 then print cix, "x"; // ie ^1
      if i <> 2 then print cix, "x^", i-1;
    end if

  next i
end fn

local fn AKS_test_for_primes
  short nx

  for nx = 1 to 10 // (0 to 9 really)
    fn coef(nx)
    print "(x-1)^" + str$(nx-1) +  " = ";
    fn show(nx)
    print
  next

  print
  print "primes (<=53): ";
  short n
  c(2) = 1 // (this manages "", which is all that call did anyway...)
  for n = 2 to 53
    if fn is_prime(n)
      print " ", n;
    end if
  next
  print
end fn

window 1,@"AKS test for Primes",fn CGRectMake(0, 0, 850, 200)

fn AKS_test_for_primes

handleevents