PageRenderTime 50ms CodeModel.GetById 21ms RepoModel.GetById 0ms app.codeStats 0ms

/cln-1.3.2/src/polynomial/elem/cl_UP_MI.h

#
C Header | 506 lines | 444 code | 35 blank | 27 comment | 96 complexity | cea3179322a68b56733904fe49200059 MD5 | raw file
Possible License(s): GPL-2.0
  1. // Univariate Polynomials over a ring of modular integers.
  2. #include "cln/GV_modinteger.h"
  3. #include "cln/modinteger.h"
  4. #include "cln/exception.h"
  5. namespace cln {
  6. // Assume a ring is a modint ring.
  7. inline cl_heap_modint_ring* TheModintRing (const cl_ring& R)
  8. { return (cl_heap_modint_ring*) R.heappointer; }
  9. // Normalize a vector: remove leading zero coefficients.
  10. // The result vector is known to have length len > 0.
  11. static inline void modint_normalize (cl_heap_modint_ring* R, cl_GV_MI& result, uintL len)
  12. {
  13. if (R->_zerop(result[len-1])) {
  14. len--;
  15. while (len > 0) {
  16. if (!R->_zerop(result[len-1]))
  17. break;
  18. len--;
  19. }
  20. var cl_GV_MI newresult = cl_GV_MI(len,R);
  21. #if 0
  22. for (var sintL i = len-1; i >= 0; i--)
  23. newresult[i] = result[i];
  24. #else
  25. cl_GV_MI::copy_elements(result,0,newresult,0,len);
  26. #endif
  27. result = newresult;
  28. }
  29. }
  30. static void modint_fprint (cl_heap_univpoly_ring* UPR, std::ostream& stream, const _cl_UP& x)
  31. {{
  32. DeclarePoly(cl_GV_MI,x);
  33. var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
  34. var sintL xlen = x.size();
  35. if (xlen == 0)
  36. fprint(stream, "0");
  37. else {
  38. var cl_string varname = get_varname(UPR);
  39. for (var sintL i = xlen-1; i >= 0; i--)
  40. if (!R->_zerop(x[i])) {
  41. if (i < xlen-1)
  42. fprint(stream, " + ");
  43. fprint(stream, "(");
  44. R->_fprint(stream, x[i]);
  45. fprint(stream, ")");
  46. if (i > 0) {
  47. fprint(stream, "*");
  48. fprint(stream, varname);
  49. if (i != 1) {
  50. fprint(stream, "^");
  51. fprintdecimal(stream, i);
  52. }
  53. }
  54. }
  55. }
  56. }}
  57. static bool modint_equal (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
  58. {{
  59. DeclarePoly(cl_GV_MI,x);
  60. DeclarePoly(cl_GV_MI,y);
  61. var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
  62. var sintL xlen = x.size();
  63. var sintL ylen = y.size();
  64. if (!(xlen == ylen))
  65. return false;
  66. for (var sintL i = xlen-1; i >= 0; i--)
  67. if (!R->_equal(x[i],y[i]))
  68. return false;
  69. return true;
  70. }}
  71. static const _cl_UP modint_zero (cl_heap_univpoly_ring* UPR)
  72. {
  73. return _cl_UP(UPR, cl_null_GV_I);
  74. }
  75. static bool modint_zerop (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
  76. {
  77. unused UPR;
  78. { DeclarePoly(cl_GV_MI,x);
  79. var sintL xlen = x.size();
  80. if (xlen == 0)
  81. return true;
  82. else
  83. return false;
  84. }}
  85. static const _cl_UP modint_plus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
  86. {{
  87. DeclarePoly(cl_GV_MI,x);
  88. DeclarePoly(cl_GV_MI,y);
  89. var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
  90. var sintL xlen = x.size();
  91. var sintL ylen = y.size();
  92. if (xlen == 0)
  93. return _cl_UP(UPR, y);
  94. if (ylen == 0)
  95. return _cl_UP(UPR, x);
  96. // Now xlen > 0, ylen > 0.
  97. if (xlen > ylen) {
  98. var cl_GV_MI result = cl_GV_MI(xlen,R);
  99. var sintL i;
  100. #if 0
  101. for (i = xlen-1; i >= ylen; i--)
  102. result[i] = x[i];
  103. #else
  104. cl_GV_MI::copy_elements(x,ylen,result,ylen,xlen-ylen);
  105. #endif
  106. for (i = ylen-1; i >= 0; i--)
  107. result[i] = R->_plus(x[i],y[i]);
  108. return _cl_UP(UPR, result);
  109. }
  110. if (xlen < ylen) {
  111. var cl_GV_MI result = cl_GV_MI(ylen,R);
  112. var sintL i;
  113. #if 0
  114. for (i = ylen-1; i >= xlen; i--)
  115. result[i] = y[i];
  116. #else
  117. cl_GV_MI::copy_elements(y,xlen,result,xlen,ylen-xlen);
  118. #endif
  119. for (i = xlen-1; i >= 0; i--)
  120. result[i] = R->_plus(x[i],y[i]);
  121. return _cl_UP(UPR, result);
  122. }
  123. // Now xlen = ylen > 0. Add and normalize simultaneously.
  124. for (var sintL i = xlen-1; i >= 0; i--) {
  125. var _cl_MI hicoeff = R->_plus(x[i],y[i]);
  126. if (!R->_zerop(hicoeff)) {
  127. var cl_GV_MI result = cl_GV_MI(i+1,R);
  128. result[i] = hicoeff;
  129. for (i-- ; i >= 0; i--)
  130. result[i] = R->_plus(x[i],y[i]);
  131. return _cl_UP(UPR, result);
  132. }
  133. }
  134. return _cl_UP(UPR, cl_null_GV_I);
  135. }}
  136. static const _cl_UP modint_uminus (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
  137. {{
  138. DeclarePoly(cl_GV_MI,x);
  139. var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
  140. var sintL xlen = x.size();
  141. if (xlen == 0)
  142. return _cl_UP(UPR, x);
  143. // Now xlen > 0.
  144. // Negate. No normalization necessary, since the degree doesn't change.
  145. var sintL i = xlen-1;
  146. var _cl_MI hicoeff = R->_uminus(x[i]);
  147. if (R->_zerop(hicoeff)) throw runtime_exception();
  148. var cl_GV_MI result = cl_GV_MI(xlen,R);
  149. result[i] = hicoeff;
  150. for (i-- ; i >= 0; i--)
  151. result[i] = R->_uminus(x[i]);
  152. return _cl_UP(UPR, result);
  153. }}
  154. static const _cl_UP modint_minus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
  155. {{
  156. DeclarePoly(cl_GV_MI,x);
  157. DeclarePoly(cl_GV_MI,y);
  158. var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
  159. var sintL xlen = x.size();
  160. var sintL ylen = y.size();
  161. if (ylen == 0)
  162. return _cl_UP(UPR, x);
  163. if (xlen == 0)
  164. return modint_uminus(UPR, _cl_UP(UPR, y));
  165. // Now xlen > 0, ylen > 0.
  166. if (xlen > ylen) {
  167. var cl_GV_MI result = cl_GV_MI(xlen,R);
  168. var sintL i;
  169. #if 0
  170. for (i = xlen-1; i >= ylen; i--)
  171. result[i] = x[i];
  172. #else
  173. cl_GV_MI::copy_elements(x,ylen,result,ylen,xlen-ylen);
  174. #endif
  175. for (i = ylen-1; i >= 0; i--)
  176. result[i] = R->_minus(x[i],y[i]);
  177. return _cl_UP(UPR, result);
  178. }
  179. if (xlen < ylen) {
  180. var cl_GV_MI result = cl_GV_MI(ylen,R);
  181. var sintL i;
  182. for (i = ylen-1; i >= xlen; i--)
  183. result[i] = R->_uminus(y[i]);
  184. for (i = xlen-1; i >= 0; i--)
  185. result[i] = R->_minus(x[i],y[i]);
  186. return _cl_UP(UPR, result);
  187. }
  188. // Now xlen = ylen > 0. Add and normalize simultaneously.
  189. for (var sintL i = xlen-1; i >= 0; i--) {
  190. var _cl_MI hicoeff = R->_minus(x[i],y[i]);
  191. if (!R->_zerop(hicoeff)) {
  192. var cl_GV_MI result = cl_GV_MI(i+1,R);
  193. result[i] = hicoeff;
  194. for (i-- ; i >= 0; i--)
  195. result[i] = R->_minus(x[i],y[i]);
  196. return _cl_UP(UPR, result);
  197. }
  198. }
  199. return _cl_UP(UPR, cl_null_GV_I);
  200. }}
  201. static const _cl_UP modint_one (cl_heap_univpoly_ring* UPR)
  202. {
  203. var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
  204. var cl_GV_MI result = cl_GV_MI(1,R);
  205. result[0] = R->_one();
  206. return _cl_UP(UPR, result);
  207. }
  208. static const _cl_UP modint_canonhom (cl_heap_univpoly_ring* UPR, const cl_I& x)
  209. {
  210. var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
  211. var cl_GV_MI result = cl_GV_MI(1,R);
  212. result[0] = R->_canonhom(x);
  213. return _cl_UP(UPR, result);
  214. }
  215. static const _cl_UP modint_mul (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
  216. {{
  217. DeclarePoly(cl_GV_MI,x);
  218. DeclarePoly(cl_GV_MI,y);
  219. var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
  220. var sintL xlen = x.size();
  221. var sintL ylen = y.size();
  222. if (xlen == 0)
  223. return _cl_UP(UPR, x);
  224. if (ylen == 0)
  225. return _cl_UP(UPR, y);
  226. // Multiply.
  227. var sintL len = xlen + ylen - 1;
  228. var cl_GV_MI result = cl_GV_MI(len,R);
  229. if (xlen < ylen) {
  230. {
  231. var sintL i = xlen-1;
  232. var _cl_MI xi = x[i];
  233. for (sintL j = ylen-1; j >= 0; j--)
  234. result[i+j] = R->_mul(xi,y[j]);
  235. }
  236. for (sintL i = xlen-2; i >= 0; i--) {
  237. var _cl_MI xi = x[i];
  238. for (sintL j = ylen-1; j > 0; j--)
  239. result[i+j] = R->_plus(result[i+j],R->_mul(xi,y[j]));
  240. /* j=0 */ result[i] = R->_mul(xi,y[0]);
  241. }
  242. } else {
  243. {
  244. var sintL j = ylen-1;
  245. var _cl_MI yj = y[j];
  246. for (sintL i = xlen-1; i >= 0; i--)
  247. result[i+j] = R->_mul(x[i],yj);
  248. }
  249. for (sintL j = ylen-2; j >= 0; j--) {
  250. var _cl_MI yj = y[j];
  251. for (sintL i = xlen-1; i > 0; i--)
  252. result[i+j] = R->_plus(result[i+j],R->_mul(x[i],yj));
  253. /* i=0 */ result[j] = R->_mul(x[0],yj);
  254. }
  255. }
  256. // Normalize (not necessary in integral domains).
  257. //modint_normalize(R,result,len);
  258. if (R->_zerop(result[len-1])) throw runtime_exception();
  259. return _cl_UP(UPR, result);
  260. }}
  261. static const _cl_UP modint_square (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
  262. {{
  263. DeclarePoly(cl_GV_MI,x);
  264. var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
  265. var sintL xlen = x.size();
  266. if (xlen == 0)
  267. return cl_UP(UPR, x);
  268. var sintL len = 2*xlen-1;
  269. var cl_GV_MI result = cl_GV_MI(len,R);
  270. if (xlen > 1) {
  271. // Loop through all 0 <= j < i <= xlen-1.
  272. {
  273. var sintL i = xlen-1;
  274. var _cl_MI xi = x[i];
  275. for (sintL j = i-1; j >= 0; j--)
  276. result[i+j] = R->_mul(xi,x[j]);
  277. }
  278. {for (sintL i = xlen-2; i >= 1; i--) {
  279. var _cl_MI xi = x[i];
  280. for (sintL j = i-1; j >= 1; j--)
  281. result[i+j] = R->_plus(result[i+j],R->_mul(xi,x[j]));
  282. /* j=0 */ result[i] = R->_mul(xi,x[0]);
  283. }}
  284. // Double.
  285. {for (sintL i = len-2; i >= 1; i--)
  286. result[i] = R->_plus(result[i],result[i]);
  287. }
  288. // Add squares.
  289. result[2*(xlen-1)] = R->_square(x[xlen-1]);
  290. for (sintL i = xlen-2; i >= 1; i--)
  291. result[2*i] = R->_plus(result[2*i],R->_square(x[i]));
  292. }
  293. result[0] = R->_square(x[0]);
  294. // Normalize (not necessary in integral domains).
  295. //modint_normalize(R,result,len);
  296. if (R->_zerop(result[len-1])) throw runtime_exception();
  297. return _cl_UP(UPR, result);
  298. }}
  299. static const _cl_UP modint_exptpos (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const cl_I& y)
  300. {
  301. var _cl_UP a = x;
  302. var cl_I b = y;
  303. while (!oddp(b)) { a = UPR->_square(a); b = b >> 1; }
  304. var _cl_UP c = a;
  305. until (b == 1)
  306. { b = b >> 1;
  307. a = UPR->_square(a);
  308. if (oddp(b)) { c = UPR->_mul(a,c); }
  309. }
  310. return c;
  311. }
  312. static const _cl_UP modint_scalmul (cl_heap_univpoly_ring* UPR, const cl_ring_element& x, const _cl_UP& y)
  313. {
  314. if (!(UPR->basering() == x.ring())) throw runtime_exception();
  315. {
  316. DeclarePoly(_cl_MI,x);
  317. DeclarePoly(cl_GV_MI,y);
  318. var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
  319. var sintL ylen = y.size();
  320. if (ylen == 0)
  321. return _cl_UP(UPR, y);
  322. if (R->_zerop(x))
  323. return _cl_UP(UPR, cl_null_GV_I);
  324. // Now ylen > 0.
  325. // No normalization necessary, since the degree doesn't change.
  326. var cl_GV_MI result = cl_GV_MI(ylen,R);
  327. for (sintL i = ylen-1; i >= 0; i--)
  328. result[i] = R->_mul(x,y[i]);
  329. return _cl_UP(UPR, result);
  330. }}
  331. static sintL modint_degree (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
  332. {
  333. unused UPR;
  334. { DeclarePoly(cl_GV_MI,x);
  335. return (sintL) x.size() - 1;
  336. }}
  337. static sintL modint_ldegree (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
  338. {{
  339. DeclarePoly(cl_GV_MI,x);
  340. var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
  341. var sintL xlen = x.size();
  342. for (sintL i = 0; i < xlen; i++) {
  343. if (!R->_zerop(x[i]))
  344. return i;
  345. }
  346. return -1;
  347. }}
  348. static const _cl_UP modint_monomial (cl_heap_univpoly_ring* UPR, const cl_ring_element& x, uintL e)
  349. {
  350. if (!(UPR->basering() == x.ring())) throw runtime_exception();
  351. { DeclarePoly(_cl_MI,x);
  352. var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
  353. if (R->_zerop(x))
  354. return _cl_UP(UPR, cl_null_GV_I);
  355. else {
  356. var sintL len = e+1;
  357. var cl_GV_MI result = cl_GV_MI(len,R);
  358. result[e] = x;
  359. return _cl_UP(UPR, result);
  360. }
  361. }}
  362. static const cl_ring_element modint_coeff (cl_heap_univpoly_ring* UPR, const _cl_UP& x, uintL index)
  363. {{
  364. DeclarePoly(cl_GV_MI,x);
  365. var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
  366. if (index < x.size())
  367. return cl_MI(R, x[index]);
  368. else
  369. return R->zero();
  370. }}
  371. static const _cl_UP modint_create (cl_heap_univpoly_ring* UPR, sintL deg)
  372. {
  373. if (deg < 0)
  374. return _cl_UP(UPR, cl_null_GV_I);
  375. else {
  376. var sintL len = deg+1;
  377. var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
  378. return _cl_UP(UPR, cl_GV_MI(len,R));
  379. }
  380. }
  381. static void modint_set_coeff (cl_heap_univpoly_ring* UPR, _cl_UP& x, uintL index, const cl_ring_element& y)
  382. {{
  383. DeclareMutablePoly(cl_GV_MI,x);
  384. if (!(UPR->basering() == y.ring())) throw runtime_exception();
  385. { DeclarePoly(_cl_MI,y);
  386. if (!(index < x.size())) throw runtime_exception();
  387. x[index] = y;
  388. }}}
  389. static void modint_finalize (cl_heap_univpoly_ring* UPR, _cl_UP& x)
  390. {{
  391. DeclareMutablePoly(cl_GV_MI,x); // NB: x is modified by reference!
  392. var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
  393. var uintL len = x.size();
  394. if (len > 0)
  395. modint_normalize(R,x,len);
  396. }}
  397. static const cl_ring_element modint_eval (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const cl_ring_element& y)
  398. {{
  399. // Method:
  400. // If x = 0, return 0.
  401. // If y = 0, return x[0].
  402. // Else compute (...(x[len-1]*y+x[len-2])*y ...)*y + x[0].
  403. DeclarePoly(cl_GV_MI,x);
  404. if (!(UPR->basering() == y.ring())) throw runtime_exception();
  405. { DeclarePoly(_cl_MI,y);
  406. var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
  407. var uintL len = x.size();
  408. if (len==0)
  409. return R->zero();
  410. if (R->_zerop(y))
  411. return cl_MI(R, x[0]);
  412. var sintL i = len-1;
  413. var _cl_MI z = x[i];
  414. for ( ; --i >= 0; )
  415. z = R->_plus(R->_mul(z,y),x[i]);
  416. return cl_MI(R, z);
  417. }}}
  418. static cl_univpoly_setops modint_setops = {
  419. modint_fprint,
  420. modint_equal
  421. };
  422. static cl_univpoly_addops modint_addops = {
  423. modint_zero,
  424. modint_zerop,
  425. modint_plus,
  426. modint_minus,
  427. modint_uminus
  428. };
  429. static cl_univpoly_mulops modint_mulops = {
  430. modint_one,
  431. modint_canonhom,
  432. modint_mul,
  433. modint_square,
  434. modint_exptpos
  435. };
  436. static cl_univpoly_modulops modint_modulops = {
  437. modint_scalmul
  438. };
  439. static cl_univpoly_polyops modint_polyops = {
  440. modint_degree,
  441. modint_ldegree,
  442. modint_monomial,
  443. modint_coeff,
  444. modint_create,
  445. modint_set_coeff,
  446. modint_finalize,
  447. modint_eval
  448. };
  449. class cl_heap_modint_univpoly_ring : public cl_heap_univpoly_ring {
  450. SUBCLASS_cl_heap_univpoly_ring()
  451. public:
  452. // Constructor.
  453. cl_heap_modint_univpoly_ring (const cl_ring& r);
  454. // Destructor.
  455. ~cl_heap_modint_univpoly_ring () {}
  456. };
  457. static void cl_heap_modint_univpoly_ring_destructor (cl_heap* pointer)
  458. {
  459. (*(cl_heap_modint_univpoly_ring*)pointer).~cl_heap_modint_univpoly_ring();
  460. }
  461. cl_class cl_class_modint_univpoly_ring = {
  462. cl_heap_modint_univpoly_ring_destructor,
  463. cl_class_flags_univpoly_ring
  464. };
  465. // Constructor.
  466. inline cl_heap_modint_univpoly_ring::cl_heap_modint_univpoly_ring (const cl_ring& r)
  467. : cl_heap_univpoly_ring (r, &modint_setops, &modint_addops, &modint_mulops, &modint_modulops, &modint_polyops)
  468. {
  469. type = &cl_class_modint_univpoly_ring;
  470. }
  471. } // namespace cln