Z Under Two
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443
  1. // Z under two
  2. // Because what's some maths code without a bad french pun :)
  3. // >>=, <<=, <<, >>, +=, +, *=, *, /, %, &=, &, ==, !=, <, >, = are implemented
  4. class zut {
  5. unsigned int* number;
  6. int size;
  7. bool irreducible;
  8. public:
  9. // Initializers -------------------------------------------------
  10. zut(unsigned int* a, int b){
  11. number=a;
  12. size=b;
  13. irreducible=false;
  14. }
  15. zut(initializer_list<unsigned int> a){
  16. number=new unsigned int[a.size()];
  17. size=a.size()*32;
  18. copy(a.begin(), a.end(), number);
  19. irreducible=false;
  20. }
  21. zut(int b) {
  22. size=b;
  23. number = new unsigned int[b/32];
  24. irreducible=false;
  25. }
  26. zut(const zut &b){
  27. size=b.size;
  28. number = new unsigned int[b.size/32];
  29. (*this)=b;
  30. irreducible=b.irreducible;
  31. }
  32. void resize(int size){
  33. this->size=size;
  34. }
  35. // Primitives ---------------------------------------------------
  36. // https://stackoverflow.com/questions/2773890/efficient-bitshifting-an-array-of-int
  37. void operator<<=(int s){
  38. for (int i = s; i > 0; i -= (i % 32)) {
  39. if (!(i % 32)) {
  40. for (int y = (size/32)-1; y > 0; y--)
  41. number[y] = number[y - 1];
  42. number[0] = 0;
  43. i -= 32;
  44. } else {
  45. for (int y = (size/32)-1; y > 0; y--)
  46. number[y] = (number[y] << (i%32)) | ((number[y - 1] >> (32-i%32)) );
  47. number[0] <<= i%32;
  48. }
  49. }
  50. }
  51. zut operator<<(int s){
  52. zut res((*this));
  53. res <<= s;
  54. return res;
  55. }
  56. void operator>>=(int s){
  57. for (int i = s; i > 0; i -= (i % 32)) {
  58. if (!(i % 32)) {
  59. for (int y = (size/32)-1; y >0; y--)
  60. number[y] = number[y - 1];
  61. number[0] = 0;
  62. i -= 32;
  63. } else {
  64. for (int y = (size / 32) - 1; y > 0; y--)
  65. number[y] = (number[y] >> (i%32)) | ((number[y - 1] << (32 - (i%32)) & 1));
  66. number[0] >>= i%32;
  67. }
  68. }
  69. }
  70. zut operator>>(int s){
  71. zut res((*this));
  72. res >>= s;
  73. return res;
  74. }
  75. // i don't think this can be made faster :p
  76. // in a z/2z field 1+1=0 and 1+0 = 1, thus it's a simple xor
  77. // also i guess you know if you are reading this code but substraction =
  78. // addition in Z/2Z so this can be(and is, ie in division) used interchangeably
  79. void operator+=(const zut& b){
  80. for (int i = 0; i < size / 32; i++){
  81. if(!(i < b.size / 32))
  82. break;
  83. number[i] ^= b.number[i];
  84. }
  85. }
  86. zut operator+(const zut& b){
  87. zut res((*this));
  88. res+=b;
  89. return res;
  90. }
  91. // https://en.wikipedia.org/wiki/Finite_field_arithmetic#C_programming_example
  92. // Ugh, bit by bit comparison, it really looks like a bottleneck to me, but
  93. // I can't find a way to make it faster...
  94. void operator*=(zut& b){
  95. zut tmp(size);
  96. tmp=(*this);
  97. (*this)=0;
  98. for (int i = 0; i < size; i++) {
  99. if(!(i < b.size))
  100. break;
  101. // this verifies if the n-th bit inside b is lit
  102. if (b.get_bit(i))
  103. (*this)+=tmp;
  104. tmp <<=1;
  105. }
  106. }
  107. zut operator*(zut& b){
  108. zut res((*this));
  109. res *= b;
  110. return res;
  111. }
  112. // /!\ primitive might fail if 2 values have different size!!
  113. zut operator/(zut& b){
  114. zut r((*this));
  115. zut q(size); q=0;
  116. int deg_a=(*this).degree();
  117. int deg_div=deg_a-b.degree();
  118. for (int i = 0; i <= deg_div; i++) {
  119. q<<=1;
  120. int y=deg_a-i;
  121. if (r.get_bit(y)) {
  122. q++; r+=(b << (deg_div-i));
  123. }
  124. }
  125. return q;
  126. }
  127. zut operator%(zut& b){
  128. zut r((*this));
  129. zut q(size); q=0;
  130. int deg_a=(*this).degree();
  131. int deg_div=deg_a-b.degree();
  132. for (int i = 0; i <= deg_div; i++) {
  133. q<<=1;
  134. int y=deg_a-i;
  135. if (r.get_bit(y)) {
  136. q++; r+=(b << (deg_div-i));
  137. }
  138. }
  139. return r;
  140. }
  141. void operator&=(int b){
  142. for (int i = 0; i < size / 32; i++)
  143. number[i] &= b;
  144. }
  145. zut operator&(int b){
  146. zut res((*this));
  147. res &= b;
  148. return res;
  149. }
  150. void operator&=(const zut& b){
  151. for (int i = 0; i < size / 32; i++){
  152. if(!(i < b.size / 32))
  153. break;
  154. number[i] &= b.number[i];
  155. }
  156. }
  157. zut operator&(const zut& b){
  158. zut res((*this));
  159. res &= b;
  160. return res;
  161. }
  162. void operator++(int a){
  163. number[0]^=1;
  164. }
  165. // Logical operators --------------------------------------------
  166. bool operator==(const zut& b){
  167. for (int i = 0; i < size / 32; i++){
  168. if(!(i < b.size / 32))
  169. break;
  170. if(number[i]!=b.number[i])
  171. return false;
  172. }
  173. return true;
  174. }
  175. bool operator==(int b){
  176. if(number[0]!=b)
  177. return false;
  178. for (int i = 1; i < size / 32; i++)
  179. if(number[i]!=0)
  180. return false;
  181. return true;
  182. }
  183. bool operator!=(const zut& b){
  184. return !((*this)==b);
  185. }
  186. bool operator!=(int b){
  187. return !((*this)==b);
  188. }
  189. bool operator<(int a){
  190. return number[0]<a;
  191. }
  192. bool operator>(int a){
  193. return number[0]>a;
  194. }
  195. void operator=(const zut& b){
  196. for (int i = 0; i < b.size / 32; i++)
  197. number[i] = b.number[i];
  198. for (int i = b.size / 32; i < size / 32; i++)
  199. number[i] = 0;
  200. }
  201. void operator=(int b){
  202. number[0] = b;
  203. for (int i = 1; i < size / 32; i++)
  204. number[i] = 0;
  205. }
  206. // Algorithms ---------------------------------------------------
  207. // should be able to do this faster...
  208. zut get_bits(int start, int end) {
  209. zut res(size); res=0;
  210. for (int i = start; i < end; i++) {
  211. if ((*this).get_bit(i))
  212. res.number[(i-start)/32] |= (1 << (i-start)%32);
  213. }
  214. return res;
  215. }
  216. bool get_bit(int bit){
  217. if (((number[(bit / 32)] & (1 << bit % 32)) >> bit % 32) & 1)
  218. return true;
  219. return false;
  220. }
  221. /* So, how does this work?
  222. * Well to begin with we are dealing with a Z/2Z field and as such have an
  223. * equation in the form of sigma(ax^n), for n up to the maximum bit size of
  224. * our current number. Knowing a is either 0 or 1 in this field, it can also
  225. * be disregarded as d0=0. This means the derivative is in
  226. * the form of sigma(nx^(n-1)). This also means that any odd power is
  227. * completely disregarded as they would be, by definition, 0.
  228. * Thus we only need to lower the power of every power by one and to
  229. * drop x^n where n is odd.
  230. */
  231. zut derivative() {
  232. return ((*this)&0xaaaaaaaa)>>1; // -> 0b101010... -> 1+3 aka mask odd numbers
  233. }
  234. // https://stackoverflow.com/questions/671815/what-is-the-fastest-most-efficient-way-to-find-the-highest-set-bit-msb-in-an-i
  235. int degree() {
  236. for (int i = (size / 32) - 1; i >= 0; i--)
  237. if (number[i])
  238. return (i * 32) + (31 - __builtin_clz(number[i]));
  239. return 0;
  240. }
  241. // mostly same as derivative: sqrt(x^n) = x^(n/2) (except 1)thus we just shift everything
  242. zut sqrt(){
  243. zut sqrt(size);
  244. zut bit(size); bit=1;
  245. for (int i=0;i<size;i++,bit<<=1)
  246. sqrt+=((bit&(*this))>>(i>>1));
  247. return sqrt;
  248. }
  249. // aren't you supposed to learn what this is in middle school? ;)
  250. zut gcd(zut b) {
  251. zut a((*this));
  252. zut mod=a%b;
  253. while(mod>2){
  254. a=b;
  255. b=mod;
  256. mod=a%b;
  257. }
  258. if(mod==0)
  259. return b;
  260. return mod;
  261. }
  262. // https://en.wikipedia.org/wiki/Factorization_of_polynomials_over_finite_fields#Square-free_factorization
  263. // ^same as berlekamp, explains it better than i ever could
  264. vector<zut> sff(vector<zut> factors={}) {
  265. zut d=(*this).derivative();
  266. // step 2 of sff(wikipedia)
  267. if(d==0){
  268. factors=(*this).sqrt().sff();
  269. factors.insert(factors.end(), factors.begin(), factors.end());
  270. }
  271. else{
  272. zut gcd=(*this).gcd(d);
  273. if(gcd==1){
  274. factors.push_back((*this));
  275. return factors;
  276. }
  277. factors=gcd.sff(factors);
  278. factors=((*this)/gcd).sff(factors);
  279. }
  280. return factors;
  281. }
  282. /*
  283. * Honestly I was about to write a whole paragraph on how this works but,
  284. * honestly, wikipedia is a much better resource than me for that. I separated
  285. * and commented briefly each subsection of the algorithm. It's not _that_
  286. * complex plus, honestly, I have no idea who I'd write that for, as I'm sure
  287. * you already knew about that before making this problem.
  288. */
  289. vector<zut> bk(){
  290. // first create the matrix Q-I
  291. vector<zut> subalgebra;
  292. vector<zut> factors;
  293. int deg = (*this).degree();
  294. zut x_2n(size*2);
  295. zut fac(size*2); fac=(*this);
  296. zut di(size); di=1;
  297. vector<zut> mat; mat.reserve(deg);
  298. for(int i=0; i<deg;i++){
  299. x_2n=1; x_2n <<=i*2;
  300. zut row(size*2);
  301. // Q = (x_2n%fac)+di), I=di
  302. row=(((x_2n%fac)+di)<<(deg+1))+di; mat.push_back(row);
  303. di<<=1;
  304. }
  305. // then compute the row reduced echelon form of the matrix
  306. int row=0;
  307. for (int j = deg*2; j>=0; j--) {
  308. int pivrow = -1;
  309. for (int i = row; i < mat.size(); i++) {
  310. if (mat[i].get_bit(j)) {
  311. pivrow = i;
  312. break;
  313. }
  314. }
  315. if (pivrow != -1 ) {
  316. swap(mat[row], mat[pivrow]);
  317. for (int i = 0; i < mat.size(); i++) {
  318. if (i != row) {
  319. if (mat[i].get_bit(j)) {
  320. mat[i]+= mat[row];
  321. }
  322. }
  323. }
  324. row++;
  325. }
  326. }
  327. // then read off berlekamp subalgebra from the null space basis
  328. bool empty=true;
  329. for(zut &row: mat) {
  330. zut row_left = row.get_bits(deg+2, (deg*2)+2);
  331. if(row_left==0){
  332. zut row_right = row.get_bits(0, deg);
  333. if(row_right!=1){
  334. // bk subalgebra takes half of the matrix thus size/2
  335. row_right.resize(size);
  336. subalgebra.push_back(row_right);
  337. empty=false;
  338. }
  339. }
  340. }
  341. // if we had no correct null basis it is irreductible
  342. if(empty){
  343. (*this).irreducible=true; factors.push_back((*this)); return factors;
  344. }
  345. // otherwise calculates irreducible factor from berlekamp subalgebra
  346. factors.push_back((*this));
  347. for(zut &vector: subalgebra) {
  348. zut test((*this).gcd(vector));
  349. test=(*this)/test;
  350. for(int i=0; i<factors.size(); i++){
  351. if((factors[i]%test)==0){
  352. if((factors[i]/test)!=1){
  353. zut help((factors[i]/test));
  354. factors[i]=test;
  355. factors.push_back(help);
  356. empty=false;
  357. break;
  358. }
  359. }
  360. }
  361. }
  362. if(empty)
  363. factors[0].irreducible=true;
  364. return factors;
  365. }
  366. // we get square free factors, then recurse with berlekamp to get irreducible
  367. // forms of all possible factors
  368. vector<zut> factorize() {
  369. vector<zut> factors=(*this).sff();
  370. bool done=false;
  371. while(!done){
  372. for(int i=0; i< factors.size(); i++) {
  373. if(!factors[i].irreducible) {
  374. vector<zut> factored = factors[i].bk();
  375. factors.erase(factors.begin()+i);
  376. factors.insert(factors.end(), factored.begin(), factored.end());
  377. break;
  378. }
  379. if(i==factors.size()-1){
  380. done=true;
  381. }
  382. }
  383. if(factors.size()==0)
  384. done=true;
  385. }
  386. return factors;
  387. }
  388. }