#include <stdio.h>
#define BIT 16
long I;
long MODULO;
long OMEGA;
unsigned long compToNum(short x, short y) {
unsigned long n = (x + y * I + MODULO) % MODULO;
return n;
}
short numToRe(unsigned long nRaw) {
unsigned long n = nRaw % MODULO;
// if (n == OMEGA) {
// return;
// }
unsigned long geta = MODULO;
unsigned long nGeta = (n < OMEGA) ? (n + geta) : n;
short reRaw = (nGeta - 1) % I;
short re = (reRaw >= (I / 2)) ? (reRaw - I) : reRaw;
return re;
}
short numToIm(unsigned long nRaw) {
unsigned long n = nRaw % MODULO;
// if (n == OMEGA) {
// return;
// }
unsigned long geta = MODULO;
unsigned long nGeta = (n < OMEGA) ? (n + geta) : n;
short imRaw = ((nGeta - 1 + I / 2) >> BIT);
short im = imRaw - I;
return im;
}
int main(void) {
I = 1 << BIT;
MODULO = I * I + 1;
OMEGA = (1 << (BIT * 2 - 1)) - (1 << (BIT - 1));
unsigned long c1;
unsigned long c2;
unsigned long c;
short re1 = 3;
short im1 = 3;
short re2 = 0;
short im2 = 5;
c1 = compToNum(re1, im1);
c2 = compToNum(re2, im2);
c = c1 * c2;
short re = numToRe(c);
short im = numToIm(c);
printf("MODULO: %ld\n", MODULO
); printf("OMEGA: %ld\n", OMEGA
); printf("c1: %ld: (%d, %d)\n", c1
, numToRe
(c1
), numToIm
(c1
)); printf("c2: %ld: (%d, %d)\n", c2
, numToRe
(c2
), numToIm
(c2
)); printf("c1 * c2: %ld: (%d, %d)", c
, re
, im
);
return 0;
}
I2luY2x1ZGUgPHN0ZGlvLmg+CgojZGVmaW5lIEJJVCAxNgoKbG9uZyBJOwpsb25nIE1PRFVMTzsKbG9uZyBPTUVHQTsKCnVuc2lnbmVkIGxvbmcgY29tcFRvTnVtKHNob3J0IHgsIHNob3J0IHkpIHsKCXVuc2lnbmVkIGxvbmcgbiA9ICh4ICsgeSAqIEkgKyBNT0RVTE8pICUgTU9EVUxPOwoKCXJldHVybiBuOwp9CgpzaG9ydCBudW1Ub1JlKHVuc2lnbmVkIGxvbmcgblJhdykgewoJdW5zaWduZWQgbG9uZyBuID0gblJhdyAlIE1PRFVMTzsKCgkvLyBpZiAobiA9PSBPTUVHQSkgewoJLy8JCXJldHVybjsKCS8vIH0KCgl1bnNpZ25lZCBsb25nIGdldGEgPSBNT0RVTE87CgoJdW5zaWduZWQgbG9uZyBuR2V0YSA9IChuIDwgT01FR0EpID8gKG4gKyBnZXRhKSA6IG47CgoJc2hvcnQgcmVSYXcgPSAobkdldGEgLSAxKSAlIEk7CgoJc2hvcnQgcmUgPSAocmVSYXcgPj0gKEkgLyAyKSkgPyAocmVSYXcgLSBJKSA6IHJlUmF3OwoKCXJldHVybiByZTsKfQoKc2hvcnQgbnVtVG9JbSh1bnNpZ25lZCBsb25nIG5SYXcpIHsKCXVuc2lnbmVkIGxvbmcgbiA9IG5SYXcgJSBNT0RVTE87CgoJLy8gaWYgKG4gPT0gT01FR0EpIHsKCS8vCQlyZXR1cm47CgkvLyB9CgoJdW5zaWduZWQgbG9uZyBnZXRhID0gTU9EVUxPOwoKCXVuc2lnbmVkIGxvbmcgbkdldGEgPSAobiA8IE9NRUdBKSA/IChuICsgZ2V0YSkgOiBuOwoKCXNob3J0IGltUmF3ID0gKChuR2V0YSAtIDEgKyBJIC8gMikgPj4gQklUKTsKCglzaG9ydCBpbSA9IGltUmF3IC0gSTsKCglyZXR1cm4gaW07Cn0KCmludCBtYWluKHZvaWQpIHsKCUkgPSAxIDw8IEJJVDsKCU1PRFVMTyA9IEkgKiBJICsgMTsKCU9NRUdBID0gKDEgPDwgKEJJVCAqIDIgLSAxKSkgLSAoMSA8PCAoQklUIC0gMSkpOwoKCXVuc2lnbmVkIGxvbmcgYzE7Cgl1bnNpZ25lZCBsb25nIGMyOwoJdW5zaWduZWQgbG9uZyBjOwoKCXNob3J0IHJlMSA9IDM7CglzaG9ydCBpbTEgPSAzOwoJc2hvcnQgcmUyID0gMDsKCXNob3J0IGltMiA9IDU7CgoJYzEgPSBjb21wVG9OdW0ocmUxLCBpbTEpOwoJYzIgPSBjb21wVG9OdW0ocmUyLCBpbTIpOwoJYyA9IGMxICogYzI7CgoJc2hvcnQgcmUgPSBudW1Ub1JlKGMpOwoJc2hvcnQgaW0gPSBudW1Ub0ltKGMpOwoKCXByaW50ZigiQklUOiAlZFxuIiwgQklUKTsKCXByaW50ZigiSTogJWxkXG4iLCBJKTsKCXByaW50ZigiTU9EVUxPOiAlbGRcbiIsIE1PRFVMTyk7CglwcmludGYoIk9NRUdBOiAlbGRcbiIsIE9NRUdBKTsKCXByaW50ZigiYzE6ICVsZDogKCVkLCAlZClcbiIsIGMxLCBudW1Ub1JlKGMxKSwgbnVtVG9JbShjMSkpOwoJcHJpbnRmKCJjMjogJWxkOiAoJWQsICVkKVxuIiwgYzIsIG51bVRvUmUoYzIpLCBudW1Ub0ltKGMyKSk7CglwcmludGYoImMxICogYzI6ICVsZDogKCVkLCAlZCkiLCBjLCByZSwgaW0pOwoKCXJldHVybiAwOwp9Cg==