#include <inttypes.h>
#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>

uint64_t *dec_to_bin(uint64_t d, uint64_t *length) {
    uint64_t *binary_form = calloc(100, sizeof(uint64_t));
    int index = 0;
    while (d != 0) {
        binary_form[index] = d % 2;
        d /= 2;
        index++;
    }

    *length = index;

    for (int i = index; i >= 0; i--) {
        printf("%ju", binary_form[i]);
    }
    printf(" index: %d\n", index);

    return binary_form;
}

uint64_t quick_pow(uint64_t *d_binary, uint64_t a, uint64_t n, uint64_t length) {
    uint64_t *powed = calloc(100, sizeof(uint64_t));

    powed[0] = a;
    for (int i = 1; i <= length; i++) {
        // powed[i] = (powed[i - 1] * powed[i - 1]) % n;
        powed[i] = (uint64_t)(((unsigned __int128)powed[i - 1] * powed[i - 1]) % n);
        // printf("powed: %ju, index: %d; ", powed[i], (i));
    }

    // check where in the binary are ones
    uint64_t multiplied = 1;
    for (int i = 0; i < length; i++) {
        if (d_binary[i] == 1) {
            // multiplied = (multiplied * powed[i]) % n;
            multiplied = (uint64_t)(((unsigned __int128)multiplied * powed[i]) % n);
        }
    }

    // printf("\nbm quick math: %ju; %ju ", multiplied, n);

    free(powed);

    return multiplied;
}

bool prime_test(uint64_t n, int a) {
    printf("\n\nprime test: %ju\n", n);
    // Miller Rabin prime test
    // choose a base: a, which should be a prime so that (n, a) = 1
    // then do 2 rounds of tests provided the first one did not fail
    // 1: a^d =k 1 mod n
    // 2: a^(d * 2^r) =k n-1 mod n
    // d = n-1 / 2^S (where S means how many time did we divide the number till we reached the first odd number)
    // S: see above
    // r = {0,... S-1}

    uint64_t d = n - 1;
    uint64_t S = 0;

    while (d % 2 == 0) {
        d = d / 2;
        S++;
        printf("%ju ", d);
    }
    printf("S: %ju", S);
    printf("\n");

    uint64_t r = S - 1; // this stores the number of elements from 0 to S-1

    // round 1
    // 1: a^d =k 1 mod n
    uint64_t length;
    uint64_t *d_binary = dec_to_bin(d, &length);
    uint64_t first_qp_res = 0;

    if ((first_qp_res = quick_pow(d_binary, a, n, length)) == 1) {
        free(d_binary);
        printf("true\n");
        return true;
    }

    // round 2
    //  2: a^(d * 2^r) =k n-1 mod n
    for (int i = 0; i <= r; i++) {
        if (first_qp_res == n - 1) {
            free(d_binary);
            printf("true\n");
            return true;
        } else {
            // first_qp_res = (first_qp_res * first_qp_res) % n;
            first_qp_res = (uint64_t)(((unsigned __int128)first_qp_res * first_qp_res) % n);
        }
    }

    printf("\nfalse\n\n");
    free(d_binary);
    return false;
}

int main() {
    prime_test(111, 5);
    prime_test(29, 2);
    prime_test(27, 2);
    prime_test(17, 2);
    prime_test(661, 2);
    prime_test(18446744073709551557UL, 2);
    prime_test(18446744073709551533UL, 3);
    return 0;
}