mpaux.c 2.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190
  1. #include "os.h"
  2. #include <mp.h>
  3. #include "dat.h"
  4. static mpdigit _mptwodata[1] = { 2 };
  5. static mpint _mptwo =
  6. {
  7. 1,
  8. 1,
  9. 1,
  10. _mptwodata,
  11. MPstatic
  12. };
  13. mpint *mptwo = &_mptwo;
  14. static mpdigit _mponedata[1] = { 1 };
  15. static mpint _mpone =
  16. {
  17. 1,
  18. 1,
  19. 1,
  20. _mponedata,
  21. MPstatic
  22. };
  23. mpint *mpone = &_mpone;
  24. static mpdigit _mpzerodata[1] = { 0 };
  25. static mpint _mpzero =
  26. {
  27. 1,
  28. 1,
  29. 0,
  30. _mpzerodata,
  31. MPstatic
  32. };
  33. mpint *mpzero = &_mpzero;
  34. static int mpmindigits = 33;
  35. // set minimum digit allocation
  36. void
  37. mpsetminbits(int n)
  38. {
  39. if(n < 0)
  40. sysfatal("mpsetminbits: n < 0");
  41. if(n == 0)
  42. n = 1;
  43. mpmindigits = DIGITS(n);
  44. }
  45. // allocate an n bit 0'd number
  46. mpint*
  47. mpnew(int n)
  48. {
  49. mpint *b;
  50. if(n < 0)
  51. sysfatal("mpsetminbits: n < 0");
  52. b = mallocz(sizeof(mpint), 1);
  53. setmalloctag(b, getcallerpc(&n));
  54. if(b == nil)
  55. sysfatal("mpnew: %r");
  56. n = DIGITS(n);
  57. if(n < mpmindigits)
  58. n = mpmindigits;
  59. b->p = (mpdigit*)mallocz(n*Dbytes, 1);
  60. if(b->p == nil)
  61. sysfatal("mpnew: %r");
  62. b->size = n;
  63. b->sign = 1;
  64. return b;
  65. }
  66. // guarantee at least n significant bits
  67. void
  68. mpbits(mpint *b, int m)
  69. {
  70. int n;
  71. n = DIGITS(m);
  72. if(b->size >= n){
  73. if(b->top >= n)
  74. return;
  75. memset(&b->p[b->top], 0, Dbytes*(n - b->top));
  76. b->top = n;
  77. return;
  78. }
  79. b->p = (mpdigit*)realloc(b->p, n*Dbytes);
  80. if(b->p == nil)
  81. sysfatal("mpbits: %r");
  82. memset(&b->p[b->top], 0, Dbytes*(n - b->top));
  83. b->size = n;
  84. b->top = n;
  85. }
  86. void
  87. mpfree(mpint *b)
  88. {
  89. if(b == nil)
  90. return;
  91. if(b->flags & MPstatic)
  92. sysfatal("freeing mp constant");
  93. memset(b->p, 0, b->size*Dbytes); // information hiding
  94. free(b->p);
  95. free(b);
  96. }
  97. void
  98. mpnorm(mpint *b)
  99. {
  100. int i;
  101. for(i = b->top-1; i >= 0; i--)
  102. if(b->p[i] != 0)
  103. break;
  104. b->top = i+1;
  105. if(b->top == 0)
  106. b->sign = 1;
  107. }
  108. mpint*
  109. mpcopy(mpint *old)
  110. {
  111. mpint *new;
  112. new = mpnew(Dbits*old->size);
  113. new->top = old->top;
  114. new->sign = old->sign;
  115. memmove(new->p, old->p, Dbytes*old->top);
  116. return new;
  117. }
  118. void
  119. mpassign(mpint *old, mpint *new)
  120. {
  121. mpbits(new, Dbits*old->top);
  122. new->sign = old->sign;
  123. new->top = old->top;
  124. memmove(new->p, old->p, Dbytes*old->top);
  125. }
  126. // number of significant bits in mantissa
  127. int
  128. mpsignif(mpint *n)
  129. {
  130. int i, j;
  131. mpdigit d;
  132. if(n->top == 0)
  133. return 0;
  134. for(i = n->top-1; i >= 0; i--){
  135. d = n->p[i];
  136. for(j = Dbits-1; j >= 0; j--){
  137. if(d & (((mpdigit)1)<<j))
  138. return i*Dbits + j + 1;
  139. }
  140. }
  141. return 0;
  142. }
  143. // k, where n = 2**k * q for odd q
  144. int
  145. mplowbits0(mpint *n)
  146. {
  147. int k, bit, digit;
  148. mpdigit d;
  149. if(n->top==0)
  150. return 0;
  151. k = 0;
  152. bit = 0;
  153. digit = 0;
  154. d = n->p[0];
  155. for(;;){
  156. if(d & (1<<bit))
  157. break;
  158. k++;
  159. bit++;
  160. if(bit==Dbits){
  161. if(++digit >= n->top)
  162. return 0;
  163. d = n->p[digit];
  164. bit = 0;
  165. }
  166. }
  167. return k;
  168. }