fpe.c 5.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240
  1. /*
  2. * File: dclxvi-20130329/fpe.c
  3. * Author: Ruben Niederhagen, Peter Schwabe
  4. * Public Domain
  5. */
  6. #include <math.h>
  7. #include <assert.h>
  8. #include "scalar.h"
  9. #include "fpe.h"
  10. #include "mul.h"
  11. extern const scalar_t bn_pminus2;
  12. extern const double bn_v;
  13. extern const double bn_v6;
  14. void fpe_short_coeffred(fpe_t rop)
  15. {
  16. mydouble carry11 = round(rop->v[11]/bn_v);
  17. rop->v[11] = remround(rop->v[11],bn_v);
  18. rop->v[0] = rop->v[0] - carry11;
  19. rop->v[3] = rop->v[3] - carry11;
  20. rop->v[6] = rop->v[6] - 4*carry11;
  21. rop->v[9] = rop->v[9] - carry11;
  22. mydouble carry0 = round(rop->v[0]/bn_v6);
  23. mydouble carry1 = round(rop->v[1]/bn_v);
  24. mydouble carry2 = round(rop->v[2]/bn_v);
  25. mydouble carry3 = round(rop->v[3]/bn_v);
  26. mydouble carry4 = round(rop->v[4]/bn_v);
  27. mydouble carry5 = round(rop->v[5]/bn_v);
  28. mydouble carry6 = round(rop->v[6]/bn_v6);
  29. mydouble carry7 = round(rop->v[7]/bn_v);
  30. mydouble carry8 = round(rop->v[8]/bn_v);
  31. mydouble carry9 = round(rop->v[9]/bn_v);
  32. mydouble carry10 = round(rop->v[10]/bn_v);
  33. rop->v[0] = remround(rop->v[0],bn_v6);
  34. rop->v[1] = remround(rop->v[1],bn_v);
  35. rop->v[2] = remround(rop->v[2],bn_v);
  36. rop->v[3] = remround(rop->v[3],bn_v);
  37. rop->v[4] = remround(rop->v[4],bn_v);
  38. rop->v[5] = remround(rop->v[5],bn_v);
  39. rop->v[6] = remround(rop->v[6],bn_v6);
  40. rop->v[7] = remround(rop->v[7],bn_v);
  41. rop->v[8] = remround(rop->v[8],bn_v);
  42. rop->v[9] = remround(rop->v[9],bn_v);
  43. rop->v[10] = remround(rop->v[10],bn_v);
  44. rop->v[1] += carry0;
  45. rop->v[2] += carry1;
  46. rop->v[3] += carry2;
  47. rop->v[4] += carry3;
  48. rop->v[5] += carry4;
  49. rop->v[6] += carry5;
  50. rop->v[7] += carry6;
  51. rop->v[8] += carry7;
  52. rop->v[9] += carry8;
  53. rop->v[10] += carry9;
  54. rop->v[11] += carry10;
  55. }
  56. // Set fpe_t rop to given value:
  57. void fpe_set(fpe_t rop, const fpe_t op)
  58. {
  59. int i;
  60. for(i=0;i<12;i++)
  61. rop->v[i] = op->v[i];
  62. }
  63. /* Communicate the fact that the fpe is reduced (and that we don't know anything more about it) */
  64. void fpe_isreduced(fpe_t rop)
  65. {
  66. setmax(rop->v[0],(long)bn_v6/2);
  67. setmax(rop->v[6],(long)bn_v6/2);
  68. setmax(rop->v[1],(long)bn_v/2);
  69. setmax(rop->v[3],(long)bn_v/2);
  70. setmax(rop->v[4],(long)bn_v/2);
  71. setmax(rop->v[7],(long)bn_v/2);
  72. setmax(rop->v[9],(long)bn_v/2);
  73. setmax(rop->v[10],(long)bn_v/2);
  74. //XXX: Change additive constant:
  75. setmax(rop->v[2],(long)bn_v/2+2331);
  76. setmax(rop->v[5],(long)bn_v/2+2331);
  77. setmax(rop->v[8],(long)bn_v/2+2331);
  78. setmax(rop->v[11],(long)bn_v/2+2331);
  79. }
  80. // Set fpe_t rop to value given in double array of length 12
  81. void fpe_set_doublearray(fpe_t rop, const mydouble op[12])
  82. {
  83. int i;
  84. for(i=0;i<12;i++)
  85. rop->v[i] = op[i];
  86. }
  87. // Set rop to one
  88. void fpe_setone(fpe_t rop)
  89. {
  90. int i;
  91. for(i=1;i<12;i++)
  92. rop->v[i] = 0.;
  93. rop->v[0] = 1;
  94. }
  95. // Set rop to zero
  96. void fpe_setzero(fpe_t rop)
  97. {
  98. int i;
  99. for(i=0;i<12;i++)
  100. rop->v[i] = 0.;
  101. }
  102. int fpe_iseq(const fpe_t op1, const fpe_t op2)
  103. {
  104. fpe_t t;
  105. fpe_sub(t, op1, op2);
  106. return fpe_iszero(t);
  107. }
  108. int fpe_isone(const fpe_t op)
  109. {
  110. fpe_t t;
  111. int i;
  112. for(i=1;i<12;i++)
  113. t->v[i] = op->v[i];
  114. t->v[0] = op->v[0] - 1.;
  115. return fpe_iszero(t);
  116. }
  117. int fpe_iszero(const fpe_t op)
  118. {
  119. fpe_t t;
  120. double d;
  121. int i;
  122. unsigned long long tr = 0;
  123. unsigned int differentbits=0;
  124. for(i=0;i<12;i++)
  125. t->v[i] = op->v[i];
  126. coeffred_round_par(t->v);
  127. //Constant-time comparison
  128. double zero = 0.;
  129. unsigned long long *zp = (unsigned long long *)&zero;;
  130. unsigned long long *tp;
  131. for(i=0;i<12;i++)
  132. {
  133. d = todouble(t->v[i]);
  134. tp = (unsigned long long *)&d;
  135. tr |= (*tp ^ *zp);
  136. }
  137. for(i=0;i<8;i++)
  138. differentbits |= i[(unsigned char*)&tr];
  139. return 1 & ((differentbits - 1) >> 8);
  140. }
  141. // Compute the negative of an fpe
  142. void fpe_neg(fpe_t rop, const fpe_t op)
  143. {
  144. int i;
  145. for(i=0;i<12;i++)
  146. rop->v[i] = -op->v[i];
  147. }
  148. // Double an fpe:
  149. void fpe_double(fpe_t rop, const fpe_t op)
  150. {
  151. int i;
  152. for(i=0;i<12;i++)
  153. rop->v[i] = op->v[i]*2;
  154. }
  155. // Double an fpe:
  156. void fpe_triple(fpe_t rop, const fpe_t op)
  157. {
  158. int i;
  159. for(i=0;i<12;i++)
  160. rop->v[i] = op->v[i]*3;
  161. }
  162. // Add two fpe, store result in rop:
  163. void fpe_add(fpe_t rop, const fpe_t op1, const fpe_t op2)
  164. {
  165. int i;
  166. for(i=0;i<12;i++)
  167. rop->v[i] = op1->v[i] + op2->v[i];
  168. }
  169. // Subtract op2 from op1, store result in rop:
  170. void fpe_sub(fpe_t rop, const fpe_t op1, const fpe_t op2)
  171. {
  172. int i;
  173. for(i=0;i<12;i++)
  174. rop->v[i] = op1->v[i] - op2->v[i];
  175. }
  176. // Multiply two fpe, store result in rop:
  177. #ifndef QHASM
  178. void fpe_mul_c(fpe_t rop, const fpe_t op1, const fpe_t op2)
  179. {
  180. mydouble h[24];
  181. polymul(h,op1->v,op2->v);
  182. degred(h);
  183. coeffred_round_par(h);
  184. int i;
  185. for (i=0;i<12;i++)
  186. rop->v[i] = h[i];
  187. }
  188. #endif
  189. // Square an fpe, store result in rop:
  190. void fpe_square(fpe_t rop, const fpe_t op)
  191. {
  192. /* Not used during pairing computation */
  193. fpe_mul(rop, op, op);
  194. }
  195. // Compute inverse of an fpe, store result in rop:
  196. void fpe_invert(fpe_t rop, const fpe_t op1)
  197. {
  198. fpe_set(rop,op1);
  199. int i;
  200. for(i = 254; i >= 0; i--)
  201. {
  202. fpe_mul(rop,rop,rop);
  203. if(scalar_getbit(bn_pminus2, i))
  204. fpe_mul(rop,rop,op1);
  205. }
  206. }
  207. // Print the element to stdout:
  208. void fpe_print(FILE * outfile, const fpe_t op)
  209. {
  210. int i;
  211. for(i=0;i<11;i++) fprintf(outfile, "%10lf, ", todouble(op->v[i]));
  212. fprintf(outfile, "%10lf", todouble(op->v[11]));
  213. }