LongArray.cs 38 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295
  1. #if !BESTHTTP_DISABLE_ALTERNATE_SSL && (!UNITY_WEBGL || UNITY_EDITOR)
  2. #pragma warning disable
  3. using System;
  4. using System.Text;
  5. using Best.HTTP.SecureProtocol.Org.BouncyCastle.Math.Raw;
  6. using Best.HTTP.SecureProtocol.Org.BouncyCastle.Utilities;
  7. namespace Best.HTTP.SecureProtocol.Org.BouncyCastle.Math.EC
  8. {
  9. internal struct LongArray
  10. {
  11. internal static bool AreAliased(ref LongArray a, ref LongArray b)
  12. {
  13. return a.m_data == b.m_data;
  14. }
  15. // TODO make m fixed for the LongArray, and hence compute T once and for all
  16. private ulong[] m_data;
  17. internal LongArray(int intLen)
  18. {
  19. m_data = new ulong[intLen];
  20. }
  21. internal LongArray(ulong[] data)
  22. {
  23. m_data = data;
  24. }
  25. internal LongArray(ulong[] data, int off, int len)
  26. {
  27. if (off == 0 && len == data.Length)
  28. {
  29. m_data = data;
  30. }
  31. else
  32. {
  33. m_data = new ulong[len];
  34. Array.Copy(data, off, m_data, 0, len);
  35. }
  36. }
  37. internal LongArray(BigInteger bigInt)
  38. {
  39. if (bigInt == null || bigInt.SignValue < 0)
  40. throw new ArgumentException("invalid F2m field value", nameof(bigInt));
  41. if (bigInt.SignValue == 0)
  42. {
  43. m_data = new ulong[1]{ 0UL };
  44. return;
  45. }
  46. byte[] barr = bigInt.ToByteArray();
  47. int barrLen = barr.Length;
  48. int barrStart = 0;
  49. if (barr[0] == 0)
  50. {
  51. // First byte is 0 to enforce highest (=sign) bit is zero.
  52. // In this case ignore barr[0].
  53. barrLen--;
  54. barrStart = 1;
  55. }
  56. int intLen = (barrLen + 7) / 8;
  57. m_data = new ulong[intLen];
  58. int iarrJ = intLen - 1;
  59. int rem = barrLen % 8 + barrStart;
  60. ulong temp = 0;
  61. int barrI = barrStart;
  62. if (barrStart < rem)
  63. {
  64. for (; barrI < rem; barrI++)
  65. {
  66. temp <<= 8;
  67. uint barrBarrI = barr[barrI];
  68. temp |= barrBarrI;
  69. }
  70. m_data[iarrJ--] = temp;
  71. }
  72. for (; iarrJ >= 0; iarrJ--)
  73. {
  74. temp = 0;
  75. for (int i = 0; i < 8; i++)
  76. {
  77. temp <<= 8;
  78. uint barrBarrI = barr[barrI++];
  79. temp |= barrBarrI;
  80. }
  81. m_data[iarrJ] = temp;
  82. }
  83. }
  84. internal void CopyTo(ulong[] z, int zOff)
  85. {
  86. Array.Copy(m_data, 0, z, zOff, m_data.Length);
  87. }
  88. internal bool IsOne()
  89. {
  90. ulong[] a = m_data;
  91. int aLen = a.Length;
  92. if (aLen < 1 || a[0] != 1UL)
  93. return false;
  94. for (int i = 1; i < aLen; ++i)
  95. {
  96. if (a[i] != 0UL)
  97. return false;
  98. }
  99. return true;
  100. }
  101. internal bool IsZero()
  102. {
  103. ulong[] a = m_data;
  104. for (int i = 0; i < a.Length; ++i)
  105. {
  106. if (a[i] != 0UL)
  107. return false;
  108. }
  109. return true;
  110. }
  111. internal int GetUsedLength()
  112. {
  113. return GetUsedLengthFrom(m_data.Length);
  114. }
  115. internal int GetUsedLengthFrom(int from)
  116. {
  117. ulong[] a = m_data;
  118. from = System.Math.Min(from, a.Length);
  119. if (from < 1)
  120. return 0;
  121. // Check if first element will act as sentinel
  122. if (a[0] != 0UL)
  123. {
  124. while (a[--from] == 0UL)
  125. {
  126. }
  127. return from + 1;
  128. }
  129. do
  130. {
  131. if (a[--from] != 0UL)
  132. {
  133. return from + 1;
  134. }
  135. }
  136. while (from > 0);
  137. return 0;
  138. }
  139. internal int Degree()
  140. {
  141. int i = m_data.Length;
  142. ulong w;
  143. do
  144. {
  145. if (i == 0)
  146. return 0;
  147. w = m_data[--i];
  148. }
  149. while (w == 0UL);
  150. return (i << 6) + BitLength(w);
  151. }
  152. private int DegreeFrom(int limit)
  153. {
  154. int i = (int)(((uint)limit + 62) >> 6);
  155. ulong w;
  156. do
  157. {
  158. if (i == 0)
  159. return 0;
  160. w = m_data[--i];
  161. }
  162. while (w == 0);
  163. return (i << 6) + BitLength(w);
  164. }
  165. private static int BitLength(ulong w)
  166. {
  167. return 64 - Longs.NumberOfLeadingZeros((long)w);
  168. }
  169. private ulong[] ResizedData(int newLen)
  170. {
  171. ulong[] newInts = new ulong[newLen];
  172. Array.Copy(m_data, 0, newInts, 0, System.Math.Min(m_data.Length, newLen));
  173. return newInts;
  174. }
  175. internal BigInteger ToBigInteger()
  176. {
  177. int usedLen = GetUsedLength();
  178. if (usedLen == 0)
  179. return BigInteger.Zero;
  180. ulong highestInt = m_data[usedLen - 1];
  181. byte[] temp = new byte[8];
  182. int barrI = 0;
  183. bool trailingZeroBytesDone = false;
  184. for (int j = 7; j >= 0; j--)
  185. {
  186. byte thisByte = (byte)(highestInt >> (8 * j));
  187. if (trailingZeroBytesDone || (thisByte != 0))
  188. {
  189. trailingZeroBytesDone = true;
  190. temp[barrI++] = thisByte;
  191. }
  192. }
  193. int barrLen = 8 * (usedLen - 1) + barrI;
  194. byte[] barr = new byte[barrLen];
  195. for (int j = 0; j < barrI; j++)
  196. {
  197. barr[j] = temp[j];
  198. }
  199. // Highest value int is done now
  200. for (int iarrJ = usedLen - 2; iarrJ >= 0; iarrJ--)
  201. {
  202. ulong mi = m_data[iarrJ];
  203. for (int j = 7; j >= 0; j--)
  204. {
  205. barr[barrI++] = (byte)(mi >> (8 * j));
  206. }
  207. }
  208. return new BigInteger(1, barr);
  209. }
  210. private static ulong ShiftUp(ulong[] x, int xOff, int count, int shift)
  211. {
  212. int shiftInv = 64 - shift;
  213. ulong prev = 0UL;
  214. for (int i = 0; i < count; ++i)
  215. {
  216. ulong next = x[xOff + i];
  217. x[xOff + i] = (next << shift) | prev;
  218. prev = next >> shiftInv;
  219. }
  220. return prev;
  221. }
  222. private static ulong ShiftUp(ulong[] x, int xOff, ulong[] z, int zOff, int count, int shift)
  223. {
  224. int shiftInv = 64 - shift;
  225. ulong prev = 0UL;
  226. for (int i = 0; i < count; ++i)
  227. {
  228. ulong next = x[xOff + i];
  229. z[zOff + i] = (next << shift) | prev;
  230. prev = next >> shiftInv;
  231. }
  232. return prev;
  233. }
  234. internal LongArray AddOne()
  235. {
  236. if (m_data.Length == 0)
  237. return new LongArray(new ulong[1]{ 1UL });
  238. int resultLen = System.Math.Max(1, GetUsedLength());
  239. ulong[] data = ResizedData(resultLen);
  240. data[0] ^= 1UL;
  241. return new LongArray(data);
  242. }
  243. private void AddShiftedByBitsSafe(LongArray other, int otherDegree, int bits)
  244. {
  245. int otherLen = (int)((uint)(otherDegree + 63) >> 6);
  246. int words = (int)((uint)bits >> 6);
  247. int shift = bits & 0x3F;
  248. if (shift == 0)
  249. {
  250. Add(m_data, words, other.m_data, 0, otherLen);
  251. return;
  252. }
  253. ulong carry = AddShiftedUp(m_data, words, other.m_data, 0, otherLen, shift);
  254. if (carry != 0UL)
  255. {
  256. m_data[otherLen + words] ^= carry;
  257. }
  258. }
  259. private static ulong AddShiftedUp(ulong[] x, int xOff, ulong[] y, int yOff, int count, int shift)
  260. {
  261. int shiftInv = 64 - shift;
  262. ulong prev = 0;
  263. for (int i = 0; i < count; ++i)
  264. {
  265. ulong next = y[yOff + i];
  266. x[xOff + i] ^= (next << shift) | prev;
  267. prev = next >> shiftInv;
  268. }
  269. return prev;
  270. }
  271. private static ulong AddShiftedDown(ulong[] x, int xOff, ulong[] y, int yOff, int count, int shift)
  272. {
  273. int shiftInv = 64 - shift;
  274. ulong prev = 0;
  275. int i = count;
  276. while (--i >= 0)
  277. {
  278. ulong next = y[yOff + i];
  279. x[xOff + i] ^= (next >> shift) | prev;
  280. prev = next << shiftInv;
  281. }
  282. return prev;
  283. }
  284. internal void AddShiftedByWords(LongArray other, int words)
  285. {
  286. int otherUsedLen = other.GetUsedLength();
  287. if (otherUsedLen == 0)
  288. return;
  289. int minLen = otherUsedLen + words;
  290. if (minLen > m_data.Length)
  291. {
  292. m_data = ResizedData(minLen);
  293. }
  294. Add(m_data, words, other.m_data, 0, otherUsedLen);
  295. }
  296. private static void Add(ulong[] x, int xOff, ulong[] y, int yOff, int count)
  297. {
  298. Nat.XorTo64(count, y, yOff, x, xOff);
  299. }
  300. private static void Add(ulong[] x, int xOff, ulong[] y, int yOff, ulong[] z, int zOff, int count)
  301. {
  302. Nat.Xor64(count, x, xOff, y, yOff, z, zOff);
  303. }
  304. private static void AddBoth(ulong[] x, int xOff, ulong[] y1, int y1Off, ulong[] y2, int y2Off, int count)
  305. {
  306. for (int i = 0; i < count; ++i)
  307. {
  308. x[xOff + i] ^= y1[y1Off + i] ^ y2[y2Off + i];
  309. }
  310. }
  311. private static void FlipWord(ulong[] buf, int off, int bit, ulong word)
  312. {
  313. int n = off + (int)((uint)bit >> 6);
  314. int shift = bit & 0x3F;
  315. if (shift == 0)
  316. {
  317. buf[n] ^= word;
  318. }
  319. else
  320. {
  321. buf[n] ^= word << shift;
  322. word = word >> (64 - shift);
  323. if (word != 0)
  324. {
  325. buf[++n] ^= word;
  326. }
  327. }
  328. }
  329. internal bool TestBitZero()
  330. {
  331. return m_data.Length > 0 && (m_data[0] & 1UL) != 0;
  332. }
  333. private static bool TestBit(ulong[] buf, int off, int n)
  334. {
  335. // theInt = n / 64
  336. int theInt = (int)((uint)n >> 6);
  337. // theBit = n % 64
  338. int theBit = n & 0x3F;
  339. ulong tester = 1UL << theBit;
  340. return (buf[off + theInt] & tester) != 0UL;
  341. }
  342. private static void FlipBit(ulong[] buf, int off, int n)
  343. {
  344. // theInt = n / 64
  345. int theInt = (int)((uint)n >> 6);
  346. // theBit = n % 64
  347. int theBit = n & 0x3F;
  348. ulong flipper = 1UL << theBit;
  349. buf[off + theInt] ^= flipper;
  350. }
  351. private static void MultiplyWord(ulong a, ulong[] b, int bLen, ulong[] c, int cOff)
  352. {
  353. if ((a & 1UL) != 0UL)
  354. {
  355. Add(c, cOff, b, 0, bLen);
  356. }
  357. int k = 1;
  358. while ((a >>= 1) != 0UL)
  359. {
  360. if ((a & 1UL) != 0UL)
  361. {
  362. ulong carry = AddShiftedUp(c, cOff, b, 0, bLen, k);
  363. if (carry != 0UL)
  364. {
  365. c[cOff + bLen] ^= carry;
  366. }
  367. }
  368. ++k;
  369. }
  370. }
  371. internal LongArray ModMultiplyLD(LongArray other, int m, int[] ks)
  372. {
  373. /*
  374. * Find out the degree of each argument and handle the zero cases
  375. */
  376. int aDeg = Degree();
  377. if (aDeg == 0)
  378. return this;
  379. int bDeg = other.Degree();
  380. if (bDeg == 0)
  381. return other;
  382. /*
  383. * Swap if necessary so that A is the smaller argument
  384. */
  385. LongArray A = this, B = other;
  386. if (aDeg > bDeg)
  387. {
  388. A = other; B = this;
  389. int tmp = aDeg; aDeg = bDeg; bDeg = tmp;
  390. }
  391. /*
  392. * Establish the word lengths of the arguments and result
  393. */
  394. int aLen = (int)((uint)(aDeg + 63) >> 6);
  395. int bLen = (int)((uint)(bDeg + 63) >> 6);
  396. int cLen = (int)((uint)(aDeg + bDeg + 62) >> 6);
  397. if (aLen == 1)
  398. {
  399. ulong a0 = A.m_data[0];
  400. if (a0 == 1UL)
  401. return B;
  402. /*
  403. * Fast path for small A, with performance dependent only on the number of set bits
  404. */
  405. ulong[] c0 = new ulong[cLen];
  406. MultiplyWord(a0, B.m_data, bLen, c0, 0);
  407. /*
  408. * Reduce the raw answer against the reduction coefficients
  409. */
  410. return ReduceResult(c0, 0, cLen, m, ks);
  411. }
  412. /*
  413. * Determine if B will get bigger during shifting
  414. */
  415. int bMax = (int)((uint)(bDeg + 7 + 63) >> 6);
  416. /*
  417. * Lookup table for the offset of each B in the tables
  418. */
  419. int[] ti = new int[16];
  420. /*
  421. * Precompute table of all 4-bit products of B
  422. */
  423. ulong[] T0 = new ulong[bMax << 4];
  424. int tOff = bMax;
  425. ti[1] = tOff;
  426. Array.Copy(B.m_data, 0, T0, tOff, bLen);
  427. for (int i = 2; i < 16; ++i)
  428. {
  429. ti[i] = (tOff += bMax);
  430. if ((i & 1) == 0)
  431. {
  432. ShiftUp(T0, (int)((uint)tOff >> 1), T0, tOff, bMax, 1);
  433. }
  434. else
  435. {
  436. Add(T0, bMax, T0, tOff - bMax, T0, tOff, bMax);
  437. }
  438. }
  439. /*
  440. * Second table with all 4-bit products of B shifted 4 bits
  441. */
  442. ulong[] T1 = new ulong[T0.Length];
  443. ShiftUp(T0, 0, T1, 0, T0.Length, 4);
  444. // shiftUp(T0, bMax, T1, bMax, tOff, 4);
  445. ulong[] a = A.m_data;
  446. ulong[] c = new ulong[cLen];
  447. uint MASK = 0xF;
  448. /*
  449. * Lopez-Dahab algorithm
  450. */
  451. for (int k = 56; k >= 0; k -= 8)
  452. {
  453. for (int j = 1; j < aLen; j += 2)
  454. {
  455. uint aVal = (uint)(a[j] >> k);
  456. uint u = aVal & MASK;
  457. uint v = (aVal >> 4) & MASK;
  458. AddBoth(c, j - 1, T0, ti[u], T1, ti[v], bMax);
  459. }
  460. ShiftUp(c, 0, cLen, 8);
  461. }
  462. for (int k = 56; k >= 0; k -= 8)
  463. {
  464. for (int j = 0; j < aLen; j += 2)
  465. {
  466. uint aVal = (uint)(a[j] >> k);
  467. uint u = aVal & MASK;
  468. uint v = (aVal >> 4) & MASK;
  469. AddBoth(c, j, T0, ti[u], T1, ti[v], bMax);
  470. }
  471. if (k > 0)
  472. {
  473. ShiftUp(c, 0, cLen, 8);
  474. }
  475. }
  476. /*
  477. * Finally the raw answer is collected, reduce it against the reduction coefficients
  478. */
  479. return ReduceResult(c, 0, cLen, m, ks);
  480. }
  481. internal LongArray ModMultiply(LongArray other, int m, int[] ks)
  482. {
  483. /*
  484. * Find out the degree of each argument and handle the zero cases
  485. */
  486. int aDeg = Degree();
  487. if (aDeg == 0)
  488. return this;
  489. int bDeg = other.Degree();
  490. if (bDeg == 0)
  491. return other;
  492. /*
  493. * Swap if necessary so that A is the smaller argument
  494. */
  495. LongArray A = this, B = other;
  496. if (aDeg > bDeg)
  497. {
  498. A = other; B = this;
  499. int tmp = aDeg; aDeg = bDeg; bDeg = tmp;
  500. }
  501. /*
  502. * Establish the word lengths of the arguments and result
  503. */
  504. int aLen = (int)((uint)(aDeg + 63) >> 6);
  505. int bLen = (int)((uint)(bDeg + 63) >> 6);
  506. int cLen = (int)((uint)(aDeg + bDeg + 62) >> 6);
  507. if (aLen == 1)
  508. {
  509. ulong a0 = A.m_data[0];
  510. if (a0 == 1UL)
  511. return B;
  512. /*
  513. * Fast path for small A, with performance dependent only on the number of set bits
  514. */
  515. ulong[] c0 = new ulong[cLen];
  516. MultiplyWord(a0, B.m_data, bLen, c0, 0);
  517. /*
  518. * Reduce the raw answer against the reduction coefficients
  519. */
  520. return ReduceResult(c0, 0, cLen, m, ks);
  521. }
  522. /*
  523. * Determine if B will get bigger during shifting
  524. */
  525. int bMax = (int)((uint)(bDeg + 7 + 63) >> 6);
  526. /*
  527. * Lookup table for the offset of each B in the tables
  528. */
  529. int[] ti = new int[16];
  530. /*
  531. * Precompute table of all 4-bit products of B
  532. */
  533. ulong[] T0 = new ulong[bMax << 4];
  534. int tOff = bMax;
  535. ti[1] = tOff;
  536. Array.Copy(B.m_data, 0, T0, tOff, bLen);
  537. for (int i = 2; i < 16; ++i)
  538. {
  539. tOff += bMax;
  540. ti[i] = tOff;
  541. if ((i & 1) == 0)
  542. {
  543. ShiftUp(T0, (int)((uint)tOff >> 1), T0, tOff, bMax, 1);
  544. }
  545. else
  546. {
  547. Add(T0, bMax, T0, tOff - bMax, T0, tOff, bMax);
  548. }
  549. }
  550. /*
  551. * Second table with all 4-bit products of B shifted 4 bits
  552. */
  553. ulong[] T1 = new ulong[T0.Length];
  554. ShiftUp(T0, 0, T1, 0, T0.Length, 4);
  555. // ShiftUp(T0, bMax, T1, bMax, tOff, 4);
  556. ulong[] a = A.m_data;
  557. ulong[] c = new ulong[cLen << 3];
  558. uint MASK = 0xF;
  559. /*
  560. * Lopez-Dahab (Modified) algorithm
  561. */
  562. for (int aPos = 0; aPos < aLen; ++aPos)
  563. {
  564. ulong aVal = a[aPos];
  565. int cOff = aPos;
  566. for (;;)
  567. {
  568. uint u = (uint)aVal & MASK; aVal >>= 4;
  569. uint v = (uint)aVal & MASK; aVal >>= 4;
  570. AddBoth(c, cOff, T0, ti[u], T1, ti[v], bMax);
  571. if (aVal == 0UL)
  572. break;
  573. cOff += cLen;
  574. }
  575. }
  576. {
  577. int cOff = c.Length;
  578. while ((cOff -= cLen) != 0)
  579. {
  580. AddShiftedUp(c, cOff - cLen, c, cOff, cLen, 8);
  581. }
  582. }
  583. /*
  584. * Finally the raw answer is collected, reduce it against the reduction coefficients
  585. */
  586. return ReduceResult(c, 0, cLen, m, ks);
  587. }
  588. //internal LongArray ModReduce(int m, int[] ks)
  589. //{
  590. // ulong[] buf = Arrays.Clone(m_data);
  591. // int rLen = ReduceInPlace(buf, 0, buf.Length, m, ks);
  592. // return new LongArray(buf, 0, rLen);
  593. //}
  594. internal LongArray Multiply(LongArray other, int m, int[] ks)
  595. {
  596. /*
  597. * Find out the degree of each argument and handle the zero cases
  598. */
  599. int aDeg = Degree();
  600. if (aDeg == 0)
  601. return this;
  602. int bDeg = other.Degree();
  603. if (bDeg == 0)
  604. return other;
  605. /*
  606. * Swap if necessary so that A is the smaller argument
  607. */
  608. LongArray A = this, B = other;
  609. if (aDeg > bDeg)
  610. {
  611. A = other; B = this;
  612. int tmp = aDeg; aDeg = bDeg; bDeg = tmp;
  613. }
  614. /*
  615. * Establish the word lengths of the arguments and result
  616. */
  617. int aLen = (int)((uint)(aDeg + 63) >> 6);
  618. int bLen = (int)((uint)(bDeg + 63) >> 6);
  619. int cLen = (int)((uint)(aDeg + bDeg + 62) >> 6);
  620. if (aLen == 1)
  621. {
  622. ulong a0 = A.m_data[0];
  623. if (a0 == 1UL)
  624. return B;
  625. /*
  626. * Fast path for small A, with performance dependent only on the number of set bits
  627. */
  628. ulong[] c0 = new ulong[cLen];
  629. MultiplyWord(a0, B.m_data, bLen, c0, 0);
  630. /*
  631. * Reduce the raw answer against the reduction coefficients
  632. */
  633. //return ReduceResult(c0, 0, cLen, m, ks);
  634. return new LongArray(c0, 0, cLen);
  635. }
  636. /*
  637. * Determine if B will get bigger during shifting
  638. */
  639. int bMax = (int)((uint)(bDeg + 7 + 63) >> 6);
  640. /*
  641. * Lookup table for the offset of each B in the tables
  642. */
  643. int[] ti = new int[16];
  644. /*
  645. * Precompute table of all 4-bit products of B
  646. */
  647. ulong[] T0 = new ulong[bMax << 4];
  648. int tOff = bMax;
  649. ti[1] = tOff;
  650. Array.Copy(B.m_data, 0, T0, tOff, bLen);
  651. for (int i = 2; i < 16; ++i)
  652. {
  653. tOff += bMax;
  654. ti[i] = tOff;
  655. if ((i & 1) == 0)
  656. {
  657. ShiftUp(T0, (int)((uint)tOff >> 1), T0, tOff, bMax, 1);
  658. }
  659. else
  660. {
  661. Add(T0, bMax, T0, tOff - bMax, T0, tOff, bMax);
  662. }
  663. }
  664. /*
  665. * Second table with all 4-bit products of B shifted 4 bits
  666. */
  667. ulong[] T1 = new ulong[T0.Length];
  668. ShiftUp(T0, 0, T1, 0, T0.Length, 4);
  669. //ShiftUp(T0, bMax, T1, bMax, tOff, 4);
  670. ulong[] a = A.m_data;
  671. ulong[] c = new ulong[cLen << 3];
  672. uint MASK = 0xF;
  673. /*
  674. * Lopez-Dahab (Modified) algorithm
  675. */
  676. for (int aPos = 0; aPos < aLen; ++aPos)
  677. {
  678. ulong aVal = a[aPos];
  679. int cOff = aPos;
  680. for (; ; )
  681. {
  682. uint u = (uint)aVal & MASK; aVal >>= 4;
  683. uint v = (uint)aVal & MASK; aVal >>= 4;
  684. AddBoth(c, cOff, T0, ti[u], T1, ti[v], bMax);
  685. if (aVal == 0UL)
  686. break;
  687. cOff += cLen;
  688. }
  689. }
  690. {
  691. int cOff = c.Length;
  692. while ((cOff -= cLen) != 0)
  693. {
  694. AddShiftedUp(c, cOff - cLen, c, cOff, cLen, 8);
  695. }
  696. }
  697. /*
  698. * Finally the raw answer is collected, reduce it against the reduction coefficients
  699. */
  700. //return ReduceResult(c, 0, cLen, m, ks);
  701. return new LongArray(c, 0, cLen);
  702. }
  703. internal void Reduce(int m, int[] ks)
  704. {
  705. ulong[] buf = m_data;
  706. int rLen = ReduceInPlace(buf, 0, buf.Length, m, ks);
  707. if (rLen < buf.Length)
  708. {
  709. m_data = new ulong[rLen];
  710. Array.Copy(buf, 0, m_data, 0, rLen);
  711. }
  712. }
  713. private static LongArray ReduceResult(ulong[] buf, int off, int len, int m, int[] ks)
  714. {
  715. int rLen = ReduceInPlace(buf, off, len, m, ks);
  716. return new LongArray(buf, off, rLen);
  717. }
  718. private static int ReduceInPlace(ulong[] buf, int off, int len, int m, int[] ks)
  719. {
  720. int mLen = (m + 63) >> 6;
  721. if (len < mLen)
  722. return len;
  723. int numBits = System.Math.Min(len << 6, (m << 1) - 1); // TODO use actual degree?
  724. int excessBits = (len << 6) - numBits;
  725. while (excessBits >= 64)
  726. {
  727. --len;
  728. excessBits -= 64;
  729. }
  730. int kLen = ks.Length, kMax = ks[kLen - 1], kNext = kLen > 1 ? ks[kLen - 2] : 0;
  731. int wordWiseLimit = System.Math.Max(m, kMax + 64);
  732. int vectorableWords = (excessBits + System.Math.Min(numBits - wordWiseLimit, m - kNext)) >> 6;
  733. if (vectorableWords > 1)
  734. {
  735. int vectorWiseWords = len - vectorableWords;
  736. ReduceVectorWise(buf, off, len, vectorWiseWords, m, ks);
  737. while (len > vectorWiseWords)
  738. {
  739. buf[off + --len] = 0L;
  740. }
  741. numBits = vectorWiseWords << 6;
  742. }
  743. if (numBits > wordWiseLimit)
  744. {
  745. ReduceWordWise(buf, off, len, wordWiseLimit, m, ks);
  746. numBits = wordWiseLimit;
  747. }
  748. if (numBits > m)
  749. {
  750. ReduceBitWise(buf, off, numBits, m, ks);
  751. }
  752. return mLen;
  753. }
  754. private static void ReduceBitWise(ulong[] buf, int off, int BitLength, int m, int[] ks)
  755. {
  756. while (--BitLength >= m)
  757. {
  758. if (TestBit(buf, off, BitLength))
  759. {
  760. ReduceBit(buf, off, BitLength, m, ks);
  761. }
  762. }
  763. }
  764. private static void ReduceBit(ulong[] buf, int off, int bit, int m, int[] ks)
  765. {
  766. FlipBit(buf, off, bit);
  767. int n = bit - m;
  768. int j = ks.Length;
  769. while (--j >= 0)
  770. {
  771. FlipBit(buf, off, ks[j] + n);
  772. }
  773. FlipBit(buf, off, n);
  774. }
  775. private static void ReduceWordWise(ulong[] buf, int off, int len, int toBit, int m, int[] ks)
  776. {
  777. int toPos = (int)((uint)toBit >> 6);
  778. while (--len > toPos)
  779. {
  780. ulong word = buf[off + len];
  781. if (word != 0)
  782. {
  783. buf[off + len] = 0UL;
  784. ReduceWord(buf, off, (len << 6), word, m, ks);
  785. }
  786. }
  787. {
  788. int partial = toBit & 0x3F;
  789. ulong word = buf[off + toPos] >> partial;
  790. if (word != 0)
  791. {
  792. buf[off + toPos] ^= word << partial;
  793. ReduceWord(buf, off, toBit, word, m, ks);
  794. }
  795. }
  796. }
  797. private static void ReduceWord(ulong[] buf, int off, int bit, ulong word, int m, int[] ks)
  798. {
  799. int offset = bit - m;
  800. int j = ks.Length;
  801. while (--j >= 0)
  802. {
  803. FlipWord(buf, off, offset + ks[j], word);
  804. }
  805. FlipWord(buf, off, offset, word);
  806. }
  807. private static void ReduceVectorWise(ulong[] buf, int off, int len, int words, int m, int[] ks)
  808. {
  809. /*
  810. * NOTE: It's important we go from highest coefficient to lowest, because for the highest
  811. * one (only) we allow the ranges to partially overlap, and therefore any changes must take
  812. * effect for the subsequent lower coefficients.
  813. */
  814. int baseBit = (words << 6) - m;
  815. int j = ks.Length;
  816. while (--j >= 0)
  817. {
  818. FlipVector(buf, off, buf, off + words, len - words, baseBit + ks[j]);
  819. }
  820. FlipVector(buf, off, buf, off + words, len - words, baseBit);
  821. }
  822. private static void FlipVector(ulong[] x, int xOff, ulong[] y, int yOff, int yLen, int bits)
  823. {
  824. xOff += (int)((uint)bits >> 6);
  825. bits &= 0x3F;
  826. if (bits == 0)
  827. {
  828. Add(x, xOff, y, yOff, yLen);
  829. }
  830. else
  831. {
  832. ulong carry = AddShiftedDown(x, xOff + 1, y, yOff, yLen, 64 - bits);
  833. x[xOff] ^= carry;
  834. }
  835. }
  836. internal LongArray ModSquare(int m, int[] ks)
  837. {
  838. int len = GetUsedLength();
  839. if (len == 0)
  840. return this;
  841. ulong[] r = new ulong[len << 1];
  842. Interleave.Expand64To128(m_data, 0, len, r, 0);
  843. return new LongArray(r, 0, ReduceInPlace(r, 0, r.Length, m, ks));
  844. }
  845. internal LongArray ModSquareN(int n, int m, int[] ks)
  846. {
  847. int len = GetUsedLength();
  848. if (len == 0)
  849. return this;
  850. int mLen = (m + 63) >> 6;
  851. ulong[] r = new ulong[mLen << 1];
  852. Array.Copy(m_data, 0, r, 0, len);
  853. while (--n >= 0)
  854. {
  855. Interleave.Expand64To128(r, 0, len, r, 0);
  856. len = ReduceInPlace(r, 0, r.Length, m, ks);
  857. }
  858. return new LongArray(r, 0, len);
  859. }
  860. internal LongArray Square(int m, int[] ks)
  861. {
  862. int len = GetUsedLength();
  863. if (len == 0)
  864. return this;
  865. ulong[] r = new ulong[len << 1];
  866. Interleave.Expand64To128(m_data, 0, len, r, 0);
  867. return new LongArray(r, 0, r.Length);
  868. }
  869. // private static LongArray ExpItohTsujii2(LongArray B, int n, int m, int[] ks)
  870. // {
  871. // LongArray t1 = B, t3 = new LongArray(new long[]{ 1L });
  872. // int scale = 1;
  873. //
  874. // int numTerms = n;
  875. // while (numTerms > 1)
  876. // {
  877. // if ((numTerms & 1) != 0)
  878. // {
  879. // t3 = t3.ModMultiply(t1, m, ks);
  880. // t1 = t1.modSquareN(scale, m, ks);
  881. // }
  882. //
  883. // LongArray t2 = t1.modSquareN(scale, m, ks);
  884. // t1 = t1.ModMultiply(t2, m, ks);
  885. // numTerms >>>= 1; scale <<= 1;
  886. // }
  887. //
  888. // return t3.ModMultiply(t1, m, ks);
  889. // }
  890. //
  891. // private static LongArray ExpItohTsujii23(LongArray B, int n, int m, int[] ks)
  892. // {
  893. // LongArray t1 = B, t3 = new LongArray(new long[]{ 1L });
  894. // int scale = 1;
  895. //
  896. // int numTerms = n;
  897. // while (numTerms > 1)
  898. // {
  899. // bool m03 = numTerms % 3 == 0;
  900. // bool m14 = !m03 && (numTerms & 1) != 0;
  901. //
  902. // if (m14)
  903. // {
  904. // t3 = t3.ModMultiply(t1, m, ks);
  905. // t1 = t1.modSquareN(scale, m, ks);
  906. // }
  907. //
  908. // LongArray t2 = t1.modSquareN(scale, m, ks);
  909. // t1 = t1.ModMultiply(t2, m, ks);
  910. //
  911. // if (m03)
  912. // {
  913. // t2 = t2.modSquareN(scale, m, ks);
  914. // t1 = t1.ModMultiply(t2, m, ks);
  915. // numTerms /= 3; scale *= 3;
  916. // }
  917. // else
  918. // {
  919. // numTerms >>>= 1; scale <<= 1;
  920. // }
  921. // }
  922. //
  923. // return t3.ModMultiply(t1, m, ks);
  924. // }
  925. //
  926. // private static LongArray ExpItohTsujii235(LongArray B, int n, int m, int[] ks)
  927. // {
  928. // LongArray t1 = B, t4 = new LongArray(new long[]{ 1L });
  929. // int scale = 1;
  930. //
  931. // int numTerms = n;
  932. // while (numTerms > 1)
  933. // {
  934. // if (numTerms % 5 == 0)
  935. // {
  936. //// t1 = ExpItohTsujii23(t1, 5, m, ks);
  937. //
  938. // LongArray t3 = t1;
  939. // t1 = t1.modSquareN(scale, m, ks);
  940. //
  941. // LongArray t2 = t1.modSquareN(scale, m, ks);
  942. // t1 = t1.ModMultiply(t2, m, ks);
  943. // t2 = t1.modSquareN(scale << 1, m, ks);
  944. // t1 = t1.ModMultiply(t2, m, ks);
  945. //
  946. // t1 = t1.ModMultiply(t3, m, ks);
  947. //
  948. // numTerms /= 5; scale *= 5;
  949. // continue;
  950. // }
  951. //
  952. // bool m03 = numTerms % 3 == 0;
  953. // bool m14 = !m03 && (numTerms & 1) != 0;
  954. //
  955. // if (m14)
  956. // {
  957. // t4 = t4.ModMultiply(t1, m, ks);
  958. // t1 = t1.modSquareN(scale, m, ks);
  959. // }
  960. //
  961. // LongArray t2 = t1.modSquareN(scale, m, ks);
  962. // t1 = t1.ModMultiply(t2, m, ks);
  963. //
  964. // if (m03)
  965. // {
  966. // t2 = t2.modSquareN(scale, m, ks);
  967. // t1 = t1.ModMultiply(t2, m, ks);
  968. // numTerms /= 3; scale *= 3;
  969. // }
  970. // else
  971. // {
  972. // numTerms >>>= 1; scale <<= 1;
  973. // }
  974. // }
  975. //
  976. // return t4.ModMultiply(t1, m, ks);
  977. // }
  978. internal LongArray ModInverse(int m, int[] ks)
  979. {
  980. /*
  981. * Fermat's Little Theorem
  982. */
  983. // LongArray A = this;
  984. // LongArray B = A.modSquare(m, ks);
  985. // LongArray R0 = B, R1 = B;
  986. // for (int i = 2; i < m; ++i)
  987. // {
  988. // R1 = R1.modSquare(m, ks);
  989. // R0 = R0.ModMultiply(R1, m, ks);
  990. // }
  991. //
  992. // return R0;
  993. /*
  994. * Itoh-Tsujii
  995. */
  996. // LongArray B = modSquare(m, ks);
  997. // switch (m)
  998. // {
  999. // case 409:
  1000. // return ExpItohTsujii23(B, m - 1, m, ks);
  1001. // case 571:
  1002. // return ExpItohTsujii235(B, m - 1, m, ks);
  1003. // case 163:
  1004. // case 233:
  1005. // case 283:
  1006. // default:
  1007. // return ExpItohTsujii2(B, m - 1, m, ks);
  1008. // }
  1009. /*
  1010. * Inversion in F2m using the extended Euclidean algorithm
  1011. *
  1012. * Input: A nonzero polynomial a(z) of degree at most m-1
  1013. * Output: a(z)^(-1) mod f(z)
  1014. */
  1015. int uzDegree = Degree();
  1016. if (uzDegree == 0)
  1017. throw new InvalidOperationException();
  1018. if (uzDegree == 1)
  1019. return this;
  1020. // u(z) := a(z)
  1021. LongArray uz = Copy();
  1022. int t = (m + 63) >> 6;
  1023. // v(z) := f(z)
  1024. LongArray vz = new LongArray(t);
  1025. ReduceBit(vz.m_data, 0, m, m, ks);
  1026. // g1(z) := 1, g2(z) := 0
  1027. LongArray g1z = new LongArray(t);
  1028. g1z.m_data[0] = 1UL;
  1029. LongArray g2z = new LongArray(t);
  1030. int[] uvDeg = new int[]{ uzDegree, m + 1 };
  1031. LongArray[] uv = new LongArray[]{ uz, vz };
  1032. int[] ggDeg = new int[]{ 1, 0 };
  1033. LongArray[] gg = new LongArray[]{ g1z, g2z };
  1034. int b = 1;
  1035. int duv1 = uvDeg[b];
  1036. int dgg1 = ggDeg[b];
  1037. int j = duv1 - uvDeg[1 - b];
  1038. for (;;)
  1039. {
  1040. if (j < 0)
  1041. {
  1042. j = -j;
  1043. uvDeg[b] = duv1;
  1044. ggDeg[b] = dgg1;
  1045. b = 1 - b;
  1046. duv1 = uvDeg[b];
  1047. dgg1 = ggDeg[b];
  1048. }
  1049. uv[b].AddShiftedByBitsSafe(uv[1 - b], uvDeg[1 - b], j);
  1050. int duv2 = uv[b].DegreeFrom(duv1);
  1051. if (duv2 == 0)
  1052. return gg[1 - b];
  1053. {
  1054. int dgg2 = ggDeg[1 - b];
  1055. gg[b].AddShiftedByBitsSafe(gg[1 - b], dgg2, j);
  1056. dgg2 += j;
  1057. if (dgg2 > dgg1)
  1058. {
  1059. dgg1 = dgg2;
  1060. }
  1061. else if (dgg2 == dgg1)
  1062. {
  1063. dgg1 = gg[b].DegreeFrom(dgg1);
  1064. }
  1065. }
  1066. j += (duv2 - duv1);
  1067. duv1 = duv2;
  1068. }
  1069. }
  1070. public override bool Equals(object obj)
  1071. {
  1072. if (obj is LongArray longArray)
  1073. return Equals(ref longArray);
  1074. return false;
  1075. }
  1076. internal bool Equals(ref LongArray other)
  1077. {
  1078. if (AreAliased(ref this, ref other))
  1079. return true;
  1080. int usedLen = GetUsedLength();
  1081. if (other.GetUsedLength() != usedLen)
  1082. return false;
  1083. for (int i = 0; i < usedLen; i++)
  1084. {
  1085. if (m_data[i] != other.m_data[i])
  1086. return false;
  1087. }
  1088. return true;
  1089. }
  1090. public override int GetHashCode()
  1091. {
  1092. int usedLen = GetUsedLength();
  1093. int hash = 1;
  1094. for (int i = 0; i < usedLen; i++)
  1095. {
  1096. ulong mi = m_data[i];
  1097. hash *= 31;
  1098. hash ^= (int)mi;
  1099. hash *= 31;
  1100. hash ^= (int)(mi >> 32);
  1101. }
  1102. return hash;
  1103. }
  1104. public LongArray Copy()
  1105. {
  1106. return new LongArray(Arrays.Clone(m_data));
  1107. }
  1108. public override string ToString()
  1109. {
  1110. int i = GetUsedLength();
  1111. if (i == 0)
  1112. return "0";
  1113. StringBuilder sb = new StringBuilder(i * 64);
  1114. sb.Append(Convert.ToString((long)m_data[--i], 2));
  1115. while (--i >= 0)
  1116. {
  1117. string s = Convert.ToString((long)m_data[i], 2);
  1118. // Add leading zeroes, except for highest significant word
  1119. int len = s.Length;
  1120. if (len < 64)
  1121. {
  1122. sb.Append('0', 64 - len);
  1123. }
  1124. sb.Append(s);
  1125. }
  1126. return sb.ToString();
  1127. }
  1128. }
  1129. }
  1130. #pragma warning restore
  1131. #endif