Add functions for primality testing

[primo-ppcg]

Related: #140

Adds the following functions for primality testing:

• (math/prime? n)
  A primality test, deterministic for all n less than 2^63.
• (math/next-prime n)
  Returns the next prime number strictly larger than n.
• (math/primes)
  A boundless prime number generator.

Also defines the following support functions in a separate C module:

• (math/powmod a b m)
  Modular exponentiation of a to the power of b mod m.
• (math/mulmod a b m)
  Modular multiplication of a and b mod m.
• (math/invmod a m)
  Modular multiplicative inverse of a mod m.
• (math/jacobi a m)
  Computes the Jacobi Symbol (a|m).

I haven't added any tests to the PR yet (now added), but it has passed this battery of tests: prime-tests.janet. This could be stripped down some before including in a test suite.

Prime testing for primes on the range 2^53 - 2^63 is about 16µs per prime on my system. Composites on average are much faster. This is down from around 600µs using a pure Janet implementation.

Possible Concerns

The C module is imported without prefix into the spork/math module. I don't know if this is the proper way to do it.

128-bit operations are used for modular multiplication where available. On Linux, I test with:

#if defined(__SIZEOF_INT128__)

which is also used in the Linux kernel, so it's probably sufficient.

For Windows I test with:

#elif defined(_WIN64) && defined(_MSC_VER)

which I can confirm works with VC 2019 and VC 2022.

The C functions use signed 64-bit integers in order to allow potentially negative arguments, which unfortunately means these functions (and therefore also the Janet functions that rely on them) don't work properly with int/u64 values exceeding 2^63. I think this could possibly be addressed, but not without a lot of code duplication in the C module.

[sogaiu]

it has passed this battery of tests: prime-tests.janet.

I gave these a try here (manually as well), and unsurprisingly, they passed :)

Also tried some random-ish invocations successfully.

My attempts were on a Linux machine. I don't have access to macos or any BSD machine -- I wonder if we can get someone to test on some of those.

Prime testing for primes on the range 2^53 - 2^63 is about 16µs per prime on my system. Composites on average are much faster. This is down from around 600µs using a pure Janet implementation.

Nice!

The C module is imported without prefix into the spork/math module.

That seems ok to me.

The C functions use signed 64-bit integers in order to allow potentially negative arguments, which unfortunately means these functions (and therefore also the Janet functions that rely on them) don't work properly with int/u64 values exceeding 2^63. I think this could possibly be addressed, but not without a lot of code duplication in the C module.

Perhaps this can be left as an exercise for the interested user down the line?

[sogaiu]

I got NetBSD running via qemu and tried out the PR and it seemed to worked fine.

Building was ok and the prime-tests.janet file gave no errors 👍

[sogaiu]

FWIW, ran the tests for e74b793 on Linux and NetBSD and things seemed fine 👍

[Techcable]

How common is testing for primes? What use case are you thinking of that makes it deserve to be in spork?

(Just curious, not trying to be combative...)

[primo-ppcg]

How common is testing for primes? What use case are you thinking of that makes it deserve to be in spork?

The most useful function is likely next-prime. A number of applications require one or more primes, often of a similar length. This of course also relies on prime?. The miller-rabin-prp? function is perhaps a bit too niche, and wouldn't need to be exported necessarily.

[sogaiu]

It was mentioned in #140 (where there was a query about interest in prime number generation) that some hash functions have initial values that depend on primes.

I ran into this when implementing some of the SHA-2 functions.

Granted, that's not directly testing for primality (as asked about above), but rather generating primes (^^;

---

I guess not everyone is really into simulating cicada populations over time :P

[primo-ppcg]

It was mentioned in #140 (where there was a query about interest in prime number generation) that> I ran into this when implementing some of the SHA-2 functions.

(take 64 (math/primes))

would definitely look cleaner, although with so few, listing them out is fine too.

[sogaiu]

(take 64 (math/primes))

would definitely look cleaner, although with so few, listing them
out is fine too.

Yes, I agree.

Using a form like the aforementioned allows one to skip the step of manually enumerating as well as verifying correctness of individual numbers and whether any had been skipped or duplicated (^^;

[sogaiu]

Tests from b828418 ran with no issues here (Linux and NetBSD) 👍

[primo-ppcg]

Tests from b828418 ran with no issues here (Linux and NetBSD) +1

Appreciated. I apologize for the force pushes, I just want to avoid a bunch of minor commits in the commit history, should this PR be merged.

[sogaiu]

I noticed today that janet itself has a primes example. Just mentioning it as there having been some additional interest in the past :)