#include <stdio.h>
#include <string.h>
#include <stdint.h>

static const uint64_t P256_P[4] = {
    0xFFFFFFFFFFFFFFFFULL, 0x00000000FFFFFFFFULL,
    0x0000000000000000ULL, 0xFFFFFFFF00000001ULL
};

static const uint64_t P256_GX[4] = {
    0xF4A13945D898C296ULL, 0x77037D812DEB33A0ULL,
    0xF8BCE6E563A440F2ULL, 0x6B17D1F2E12C4247ULL
};

static void bn_copy(uint64_t *a, const uint64_t *b)
{ a[0]=b[0]; a[1]=b[1]; a[2]=b[2]; a[3]=b[3]; }

static int bn_lt(const uint64_t *a, const uint64_t *b)
{
    for(int i=3;i>=0;i--) { if(a[i]<b[i]) return 1; if(a[i]>b[i]) return 0; }
    return 0;
}

static uint64_t bn_sub(uint64_t *a, const uint64_t *b)
{
    uint64_t borrow=0;
    for(int i=0;i<4;i++){
        uint64_t diff = a[i] - b[i] - borrow;
        uint64_t new_borrow = (a[i] < b[i]) ? 1 : ((a[i] == b[i]) ? borrow : 0);
        borrow = new_borrow;
        a[i] = diff;
    }
    return borrow;
}

static void bn_mod_p256(uint64_t *a)
{
    while(bn_lt(P256_P,a)){bn_sub(a,P256_P);}
}

static uint64_t mul64(uint64_t a, uint64_t b, uint64_t *lo)
{
    uint32_t a_lo = (uint32_t)a, a_hi = (uint32_t)(a >> 32);
    uint32_t b_lo = (uint32_t)b, b_hi = (uint32_t)(b >> 32);
    uint64_t ll = (uint64_t)a_lo * b_lo;
    uint64_t lh = (uint64_t)a_lo * b_hi;
    uint64_t hl = (uint64_t)a_hi * b_lo;
    uint64_t hh = (uint64_t)a_hi * b_hi;
    uint64_t mid = (ll >> 32) + (uint32_t)lh + (uint32_t)hl;
    *lo = (mid << 32) | (uint32_t)ll;
    return hh + (mid >> 32) + (lh >> 32) + (hl >> 32);
}

static void print_limbs(const char *label, uint64_t *v, int n) {
    printf("%s: ", label);
    for(int i=n-1;i>=0;i--) printf("%016lx", v[i]);
    printf("\n");
}

int main() {
    uint64_t a[4], b[4];
    bn_copy(a, P256_GX);
    bn_copy(b, P256_GX);
    
    /* Full schoolbook product */
    uint64_t tmp[8]={0};
    for(int i=0;i<4;i++){
        uint64_t carry=0;
        for(int j=0;j<4;j++){
            uint64_t lo;
            uint64_t hi = mul64(a[i], b[j], &lo);
            uint64_t sum = tmp[i+j] + lo;
            uint64_t c1 = (sum < lo) ? 1 : 0;
            sum += carry;
            c1 += (sum < carry) ? 1 : 0;
            tmp[i+j] = sum;
            carry = hi + c1;
        }
        tmp[i+4]=carry;
    }
    
    printf("Full 512-bit product of GX*GX:\n");
    print_limbs("tmp", tmp, 8);
    
    /* Convert to 32-bit words */
    uint32_t t[16];
    for(int i=0;i<4;i++){
        t[2*i]   = (uint32_t)(tmp[i]);
        t[2*i+1] = (uint32_t)(tmp[i]>>32);
    }
    for(int i=0;i<4;i++){
        t[8+2*i]   = (uint32_t)(tmp[4+i]);
        t[8+2*i+1] = (uint32_t)(tmp[4+i]>>32);
    }
    
    printf("\n32-bit words:\n");
    for(int i=15;i>=0;i--) printf("%08x", t[i]); printf("\n");
    
    /* Verify using Python computation */
    printf("\nExpected product from Python:\n");
    printf("(awaiting input)\n");
    
    /* Test with simpler numbers to isolate the issue */
    printf("\n--- Testing with known values ---\n");
    {
        uint64_t small_a[4] = {5, 0, 0, 0};
        uint64_t small_b[4] = {7, 0, 0, 0};
        
        uint64_t small_tmp[8]={0};
        for(int i=0;i<4;i++){
            uint64_t carry=0;
            for(int j=0;j<4;j++){
                uint64_t lo;
                uint64_t hi = mul64(small_a[i], small_b[j], &lo);
                uint64_t sum = small_tmp[i+j] + lo;
                uint64_t c1 = (sum < lo) ? 1 : 0;
                sum += carry;
                c1 |= (sum < carry) ? 1 : 0;
                small_tmp[i+j] = sum;
                carry = hi + c1;
            }
            small_tmp[i+4]=carry;
        }
        printf("5*7 product:"); print_limbs("", small_tmp, 8);
    }
    
    /* Test: Gx * 1 */
    {
        uint64_t one[4] = {1, 0, 0, 0};
        uint64_t one_tmp[8]={0};
        for(int i=0;i<4;i++){
            uint64_t carry=0;
            for(int j=0;j<4;j++){
                uint64_t lo;
                uint64_t hi = mul64(P256_GX[i], one[j], &lo);
                uint64_t sum = one_tmp[i+j] + lo;
                uint64_t c1 = (sum < lo) ? 1 : 0;
                sum += carry;
                c1 |= (sum < carry) ? 1 : 0;
                one_tmp[i+j] = sum;
                carry = hi + c1;
            }
            one_tmp[i+4]=carry;
        }
        printf("GX*1 product:"); print_limbs("", one_tmp, 8);
    }
    
    return 0;
}
