mul.c 4.8 KB

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