Find Involutory Exponents

find_exp.c

/**
Find involutory exponents for a modulus like a Mersenne prime: 2^31-1
Uses brute force
Very fast for small numbers, very slow for anything more than 16-bits.

gcc -O2 find_exp.c -ofind_exp -lcrypto
*/
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <openssl/bn.h>

void 
find_involutory_exponents(BIGNUM *p) {
    BIGNUM *m = NULL, *e = NULL, *ee = NULL;
    BN_CTX *ctx = NULL;
    uint32_t cnt = 0;
    
    m = BN_new();
    e = BN_new();
    ee = BN_new();
    ctx = BN_CTX_new();
    
    // m = p - 1
    BN_sub(m, p, BN_value_one());
    
    // Loop over all potential e values from 1 to m-1
    for (BN_set_word(e, 1); BN_cmp(e, m) < 0; BN_add(e, e, BN_value_one())) {
        // compute e * e % m
        BN_mod_sqr(ee, e, m, ctx);
        
        // Check if ee == 1 % m
        if (BN_cmp(ee, BN_value_one()) == 0) {
            //if (cnt && !(cnt & 7)) putchar('\n');
            printf("%d : %s\n", (cnt + 1), BN_bn2hex(e));
            cnt++;
        }
    }
    
    printf("\n");
    
    // Clean up
    BN_free(m);
    BN_free(e);
    BN_free(ee);
    BN_CTX_free(ctx);
}

int 
main(int argc, char *argv[]) {
    puts("\nFind involutory Exponents for Modular Exponentiation.\n");
    
    if (argc != 2) {
        printf("find_invol <modulus>\n");
        return 0;
    }
    
    BIGNUM *p = NULL;
    
    
    do {
        p = BN_new();
        
        if (!p) {
            printf("BN_new() failed...\n");
            break;
        }
        
        // try as decimal first
        int result = BN_dec2bn(&p, argv[1]);
        
        if (!result) {
            // then hexadecimal
            result = BN_hex2bn(&p, argv[1]);
            if (!result) {
                printf("Can't parse input.\n");
                break;
            }
        }
        
        // try to find involutory exponents (can take along time)
        find_involutory_exponents(p);
    } while(0);
    
    // Clean up
    if (p) BN_free(p);
    return 0;
}

find_exp2.c

/**
Find involutory exponents for a modulus like a Mersenne prime: 2^31-1
Uses factorisation and the Chinese Remainder Theorem to find valid exponents.
Very fast for small numbers...

gcc -O2 find_exp.c -ofind_exp -lcrypto
*/
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <openssl/bn.h>

// Structure to hold prime factors and their exponents
typedef struct {
    BIGNUM *p;      // Prime factor
    int exp;        // Exponent
    BIGNUM *mod;    // p^exp
    int num_sols;   // Number of solutions modulo p^exp
    BIGNUM **sols;  // Solutions modulo p^exp
} PrimeFactor;

// Function to perform trial division factoring
int factorize(BIGNUM *m, PrimeFactor **factors, int *num_factors) {
    BN_CTX *ctx = BN_CTX_new();
    BIGNUM *i = BN_new();
    const BIGNUM *one = BN_value_one();
    BIGNUM *zero = BN_new();
    BIGNUM *remainder = BN_new();
    BIGNUM *m_copy = BN_new();

    BN_copy(m_copy, m);
    BN_zero(zero);
    *num_factors = 0;

    // Start from 2
    BN_set_word(i, 2);

    while (BN_cmp(i, m_copy) <= 0) {
        int exp = 0;

        // Check if i divides m_copy
        while (1) {
            BN_div(NULL, remainder, m_copy, i, ctx);
            if (BN_cmp(remainder, zero) == 0) {
                // i divides m_copy
                exp++;
                BN_div(m_copy, NULL, m_copy, i, ctx);
            } else {
                break;
            }
        }

        if (exp > 0) {
            // Store the prime factor and its exponent
            factors[*num_factors] = malloc(sizeof(PrimeFactor));
            factors[*num_factors]->p = BN_dup(i);
            factors[*num_factors]->exp = exp;
            factors[*num_factors]->mod = BN_new();
            BN_exp(factors[*num_factors]->mod, i, BN_value_one(), ctx);
            BN_copy(factors[*num_factors]->mod, i);
            for (int j = 1; j < exp; j++) {
                BN_mul(factors[*num_factors]->mod, factors[*num_factors]->mod, i, ctx);
            }
            (*num_factors)++;
        }

        BN_add(i, i, one);
    }

    if (BN_cmp(m_copy, one) > 0) {
        // m_copy is a prime number larger than sqrt(m)
        factors[*num_factors] = malloc(sizeof(PrimeFactor));
        factors[*num_factors]->p = BN_dup(m_copy);
        factors[*num_factors]->exp = 1;
        factors[*num_factors]->mod = BN_dup(m_copy);
        (*num_factors)++;
    }

    BN_free(i);
    BN_free(zero);
    BN_free(remainder);
    BN_free(m_copy);
    BN_CTX_free(ctx);

    return 0;
}

// Function to compute solutions modulo p^e
void compute_solutions_mod_pe(PrimeFactor *factor) {
    BN_CTX *ctx = BN_CTX_new();
    const BIGNUM *one = BN_value_one();
    BIGNUM *neg_one = BN_new();
    BIGNUM *tmp = BN_new();
    BN_sub(neg_one, factor->mod, one); // -1 mod p^e

    if (BN_cmp(factor->p, BN_value_one()) == 0) {
        // Skip p = 1
        factor->num_sols = 0;
        factor->sols = NULL;
    } else if (BN_is_word(factor->p, 2)) {
        // p = 2
        if (factor->exp == 1) {
            // For 2^1, only x ≡ 1 mod 2
            factor->num_sols = 1;
            factor->sols = malloc(sizeof(BIGNUM *) * 1);
            factor->sols[0] = BN_new();
            BN_one(factor->sols[0]);
        } else if (factor->exp == 2) {
            // For 2^2, x ≡ 1 mod 4 and x ≡ 3 mod 4
            factor->num_sols = 2;
            factor->sols = malloc(sizeof(BIGNUM *) * 2);
            factor->sols[0] = BN_new();
            factor->sols[1] = BN_new();
            BN_one(factor->sols[0]);
            BN_sub(factor->sols[1], factor->mod, BN_value_one()); // x ≡ -1 mod 4
        } else {
            // For 2^e where e >= 3, four solutions
            factor->num_sols = 4;
            factor->sols = malloc(sizeof(BIGNUM *) * 4);
            // x ≡ 1 mod 2^e
            factor->sols[0] = BN_new();
            BN_one(factor->sols[0]);
            // x ≡ -1 mod 2^e
            factor->sols[1] = BN_new();
            BN_copy(factor->sols[1], neg_one);
            // x ≡ 1 + 2^{e-1} mod 2^e
            factor->sols[2] = BN_new();
            BN_zero(tmp);
            BN_set_bit(tmp, factor->exp - 1); // tmp = 2^{e-1}
            BN_add(factor->sols[2], BN_value_one(), tmp);
            // x ≡ -1 + 2^{e-1} mod 2^e
            factor->sols[3] = BN_new();
            BN_add(factor->sols[3], neg_one, tmp);
        }
    } else {
        // For odd primes, solutions are x ≡ ±1 mod p^e
        factor->num_sols = 2;
        factor->sols = malloc(sizeof(BIGNUM *) * 2);
        // x ≡ 1 mod p^e
        factor->sols[0] = BN_new();
        BN_one(factor->sols[0]);
        // x ≡ -1 mod p^e
        factor->sols[1] = BN_new();
        BN_copy(factor->sols[1], neg_one);
    }

    BN_free(neg_one);
    BN_free(tmp);
    BN_CTX_free(ctx);
}

// Comparator function for qsort to sort BIGNUM pointers in ascending order
int compare_bn(const void *a, const void *b) {
    const BIGNUM *bn_a = *(const BIGNUM **)a;
    const BIGNUM *bn_b = *(const BIGNUM **)b;
    return BN_cmp(bn_a, bn_b);
}

// Function to compute solutions modulo m using CRT
void compute_solutions(PrimeFactor **factors, int num_factors, BIGNUM ***solutions, int *num_solutions, BIGNUM *m) {
    BN_CTX *ctx = BN_CTX_new();
    BIGNUM *M = BN_new();
    BN_copy(M, m);

    // Compute the number of total solutions
    int total_solutions = 1;
    for (int i = 0; i < num_factors; i++) {
        compute_solutions_mod_pe(factors[i]);
        total_solutions *= factors[i]->num_sols;
    }
    *num_solutions = total_solutions;

    // Allocate array for solutions
    *solutions = malloc(sizeof(BIGNUM *) * total_solutions);

    // Prepare arrays for CRT
    BIGNUM **moduli = malloc(sizeof(BIGNUM *) * num_factors);
    BIGNUM **Mi = malloc(sizeof(BIGNUM *) * num_factors);
    BIGNUM **invMi = malloc(sizeof(BIGNUM *) * num_factors);

    // Compute Mi = M / moduli[i] and inverse of Mi modulo moduli[i]
    for (int i = 0; i < num_factors; i++) {
        moduli[i] = BN_dup(factors[i]->mod);
    }

    // Generate all combinations of residues
    int sol_index = 0;
    int *indices = calloc(num_factors, sizeof(int));
    while (sol_index < total_solutions) {
        BIGNUM *sum = BN_new();
        BN_zero(sum);

        for (int i = 0; i < num_factors; i++) {
            Mi[i] = BN_new();
            BN_div(Mi[i], NULL, M, moduli[i], ctx);

            invMi[i] = BN_mod_inverse(NULL, Mi[i], moduli[i], ctx);
            if (!invMi[i]) {
                printf("Error computing inverse modulo.\n");
                exit(1);
            }

            BIGNUM *term = BN_new();
            // term = residue_i * invMi_i * Mi_i
            BN_mul(term, factors[i]->sols[indices[i]], invMi[i], ctx);
            BN_mul(term, term, Mi[i], ctx);

            BN_add(sum, sum, term);

            BN_free(term);
            BN_free(Mi[i]);
            BN_free(invMi[i]);
        }

        // Solution modulo M
        (*solutions)[sol_index] = BN_new();
        BN_mod((*solutions)[sol_index], sum, M, ctx);

        BN_free(sum);

        // Increment indices
        indices[0]++;
        for (int i = 0; i < num_factors; i++) {
            if (indices[i] >= factors[i]->num_sols && i + 1 < num_factors) {
                indices[i] = 0;
                indices[i + 1]++;
            }
        }

        sol_index++;
    }

    // Sort the solutions in ascending order
    qsort(*solutions, total_solutions, sizeof(BIGNUM *), compare_bn);

    // Free allocated memory
    for (int i = 0; i < num_factors; i++) {
        BN_free(moduli[i]);
    }
    free(moduli);
    free(Mi);
    free(invMi);
    free(indices);
    BN_free(M);
    BN_CTX_free(ctx);
}

void free_factors(PrimeFactor **factors, int num_factors) {
    for (int i = 0; i < num_factors; i++) {
        BN_free(factors[i]->p);
        BN_free(factors[i]->mod);
        for (int j = 0; j < factors[i]->num_sols; j++) {
            BN_free(factors[i]->sols[j]);
        }
        free(factors[i]->sols);
        free(factors[i]);
    }
}

void find_involutory_exponents(BIGNUM *p) {
    BIGNUM *m = BN_new();
    BN_CTX *ctx = BN_CTX_new();
    BN_sub(m, p, BN_value_one()); // m = p - 1

    PrimeFactor *factors[64];
    int num_factors = 0;

    // Factorize m
    factorize(m, factors, &num_factors);

    // Compute solutions using CRT
    BIGNUM **solutions;
    int num_solutions = 0;
    compute_solutions(factors, num_factors, &solutions, &num_solutions, m);

    // Print solutions
    printf("involutory exponents modulo ");
    char *m_str = BN_bn2dec(m);
    printf("%s:\n", m_str);
    OPENSSL_free(m_str);

    for (int i = 0; i < num_solutions; i++) {
        char *sol_str = BN_bn2hex(solutions[i]);
        printf("%d : %s\n", (i + 1), sol_str);
        OPENSSL_free(sol_str);
        BN_free(solutions[i]);
    }
    free(solutions);

    // Clean up
    free_factors(factors, num_factors);
    BN_free(m);
    BN_CTX_free(ctx);
}

int main(int argc, char *argv[]) {
    puts("\nFind involutory Exponents for Modular Exponentiation.\n");

    if (argc != 2) {
        printf("Usage: find_invol <modulus>\n");
        return 0;
    }

    BIGNUM *p = NULL;

    do {
        p = BN_new();

        if (!p) {
            printf("BN_new() failed...\n");
            break;
        }

        // Try as decimal first
        int result = BN_dec2bn(&p, argv[1]);

        if (!result) {
            // Then hexadecimal
            result = BN_hex2bn(&p, argv[1]);
            if (!result) {
                printf("Can't parse input.\n");
                break;
            }
        }

        // Find involutory exponents
        find_involutory_exponents(p);
    } while (0);

    // Clean up
    if (p) BN_free(p);
    return 0;
}