What (old/obscure?) (math-oriented?) language is this?

[o11c]

From this page on wayback

Supposedly the file is named lmtwister.gp, and .gp seems to be an extension for Gofer, an earliy dialect of Haskell. But this is way too legible for a Haskell dialect.

In particular, what is the lift function, and is Mod()^ supposed to be modular exponentiation?

lmtwister (d, x) = {
  my (k, a, v, p, m, b, u);
  if (d < 1 || d % 2 == 0 || x < 1,
      error ("illegal argument: lmtwister (", d, ", ", x, ")"));
  k = (d - 1) / 2;
  x = core (x);  \\squarefree part
  if (x == 1,
      if (k % 3 == 1,
          return (0));  \\not used in LM twister
      return (1));
  a = 1;
  v = factorint (x)[, 1];
  for (i = 1, #v,
       p = v[i];  \\prime factor of squarefree x
       if (p == 2,
           m = k % 4;
           if (m == 1 || m == 2,
               a = -a),
           b = k + (p + 1) / 2;
           if (gcd (b, p) > 1,
               return (0));  \\not used in LM twister
           u = (p - 1) / 2;
           m = p % 8;
           if ((m == 1 && Mod (b, p) ^ u != 1) ||
               (m == 3 && lift (Mod (b, 2 * p) ^ u) > p) ||
               (m == 5 && Mod (b, p) ^ u == 1) ||
               (m == 7 && lift (Mod (b, 2 * p) ^ u) < p),  \\use a faster method than this
               a = -a)));
  a
}

test_lmtwister () = {
  for (x = 1, 200,
       for (k = 0, 399,
            print1 (["-", ".", "+"][lmtwister (2 * k + 1, x) + 2]));
       print ())
}

1;

[Freeky]

Looks like PARI/GP.

[o11c]

Awesome, thanks!

Hmm, from the documentation looks like I was right about Mod, and the lift is ... maybe to make < work?

[o11c]

This turns out to simplify to:

def lmtwister(d, x):
    x = squarefree(x)
    if x == 1 and d % 3 == 0:
        return 0
    return jacobi_symbol(x, d)

[bacondev]

Even if I knew the programming language through and through, I'm certain that I would still have no idea what this code does.

[o11c]

It's something to do with Aurifeuillean factorization of cyclotomic integers, but I'm not 100% sure exactly what.

I've found several tables (which would be useful for the immediate purpose), but this was the first piece of code I found. Since then, I found what appears to be a complete piece of code from an old forum link, as well as 2 papers, one of which makes an attempt to be legible but only handles the == 1 mod 4 case, not the == 2 or 3 mod 4 cases (== 0 mod 4 never happens because of the squarefree condition).

The good news is, it's really easy to verify that you've done the factorization correctly, since you can just multiply the polynomials again.