prime power sums

$ cat PrimePowerSums.icn
procedure main(args)

  if *args == 0
  then args := [ 5000 ]
  else args[1] := integer(args[1]) |
         stop("Expecting an integer argument, not \"" || args[1] || "\".")

  every write(prime_power_sums(args[1]))

  write(&errout, "max:  ", args[1], ", time:  ", time_fmt(&time))

  end # main


procedure partition(lst, n)

  # Return the index i into the given ordered list such that
  # lst[1:i] < n <= lst[i+1:].

  local l, r, m

  l := 1
  r := *lst + 1

  # lst[1:l-1] < n and n <= lst[r:n]

  while l < r do {
    m := (l + r)/2
    if lst[m] < n
    then l := m + 1
    else r := m
    }

  return l - 1

  end # partition


procedure prime_power_sums(max)

  # Generate all the prime power sums less then the given value.  A prime power
  # sum is a value v with the property that p = p1^2 + p2^3 + p3^4, where p1,
  # p2, and p3 are primes.

  local p2, p3, p4, p2s, p3c, p4q
  local pps, ppps
  local filter
  local max3, max2

  filter := table(0)

  every p4 := primes_at_most(max) do {

    p4q := p4*p4*p4*p4
    max3 := max - p4q
    if max3 < 4 then break

    every p3 := primes_at_most(max3) do {

      p3c := p3*p3*p3
      ppps := p4q + p3c
      max2 := max - ppps
      if max2 < 4 then break

      every p2 := primes_at_most(max2) do {

        p2s := p2*p2
        if p2s > max2 then break

        pps := ppps + p2s
        if filter[pps] == 0 then {
          filter[pps] := 1
          suspend pps
          }
        }
      }
    }

  end # prime_power_sums


procedure primes_at_most(max)

  # Generate the sequence of primes (in increasing order) that are at most the
  # given value.

  local i

  static primes, primes_upper_bound
  initial {
    primes := []
    primes_upper_bound := 0
    }

  if primes_upper_bound < max then {
    primes := primes_upto(max)
    primes_upper_bound := max
    }

  every suspend !primes[1:partition(primes, max)]

  end # primes_at_most


procedure primes_upto(max)

  # Return a list of primes (in ascending order) no greater than the given value.

  local sieve, i, primes

  sieve := list(max, 1)

  every i := 2 to max do
    if sieve[i] == 1 then
      every sieve[i*2 to max by i] := 0

  primes := []
  every i := 2 to max do
    if sieve[i] == 1 then
      put(primes, i)

  return primes

  end # primes_upto


procedure time_fmt(t)

  # Return the given time in ms as a conviently formatted string.

  if t < 1000 then
    return string(t) || " ms"

  t := integer(t/100.0 + 0.5)/10.0

  if t < 60 then
    return string(t) || " sec"

  t := integer((t/60.0)*10 + 0.5)/10.0

  return string(t) || " min"

  end # time_fmt


$ icont PrimePowerSums.icn
Translating:
PrimePowerSums.icn:
  main
  partition
  prime_power_sums
  primes_at_most
  primes_upto
  time_fmt
No errors
Linking:

$ PrimePowerSums 50
28
33
49
47
max:  50, time:  0 ms

$ for i in 2 3 4 5 6 ; do PrimePowerSums 5e$i | wc -l ; done
max:  500, time:  0 ms
53
max:  5000, time:  10 ms
395
max:  50000, time:  110 ms
2579
max:  500000, time:  1.2 sec
18899
max:  5000000, time:  13.4 sec
138932

$ PrimePowerSums 5e7 > /dev/null
max:  50000000, time:  3.1 min

$