cint university / 301 / ced, the cint editor

301-01: when 64 bits are not enough

This is lesson 301-01 of CU-301.1, checked against cint 0.2.0.

You will watch a product of two 64-bit numbers outgrow I64, compute it exactly in I128 with mul_full and meet the integer types that go up to 1024 bits.

the program

Fermat's test checks whether a number m could be prime. If m is an odd prime, then 2 to the power m - 1, divided by m, leaves a remainder of 1. A remainder other than 1 proves that m is not prime: it is composite, a product of smaller numbers. The power is far too large to compute, so the program keeps only the remainder after each multiplication:

// Fermat's test: if m is prime, then 2 to the power m - 1,
// divided by m, leaves 1. A remainder other than 1 proves m composite.

// a * b % m, for a and b below m.
I64 mul_mod(I64 a, I64 b, I64 m) {
    return a * b % m;
}

// base to the power exp, modulo m, by repeated squaring.
I64 pow_mod(I64 base, I64 exp, I64 m) {
    I64 result = 1;
    I64 b = base % m;
    I64 e = exp;
    while (e > 0) {
        if (e % 2 == 1) {
            result = mul_mod(result, b, m);
        }
        b = mul_mod(b, b, m);
        e /= 2;
    }
    return result;
}

void main() {
    I64[3] tests = [1_000_000_007, 2_305_843_009_213_693_951, 2_305_843_009_213_693_953];
    for i in 0..3 {
        I64 m = tests[i];
        "{m:>19}: remainder {pow_mod(2, m - 1, m)}\n";
    }
}

Run this program in ced

The first number, 1000000007, is a prime often used in programming contests. The second is 2⁶¹ - 1, which is also prime. The third is 2⁶¹ + 1, which 3 divides.

run it

$ cint run fermat.ci
         1000000007: remainder 1
fault[E_OVERFLOW]: the exact result is outside the type's range
  --> fermat.ci:6:14
   |
 6 |     return a * b % m;
   |              ^ mul.checked.i64
   |
   = operation: mul.checked.i64
   = operand-count: 2
   = operand: I64 4294967296
   = operand: I64 4294967296
   = exact:   18446744073709551616
   = limit:   I64 9223372036854775807
   = revision: c9b0d558c225959442c702c2a4e5015ffb1547a735342e30e49711666e8e56f3
   = source-map: none
   = address: none
   = stack-depth: 2
   = called from: fermat.ci:28:30
   = called from: fermat.ci:18:13
   = fuel-consumed: 97 (not part of the fault record)
   = state: the program stopped; stdout holds what it wrote before the fault

The first test passes. For 1000000007, every product of two numbers below it is about 10¹⁸ at most, far below I64.max. The second test squares b again and again: 2, 4, 16, 256, 65536, 4294967296. The next square is 2⁶⁴, and the fault record shows it as the exact result. A number below 2⁶¹ times another can need up to 122 bits, and I64 has 63 bits for the value and one for the sign.

the fix: an exact 128-bit product

mul_full multiplies two I64 values and returns the exact product as an I128, which can never overflow. The remainder is then small enough to convert back:

// a * b % m, for a and b below m, through the exact 128-bit product.
I64 mul_mod(I64 a, I64 b, I64 m) {
    I128 product = mul_full(a, b);
    return (product % (m as I128)) as I64;
}

Run this program in ced

Nothing else in the program changes:

$ cint run fermat.ci
         1000000007: remainder 1
2305843009213693951: remainder 1
2305843009213693953: remainder 4

Both primes leave 1. The remainder 4 proves that 2305843009213693953 is not prime, as expected, since 3 divides it. A remainder of 1 does not by itself prove a number prime, because some composite numbers also leave 1. The smallest is 341, which is 11 times 31. Programs that need certainty use a stronger test, such as Miller-Rabin with a fixed set of bases, and it needs this same exact product.

the wide integer types

cint 0.2.0 has signed integers in eight sizes, from I8 up. This table shows I64 and the wider ones. Each wider type holds every value of the narrower ones:

Type	Largest value	Digits in the largest value
`I64`	9223372036854775807	19
`I128`	170141183460469231731687303715884105727	39
`I256`	5789604461865809771...	77
`I512`	6703903964971298549...	154
`I1024`	8988465674311579538...	308

Each T.max above was printed by cint. I128 to I1024 are called the wide types. They take the same operators as I64, and every plain operator on them is checked the same way.

cint 0.2.0 also has the unsigned wide types U128, U256 and U512. mul_full of two U64 values gives a U128:

void main() {
    U64 a = U64.max;
    U128 square = mul_full(a, a);
    "{a} squared is {square}\n";
    "U128 runs up to {U128.max}\n";
}

Run this program in ced

18446744073709551615 squared is 340282366920938463426481119284349108225
U128 runs up to 340282366920938463463374607431768211455

widths never mix silently

An I64 value does not become an I128 on its own. Each of these lines is refused with error C2001, type mismatch:

I128 product = a * a;      // a * a is an I64 product, and it may already have overflowed
"{product % m}\n";         // product is I128 and m is I64

Run this program in ced

Convert a value with as to widen it. (a as I128) * (a as I128) computes the same exact product as mul_full(a, a). Widening never fails, since every I64 value is an I128 value.

Narrowing can fail, and as checks it. Converting an I128 that holds 2⁶³ to I64 faults:

fault[E_NARROW]: the value does not fit the destination type
  --> narrow.ci:5:11
   |
 5 |     "{big as I64}\n";
   |           ^ as.checked.i128.i64
   |
   = operation: as.checked.i128.i64
   = operand-count: 1
   = operand: I128 9223372036854775808
   = exact:   9223372036854775808
   = limit:   I64 9223372036854775807
   ...

In mul_mod, (product % (m as I128)) as I64 cannot fault, because a remainder after dividing by m is always less than m, and m is an I64.

a wider type moves the edge

I128 has an edge of its own. Doubling 1 in an I128 faults at the 127th doubling:

fault[E_OVERFLOW]: the exact result is outside the type's range
  --> over128.ci:4:16
   |
 4 |         grains *= 2;
   |                ^ mul.checked.i128
   |
   = operation: mul.checked.i128
   = operand-count: 2
   = operand: I128 85070591730234615865843651857942052864
   = operand: I128 2
   = exact:   170141183460469231731687303715884105728
   = limit:   I128 170141183460469231731687303715884105727
   ...

Pick the narrowest type that holds every value the program can reach, as mul_mod does: its inputs stay in I64, and only the product needs I128. A wide value takes more memory and each operation on it does more work.

refused on wide types in 0.2.0

A few forms that work on I64 are refused on the wide types in cint 0.2.0, with error C9102, not yet supported:

A switch cannot test a wide value.
as?, the conversion that returns an error instead of faulting, cannot convert from a wide type.
The built-ins that return an error instead of faulting, such as add_result, do not take wide values.
muldiv cannot give a wide result from narrow operands.

In place of a switch, compare with if and else if. In place of the others, convert a value that is known to fit to I64 first, or widen the operands with as.

elsewhere

The same test was written in C with uint64_t and a * b % m. It was built with GCC 13.3 and with Clang 18.1, using -O2 -Wall -Wextra on Ubuntu 24.04 for x86-64. Neither gave a warning, and both builds printed:

         1000000007: remainder 1
2305843009213693951: remainder 0
2305843009213693953: remainder 0

Unsigned arithmetic in C wraps without an error, so the products lost their high bits. A remainder of 0 would mean m divides a power of 2, which no odd number above 1 does. The program called a prime composite, and nothing signaled it. Computing the product as unsigned __int128, an extension that GCC and Clang provide on 64-bit targets, printed 1, 1 and 4, as cint does.

__int128 is signed and unchecked. A C program that doubled an __int128 128 times, built with the same compilers, printed that the result was zero at -O0 and positive at -O2, with no warning at either. Signed overflow in C is undefined behavior, so the C standard places no requirement on what the program does. C23 adds _BitInt(N) for integers of a chosen width. Clang 18.1 built a _BitInt(256) program, and GCC 13.3 did not recognize _BitInt. GCC's release notes list it as new in GCC 14.

Python's int has no fixed size. pow(2, m - 1, m) in Python 3.13 on the same machine printed 1, 1 and 4.

walk-through

mul_full(a, b) takes two values of the same type and returns their exact product in the type twice as wide: I64 to I128, I128 to I256, I256 to I512 and I512 to I1024.

m as I128 widens m so that % has two operands of one type.

e % 2 == 1 tests the lowest bit of the exponent. Repeated squaring multiplies result by b for each 1 bit, so about 120 multiplications cover an exponent near 2⁶¹.

I64[3] tests is declared inside main, since an array cannot be declared at the top level of a program in cint 0.2.0.

try this

Add 9223372036854775783, the largest prime below 2⁶³, to tests in the version with mul_full. The array and the loop then need a length of 4. Does the product still fit, and what does the new line print?
Replace mul_full(a, b) with (a as I128) * (b as I128). Does any output change?
The first version is safe for every m up to some size. What is the largest such m, and why?

answer to 3

Every product in the program multiplies two numbers below m, so the largest is (m - 1)². It fits in an `I64` while m - 1 is at most 3037000499, the integer square root of `I64.max`. So the first version is safe for every m up to 3037000500.

reference

integer types, wide integers, mul_full, conversions, checked arithmetic

Next: 301-04, big exact numbers.