Problem. A programme needs base**exponent % modulus, and a direct loop multiplies once for every exponent step. Baseline complexity. The shown baseline performs O(e) modular multiplications for a non-negative exponent e. Solution. Square-and-multiply reads the binary digits of e and reuses squared factors. Measured result. With e = 500,000, the executed comparison reduced the shown loop from 500,000 to 25 modular multiplications; its best timed sample was 12,126.90 times faster on this host.[1]
Reproduce the result
Complete code
from __future__ import annotationsimport gcimport jsonimport platformimport sysimport timeitimport tracemallocMODULUS = 1_000_000_007BASE = 987_654_321EXPONENT = 500_000REPEATS = 5CALLS_PER_SAMPLE = 3def repeated_multiply(base: int, exponent: int, modulus: int) -> int: """Compute base**exponent % modulus with one multiply per exponent step.""" if exponent < 0 or modulus == 0: raise ValueError("exponent must be non-negative and modulus must be non-zero") value = 1 % modulus base %= modulus for _ in range(exponent): value = (value * base) % modulus return valuedef square_and_multiply(base: int, exponent: int, modulus: int) -> int: """Compute base**exponent % modulus by consuming the exponent's binary digits.""" if exponent < 0 or modulus == 0: raise ValueError("exponent must be non-negative and modulus must be non-zero") value = 1 % modulus factor = base % modulus while exponent: if exponent & 1: value = (value * factor) % modulus exponent >>= 1 if exponent: factor = (factor * factor) % modulus return valuedef operation_counts(exponent: int) -> dict[str, int]: """Count modular multiplications in the two implementations.""" if exponent == 0: return {"repeated": 0, "square_and_multiply": 0} repeated = exponent squares = exponent.bit_length() - 1 selected_multiplications = exponent.bit_count() return { "repeated": repeated, "square_and_multiply": squares + selected_multiplications, }def peak_python_allocations(function) -> int: gc.collect() tracemalloc.start() function(BASE, EXPONENT, MODULUS) _, peak = tracemalloc.get_traced_memory() tracemalloc.stop() return peakdef samples(function) -> list[float]: timer = timeit.Timer(lambda: function(BASE, EXPONENT, MODULUS)) return timer.repeat(repeat=REPEATS, number=CALLS_PER_SAMPLE)def fixed_seconds(values: list[float]) -> list[str]: """Render samples in fixed-point seconds for direct visual comparison.""" return [f"{value:.9f}" for value in values]def main() -> None: cases = [(2, 0), (2, 1), (7, 13), (BASE, EXPONENT), (-5, 31)] correctness = [ repeated_multiply(base, exponent, MODULUS) == square_and_multiply(base, exponent, MODULUS) == pow(base, exponent, MODULUS) for base, exponent in cases ] if not all(correctness): raise AssertionError("implementations disagree") invalid_input = [] for args in ((2, -1, MODULUS), (2, 4, 0)): try: square_and_multiply(*args) except ValueError: invalid_input.append(True) if invalid_input != [True, True]: raise AssertionError("invalid-input checks failed") repeated_samples = samples(repeated_multiply) binary_samples = samples(square_and_multiply) counts = operation_counts(EXPONENT) result = { "python": sys.version.split()[0], "implementation": platform.python_implementation(), "platform": platform.platform(), "machine": platform.machine(), "base": BASE, "exponent": EXPONENT, "modulus": MODULUS, "repeat": REPEATS, "calls_per_sample": CALLS_PER_SAMPLE, "correctness_cases": len(cases), "all_results_equal_to_builtin_pow": all(correctness), "invalid_input_checks": "negative exponent and zero modulus raise ValueError", "multiplications_repeated": counts["repeated"], "multiplications_square_and_multiply": counts["square_and_multiply"], "timing_sample_format": "fixed-point decimal strings in seconds (9 decimal places)", "repeated_samples_seconds": fixed_seconds(repeated_samples), "square_and_multiply_samples_seconds": fixed_seconds(binary_samples), "best_repeated_seconds_per_call": f"{min(repeated_samples) / CALLS_PER_SAMPLE:.9f}", "best_square_and_multiply_seconds_per_call": f"{min(binary_samples) / CALLS_PER_SAMPLE:.9f}", "best_ratio_repeated_over_square_and_multiply": ( min(repeated_samples) / min(binary_samples) ), "tracemalloc_peak_repeated_bytes": peak_python_allocations(repeated_multiply), "tracemalloc_peak_square_and_multiply_bytes": peak_python_allocations(square_and_multiply), } print(json.dumps(result, indent=2, sort_keys=True))if __name__ == "__main__": main()
Output
{ "all_results_equal_to_builtin_pow": true, "base": 987654321, "best_ratio_repeated_over_square_and_multiply": 12126.897548934625, "best_repeated_seconds_per_call": "0.114004999", "best_square_and_multiply_seconds_per_call": "0.000009401", "calls_per_sample": 3, "correctness_cases": 5, "exponent": 500000, "implementation": "CPython", "invalid_input_checks": "negative exponent and zero modulus raise ValueError", "machine": "aarch64", "modulus": 1000000007, "multiplications_repeated": 500000, "multiplications_square_and_multiply": 25, "platform": "Linux-6.12.47+rpt-rpi-v8-aarch64-with-glibc2.41", "python": "3.13.5", "repeat": 5, "repeated_samples_seconds": [ "0.352864558", "0.348574345", "0.345503772", "0.342014997", "0.342542642" ], "square_and_multiply_samples_seconds": [ "0.000042908", "0.000029407", "0.000028278", "0.000028204", "0.000028203" ], "timing_sample_format": "fixed-point decimal strings in seconds (9 decimal places)", "tracemalloc_peak_repeated_bytes": 200, "tracemalloc_peak_square_and_multiply_bytes": 160}
Environment
The program ran with CPython 3.13.5 on Linux 6.12.47+rpt-rpi-v8, aarch64, with glibc 2.41. It uses only the standard library. Python documents timeit as a tool for timing small code snippets; its default timer is time.perf_counter().[3] The allocation figures come from tracemalloc, which traces memory blocks allocated by Python rather than total process memory.[4]
Methodology
Before timing, the program compared both implementations and Python’s three-argument pow across five fixed cases, including exponent zero and a negative base. It separately checked that a negative exponent and a zero modulus raise ValueError in the teaching functions. Each timing sample ran three calls; there were five samples per function. Each sample array is stored as fixed-point decimal strings in seconds, with nine places, so both implementations use the same human-readable notation. The timeit timing path temporarily disables garbage collection by default, so the recorded samples do not include GC work during the timed loops.[3] The input was fixed at base 987,654,321, exponent 500,000 and modulus 1,000,000,007. The reported ratio uses the fastest total from each five-sample group, divided by the same three calls. tracemalloc starts immediately before one call and stops when it returns. No random input, native RSS measurement or claim about other machines is involved.
The problem
Modular exponentiation appears in public-key encryption and digital signatures, where a calculation is reduced modulo a positive integer.[1] The mathematical task is compact: compute a raised to e, then keep the remainder after division by m. The first implementation follows the definition with admirable patience. It starts at one and multiplies by a, modulo m, exactly e times.
That is acceptable for a small exponent. It becomes an awkward loop counter when e reaches hundreds of thousands. Keeping the modulus in every step prevents the intermediate value from becoming a**e, but it does not remove the repeated multiplications.
Baseline complexity
For the exact repeated_multiply function, the loop body runs e times. Its multiplication count is therefore O(e), while its explicit Python state remains O(1). Those symbols describe the count of modular multiplications, not a promise that every multiplication costs the same amount of wall-clock time. Integer multiplication and remainder work depend on operand size; Python call overhead, CPU frequency and allocator state also exist outside the tidy notation.
The comparison uses a modulus near one billion, so the residues remain bounded by that modulus. It records 500,000 modular multiplications for the direct loop. The tracemalloc peak was 200 bytes for one measured call, but that narrow figure is Python-traced allocation during the selected interval, not resident memory for the interpreter.
The solution
The exponent has a binary representation. Repeated squaring builds factors for powers of two: a, a², a⁴, a⁸, and so on. When a bit of the exponent is one, the current factor belongs in the answer. The Handbook of Applied Cryptography describes general exponentiation techniques as repeated square-and-multiply algorithms and gives a right-to-left binary form.[1]
The shown square_and_multiply function keeps two values. factor starts as a mod m and is squared after each remaining bit. value starts as one and is multiplied by factor only when the current low bit is set. Shifting the exponent right discards that bit. When the exponent reaches zero, value is the required residue.
For a positive exponent with b binary digits and p set bits, this version makes (b - 1) + p modular multiplications. That is O(log e) because both b and p are bounded by the number of bits. For 500,000, the program found 19 bits and 7 set bits, giving 18 squarings plus 7 selected multiplications. The 25 is an operation count for this loop. It is not a count of CPU instructions.
Python’s pow(base, exp, mod) returns the modular power and the documentation says the three-argument form is more efficient than calculating pow(base, exp) % mod.[2] Use that built-in for ordinary Python work. The two teaching functions exist to make the work count visible rather than to replace the interpreter’s implementation.
What the measurement shows
The functions agreed with each other and with pow in all five fixed correctness cases. On this aarch64 host, the fastest direct-loop sample was 0.114004999 seconds per call. The fastest square-and-multiply sample was 0.000009401 seconds per call. Dividing the first best sample by the second gives 12,126.90 for this input and runtime.
The raw sample groups are in the Output block. Their absolute values should travel with the machine rather than being treated as a portable score. The operation gap is less fragile: the same exponent makes the shown baseline perform 500,000 modular multiplications and the shown binary loop 25. Different exponents have different bit patterns, and different moduli change the cost of each arithmetic operation.
The allocation peaks were 200 bytes for repeated multiplication and 160 bytes for square-and-multiply within the deliberately narrow tracing boundary. That 40-byte difference is too small to support a general memory claim. The useful result is the measured removal of repeated arithmetic for this one workload.
Where it stops helping
For a tiny exponent, the plain loop is often easier to inspect and its extra work may be irrelevant. The teaching version also rejects negative exponents, while Python’s three-argument pow documents conditions for modular inverses when the exponent is negative.[2] A real cryptographic implementation has further requirements, including a reviewed implementation and an appropriate side-channel model; this article measures neither.
Square-and-multiply saves multiplication count when an exponent is known and non-negative. It does not make arbitrary arithmetic free, and it does not settle which implementation wins for every interpreter, modulus or workload. It does give the exponent’s bits a proper job, which is more than can be said for a loop that counts all the way to 500,000.
Sources
[1] Menezes, van Oorschot and Vanstone, Handbook of Applied Cryptography, chapter 14
read0nly 