dmisc.c 4.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224
  1. /****************************************************************
  2. The author of this software is David M. Gay.
  3. Copyright (C) 1998 by Lucent Technologies
  4. All Rights Reserved
  5. Permission to use, copy, modify, and distribute this software and
  6. its documentation for any purpose and without fee is hereby
  7. granted, provided that the above copyright notice appear in all
  8. copies and that both that the copyright notice and this
  9. permission notice and warranty disclaimer appear in supporting
  10. documentation, and that the name of Lucent or any of its entities
  11. not be used in advertising or publicity pertaining to
  12. distribution of the software without specific, written prior
  13. permission.
  14. LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
  15. INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
  16. IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
  17. SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  18. WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
  19. IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
  20. ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
  21. THIS SOFTWARE.
  22. ****************************************************************/
  23. /* Please send bug reports to David M. Gay (dmg at acm dot org,
  24. * with " at " changed at "@" and " dot " changed to "."). */
  25. #include "gdtoaimp.h"
  26. #ifndef MULTIPLE_THREADS
  27. char *dtoa_result;
  28. #endif
  29. char *
  30. #ifdef KR_headers
  31. rv_alloc(i) int i;
  32. #else
  33. rv_alloc(int i)
  34. #endif
  35. {
  36. int j, k, *r;
  37. j = sizeof(ULong);
  38. for(k = 0;
  39. sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i;
  40. j <<= 1)
  41. k++;
  42. r = (int*)Balloc(k);
  43. if (r == NULL)
  44. return (
  45. #ifndef MULTIPLE_THREADS
  46. dtoa_result =
  47. #endif
  48. NULL);
  49. *r = k;
  50. return
  51. #ifndef MULTIPLE_THREADS
  52. dtoa_result =
  53. #endif
  54. (char *)(r+1);
  55. }
  56. char *
  57. #ifdef KR_headers
  58. nrv_alloc(s, rve, n) char *s, **rve; int n;
  59. #else
  60. nrv_alloc(char *s, char **rve, int n)
  61. #endif
  62. {
  63. char *rv, *t;
  64. t = rv = rv_alloc(n);
  65. if (t == NULL)
  66. return (NULL);
  67. while((*t = *s++) !=0)
  68. t++;
  69. if (rve)
  70. *rve = t;
  71. return rv;
  72. }
  73. /* freedtoa(s) must be used to free values s returned by dtoa
  74. * when MULTIPLE_THREADS is #defined. It should be used in all cases,
  75. * but for consistency with earlier versions of dtoa, it is optional
  76. * when MULTIPLE_THREADS is not defined.
  77. */
  78. void
  79. #ifdef KR_headers
  80. freedtoa(s) char *s;
  81. #else
  82. freedtoa(char *s)
  83. #endif
  84. {
  85. Bigint *b = (Bigint *)((int *)s - 1);
  86. b->maxwds = 1 << (b->k = *(int*)b);
  87. Bfree(b);
  88. #ifndef MULTIPLE_THREADS
  89. if (s == dtoa_result)
  90. dtoa_result = 0;
  91. #endif
  92. }
  93. int
  94. quorem
  95. #ifdef KR_headers
  96. (b, S) Bigint *b, *S;
  97. #else
  98. (Bigint *b, Bigint *S)
  99. #endif
  100. {
  101. int n;
  102. ULong *bx, *bxe, q, *sx, *sxe;
  103. #ifdef ULLong
  104. ULLong borrow, carry, y, ys;
  105. #else
  106. ULong borrow, carry, y, ys;
  107. #ifdef Pack_32
  108. ULong si, z, zs;
  109. #endif
  110. #endif
  111. n = S->wds;
  112. #ifdef DEBUG
  113. /*debug*/ if (b->wds > n)
  114. /*debug*/ Bug("oversize b in quorem");
  115. #endif
  116. if (b->wds < n)
  117. return 0;
  118. sx = S->x;
  119. sxe = sx + --n;
  120. bx = b->x;
  121. bxe = bx + n;
  122. q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
  123. #ifdef DEBUG
  124. /*debug*/ if (q > 9)
  125. /*debug*/ Bug("oversized quotient in quorem");
  126. #endif
  127. if (q) {
  128. borrow = 0;
  129. carry = 0;
  130. do {
  131. #ifdef ULLong
  132. ys = *sx++ * (ULLong)q + carry;
  133. carry = ys >> 32;
  134. y = *bx - (ys & 0xffffffffUL) - borrow;
  135. borrow = y >> 32 & 1UL;
  136. *bx++ = y & 0xffffffffUL;
  137. #else
  138. #ifdef Pack_32
  139. si = *sx++;
  140. ys = (si & 0xffff) * q + carry;
  141. zs = (si >> 16) * q + (ys >> 16);
  142. carry = zs >> 16;
  143. y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
  144. borrow = (y & 0x10000) >> 16;
  145. z = (*bx >> 16) - (zs & 0xffff) - borrow;
  146. borrow = (z & 0x10000) >> 16;
  147. Storeinc(bx, z, y);
  148. #else
  149. ys = *sx++ * q + carry;
  150. carry = ys >> 16;
  151. y = *bx - (ys & 0xffff) - borrow;
  152. borrow = (y & 0x10000) >> 16;
  153. *bx++ = y & 0xffff;
  154. #endif
  155. #endif
  156. }
  157. while(sx <= sxe);
  158. if (!*bxe) {
  159. bx = b->x;
  160. while(--bxe > bx && !*bxe)
  161. --n;
  162. b->wds = n;
  163. }
  164. }
  165. if (cmp(b, S) >= 0) {
  166. q++;
  167. borrow = 0;
  168. carry = 0;
  169. bx = b->x;
  170. sx = S->x;
  171. do {
  172. #ifdef ULLong
  173. ys = *sx++ + carry;
  174. carry = ys >> 32;
  175. y = *bx - (ys & 0xffffffffUL) - borrow;
  176. borrow = y >> 32 & 1UL;
  177. *bx++ = y & 0xffffffffUL;
  178. #else
  179. #ifdef Pack_32
  180. si = *sx++;
  181. ys = (si & 0xffff) + carry;
  182. zs = (si >> 16) + (ys >> 16);
  183. carry = zs >> 16;
  184. y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
  185. borrow = (y & 0x10000) >> 16;
  186. z = (*bx >> 16) - (zs & 0xffff) - borrow;
  187. borrow = (z & 0x10000) >> 16;
  188. Storeinc(bx, z, y);
  189. #else
  190. ys = *sx++ + carry;
  191. carry = ys >> 16;
  192. y = *bx - (ys & 0xffff) - borrow;
  193. borrow = (y & 0x10000) >> 16;
  194. *bx++ = y & 0xffff;
  195. #endif
  196. #endif
  197. }
  198. while(sx <= sxe);
  199. bx = b->x;
  200. bxe = bx + n;
  201. if (!*bxe) {
  202. while(--bxe > bx && !*bxe)
  203. --n;
  204. b->wds = n;
  205. }
  206. }
  207. return q;
  208. }