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

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 const uint64_t GX[4] = {
    0xF4A13945D898C296ULL, 0x77037D812DEB33A0ULL,
    0xF8BCE6E563A440F2ULL, 0x6B17D1F2E12C4247ULL
};

int main() {
    uint64_t a[4], b[4];
    for(int i=0;i<4;i++) a[i]=GX[i];
    for(int i=0;i<4;i++) b[i]=GX[i];

    uint64_t tmp[8]={0};
    for(int i=0;i<4;i++){
        printf("\n--- i=%d (a[%d]=0x%016lx) ---\n", i, i, a[i]);
        uint64_t carry=0;
        for(int j=0;j<4;j++){
            uint64_t lo;
            uint64_t hi = mul64(a[i], b[j], &lo);
            printf("  j=%d: mul64(a[%d], b[%d], &lo) => hi=0x%016lx lo=0x%016lx\n", j, i, j, hi, lo);
            printf("    tmp[%d] before: 0x%016lx; carry before: 0x%016lx\n", i+j, tmp[i+j], carry);
            
            uint64_t sum = tmp[i+j] + lo;
            uint64_t c1 = (sum < lo) ? 1 : 0;
            printf("    sum1 (tmp+lo) = 0x%016lx, c1=%lu\n", sum, c1);
            
            sum += carry;
            c1 |= (sum < carry) ? 1 : 0;
            printf("    sum2 (+carry) = 0x%016lx, c1=%lu\n", sum, c1);
            
            tmp[i+j] = sum;
            carry = hi + c1;
            printf("    tmp[%d] after: 0x%016lx, new carry=0x%016lx\n", i+j, tmp[i+j], carry);
        }
        tmp[i+4]=carry;
        printf("  tmp[%d]=carry=0x%016lx\n", i+4, tmp[i+4]);
        
        printf("  Current tmp: ");
        for(int k=7;k>=0;k--) printf("%016lx", tmp[k]); printf("\n");
    }
    
    printf("\n--- Final tmp ---\n");
    for(int k=7;k>=0;k--) printf("tmp[%d] = 0x%016lx\n", k, tmp[k]);
    
    /* Expected values from Python */
    uint64_t expected[8] = {
        0x30edd460c607afe4ULL, 0x1ffd1a5bd36cd015ULL,
        0xc89de54fe3a13ecaULL, 0xb8d5b1119b64b5c2ULL,
        0x70387d127a032f22ULL, 0xd92ade46207ebdf8ULL,
        0xd4d2435551cd0b7cULL, 0x2cccebb86e0def74ULL
    };
    printf("\n--- Expected ---\n");
    for(int k=7;k>=0;k--) printf("tmp[%d] = 0x%016lx\n", k, expected[k]);
    
    printf("\n--- Mismatches ---\n");
    for(int k=0;k<8;k++) {
        if(tmp[k] != expected[k]) {
            printf("tmp[%d]: got 0x%016lx, expected 0x%016lx (diff = %ld)\n", 
                   k, tmp[k], expected[k], (int64_t)(tmp[k] - expected[k]));
        }
    }
    
    return 0;
}
