Support larger prime power bases in create_galois_field (issue #61)

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: lhs
 Title: Latin Hypercube Samples
-Version: 1.3.0
+Version: 1.4.0
 Authors@R: 
     person(given = "Rob",
            family = "Carnell",

NEWS

@@ -140,3 +140,10 @@ Changes in version 1.3.0 (2026-04-19)
 
 - Added ability to interrupt from long running C++ code processes
 - Updated R dependency to 4.0.0 and C++ to 17
+
+Changes in version 1.4.0 (2026-06-06)
+
+- Extended the Galois field power-cycle table (src/xtn.h) to cover prime power
+  bases up to 127, so create_galois_field() now supports GF(53^2), GF(59^2),
+  and other larger prime powers (issue #61).  Restored the generator script
+  etc/CreatePowerCycle.R used to build the table.

etc/CreatePowerCycle.R

@@ -0,0 +1,175 @@
+# Copyright 2026 Robert Carnell
+#
+# Generator for src/xtn.h
+#
+# src/xtn.h contains the "power cycle" representation of the Galois fields
+# GF(p^n) that are used by GaloisField.cpp.  For each prime power p^n the file
+# records a vector `xtn` of length n that defines a primitive (and therefore
+# irreducible) polynomial
+#
+#     x^n = xtn[1] + xtn[2] x + ... + xtn[n] x^(n-1)   (mod p)
+#
+# Because x is a primitive element, its powers cycle through every non-zero
+# element of the field (see Table B.3, "Power cycle II", p. 316 of Bailey,
+# "Design of Comparative Experiments").  GaloisField.cpp only requires the
+# polynomial to be irreducible, but primitive polynomials are used to stay
+# faithful to the original tables.
+#
+# Exponents n >= 2 are emitted while p^n < 1e9 (the same threshold used to
+# generate the original table).  The historical entries for p <= 47 are kept
+# as-is in src/xtn.h; this script was used to generate the appended entries for
+# 53 <= p <= 127 (issue #61).  Re-run with a larger `max_base` to extend the
+# table further, e.g. for applications such as covering arrays.
+#
+# Usage (from the package root):
+#   source("etc/CreatePowerCycle.R")
+#   writeLines(make_power_cycle_block(min_base = 53, max_base = 127))
+
+is_prime <- function(m) {
+  if (m < 2) return(FALSE)
+  if (m %% 2 == 0) return(m == 2)
+  i <- 3
+  while (i * i <= m) {
+    if (m %% i == 0) return(FALSE)
+    i <- i + 2
+  }
+  TRUE
+}
+
+prime_factors <- function(m) {
+  f <- numeric(0)
+  d <- 2
+  while (d * d <= m) {
+    while (m %% d == 0) {
+      f <- c(f, d)
+      m <- m / d
+    }
+    d <- d + 1
+  }
+  if (m > 1) f <- c(f, m)
+  unique(f)
+}
+
+primitive_roots <- function(p) {
+  if (p == 2) return(1)
+  order <- p - 1
+  fs <- prime_factors(order)
+  res <- integer(0)
+  for (g in 2:(p - 1)) {
+    # use modular exponentiation to avoid overflow for larger p
+    is_pr <- all(vapply(fs, function(q) powmod_int(g, order / q, p) != 1, logical(1)))
+    if (is_pr) res <- c(res, g)
+  }
+  res
+}
+
+# modular exponentiation of an integer base (avoids overflow)
+powmod_int <- function(base, e, mod) {
+  result <- 1
+  base <- base %% mod
+  while (e > 0) {
+    if (e %% 2 == 1) result <- (result * base) %% mod
+    e <- e %/% 2
+    if (e > 0) base <- (base * base) %% mod
+  }
+  result
+}
+
+# multiply two polynomials (length n, index 1 == coefficient on x^0) modulo the
+# field relation x^n = sum(cmod[j] x^(j-1))
+poly_mul <- function(a, b, p, n, cmod) {
+  long <- numeric(2 * n - 1)
+  for (i in 1:n) {
+    if (a[i] == 0) next
+    for (j in 1:n) {
+      if (b[j] != 0) long[i + j - 1] <- long[i + j - 1] + a[i] * b[j]
+    }
+  }
+  if (n > 1) {
+    for (i in (2 * n - 1):(n + 1)) {
+      v <- long[i] %% p
+      if (v != 0) {
+        for (j in 1:n) {
+          if (cmod[j] != 0) long[i - n + j] <- long[i - n + j] + v * cmod[j]
+        }
+      }
+      long[i] <- 0
+    }
+  }
+  long[1:n] %% p
+}
+
+poly_pow <- function(base, e, p, n, cmod) {
+  result <- c(1, numeric(n - 1))
+  b <- base
+  while (e > 0) {
+    if (e %% 2 == 1) result <- poly_mul(result, b, p, n, cmod)
+    e <- e %/% 2
+    if (e > 0) b <- poly_mul(b, b, p, n, cmod)
+  }
+  result
+}
+
+is_primitive <- function(cmod, p, n, order, fs) {
+  x <- c(0, 1, numeric(n - 2))
+  one <- c(1, numeric(n - 1))
+  if (!all(poly_pow(x, order, p, n, cmod) == one)) return(FALSE)
+  for (d in fs) {
+    if (all(poly_pow(x, order / d, p, n, cmod) == one)) return(FALSE)
+  }
+  TRUE
+}
+
+# the smallest primitive polynomial for GF(p^n); the constant term of a
+# primitive polynomial must be (-1)^(n+1) times a primitive root of GF(p), so we
+# only search those constant terms
+find_primitive <- function(p, n) {
+  order <- p^n - 1
+  fs <- prime_factors(order)
+  sign <- if (n %% 2 == 1) 1 else -1
+  c0set <- sort(unique((sign * primitive_roots(p)) %% p))
+  c0set <- c0set[c0set != 0]
+  rest_grid <- if (n == 1) NULL else as.matrix(expand.grid(rep(list(0:(p - 1)), n - 1)))
+  for (c0 in c0set) {
+    if (n == 1) {
+      if (is_primitive(c0, p, n, order, fs)) return(c0)
+    } else {
+      for (r in seq_len(nrow(rest_grid))) {
+        # rev() so the highest-order coefficient varies fastest, matching the
+        # search order used to generate the committed entries
+        cmod <- c(c0, rev(rest_grid[r, ]))
+        if (is_primitive(cmod, p, n, order, fs)) return(as.numeric(cmod))
+      }
+    }
+  }
+  stop(sprintf("no primitive polynomial found for GF(%d^%d)", p, n))
+}
+
+# build the C++ else-if blocks for all prime bases in [min_base, max_base]
+make_power_cycle_block <- function(min_base, max_base, limit = 1e9) {
+  primes <- Filter(is_prime, min_base:max_base)
+  out <- character(0)
+  for (p in primes) {
+    n <- 2
+    while (p^n < limit) {
+      cmod <- find_primitive(p, n)
+      q <- p^n
+      terms <- character(0)
+      for (j in seq_len(n)) {
+        cj <- cmod[j]
+        if (cj == 0) next
+        if (j == 1) terms <- c(terms, as.character(cj))
+        else if (j == 2) terms <- c(terms, paste0(if (cj == 1) "" else cj, "x"))
+        else terms <- c(terms, paste0(if (cj == 1) "" else cj, "x^", j - 1))
+      }
+      rhs <- if (length(terms)) paste(terms, collapse = " + ") else "0"
+      out <- c(out,
+               sprintf("// GF(%d^%d) = GF(%d)", p, n, q),
+               sprintf("// x^%d = %s", n, rhs),
+               sprintf("else if (q == primes::ipow(%d,%d)) xtn = {%s};",
+                       p, n, paste(cmod, collapse = ", ")))
+      n <- n + 1
+    }
+  }
+  out
+}

src/xtn.h

@@ -1,6 +1,6 @@
 /*
  * NOTE: This file should be excluded from the Doxygen build
- ****** COMPUTER GENERATED etc/CreatePowerCylce.R ******
+ ****** COMPUTER GENERATED etc/CreatePowerCycle.R ******
  *
  * file xtn.h
  * author Robert Carnell
@@ -385,4 +385,165 @@ else if (q == primes::ipow(47,4)) xtn = {42, 0, 0, 46};
 // GF(47^5) = GF(229345007)
 // x^5 = 5 + 46x
 else if (q == primes::ipow(47,5)) xtn = {5, 46, 0, 0, 0};
+// ---------------------------------------------------------------------------
+// Bases for prime powers with base p in (47, 127] added to address issue #61
+// (https://github.com/bertcarnell/lhs/issues/61).  As with the entries above,
+// each xtn defines a primitive (hence irreducible) polynomial
+//   x^n = xtn[0] + xtn[1] x + ... + xtn[n-1] x^(n-1)  (mod p)
+// and exponents are included while p^n < 1e9.  See etc/CreatePowerCycle.R.
+// ---------------------------------------------------------------------------
+// GF(53^2) = GF(2809)
+// x^2 = 2 + 5x
+else if (q == primes::ipow(53,2)) xtn = {2, 5};
+// GF(53^3) = GF(148877)
+// x^3 = 2 + x^2
+else if (q == primes::ipow(53,3)) xtn = {2, 0, 1};
+// GF(53^4) = GF(7890481)
+// x^4 = 2 + x^2 + x^3
+else if (q == primes::ipow(53,4)) xtn = {2, 0, 1, 1};
+// GF(53^5) = GF(418195493)
+// x^5 = 2 + 2x^4
+else if (q == primes::ipow(53,5)) xtn = {2, 0, 0, 0, 2};
+// GF(59^2) = GF(3481)
+// x^2 = 3 + x
+else if (q == primes::ipow(59,2)) xtn = {3, 1};
+// GF(59^3) = GF(205379)
+// x^3 = 2 + 4x^2
+else if (q == primes::ipow(59,3)) xtn = {2, 0, 4};
+// GF(59^4) = GF(12117361)
+// x^4 = 3 + 19x^3
+else if (q == primes::ipow(59,4)) xtn = {3, 0, 0, 19};
+// GF(59^5) = GF(714924299)
+// x^5 = 2 + 3x^4
+else if (q == primes::ipow(59,5)) xtn = {2, 0, 0, 0, 3};
+// GF(61^2) = GF(3721)
+// x^2 = 2 + 3x
+else if (q == primes::ipow(61,2)) xtn = {2, 3};
+// GF(61^3) = GF(226981)
+// x^3 = 2 + 3x^2
+else if (q == primes::ipow(61,3)) xtn = {2, 0, 3};
+// GF(61^4) = GF(13845841)
+// x^4 = 2 + 2x^3
+else if (q == primes::ipow(61,4)) xtn = {2, 0, 0, 2};
+// GF(61^5) = GF(844596301)
+// x^5 = 2 + 8x^4
+else if (q == primes::ipow(61,5)) xtn = {2, 0, 0, 0, 8};
+// GF(67^2) = GF(4489)
+// x^2 = 4 + 2x
+else if (q == primes::ipow(67,2)) xtn = {4, 2};
+// GF(67^3) = GF(300763)
+// x^3 = 2 + 2x^2
+else if (q == primes::ipow(67,3)) xtn = {2, 0, 2};
+// GF(67^4) = GF(20151121)
+// x^4 = 4 + 5x^3
+else if (q == primes::ipow(67,4)) xtn = {4, 0, 0, 5};
+// GF(71^2) = GF(5041)
+// x^2 = 2 + 3x
+else if (q == primes::ipow(71,2)) xtn = {2, 3};
+// GF(71^3) = GF(357911)
+// x^3 = 7 + 4x^2
+else if (q == primes::ipow(71,3)) xtn = {7, 0, 4};
+// GF(71^4) = GF(25411681)
+// x^4 = 2 + 3x^3
+else if (q == primes::ipow(71,4)) xtn = {2, 0, 0, 3};
+// GF(73^2) = GF(5329)
+// x^2 = 5 + x
+else if (q == primes::ipow(73,2)) xtn = {5, 1};
+// GF(73^3) = GF(389017)
+// x^3 = 5 + x^2
+else if (q == primes::ipow(73,3)) xtn = {5, 0, 1};
+// GF(73^4) = GF(28398241)
+// x^4 = 5 + 5x^3
+else if (q == primes::ipow(73,4)) xtn = {5, 0, 0, 5};
+// GF(79^2) = GF(6241)
+// x^2 = 2 + 2x
+else if (q == primes::ipow(79,2)) xtn = {2, 2};
+// GF(79^3) = GF(493039)
+// x^3 = 3 + 3x^2
+else if (q == primes::ipow(79,3)) xtn = {3, 0, 3};
+// GF(79^4) = GF(38950081)
+// x^4 = 2 + 10x^3
+else if (q == primes::ipow(79,4)) xtn = {2, 0, 0, 10};
+// GF(83^2) = GF(6889)
+// x^2 = 3 + 12x
+else if (q == primes::ipow(83,2)) xtn = {3, 12};
+// GF(83^3) = GF(571787)
+// x^3 = 2 + 2x^2
+else if (q == primes::ipow(83,3)) xtn = {2, 0, 2};
+// GF(83^4) = GF(47458321)
+// x^4 = 3 + 16x^3
+else if (q == primes::ipow(83,4)) xtn = {3, 0, 0, 16};
+// GF(89^2) = GF(7921)
+// x^2 = 3 + 5x
+else if (q == primes::ipow(89,2)) xtn = {3, 5};
+// GF(89^3) = GF(704969)
+// x^3 = 3 + 4x^2
+else if (q == primes::ipow(89,3)) xtn = {3, 0, 4};
+// GF(89^4) = GF(62742241)
+// x^4 = 3 + 14x^3
+else if (q == primes::ipow(89,4)) xtn = {3, 0, 0, 14};
+// GF(97^2) = GF(9409)
+// x^2 = 5 + 3x
+else if (q == primes::ipow(97,2)) xtn = {5, 3};
+// GF(97^3) = GF(912673)
+// x^3 = 5 + x^2
+else if (q == primes::ipow(97,3)) xtn = {5, 0, 1};
+// GF(97^4) = GF(88529281)
+// x^4 = 5 + 14x^3
+else if (q == primes::ipow(97,4)) xtn = {5, 0, 0, 14};
+// GF(101^2) = GF(10201)
+// x^2 = 2 + 6x
+else if (q == primes::ipow(101,2)) xtn = {2, 6};
+// GF(101^3) = GF(1030301)
+// x^3 = 2 + 2x^2
+else if (q == primes::ipow(101,3)) xtn = {2, 0, 2};
+// GF(101^4) = GF(104060401)
+// x^4 = 2 + 3x^3
+else if (q == primes::ipow(101,4)) xtn = {2, 0, 0, 3};
+// GF(103^2) = GF(10609)
+// x^2 = 2 + 2x
+else if (q == primes::ipow(103,2)) xtn = {2, 2};
+// GF(103^3) = GF(1092727)
+// x^3 = 5 + 4x^2
+else if (q == primes::ipow(103,3)) xtn = {5, 0, 4};
+// GF(103^4) = GF(112550881)
+// x^4 = 2 + 2x^3
+else if (q == primes::ipow(103,4)) xtn = {2, 0, 0, 2};
+// GF(107^2) = GF(11449)
+// x^2 = 3 + 3x
+else if (q == primes::ipow(107,2)) xtn = {3, 3};
+// GF(107^3) = GF(1225043)
+// x^3 = 2 + 5x^2
+else if (q == primes::ipow(107,3)) xtn = {2, 0, 5};
+// GF(107^4) = GF(131079601)
+// x^4 = 3 + 27x^3
+else if (q == primes::ipow(107,4)) xtn = {3, 0, 0, 27};
+// GF(109^2) = GF(11881)
+// x^2 = 6 + 3x
+else if (q == primes::ipow(109,2)) xtn = {6, 3};
+// GF(109^3) = GF(1295029)
+// x^3 = 6 + x^2
+else if (q == primes::ipow(109,3)) xtn = {6, 0, 1};
+// GF(109^4) = GF(141158161)
+// x^4 = 6 + x^3
+else if (q == primes::ipow(109,4)) xtn = {6, 0, 0, 1};
+// GF(113^2) = GF(12769)
+// x^2 = 3 + 3x
+else if (q == primes::ipow(113,2)) xtn = {3, 3};
+// GF(113^3) = GF(1442897)
+// x^3 = 3 + 3x^2
+else if (q == primes::ipow(113,3)) xtn = {3, 0, 3};
+// GF(113^4) = GF(163047361)
+// x^4 = 3 + 3x^3
+else if (q == primes::ipow(113,4)) xtn = {3, 0, 0, 3};
+// GF(127^2) = GF(16129)
+// x^2 = 9 + 2x
+else if (q == primes::ipow(127,2)) xtn = {9, 2};
+// GF(127^3) = GF(2048383)
+// x^3 = 3 + 3x^2
+else if (q == primes::ipow(127,3)) xtn = {3, 0, 3};
+// GF(127^4) = GF(260144641)
+// x^4 = 9 + 2x^3
+else if (q == primes::ipow(127,4)) xtn = {9, 0, 0, 2};
+
 #endif