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

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

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);}
}

/* The correct 512-bit product of GX*GX from Python */
static const uint64_t correct_product[8] = {
    0x30edd460c607afe4ULL, 0x1ffd1a5bd36cd015ULL,
    0xc89de54fe3a13ecaULL, 0xb8d5b1119b64b5c2ULL,
    0x70387d127a032f22ULL, 0xd92ade46207ebdf8ULL,
    0xd4d2435551cd0b7cULL, 0x2cccebb86e0def74ULL
};

static void print_limbs32(const char *label, int64_t *v, int n) {
    printf("%s: ", label);
    for(int i=n-1;i>=0;i--) {
        if(v[i] < 0) printf("-%08llx", (unsigned long long)(-v[i]));
        else printf("+%08llx", (unsigned long long)v[i]);
        if(i>0) printf(" ");
    }
    printf("\n");
}

static void print_hex32(const char *label, uint32_t *v, int n) {
    printf("%s: ", label);
    for(int i=n-1;i>=0;i--) printf("%08x", v[i]);
    printf("\n");
}

int main() {
    uint64_t tmp[8];
    for(int i=0;i<8;i++) tmp[i] = correct_product[i];
    
    /* 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("Input product:\n");
    print_hex32("t", t, 16);
    
    /* NIST P-256 reduction */
    int64_t s[15];
    for(int i=0;i<15;i++) s[i]=0;
    for(int i=0;i<8;i++) s[i]=t[i];
    
    printf("\nInitial s[0..7] = low words\n");
    
    /* First fold: high words t[8..15] */
    for(int k=8;k<16;k++){
        printf("Fold t[%d]=0x%08x: s[%d]+=, s[%d]-=, s[%d]-=, s[%d]+=\n",
               k, t[k], k-8, k-5, k-2, k-1);
        s[k-8]+=t[k]; s[k-5]-=t[k]; s[k-2]-=t[k]; s[k-1]+=t[k];
    }
    
    printf("\nAfter first fold:\n");
    for(int i=0;i<15;i++) printf("  s[%d] = 0x%016llx (%lld)\n", i, (unsigned long long)s[i], (long long)s[i]);
    
    /* Second fold */
    for(int k=8;k<15;k++){
        if(s[k]==0) continue;
        int64_t v=s[k]; s[k]=0;
        printf("\nFold s[%d]=%lld: s[%d]+=, s[%d]-=, s[%d]-=, s[%d]+=\n", k, (long long)v, k-8, k-5, k-2, k-1);
        s[k-8]+=v; s[k-5]-=v; s[k-2]-=v; s[k-1]+=v;
    }
    
    /* Third fold */
    for(int k=8;k<15;k++){
        if(s[k]==0) continue;
        int64_t v=s[k]; s[k]=0;
        printf("\nFold2 s[%d]=%lld: s[%d]+=, s[%d]-=, s[%d]-=, s[%d]+=\n", k, (long long)v, k-8, k-5, k-2, k-1);
        s[k-8]+=v; s[k-5]-=v; s[k-2]-=v; s[k-1]+=v;
    }
    
    printf("\nAfter folds:\n");
    for(int i=0;i<15;i++) printf("  s[%d] = 0x%016llx (%lld)\n", i, (unsigned long long)s[i], (long long)s[i]);
    
    /* Carry propagation */
    for(int i=0;i<14;i++){
        int64_t carry=s[i]>>32;
        s[i]&=0xFFFFFFFF;
        s[i+1]+=carry;
        printf("Propagate s[%d]: carry=%lld, s[%d]=0x%08llx, s[%d]=0x%016llx\n", 
               i, (long long)carry, i, (unsigned long long)s[i], i+1, (unsigned long long)s[i+1]);
    }
    
    printf("\nAfter carry propagation:\n");
    for(int i=0;i<15;i++) printf("  s[%d] = 0x%016llx (%lld)\n", i, (unsigned long long)s[i], (long long)s[i]);
    
    /* Fold high words from carry propagation */
    for(int k=8;k<15;k++){
        if(s[k]==0) continue;
        int64_t v=s[k]; s[k]=0;
        printf("Fold3 s[%d]=%lld: s[%d]+=, s[%d]-=, s[%d]-=, s[%d]+=\n", k, (long long)v, k-8, k-5, k-2, k-1);
        s[k-8]+=v; s[k-5]-=v; s[k-2]-=v; s[k-1]+=v;
    }
    
    /* Final carry propagation */
    for(int i=0;i<14;i++){
        int64_t carry=s[i]>>32;
        s[i]&=0xFFFFFFFF;
        s[i+1]+=carry;
    }
    /* Fold any remaining high words after final carry propagation */
    for(int k=8;k<15;k++){
        if(s[k]==0) continue;
        int64_t v=s[k]; s[k]=0;
        s[k-8]+=v; s[k-5]-=v; s[k-2]-=v; s[k-1]+=v;
    }
    /* Final propagation after all folds */
    for(int i=0;i<14;i++){
        int64_t carry=s[i]>>32;
        s[i]&=0xFFFFFFFF;
        s[i+1]+=carry;
    }

    printf("\nFinal s[0..7]:\n");
    for(int i=0;i<8;i++) printf("  s[%d] = %08llx\n", i, (unsigned long long)(s[i] & 0xFFFFFFFF));
    
    uint64_t result[4];
    for(int i=0;i<4;i++)
        result[i] = (uint64_t)s[2*i] | ((uint64_t)s[2*i+1] << 32);
    
    bn_mod_p256(result);
    
    printf("\nFinal result:\n");
    for(int i=3;i>=0;i--) printf("  result[%d] = 0x%016lx\n", i, (unsigned long)result[i]);
    
    uint64_t expected[4] = {0x002ae56c426b3f8c, 0x33b699495d694dd1, 0x81819a5e0e3690d8, 0x98f6b84d29bef2b2};
    printf("\nExpected:\n");
    for(int i=3;i>=0;i--) printf("  expected[%d] = 0x%016lx\n", i, (unsigned long)expected[i]);
    
    printf("\nMatch: %s\n", memcmp(result, expected, 32) == 0 ? "PASS" : "FAIL");
    
    return 0;
}
