bid128_minmax.c 38 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095
  1. /* Copyright (C) 2007-2022 Free Software Foundation, Inc.
  2. This file is part of GCC.
  3. GCC is free software; you can redistribute it and/or modify it under
  4. the terms of the GNU General Public License as published by the Free
  5. Software Foundation; either version 3, or (at your option) any later
  6. version.
  7. GCC is distributed in the hope that it will be useful, but WITHOUT ANY
  8. WARRANTY; without even the implied warranty of MERCHANTABILITY or
  9. FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
  10. for more details.
  11. Under Section 7 of GPL version 3, you are granted additional
  12. permissions described in the GCC Runtime Library Exception, version
  13. 3.1, as published by the Free Software Foundation.
  14. You should have received a copy of the GNU General Public License and
  15. a copy of the GCC Runtime Library Exception along with this program;
  16. see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
  17. <http://www.gnu.org/licenses/>. */
  18. #define BID_128RES
  19. #include "bid_internal.h"
  20. /*****************************************************************************
  21. * BID128 minimum number
  22. *****************************************************************************/
  23. #if DECIMAL_CALL_BY_REFERENCE
  24. void
  25. bid128_minnum (UINT128 * pres, UINT128 * px,
  26. UINT128 * py _EXC_FLAGS_PARAM) {
  27. UINT128 x = *px;
  28. UINT128 y = *py;
  29. #else
  30. UINT128
  31. bid128_minnum (UINT128 x, UINT128 y _EXC_FLAGS_PARAM) {
  32. #endif
  33. UINT128 res;
  34. int exp_x, exp_y;
  35. int diff;
  36. UINT128 sig_x, sig_y;
  37. UINT192 sig_n_prime192;
  38. UINT256 sig_n_prime256;
  39. char x_is_zero = 0, y_is_zero = 0;
  40. BID_SWAP128 (x);
  41. BID_SWAP128 (y);
  42. // check for non-canonical x
  43. if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
  44. x.w[1] = x.w[1] & 0xfe003fffffffffffull; // clear out G[6]-G[16]
  45. // check for non-canonical NaN payload
  46. if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
  47. (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
  48. (x.w[0] > 0x38c15b09ffffffffull))) {
  49. x.w[1] = x.w[1] & 0xffffc00000000000ull;
  50. x.w[0] = 0x0ull;
  51. }
  52. } else if ((x.w[1] & MASK_ANY_INF) == MASK_INF) { // x = inf
  53. x.w[1] = x.w[1] & (MASK_SIGN | MASK_INF);
  54. x.w[0] = 0x0ull;
  55. } else { // x is not special
  56. // check for non-canonical values - treated as zero
  57. if ((x.w[1] & MASK_STEERING_BITS) == MASK_STEERING_BITS) { // G0_G1=11
  58. // non-canonical
  59. x.w[1] = (x.w[1] & MASK_SIGN) | ((x.w[1] << 2) & MASK_EXP);
  60. x.w[0] = 0x0ull;
  61. } else { // G0_G1 != 11
  62. if ((x.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
  63. ((x.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull
  64. && x.w[0] > 0x378d8e63ffffffffull)) {
  65. // x is non-canonical if coefficient is larger than 10^34 -1
  66. x.w[1] = (x.w[1] & MASK_SIGN) | (x.w[1] & MASK_EXP);
  67. x.w[0] = 0x0ull;
  68. } else { // canonical
  69. ;
  70. }
  71. }
  72. }
  73. // check for non-canonical y
  74. if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NAN
  75. y.w[1] = y.w[1] & 0xfe003fffffffffffull; // clear out G[6]-G[16]
  76. // check for non-canonical NaN payload
  77. if (((y.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
  78. (((y.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
  79. (y.w[0] > 0x38c15b09ffffffffull))) {
  80. y.w[1] = y.w[1] & 0xffffc00000000000ull;
  81. y.w[0] = 0x0ull;
  82. }
  83. } else if ((y.w[1] & MASK_ANY_INF) == MASK_INF) { // y = inf
  84. y.w[1] = y.w[1] & (MASK_SIGN | MASK_INF);
  85. y.w[0] = 0x0ull;
  86. } else { // y is not special
  87. // check for non-canonical values - treated as zero
  88. if ((y.w[1] & MASK_STEERING_BITS) == MASK_STEERING_BITS) { // G0_G1=11
  89. // non-canonical
  90. y.w[1] = (y.w[1] & MASK_SIGN) | ((y.w[1] << 2) & MASK_EXP);
  91. y.w[0] = 0x0ull;
  92. } else { // G0_G1 != 11
  93. if ((y.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
  94. ((y.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull
  95. && y.w[0] > 0x378d8e63ffffffffull)) {
  96. // y is non-canonical if coefficient is larger than 10^34 -1
  97. y.w[1] = (y.w[1] & MASK_SIGN) | (y.w[1] & MASK_EXP);
  98. y.w[0] = 0x0ull;
  99. } else { // canonical
  100. ;
  101. }
  102. }
  103. }
  104. // NaN (CASE1)
  105. if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
  106. if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNaN
  107. // if x is SNAN, then return quiet (x)
  108. *pfpsf |= INVALID_EXCEPTION; // set exception if SNaN
  109. x.w[1] = x.w[1] & 0xfdffffffffffffffull; // quietize x
  110. res = x;
  111. } else { // x is QNaN
  112. if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NAN
  113. if ((y.w[1] & MASK_SNAN) == MASK_SNAN) { // y is SNAN
  114. *pfpsf |= INVALID_EXCEPTION; // set invalid flag
  115. }
  116. res = x;
  117. } else {
  118. res = y;
  119. }
  120. }
  121. BID_RETURN (res);
  122. } else if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NaN, but x is not
  123. if ((y.w[1] & MASK_SNAN) == MASK_SNAN) {
  124. *pfpsf |= INVALID_EXCEPTION; // set exception if SNaN
  125. y.w[1] = y.w[1] & 0xfdffffffffffffffull; // quietize y
  126. res = y;
  127. } else {
  128. // will return x (which is not NaN)
  129. res = x;
  130. }
  131. BID_RETURN (res);
  132. }
  133. // SIMPLE (CASE2)
  134. // if all the bits are the same, these numbers are equal (not Greater).
  135. if (x.w[0] == y.w[0] && x.w[1] == y.w[1]) {
  136. res = x;
  137. BID_RETURN (res);
  138. }
  139. // INFINITY (CASE3)
  140. if ((x.w[1] & MASK_INF) == MASK_INF) {
  141. // if x is neg infinity, there is no way it is greater than y, return 0
  142. res = (((x.w[1] & MASK_SIGN) == MASK_SIGN)) ? x : y;
  143. BID_RETURN (res);
  144. } else if ((y.w[1] & MASK_INF) == MASK_INF) {
  145. // x is finite, so if y is positive infinity, then x is less, return 0
  146. // if y is negative infinity, then x is greater, return 1
  147. res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? y : x;
  148. BID_RETURN (res);
  149. }
  150. // CONVERT X
  151. sig_x.w[1] = x.w[1] & 0x0001ffffffffffffull;
  152. sig_x.w[0] = x.w[0];
  153. exp_x = (x.w[1] >> 49) & 0x000000000003fffull;
  154. // CONVERT Y
  155. exp_y = (y.w[1] >> 49) & 0x0000000000003fffull;
  156. sig_y.w[1] = y.w[1] & 0x0001ffffffffffffull;
  157. sig_y.w[0] = y.w[0];
  158. // ZERO (CASE4)
  159. // some properties:
  160. // (+ZERO == -ZERO) => therefore ignore the sign
  161. // (ZERO x 10^A == ZERO x 10^B) for any valid A, B => ignore the exponent
  162. // field
  163. // (Any non-canonical # is considered 0)
  164. if ((sig_x.w[1] == 0) && (sig_x.w[0] == 0)) {
  165. x_is_zero = 1;
  166. }
  167. if ((sig_y.w[1] == 0) && (sig_y.w[0] == 0)) {
  168. y_is_zero = 1;
  169. }
  170. if (x_is_zero && y_is_zero) {
  171. // if both numbers are zero, neither is greater => return either number
  172. res = x;
  173. BID_RETURN (res);
  174. } else if (x_is_zero) {
  175. // is x is zero, it is greater if Y is negative
  176. res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? y : x;
  177. BID_RETURN (res);
  178. } else if (y_is_zero) {
  179. // is y is zero, X is greater if it is positive
  180. res = ((x.w[1] & MASK_SIGN) != MASK_SIGN) ? y : x;
  181. BID_RETURN (res);
  182. }
  183. // OPPOSITE SIGN (CASE5)
  184. // now, if the sign bits differ, x is greater if y is negative
  185. if (((x.w[1] ^ y.w[1]) & MASK_SIGN) == MASK_SIGN) {
  186. res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? y : x;
  187. BID_RETURN (res);
  188. }
  189. // REDUNDANT REPRESENTATIONS (CASE6)
  190. // if exponents are the same, then we have a simple comparison of
  191. // the significands
  192. if (exp_y == exp_x) {
  193. res = (((sig_x.w[1] > sig_y.w[1])
  194. || (sig_x.w[1] == sig_y.w[1]
  195. && sig_x.w[0] >= sig_y.w[0])) ^ ((x.w[1] & MASK_SIGN) ==
  196. MASK_SIGN)) ? y : x;
  197. BID_RETURN (res);
  198. }
  199. // if both components are either bigger or smaller, it is clear what
  200. // needs to be done
  201. if (sig_x.w[1] >= sig_y.w[1] && sig_x.w[0] >= sig_y.w[0]
  202. && exp_x > exp_y) {
  203. res = ((x.w[1] & MASK_SIGN) != MASK_SIGN) ? y : x;
  204. BID_RETURN (res);
  205. }
  206. if (sig_x.w[1] <= sig_y.w[1] && sig_x.w[0] <= sig_y.w[0]
  207. && exp_x < exp_y) {
  208. res = ((x.w[1] & MASK_SIGN) == MASK_SIGN) ? y : x;
  209. BID_RETURN (res);
  210. }
  211. diff = exp_x - exp_y;
  212. // if |exp_x - exp_y| < 33, it comes down to the compensated significand
  213. if (diff > 0) { // to simplify the loop below,
  214. // if exp_x is 33 greater than exp_y, no need for compensation
  215. if (diff > 33) {
  216. // difference cannot be greater than 10^33
  217. res = ((x.w[1] & MASK_SIGN) != MASK_SIGN) ? y : x;
  218. BID_RETURN (res);
  219. }
  220. if (diff > 19) { //128 by 128 bit multiply -> 256 bits
  221. __mul_128x128_to_256 (sig_n_prime256, sig_x, ten2k128[diff - 20]);
  222. // if postitive, return whichever significand is larger
  223. // (converse if negative)
  224. res = ((((sig_n_prime256.w[3] > 0) || sig_n_prime256.w[2] > 0)
  225. || (sig_n_prime256.w[1] > sig_y.w[1])
  226. || (sig_n_prime256.w[1] == sig_y.w[1]
  227. && sig_n_prime256.w[0] >
  228. sig_y.w[0])) ^ ((y.w[1] & MASK_SIGN) ==
  229. MASK_SIGN)) ? y : x;
  230. BID_RETURN (res);
  231. }
  232. __mul_64x128_to_192 (sig_n_prime192, ten2k64[diff], sig_x);
  233. // if postitive, return whichever significand is larger
  234. // (converse if negative)
  235. res =
  236. (((sig_n_prime192.w[2] > 0) || (sig_n_prime192.w[1] > sig_y.w[1])
  237. || (sig_n_prime192.w[1] == sig_y.w[1]
  238. && sig_n_prime192.w[0] >
  239. sig_y.w[0])) ^ ((y.w[1] & MASK_SIGN) == MASK_SIGN)) ? y : x;
  240. BID_RETURN (res);
  241. }
  242. diff = exp_y - exp_x;
  243. // if exp_x is 33 less than exp_y, no need for compensation
  244. if (diff > 33) {
  245. res = ((x.w[1] & MASK_SIGN) == MASK_SIGN) ? y : x;
  246. BID_RETURN (res);
  247. }
  248. if (diff > 19) { //128 by 128 bit multiply -> 256 bits
  249. // adjust the y significand upwards
  250. __mul_128x128_to_256 (sig_n_prime256, sig_y, ten2k128[diff - 20]);
  251. // if postitive, return whichever significand is larger
  252. // (converse if negative)
  253. res =
  254. ((sig_n_prime256.w[3] != 0 || sig_n_prime256.w[2] != 0
  255. || (sig_n_prime256.w[1] > sig_x.w[1]
  256. || (sig_n_prime256.w[1] == sig_x.w[1]
  257. && sig_n_prime256.w[0] >
  258. sig_x.w[0]))) ^ ((x.w[1] & MASK_SIGN) ==
  259. MASK_SIGN)) ? x : y;
  260. BID_RETURN (res);
  261. }
  262. // adjust the y significand upwards
  263. __mul_64x128_to_192 (sig_n_prime192, ten2k64[diff], sig_y);
  264. // if postitive, return whichever significand is larger (converse if negative)
  265. res =
  266. ((sig_n_prime192.w[2] != 0
  267. || (sig_n_prime192.w[1] > sig_x.w[1]
  268. || (sig_n_prime192.w[1] == sig_x.w[1]
  269. && sig_n_prime192.w[0] > sig_x.w[0])))
  270. ^ ((y.w[1] & MASK_SIGN) == MASK_SIGN)) ? x : y;
  271. BID_RETURN (res);
  272. }
  273. /*****************************************************************************
  274. * BID128 minimum magnitude function - returns greater of two numbers
  275. *****************************************************************************/
  276. #if DECIMAL_CALL_BY_REFERENCE
  277. void
  278. bid128_minnum_mag (UINT128 * pres, UINT128 * px,
  279. UINT128 * py _EXC_FLAGS_PARAM) {
  280. UINT128 x = *px;
  281. UINT128 y = *py;
  282. #else
  283. UINT128
  284. bid128_minnum_mag (UINT128 x, UINT128 y _EXC_FLAGS_PARAM) {
  285. #endif
  286. UINT128 res;
  287. int exp_x, exp_y;
  288. int diff;
  289. UINT128 sig_x, sig_y;
  290. UINT192 sig_n_prime192;
  291. UINT256 sig_n_prime256;
  292. BID_SWAP128 (x);
  293. BID_SWAP128 (y);
  294. // check for non-canonical x
  295. if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
  296. x.w[1] = x.w[1] & 0xfe003fffffffffffull; // clear out G[6]-G[16]
  297. // check for non-canonical NaN payload
  298. if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
  299. (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
  300. (x.w[0] > 0x38c15b09ffffffffull))) {
  301. x.w[1] = x.w[1] & 0xffffc00000000000ull;
  302. x.w[0] = 0x0ull;
  303. }
  304. } else if ((x.w[1] & MASK_ANY_INF) == MASK_INF) { // x = inf
  305. x.w[1] = x.w[1] & (MASK_SIGN | MASK_INF);
  306. x.w[0] = 0x0ull;
  307. } else { // x is not special
  308. // check for non-canonical values - treated as zero
  309. if ((x.w[1] & MASK_STEERING_BITS) == MASK_STEERING_BITS) { // G0_G1=11
  310. // non-canonical
  311. x.w[1] = (x.w[1] & MASK_SIGN) | ((x.w[1] << 2) & MASK_EXP);
  312. x.w[0] = 0x0ull;
  313. } else { // G0_G1 != 11
  314. if ((x.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
  315. ((x.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull
  316. && x.w[0] > 0x378d8e63ffffffffull)) {
  317. // x is non-canonical if coefficient is larger than 10^34 -1
  318. x.w[1] = (x.w[1] & MASK_SIGN) | (x.w[1] & MASK_EXP);
  319. x.w[0] = 0x0ull;
  320. } else { // canonical
  321. ;
  322. }
  323. }
  324. }
  325. // check for non-canonical y
  326. if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NAN
  327. y.w[1] = y.w[1] & 0xfe003fffffffffffull; // clear out G[6]-G[16]
  328. // check for non-canonical NaN payload
  329. if (((y.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
  330. (((y.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
  331. (y.w[0] > 0x38c15b09ffffffffull))) {
  332. y.w[1] = y.w[1] & 0xffffc00000000000ull;
  333. y.w[0] = 0x0ull;
  334. }
  335. } else if ((y.w[1] & MASK_ANY_INF) == MASK_INF) { // y = inf
  336. y.w[1] = y.w[1] & (MASK_SIGN | MASK_INF);
  337. y.w[0] = 0x0ull;
  338. } else { // y is not special
  339. // check for non-canonical values - treated as zero
  340. if ((y.w[1] & MASK_STEERING_BITS) == MASK_STEERING_BITS) { // G0_G1=11
  341. // non-canonical
  342. y.w[1] = (y.w[1] & MASK_SIGN) | ((y.w[1] << 2) & MASK_EXP);
  343. y.w[0] = 0x0ull;
  344. } else { // G0_G1 != 11
  345. if ((y.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
  346. ((y.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull
  347. && y.w[0] > 0x378d8e63ffffffffull)) {
  348. // y is non-canonical if coefficient is larger than 10^34 -1
  349. y.w[1] = (y.w[1] & MASK_SIGN) | (y.w[1] & MASK_EXP);
  350. y.w[0] = 0x0ull;
  351. } else { // canonical
  352. ;
  353. }
  354. }
  355. }
  356. // NaN (CASE1)
  357. if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
  358. if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNaN
  359. // if x is SNAN, then return quiet (x)
  360. *pfpsf |= INVALID_EXCEPTION; // set exception if SNaN
  361. x.w[1] = x.w[1] & 0xfdffffffffffffffull; // quietize x
  362. res = x;
  363. } else { // x is QNaN
  364. if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NAN
  365. if ((y.w[1] & MASK_SNAN) == MASK_SNAN) { // y is SNAN
  366. *pfpsf |= INVALID_EXCEPTION; // set invalid flag
  367. }
  368. res = x;
  369. } else {
  370. res = y;
  371. }
  372. }
  373. BID_RETURN (res);
  374. } else if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NaN, but x is not
  375. if ((y.w[1] & MASK_SNAN) == MASK_SNAN) {
  376. *pfpsf |= INVALID_EXCEPTION; // set exception if SNaN
  377. y.w[1] = y.w[1] & 0xfdffffffffffffffull; // quietize y
  378. res = y;
  379. } else {
  380. // will return x (which is not NaN)
  381. res = x;
  382. }
  383. BID_RETURN (res);
  384. }
  385. // SIMPLE (CASE2)
  386. // if all the bits are the same, these numbers are equal (not Greater).
  387. if (x.w[0] == y.w[0] && x.w[1] == y.w[1]) {
  388. res = y;
  389. BID_RETURN (res);
  390. }
  391. // INFINITY (CASE3)
  392. if ((x.w[1] & MASK_INF) == MASK_INF) {
  393. // if x infinity, it has maximum magnitude.
  394. // Check if magnitudes are equal. If x is negative, return it.
  395. res = ((x.w[1] & MASK_SIGN) == MASK_SIGN
  396. && (y.w[1] & MASK_INF) == MASK_INF) ? x : y;
  397. BID_RETURN (res);
  398. } else if ((y.w[1] & MASK_INF) == MASK_INF) {
  399. // x is finite, so if y is infinity, then x is less in magnitude
  400. res = x;
  401. BID_RETURN (res);
  402. }
  403. // CONVERT X
  404. sig_x.w[1] = x.w[1] & 0x0001ffffffffffffull;
  405. sig_x.w[0] = x.w[0];
  406. exp_x = (x.w[1] >> 49) & 0x000000000003fffull;
  407. // CONVERT Y
  408. exp_y = (y.w[1] >> 49) & 0x0000000000003fffull;
  409. sig_y.w[1] = y.w[1] & 0x0001ffffffffffffull;
  410. sig_y.w[0] = y.w[0];
  411. // ZERO (CASE4)
  412. // some properties:
  413. // (+ZERO == -ZERO) => therefore ignore the sign
  414. // (ZERO x 10^A == ZERO x 10^B) for any valid A, B =>
  415. // therefore ignore the exponent field
  416. // (Any non-canonical # is considered 0)
  417. if ((sig_x.w[1] == 0) && (sig_x.w[0] == 0)) {
  418. res = x;
  419. BID_RETURN (res);
  420. }
  421. if ((sig_y.w[1] == 0) && (sig_y.w[0] == 0)) {
  422. res = y;
  423. BID_RETURN (res);
  424. }
  425. // REDUNDANT REPRESENTATIONS (CASE6)
  426. // check if exponents are the same and significands are the same
  427. if (exp_y == exp_x && sig_x.w[1] == sig_y.w[1]
  428. && sig_x.w[0] == sig_y.w[0]) {
  429. if (x.w[1] & 0x8000000000000000ull) { // x is negative
  430. res = x;
  431. BID_RETURN (res);
  432. } else {
  433. res = y;
  434. BID_RETURN (res);
  435. }
  436. } else if (((sig_x.w[1] > sig_y.w[1] || (sig_x.w[1] == sig_y.w[1]
  437. && sig_x.w[0] > sig_y.w[0]))
  438. && exp_x == exp_y)
  439. || ((sig_x.w[1] > sig_y.w[1]
  440. || (sig_x.w[1] == sig_y.w[1]
  441. && sig_x.w[0] >= sig_y.w[0]))
  442. && exp_x > exp_y)) {
  443. // if both components are either bigger or smaller, it is clear what
  444. // needs to be done; also if the magnitudes are equal
  445. res = y;
  446. BID_RETURN (res);
  447. } else if (((sig_y.w[1] > sig_x.w[1] || (sig_y.w[1] == sig_x.w[1]
  448. && sig_y.w[0] > sig_x.w[0]))
  449. && exp_y == exp_x)
  450. || ((sig_y.w[1] > sig_x.w[1]
  451. || (sig_y.w[1] == sig_x.w[1]
  452. && sig_y.w[0] >= sig_x.w[0]))
  453. && exp_y > exp_x)) {
  454. res = x;
  455. BID_RETURN (res);
  456. } else {
  457. ; // continue
  458. }
  459. diff = exp_x - exp_y;
  460. // if |exp_x - exp_y| < 33, it comes down to the compensated significand
  461. if (diff > 0) { // to simplify the loop below,
  462. // if exp_x is 33 greater than exp_y, no need for compensation
  463. if (diff > 33) {
  464. res = y; // difference cannot be greater than 10^33
  465. BID_RETURN (res);
  466. }
  467. if (diff > 19) { //128 by 128 bit multiply -> 256 bits
  468. __mul_128x128_to_256 (sig_n_prime256, sig_x, ten2k128[diff - 20]);
  469. // if positive, return whichever significand is larger
  470. // (converse if negative)
  471. if (sig_n_prime256.w[3] == 0 && (sig_n_prime256.w[2] == 0)
  472. && sig_n_prime256.w[1] == sig_y.w[1]
  473. && (sig_n_prime256.w[0] == sig_y.w[0])) {
  474. res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? y : x; // if equal
  475. BID_RETURN (res);
  476. }
  477. res = (((sig_n_prime256.w[3] > 0) || sig_n_prime256.w[2] > 0)
  478. || (sig_n_prime256.w[1] > sig_y.w[1])
  479. || (sig_n_prime256.w[1] == sig_y.w[1]
  480. && sig_n_prime256.w[0] > sig_y.w[0])) ? y : x;
  481. BID_RETURN (res);
  482. }
  483. __mul_64x128_to_192 (sig_n_prime192, ten2k64[diff], sig_x);
  484. // if positive, return whichever significand is larger
  485. // (converse if negative)
  486. if ((sig_n_prime192.w[2] == 0) && sig_n_prime192.w[1] == sig_y.w[1]
  487. && (sig_n_prime192.w[0] == sig_y.w[0])) {
  488. // if = in magnitude, return +, (if possible)
  489. res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? y : x;
  490. BID_RETURN (res);
  491. }
  492. res = ((sig_n_prime192.w[2] > 0)
  493. || (sig_n_prime192.w[1] > sig_y.w[1])
  494. || (sig_n_prime192.w[1] == sig_y.w[1]
  495. && sig_n_prime192.w[0] > sig_y.w[0])) ? y : x;
  496. BID_RETURN (res);
  497. }
  498. diff = exp_y - exp_x;
  499. // if exp_x is 33 less than exp_y, no need for compensation
  500. if (diff > 33) {
  501. res = x;
  502. BID_RETURN (res);
  503. }
  504. if (diff > 19) { //128 by 128 bit multiply -> 256 bits
  505. // adjust the y significand upwards
  506. __mul_128x128_to_256 (sig_n_prime256, sig_y, ten2k128[diff - 20]);
  507. // if positive, return whichever significand is larger
  508. // (converse if negative)
  509. if (sig_n_prime256.w[3] == 0 && (sig_n_prime256.w[2] == 0)
  510. && sig_n_prime256.w[1] == sig_x.w[1]
  511. && (sig_n_prime256.w[0] == sig_x.w[0])) {
  512. // if = in magnitude, return +, (if possible)
  513. res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? y : x;
  514. BID_RETURN (res);
  515. }
  516. res = (sig_n_prime256.w[3] == 0 && sig_n_prime256.w[2] == 0
  517. && (sig_n_prime256.w[1] < sig_x.w[1]
  518. || (sig_n_prime256.w[1] == sig_x.w[1]
  519. && sig_n_prime256.w[0] < sig_x.w[0]))) ? y : x;
  520. BID_RETURN (res);
  521. }
  522. // adjust the y significand upwards
  523. __mul_64x128_to_192 (sig_n_prime192, ten2k64[diff], sig_y);
  524. // if positive, return whichever significand is larger (converse if negative)
  525. if ((sig_n_prime192.w[2] == 0) && sig_n_prime192.w[1] == sig_x.w[1]
  526. && (sig_n_prime192.w[0] == sig_x.w[0])) {
  527. // if = in magnitude, return +, if possible)
  528. res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? y : x;
  529. BID_RETURN (res);
  530. }
  531. res = (sig_n_prime192.w[2] == 0
  532. && (sig_n_prime192.w[1] < sig_x.w[1]
  533. || (sig_n_prime192.w[1] == sig_x.w[1]
  534. && sig_n_prime192.w[0] < sig_x.w[0]))) ? y : x;
  535. BID_RETURN (res);
  536. }
  537. /*****************************************************************************
  538. * BID128 maximum function - returns greater of two numbers
  539. *****************************************************************************/
  540. #if DECIMAL_CALL_BY_REFERENCE
  541. void
  542. bid128_maxnum (UINT128 * pres, UINT128 * px,
  543. UINT128 * py _EXC_FLAGS_PARAM) {
  544. UINT128 x = *px;
  545. UINT128 y = *py;
  546. #else
  547. UINT128
  548. bid128_maxnum (UINT128 x, UINT128 y _EXC_FLAGS_PARAM) {
  549. #endif
  550. UINT128 res;
  551. int exp_x, exp_y;
  552. int diff;
  553. UINT128 sig_x, sig_y;
  554. UINT192 sig_n_prime192;
  555. UINT256 sig_n_prime256;
  556. char x_is_zero = 0, y_is_zero = 0;
  557. BID_SWAP128 (x);
  558. BID_SWAP128 (y);
  559. // check for non-canonical x
  560. if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
  561. x.w[1] = x.w[1] & 0xfe003fffffffffffull; // clear out G[6]-G[16]
  562. // check for non-canonical NaN payload
  563. if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
  564. (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
  565. (x.w[0] > 0x38c15b09ffffffffull))) {
  566. x.w[1] = x.w[1] & 0xffffc00000000000ull;
  567. x.w[0] = 0x0ull;
  568. }
  569. } else if ((x.w[1] & MASK_ANY_INF) == MASK_INF) { // x = inf
  570. x.w[1] = x.w[1] & (MASK_SIGN | MASK_INF);
  571. x.w[0] = 0x0ull;
  572. } else { // x is not special
  573. // check for non-canonical values - treated as zero
  574. if ((x.w[1] & MASK_STEERING_BITS) == MASK_STEERING_BITS) { // G0_G1=11
  575. // non-canonical
  576. x.w[1] = (x.w[1] & MASK_SIGN) | ((x.w[1] << 2) & MASK_EXP);
  577. x.w[0] = 0x0ull;
  578. } else { // G0_G1 != 11
  579. if ((x.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
  580. ((x.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull
  581. && x.w[0] > 0x378d8e63ffffffffull)) {
  582. // x is non-canonical if coefficient is larger than 10^34 -1
  583. x.w[1] = (x.w[1] & MASK_SIGN) | (x.w[1] & MASK_EXP);
  584. x.w[0] = 0x0ull;
  585. } else { // canonical
  586. ;
  587. }
  588. }
  589. }
  590. // check for non-canonical y
  591. if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NAN
  592. y.w[1] = y.w[1] & 0xfe003fffffffffffull; // clear out G[6]-G[16]
  593. // check for non-canonical NaN payload
  594. if (((y.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
  595. (((y.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
  596. (y.w[0] > 0x38c15b09ffffffffull))) {
  597. y.w[1] = y.w[1] & 0xffffc00000000000ull;
  598. y.w[0] = 0x0ull;
  599. }
  600. } else if ((y.w[1] & MASK_ANY_INF) == MASK_INF) { // y = inf
  601. y.w[1] = y.w[1] & (MASK_SIGN | MASK_INF);
  602. y.w[0] = 0x0ull;
  603. } else { // y is not special
  604. // check for non-canonical values - treated as zero
  605. if ((y.w[1] & MASK_STEERING_BITS) == MASK_STEERING_BITS) { // G0_G1=11
  606. // non-canonical
  607. y.w[1] = (y.w[1] & MASK_SIGN) | ((y.w[1] << 2) & MASK_EXP);
  608. y.w[0] = 0x0ull;
  609. } else { // G0_G1 != 11
  610. if ((y.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
  611. ((y.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull
  612. && y.w[0] > 0x378d8e63ffffffffull)) {
  613. // y is non-canonical if coefficient is larger than 10^34 -1
  614. y.w[1] = (y.w[1] & MASK_SIGN) | (y.w[1] & MASK_EXP);
  615. y.w[0] = 0x0ull;
  616. } else { // canonical
  617. ;
  618. }
  619. }
  620. }
  621. // NaN (CASE1)
  622. if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
  623. if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNaN
  624. // if x is SNAN, then return quiet (x)
  625. *pfpsf |= INVALID_EXCEPTION; // set exception if SNaN
  626. x.w[1] = x.w[1] & 0xfdffffffffffffffull; // quietize x
  627. res = x;
  628. } else { // x is QNaN
  629. if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NAN
  630. if ((y.w[1] & MASK_SNAN) == MASK_SNAN) { // y is SNAN
  631. *pfpsf |= INVALID_EXCEPTION; // set invalid flag
  632. }
  633. res = x;
  634. } else {
  635. res = y;
  636. }
  637. }
  638. BID_RETURN (res);
  639. } else if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NaN, but x is not
  640. if ((y.w[1] & MASK_SNAN) == MASK_SNAN) {
  641. *pfpsf |= INVALID_EXCEPTION; // set exception if SNaN
  642. y.w[1] = y.w[1] & 0xfdffffffffffffffull; // quietize y
  643. res = y;
  644. } else {
  645. // will return x (which is not NaN)
  646. res = x;
  647. }
  648. BID_RETURN (res);
  649. }
  650. // SIMPLE (CASE2)
  651. // if all the bits are the same, these numbers are equal (not Greater).
  652. if (x.w[0] == y.w[0] && x.w[1] == y.w[1]) {
  653. res = x;
  654. BID_RETURN (res);
  655. }
  656. // INFINITY (CASE3)
  657. if ((x.w[1] & MASK_INF) == MASK_INF) {
  658. res = ((x.w[1] & MASK_SIGN) == MASK_SIGN) ? y : x;
  659. BID_RETURN (res);
  660. } else if ((y.w[1] & MASK_INF) == MASK_INF) {
  661. // x is finite, so if y is positive infinity, then x is less, return 0
  662. // if y is negative infinity, then x is greater, return 1
  663. res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? x : y;
  664. BID_RETURN (res);
  665. }
  666. // CONVERT X
  667. sig_x.w[1] = x.w[1] & 0x0001ffffffffffffull;
  668. sig_x.w[0] = x.w[0];
  669. exp_x = (x.w[1] >> 49) & 0x000000000003fffull;
  670. // CONVERT Y
  671. exp_y = (y.w[1] >> 49) & 0x0000000000003fffull;
  672. sig_y.w[1] = y.w[1] & 0x0001ffffffffffffull;
  673. sig_y.w[0] = y.w[0];
  674. // ZERO (CASE4)
  675. // some properties:
  676. // (+ZERO == -ZERO) => therefore ignore the sign
  677. // (ZERO x 10^A == ZERO x 10^B) for any valid A, B =>
  678. // therefore ignore the exponent field
  679. // (Any non-canonical # is considered 0)
  680. if ((sig_x.w[1] == 0) && (sig_x.w[0] == 0)) {
  681. x_is_zero = 1;
  682. }
  683. if ((sig_y.w[1] == 0) && (sig_y.w[0] == 0)) {
  684. y_is_zero = 1;
  685. }
  686. if (x_is_zero && y_is_zero) {
  687. // if both numbers are zero, neither is greater => return either number
  688. res = x;
  689. BID_RETURN (res);
  690. } else if (x_is_zero) {
  691. // is x is zero, it is greater if Y is negative
  692. res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? x : y;
  693. BID_RETURN (res);
  694. } else if (y_is_zero) {
  695. // is y is zero, X is greater if it is positive
  696. res = ((x.w[1] & MASK_SIGN) != MASK_SIGN) ? x : y;
  697. BID_RETURN (res);
  698. }
  699. // OPPOSITE SIGN (CASE5)
  700. // now, if the sign bits differ, x is greater if y is negative
  701. if (((x.w[1] ^ y.w[1]) & MASK_SIGN) == MASK_SIGN) {
  702. res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? x : y;
  703. BID_RETURN (res);
  704. }
  705. // REDUNDANT REPRESENTATIONS (CASE6)
  706. // if exponents are the same, then we have a simple comparison of
  707. // the significands
  708. if (exp_y == exp_x) {
  709. res = (((sig_x.w[1] > sig_y.w[1]) || (sig_x.w[1] == sig_y.w[1] &&
  710. sig_x.w[0] >= sig_y.w[0])) ^
  711. ((x.w[1] & MASK_SIGN) == MASK_SIGN)) ? x : y;
  712. BID_RETURN (res);
  713. }
  714. // if both components are either bigger or smaller, it is clear what
  715. // needs to be done
  716. if ((sig_x.w[1] > sig_y.w[1]
  717. || (sig_x.w[1] == sig_y.w[1] && sig_x.w[0] > sig_y.w[0]))
  718. && exp_x >= exp_y) {
  719. res = ((x.w[1] & MASK_SIGN) != MASK_SIGN) ? x : y;
  720. BID_RETURN (res);
  721. }
  722. if ((sig_x.w[1] < sig_y.w[1]
  723. || (sig_x.w[1] == sig_y.w[1] && sig_x.w[0] < sig_y.w[0]))
  724. && exp_x <= exp_y) {
  725. res = ((x.w[1] & MASK_SIGN) == MASK_SIGN) ? x : y;
  726. BID_RETURN (res);
  727. }
  728. diff = exp_x - exp_y;
  729. // if |exp_x - exp_y| < 33, it comes down to the compensated significand
  730. if (diff > 0) { // to simplify the loop below,
  731. // if exp_x is 33 greater than exp_y, no need for compensation
  732. if (diff > 33) {
  733. // difference cannot be greater than 10^33
  734. res = ((x.w[1] & MASK_SIGN) != MASK_SIGN) ? x : y;
  735. BID_RETURN (res);
  736. }
  737. if (diff > 19) { //128 by 128 bit multiply -> 256 bits
  738. __mul_128x128_to_256 (sig_n_prime256, sig_x, ten2k128[diff - 20]);
  739. // if postitive, return whichever significand is larger
  740. // (converse if negative)
  741. res = ((((sig_n_prime256.w[3] > 0) || sig_n_prime256.w[2] > 0)
  742. || (sig_n_prime256.w[1] > sig_y.w[1])
  743. || (sig_n_prime256.w[1] == sig_y.w[1]
  744. && sig_n_prime256.w[0] >
  745. sig_y.w[0])) ^ ((y.w[1] & MASK_SIGN) ==
  746. MASK_SIGN)) ? x : y;
  747. BID_RETURN (res);
  748. }
  749. __mul_64x128_to_192 (sig_n_prime192, ten2k64[diff], sig_x);
  750. // if postitive, return whichever significand is larger
  751. // (converse if negative)
  752. res =
  753. (((sig_n_prime192.w[2] > 0) || (sig_n_prime192.w[1] > sig_y.w[1])
  754. || (sig_n_prime192.w[1] == sig_y.w[1]
  755. && sig_n_prime192.w[0] >
  756. sig_y.w[0])) ^ ((y.w[1] & MASK_SIGN) == MASK_SIGN)) ? x : y;
  757. BID_RETURN (res);
  758. }
  759. diff = exp_y - exp_x;
  760. // if exp_x is 33 less than exp_y, no need for compensation
  761. if (diff > 33) {
  762. res = ((x.w[1] & MASK_SIGN) == MASK_SIGN) ? x : y;
  763. BID_RETURN (res);
  764. }
  765. if (diff > 19) { //128 by 128 bit multiply -> 256 bits
  766. // adjust the y significand upwards
  767. __mul_128x128_to_256 (sig_n_prime256, sig_y, ten2k128[diff - 20]);
  768. // if postitive, return whichever significand is larger
  769. // (converse if negative)
  770. res =
  771. ((sig_n_prime256.w[3] != 0 || sig_n_prime256.w[2] != 0
  772. || (sig_n_prime256.w[1] > sig_x.w[1]
  773. || (sig_n_prime256.w[1] == sig_x.w[1]
  774. && sig_n_prime256.w[0] >
  775. sig_x.w[0]))) ^ ((x.w[1] & MASK_SIGN) !=
  776. MASK_SIGN)) ? x : y;
  777. BID_RETURN (res);
  778. }
  779. // adjust the y significand upwards
  780. __mul_64x128_to_192 (sig_n_prime192, ten2k64[diff], sig_y);
  781. // if postitive, return whichever significand is larger (converse if negative)
  782. res =
  783. ((sig_n_prime192.w[2] != 0
  784. || (sig_n_prime192.w[1] > sig_x.w[1]
  785. || (sig_n_prime192.w[1] == sig_x.w[1]
  786. && sig_n_prime192.w[0] >
  787. sig_x.w[0]))) ^ ((y.w[1] & MASK_SIGN) !=
  788. MASK_SIGN)) ? x : y;
  789. BID_RETURN (res);
  790. }
  791. /*****************************************************************************
  792. * BID128 maximum magnitude function - returns greater of two numbers
  793. *****************************************************************************/
  794. #if DECIMAL_CALL_BY_REFERENCE
  795. void
  796. bid128_maxnum_mag (UINT128 * pres, UINT128 * px,
  797. UINT128 * py _EXC_FLAGS_PARAM) {
  798. UINT128 x = *px;
  799. UINT128 y = *py;
  800. #else
  801. UINT128
  802. bid128_maxnum_mag (UINT128 x, UINT128 y _EXC_FLAGS_PARAM) {
  803. #endif
  804. UINT128 res;
  805. int exp_x, exp_y;
  806. int diff;
  807. UINT128 sig_x, sig_y;
  808. UINT192 sig_n_prime192;
  809. UINT256 sig_n_prime256;
  810. BID_SWAP128 (x);
  811. BID_SWAP128 (y);
  812. // check for non-canonical x
  813. if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
  814. x.w[1] = x.w[1] & 0xfe003fffffffffffull; // clear out G[6]-G[16]
  815. // check for non-canonical NaN payload
  816. if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
  817. (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
  818. (x.w[0] > 0x38c15b09ffffffffull))) {
  819. x.w[1] = x.w[1] & 0xffffc00000000000ull;
  820. x.w[0] = 0x0ull;
  821. }
  822. } else if ((x.w[1] & MASK_ANY_INF) == MASK_INF) { // x = inf
  823. x.w[1] = x.w[1] & (MASK_SIGN | MASK_INF);
  824. x.w[0] = 0x0ull;
  825. } else { // x is not special
  826. // check for non-canonical values - treated as zero
  827. if ((x.w[1] & MASK_STEERING_BITS) == MASK_STEERING_BITS) { // G0_G1=11
  828. // non-canonical
  829. x.w[1] = (x.w[1] & MASK_SIGN) | ((x.w[1] << 2) & MASK_EXP);
  830. x.w[0] = 0x0ull;
  831. } else { // G0_G1 != 11
  832. if ((x.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
  833. ((x.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull
  834. && x.w[0] > 0x378d8e63ffffffffull)) {
  835. // x is non-canonical if coefficient is larger than 10^34 -1
  836. x.w[1] = (x.w[1] & MASK_SIGN) | (x.w[1] & MASK_EXP);
  837. x.w[0] = 0x0ull;
  838. } else { // canonical
  839. ;
  840. }
  841. }
  842. }
  843. // check for non-canonical y
  844. if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NAN
  845. y.w[1] = y.w[1] & 0xfe003fffffffffffull; // clear out G[6]-G[16]
  846. // check for non-canonical NaN payload
  847. if (((y.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
  848. (((y.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
  849. (y.w[0] > 0x38c15b09ffffffffull))) {
  850. y.w[1] = y.w[1] & 0xffffc00000000000ull;
  851. y.w[0] = 0x0ull;
  852. }
  853. } else if ((y.w[1] & MASK_ANY_INF) == MASK_INF) { // y = inf
  854. y.w[1] = y.w[1] & (MASK_SIGN | MASK_INF);
  855. y.w[0] = 0x0ull;
  856. } else { // y is not special
  857. // check for non-canonical values - treated as zero
  858. if ((y.w[1] & MASK_STEERING_BITS) == MASK_STEERING_BITS) { // G0_G1=11
  859. // non-canonical
  860. y.w[1] = (y.w[1] & MASK_SIGN) | ((y.w[1] << 2) & MASK_EXP);
  861. y.w[0] = 0x0ull;
  862. } else { // G0_G1 != 11
  863. if ((y.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
  864. ((y.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull &&
  865. y.w[0] > 0x378d8e63ffffffffull)) {
  866. // y is non-canonical if coefficient is larger than 10^34 -1
  867. y.w[1] = (y.w[1] & MASK_SIGN) | (y.w[1] & MASK_EXP);
  868. y.w[0] = 0x0ull;
  869. } else { // canonical
  870. ;
  871. }
  872. }
  873. }
  874. // NaN (CASE1)
  875. if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
  876. if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNaN
  877. // if x is SNAN, then return quiet (x)
  878. *pfpsf |= INVALID_EXCEPTION; // set exception if SNaN
  879. x.w[1] = x.w[1] & 0xfdffffffffffffffull; // quietize x
  880. res = x;
  881. } else { // x is QNaN
  882. if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NAN
  883. if ((y.w[1] & MASK_SNAN) == MASK_SNAN) { // y is SNAN
  884. *pfpsf |= INVALID_EXCEPTION; // set invalid flag
  885. }
  886. res = x;
  887. } else {
  888. res = y;
  889. }
  890. }
  891. BID_RETURN (res);
  892. } else if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NaN, but x is not
  893. if ((y.w[1] & MASK_SNAN) == MASK_SNAN) {
  894. *pfpsf |= INVALID_EXCEPTION; // set exception if SNaN
  895. y.w[1] = y.w[1] & 0xfdffffffffffffffull; // quietize y
  896. res = y;
  897. } else {
  898. // will return x (which is not NaN)
  899. res = x;
  900. }
  901. BID_RETURN (res);
  902. }
  903. // SIMPLE (CASE2)
  904. // if all the bits are the same, these numbers are equal (not Greater).
  905. if (x.w[0] == y.w[0] && x.w[1] == y.w[1]) {
  906. res = y;
  907. BID_RETURN (res);
  908. }
  909. // INFINITY (CASE3)
  910. if ((x.w[1] & MASK_INF) == MASK_INF) {
  911. // if x infinity, it has maximum magnitude
  912. res = ((x.w[1] & MASK_SIGN) == MASK_SIGN
  913. && (y.w[1] & MASK_INF) == MASK_INF) ? y : x;
  914. BID_RETURN (res);
  915. } else if ((y.w[1] & MASK_INF) == MASK_INF) {
  916. // x is finite, so if y is positive infinity, then x is less, return 0
  917. // if y is negative infinity, then x is greater, return 1
  918. res = y;
  919. BID_RETURN (res);
  920. }
  921. // CONVERT X
  922. sig_x.w[1] = x.w[1] & 0x0001ffffffffffffull;
  923. sig_x.w[0] = x.w[0];
  924. exp_x = (x.w[1] >> 49) & 0x000000000003fffull;
  925. // CONVERT Y
  926. exp_y = (y.w[1] >> 49) & 0x0000000000003fffull;
  927. sig_y.w[1] = y.w[1] & 0x0001ffffffffffffull;
  928. sig_y.w[0] = y.w[0];
  929. // ZERO (CASE4)
  930. // some properties:
  931. // (+ZERO == -ZERO) => therefore ignore the sign
  932. // (ZERO x 10^A == ZERO x 10^B) for any valid A, B =>
  933. // therefore ignore the exponent field
  934. // (Any non-canonical # is considered 0)
  935. if ((sig_x.w[1] == 0) && (sig_x.w[0] == 0)) {
  936. res = y;
  937. BID_RETURN (res);
  938. }
  939. if ((sig_y.w[1] == 0) && (sig_y.w[0] == 0)) {
  940. res = x;
  941. BID_RETURN (res);
  942. }
  943. // REDUNDANT REPRESENTATIONS (CASE6)
  944. if (exp_y == exp_x && sig_x.w[1] == sig_y.w[1]
  945. && sig_x.w[0] == sig_y.w[0]) {
  946. // check if exponents are the same and significands are the same
  947. if (x.w[1] & 0x8000000000000000ull) { // x is negative
  948. res = y;
  949. BID_RETURN (res);
  950. } else {
  951. res = x;
  952. BID_RETURN (res);
  953. }
  954. } else if (((sig_x.w[1] > sig_y.w[1] || (sig_x.w[1] == sig_y.w[1]
  955. && sig_x.w[0] > sig_y.w[0]))
  956. && exp_x == exp_y)
  957. || ((sig_x.w[1] > sig_y.w[1]
  958. || (sig_x.w[1] == sig_y.w[1]
  959. && sig_x.w[0] >= sig_y.w[0]))
  960. && exp_x > exp_y)) {
  961. // if both components are either bigger or smaller, it is clear what
  962. // needs to be done; also if the magnitudes are equal
  963. res = x;
  964. BID_RETURN (res);
  965. } else if (((sig_y.w[1] > sig_x.w[1] || (sig_y.w[1] == sig_x.w[1]
  966. && sig_y.w[0] > sig_x.w[0]))
  967. && exp_y == exp_x)
  968. || ((sig_y.w[1] > sig_x.w[1]
  969. || (sig_y.w[1] == sig_x.w[1]
  970. && sig_y.w[0] >= sig_x.w[0]))
  971. && exp_y > exp_x)) {
  972. res = y;
  973. BID_RETURN (res);
  974. } else {
  975. ; // continue
  976. }
  977. diff = exp_x - exp_y;
  978. // if |exp_x - exp_y| < 33, it comes down to the compensated significand
  979. if (diff > 0) { // to simplify the loop below,
  980. // if exp_x is 33 greater than exp_y, no need for compensation
  981. if (diff > 33) {
  982. res = x; // difference cannot be greater than 10^33
  983. BID_RETURN (res);
  984. }
  985. if (diff > 19) { //128 by 128 bit multiply -> 256 bits
  986. __mul_128x128_to_256 (sig_n_prime256, sig_x, ten2k128[diff - 20]);
  987. // if postitive, return whichever significand is larger
  988. // (converse if negative)
  989. if (sig_n_prime256.w[3] == 0 && (sig_n_prime256.w[2] == 0)
  990. && sig_n_prime256.w[1] == sig_y.w[1]
  991. && (sig_n_prime256.w[0] == sig_y.w[0])) {
  992. res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? x : y; // if equal
  993. BID_RETURN (res);
  994. }
  995. res = (((sig_n_prime256.w[3] > 0) || sig_n_prime256.w[2] > 0)
  996. || (sig_n_prime256.w[1] > sig_y.w[1])
  997. || (sig_n_prime256.w[1] == sig_y.w[1]
  998. && sig_n_prime256.w[0] > sig_y.w[0])) ? x : y;
  999. BID_RETURN (res);
  1000. }
  1001. __mul_64x128_to_192 (sig_n_prime192, ten2k64[diff], sig_x);
  1002. // if postitive, return whichever significand is larger (converse if negative)
  1003. if ((sig_n_prime192.w[2] == 0) && sig_n_prime192.w[1] == sig_y.w[1]
  1004. && (sig_n_prime192.w[0] == sig_y.w[0])) {
  1005. // if equal, return positive magnitude
  1006. res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? x : y;
  1007. BID_RETURN (res);
  1008. }
  1009. res = ((sig_n_prime192.w[2] > 0)
  1010. || (sig_n_prime192.w[1] > sig_y.w[1])
  1011. || (sig_n_prime192.w[1] == sig_y.w[1]
  1012. && sig_n_prime192.w[0] > sig_y.w[0])) ? x : y;
  1013. BID_RETURN (res);
  1014. }
  1015. diff = exp_y - exp_x;
  1016. // if exp_x is 33 less than exp_y, no need for compensation
  1017. if (diff > 33) {
  1018. res = y;
  1019. BID_RETURN (res);
  1020. }
  1021. if (diff > 19) { //128 by 128 bit multiply -> 256 bits
  1022. // adjust the y significand upwards
  1023. __mul_128x128_to_256 (sig_n_prime256, sig_y, ten2k128[diff - 20]);
  1024. // if postitive, return whichever significand is larger
  1025. // (converse if negative)
  1026. if (sig_n_prime256.w[3] == 0 && (sig_n_prime256.w[2] == 0)
  1027. && sig_n_prime256.w[1] == sig_x.w[1]
  1028. && (sig_n_prime256.w[0] == sig_x.w[0])) {
  1029. // if equal, return positive (if possible)
  1030. res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? x : y;
  1031. BID_RETURN (res);
  1032. }
  1033. res = (sig_n_prime256.w[3] == 0 && sig_n_prime256.w[2] == 0
  1034. && (sig_n_prime256.w[1] < sig_x.w[1]
  1035. || (sig_n_prime256.w[1] == sig_x.w[1]
  1036. && sig_n_prime256.w[0] < sig_x.w[0]))) ? x : y;
  1037. BID_RETURN (res);
  1038. }
  1039. // adjust the y significand upwards
  1040. __mul_64x128_to_192 (sig_n_prime192, ten2k64[diff], sig_y);
  1041. // if postitive, return whichever significand is larger (converse if negative)
  1042. if ((sig_n_prime192.w[2] == 0) && sig_n_prime192.w[1] == sig_x.w[1]
  1043. && (sig_n_prime192.w[0] == sig_x.w[0])) {
  1044. // if equal, return positive (if possible)
  1045. res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? x : y;
  1046. BID_RETURN (res);
  1047. }
  1048. res = (sig_n_prime192.w[2] == 0
  1049. && (sig_n_prime192.w[1] < sig_x.w[1]
  1050. || (sig_n_prime192.w[1] == sig_x.w[1]
  1051. && sig_n_prime192.w[0] < sig_x.w[0]))) ? x : y;
  1052. BID_RETURN (res);
  1053. }