mul.c 4.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124
  1. /*
  2. * File: dclxvi-20130329/mul.c
  3. * Author: Ruben Niederhagen, Peter Schwabe
  4. * Public Domain
  5. */
  6. #include <math.h>
  7. #include "mul.h"
  8. #include "mydouble.h"
  9. extern const double bn_v;
  10. extern const double bn_v6;
  11. void polymul(mydouble *h, const mydouble *f, const mydouble *g)
  12. {
  13. mydouble t[24];
  14. t[0] = f[0]*g[0];
  15. t[1] = f[0]*g[1] + f[1]*g[0];
  16. t[2] = 6*f[1]*g[1] + (f[0]*g[2] + f[2]*g[0]);
  17. t[3] = (f[1]*g[2] + f[2]*g[1])*6 + (f[0]*g[3] + f[3]*g[0]);
  18. t[4] = (f[1]*g[3] + f[2]*g[2] + f[3]*g[1])*6 + (f[0]*g[4] + f[4]*g[0]);
  19. t[5] = (f[1]*g[4] + f[2]*g[3] + f[3]*g[2] + f[4]*g[1])*6 + (f[0]*g[5] + f[5]*g[0]);
  20. t[6] = (f[1]*g[5] + f[2]*g[4] + f[3]*g[3] + f[4]*g[2] + f[5]*g[1])*6 + f[0]*g[6] + f[6]*g[0];
  21. t[7] = (f[0]*g[7] + f[1]*g[6] + f[2]*g[5] + f[3]*g[4] + f[4]*g[3] + f[5]*g[2] + f[6]*g[1] + f[7]*g[0]);
  22. t[8] = (f[1]*g[7] + f[7]*g[1])*6 + (f[0]*g[8] + f[2]*g[6] + f[3]*g[5] + f[4]*g[4] + f[5]*g[3] + f[6]*g[2] + f[8]*g[0]);
  23. t[9] = (f[1]*g[8] + f[2]*g[7] + f[7]*g[2] + f[8]*g[1])*6 + (f[0]*g[9] + f[3]*g[6] + f[4]*g[5] + f[5]*g[4] + f[6]*g[3] + f[9]*g[0]);
  24. t[10] = (f[1]*g[9] + f[2]*g[8] + f[3]*g[7] + f[7]*g[3] + f[8]*g[2] + f[9]*g[1])*6 + (f[0]*g[10] + f[4]*g[6] + f[5]*g[5] + f[6]*g[4] + f[10]*g[0]);
  25. t[11] = (f[1]*g[10] + f[2]*g[9] + f[3]*g[8] + f[4]*g[7] + f[7]*g[4] + f[8]*g[3] + f[9]*g[2] + f[10]*g[1])*6 + (f[0]*g[11] + f[5]*g[6] + f[6]*g[5] + f[11]*g[0]);
  26. t[12] = (f[1]*g[11] + f[2]*g[10] + f[3]*g[9] + f[4]*g[8] + f[5]*g[7] + f[7]*g[5] + f[8]*g[4] + f[9]*g[3] + f[10]*g[2] + f[11]*g[1])*6 + f[6]*g[6];
  27. t[13] = (f[2]*g[11] + f[3]*g[10] + f[4]*g[9] + f[5]*g[8] + f[6]*g[7] + f[7]*g[6] + f[8]*g[5] + f[9]*g[4] + f[10]*g[3] + f[11]*g[2]);
  28. t[14] = f[7]*g[7]*6 + (f[3]*g[11] + f[4]*g[10] + f[5]*g[9] + f[6]*g[8] + f[8]*g[6] + f[9]*g[5] + f[10]*g[4] + f[11]*g[3]);
  29. t[15] = (f[7]*g[8] + f[8]*g[7])*6 + (f[4]*g[11] + f[5]*g[10] + f[6]*g[9] + f[9]*g[6] + f[10]*g[5] + f[11]*g[4]);
  30. t[16] = (f[7]*g[9] + f[8]*g[8] + f[9]*g[7])*6 + (f[5]*g[11] + f[6]*g[10] + f[10]*g[6] + f[11]*g[5]);
  31. t[17] = (f[7]*g[10] + f[8]*g[9] + f[9]*g[8] + f[10]*g[7])*6 + (f[6]*g[11] + f[11]*g[6]);
  32. t[18] = (f[7]*g[11] + f[8]*g[10] + f[9]*g[9] + f[10]*g[8] + f[11]*g[7])*6;
  33. t[19] = (f[8]*g[11] + f[9]*g[10] + f[10]*g[9] + f[11]*g[8]);
  34. t[20] = (f[9]*g[11] + f[10]*g[10] + f[11]*g[9]);
  35. t[21] = (f[10]*g[11] + f[11]*g[10]);
  36. t[22] = f[11]*g[11];
  37. int i;
  38. for(i=0;i<23;i++)
  39. h[i]=t[i];
  40. }
  41. void degred(mydouble *h)
  42. {
  43. h[0] = h[0] - h[12] + 6*h[15] - 2*h[18] - 6*h[21];
  44. h[1] = h[1] - h[13] + h[16] - 2*h[19] - h[22];
  45. h[2] = h[2] - h[14] + h[17] - 2*h[20];
  46. h[3] = h[3] - h[12] + 5*h[15] - h[18] - 8*h[21];
  47. h[4] = h[4] - 6*h[13] + 5*h[16] - 6*h[19] - 8*h[22];
  48. h[5] = h[5] - 6*h[14] + 5*h[17] - 6*h[20];
  49. h[6] = h[6] - 4*h[12] + 18*h[15] - 3*h[18];
  50. h[6] -= 30*h[21];
  51. h[7] = h[7] - 4*h[13] + 3*h[16] - 3*h[19] - 5*h[22];
  52. h[8] = h[8] - 4*h[14] + 3*h[17] - 3*h[20];
  53. h[9] = h[9] - h[12] + 2*h[15] + h[18] - 9*h[21];
  54. h[10] = h[10] - 6*h[13] + 2*h[16] + 6*h[19] - 9*h[22];
  55. h[11] = h[11] - 6*h[14] + 2*h[17] + 6*h[20];
  56. }
  57. void coeffred_round_par(mydouble *h)
  58. {
  59. mydouble carry = 0;
  60. carry = round(h[1]/bn_v);
  61. h[1] = remround(h[1],bn_v);
  62. h[2] += carry;
  63. carry = round(h[4]/bn_v);
  64. h[4] = remround(h[4],bn_v);
  65. h[5] += carry;
  66. carry = round(h[7]/bn_v);
  67. h[7] = remround(h[7],bn_v);
  68. h[8] += carry;
  69. carry = round(h[10]/bn_v);
  70. h[10] = remround(h[10],bn_v);
  71. h[11] += carry;
  72. carry = round(h[2]/bn_v);
  73. h[2] = remround(h[2],bn_v);
  74. h[3] += carry;
  75. carry = round(h[5]/bn_v);
  76. h[5] = remround(h[5],bn_v);
  77. h[6] += carry;
  78. carry = round(h[8]/bn_v);
  79. h[8] = remround(h[8],bn_v);
  80. h[9] += carry;
  81. carry = round(h[11]/bn_v);
  82. h[11] = remround(h[11],bn_v);
  83. h[0] = h[0] - carry;
  84. h[3] = h[3] - carry;
  85. h[6] = h[6] - 4*carry;
  86. h[9] = h[9] - carry;
  87. carry = round(h[0]/bn_v6); // h0 = 2^53 - 1
  88. h[0] = remround(h[0],bn_v6); // carry = (2^53-1)/6v = 763549741
  89. h[1] += carry; // h1 = v+763549741 = 765515821
  90. carry = round(h[3]/bn_v); // h3 = 2^53 - 1
  91. h[3] = remround(h[3],bn_v); // carry = (2^53-1)/v = 4581298449
  92. h[4] += carry; // h4 = v + 4581298449 = 4583264529
  93. carry = round(h[6]/bn_v6);
  94. h[6] = remround(h[6],bn_v6);
  95. h[7] += carry;
  96. carry = round(h[9]/bn_v);
  97. h[9] = remround(h[9],bn_v);
  98. h[10] += carry;
  99. carry = round(h[1]/bn_v); // carry = 765515821/v = 389
  100. h[1] = remround(h[1],bn_v);
  101. h[2] += carry;
  102. carry = round(h[4]/bn_v); // carry = 4583264529/v = 2331
  103. h[4] = remround(h[4],bn_v);
  104. h[5] += carry;
  105. carry = round(h[7]/bn_v);
  106. h[7] = remround(h[7],bn_v);
  107. h[8] += carry;
  108. carry = round(h[10]/bn_v);
  109. h[10] = remround(h[10],bn_v);
  110. h[11] += carry;
  111. }