ieee754.c 8.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325
  1. /*
  2. ** 2013-04-17
  3. **
  4. ** The author disclaims copyright to this source code. In place of
  5. ** a legal notice, here is a blessing:
  6. **
  7. ** May you do good and not evil.
  8. ** May you find forgiveness for yourself and forgive others.
  9. ** May you share freely, never taking more than you give.
  10. **
  11. ******************************************************************************
  12. **
  13. ** This SQLite extension implements functions for the exact display
  14. ** and input of IEEE754 Binary64 floating-point numbers.
  15. **
  16. ** ieee754(X)
  17. ** ieee754(Y,Z)
  18. **
  19. ** In the first form, the value X should be a floating-point number.
  20. ** The function will return a string of the form 'ieee754(Y,Z)' where
  21. ** Y and Z are integers such that X==Y*pow(2,Z).
  22. **
  23. ** In the second form, Y and Z are integers which are the mantissa and
  24. ** base-2 exponent of a new floating point number. The function returns
  25. ** a floating-point value equal to Y*pow(2,Z).
  26. **
  27. ** Examples:
  28. **
  29. ** ieee754(2.0) -> 'ieee754(2,0)'
  30. ** ieee754(45.25) -> 'ieee754(181,-2)'
  31. ** ieee754(2, 0) -> 2.0
  32. ** ieee754(181, -2) -> 45.25
  33. **
  34. ** Two additional functions break apart the one-argument ieee754()
  35. ** result into separate integer values:
  36. **
  37. ** ieee754_mantissa(45.25) -> 181
  38. ** ieee754_exponent(45.25) -> -2
  39. **
  40. ** These functions convert binary64 numbers into blobs and back again.
  41. **
  42. ** ieee754_from_blob(x'3ff0000000000000') -> 1.0
  43. ** ieee754_to_blob(1.0) -> x'3ff0000000000000'
  44. **
  45. ** In all single-argument functions, if the argument is an 8-byte blob
  46. ** then that blob is interpreted as a big-endian binary64 value.
  47. **
  48. **
  49. ** EXACT DECIMAL REPRESENTATION OF BINARY64 VALUES
  50. ** -----------------------------------------------
  51. **
  52. ** This extension in combination with the separate 'decimal' extension
  53. ** can be used to compute the exact decimal representation of binary64
  54. ** values. To begin, first compute a table of exponent values:
  55. **
  56. ** CREATE TABLE pow2(x INTEGER PRIMARY KEY, v TEXT);
  57. ** WITH RECURSIVE c(x,v) AS (
  58. ** VALUES(0,'1')
  59. ** UNION ALL
  60. ** SELECT x+1, decimal_mul(v,'2') FROM c WHERE x+1<=971
  61. ** ) INSERT INTO pow2(x,v) SELECT x, v FROM c;
  62. ** WITH RECURSIVE c(x,v) AS (
  63. ** VALUES(-1,'0.5')
  64. ** UNION ALL
  65. ** SELECT x-1, decimal_mul(v,'0.5') FROM c WHERE x-1>=-1075
  66. ** ) INSERT INTO pow2(x,v) SELECT x, v FROM c;
  67. **
  68. ** Then, to compute the exact decimal representation of a floating
  69. ** point value (the value 47.49 is used in the example) do:
  70. **
  71. ** WITH c(n) AS (VALUES(47.49))
  72. ** ---------------^^^^^---- Replace with whatever you want
  73. ** SELECT decimal_mul(ieee754_mantissa(c.n),pow2.v)
  74. ** FROM pow2, c WHERE pow2.x=ieee754_exponent(c.n);
  75. **
  76. ** Here is a query to show various boundry values for the binary64
  77. ** number format:
  78. **
  79. ** WITH c(name,bin) AS (VALUES
  80. ** ('minimum positive value', x'0000000000000001'),
  81. ** ('maximum subnormal value', x'000fffffffffffff'),
  82. ** ('mininum positive nornal value', x'0010000000000000'),
  83. ** ('maximum value', x'7fefffffffffffff'))
  84. ** SELECT c.name, decimal_mul(ieee754_mantissa(c.bin),pow2.v)
  85. ** FROM pow2, c WHERE pow2.x=ieee754_exponent(c.bin);
  86. **
  87. */
  88. #include "sqlite3ext.h"
  89. SQLITE_EXTENSION_INIT1
  90. #include <assert.h>
  91. #include <string.h>
  92. /* Mark a function parameter as unused, to suppress nuisance compiler
  93. ** warnings. */
  94. #ifndef UNUSED_PARAMETER
  95. # define UNUSED_PARAMETER(X) (void)(X)
  96. #endif
  97. /*
  98. ** Implementation of the ieee754() function
  99. */
  100. static void ieee754func(
  101. sqlite3_context *context,
  102. int argc,
  103. sqlite3_value **argv
  104. ){
  105. if( argc==1 ){
  106. sqlite3_int64 m, a;
  107. double r;
  108. int e;
  109. int isNeg;
  110. char zResult[100];
  111. assert( sizeof(m)==sizeof(r) );
  112. if( sqlite3_value_type(argv[0])==SQLITE_BLOB
  113. && sqlite3_value_bytes(argv[0])==sizeof(r)
  114. ){
  115. const unsigned char *x = sqlite3_value_blob(argv[0]);
  116. unsigned int i;
  117. sqlite3_uint64 v = 0;
  118. for(i=0; i<sizeof(r); i++){
  119. v = (v<<8) | x[i];
  120. }
  121. memcpy(&r, &v, sizeof(r));
  122. }else{
  123. r = sqlite3_value_double(argv[0]);
  124. }
  125. if( r<0.0 ){
  126. isNeg = 1;
  127. r = -r;
  128. }else{
  129. isNeg = 0;
  130. }
  131. memcpy(&a,&r,sizeof(a));
  132. if( a==0 ){
  133. e = 0;
  134. m = 0;
  135. }else{
  136. e = a>>52;
  137. m = a & ((((sqlite3_int64)1)<<52)-1);
  138. if( e==0 ){
  139. m <<= 1;
  140. }else{
  141. m |= ((sqlite3_int64)1)<<52;
  142. }
  143. while( e<1075 && m>0 && (m&1)==0 ){
  144. m >>= 1;
  145. e++;
  146. }
  147. if( isNeg ) m = -m;
  148. }
  149. switch( *(int*)sqlite3_user_data(context) ){
  150. case 0:
  151. sqlite3_snprintf(sizeof(zResult), zResult, "ieee754(%lld,%d)",
  152. m, e-1075);
  153. sqlite3_result_text(context, zResult, -1, SQLITE_TRANSIENT);
  154. break;
  155. case 1:
  156. sqlite3_result_int64(context, m);
  157. break;
  158. case 2:
  159. sqlite3_result_int(context, e-1075);
  160. break;
  161. }
  162. }else{
  163. sqlite3_int64 m, e, a;
  164. double r;
  165. int isNeg = 0;
  166. m = sqlite3_value_int64(argv[0]);
  167. e = sqlite3_value_int64(argv[1]);
  168. /* Limit the range of e. Ticket 22dea1cfdb9151e4 2021-03-02 */
  169. if( e>10000 ){
  170. e = 10000;
  171. }else if( e<-10000 ){
  172. e = -10000;
  173. }
  174. if( m<0 ){
  175. isNeg = 1;
  176. m = -m;
  177. if( m<0 ) return;
  178. }else if( m==0 && e>-1000 && e<1000 ){
  179. sqlite3_result_double(context, 0.0);
  180. return;
  181. }
  182. while( (m>>32)&0xffe00000 ){
  183. m >>= 1;
  184. e++;
  185. }
  186. while( m!=0 && ((m>>32)&0xfff00000)==0 ){
  187. m <<= 1;
  188. e--;
  189. }
  190. e += 1075;
  191. if( e<=0 ){
  192. /* Subnormal */
  193. if( 1-e >= 64 ){
  194. m = 0;
  195. }else{
  196. m >>= 1-e;
  197. }
  198. e = 0;
  199. }else if( e>0x7ff ){
  200. e = 0x7ff;
  201. }
  202. a = m & ((((sqlite3_int64)1)<<52)-1);
  203. a |= e<<52;
  204. if( isNeg ) a |= ((sqlite3_uint64)1)<<63;
  205. memcpy(&r, &a, sizeof(r));
  206. sqlite3_result_double(context, r);
  207. }
  208. }
  209. /*
  210. ** Functions to convert between blobs and floats.
  211. */
  212. static void ieee754func_from_blob(
  213. sqlite3_context *context,
  214. int argc,
  215. sqlite3_value **argv
  216. ){
  217. UNUSED_PARAMETER(argc);
  218. if( sqlite3_value_type(argv[0])==SQLITE_BLOB
  219. && sqlite3_value_bytes(argv[0])==sizeof(double)
  220. ){
  221. double r;
  222. const unsigned char *x = sqlite3_value_blob(argv[0]);
  223. unsigned int i;
  224. sqlite3_uint64 v = 0;
  225. for(i=0; i<sizeof(r); i++){
  226. v = (v<<8) | x[i];
  227. }
  228. memcpy(&r, &v, sizeof(r));
  229. sqlite3_result_double(context, r);
  230. }
  231. }
  232. static void ieee754func_to_blob(
  233. sqlite3_context *context,
  234. int argc,
  235. sqlite3_value **argv
  236. ){
  237. UNUSED_PARAMETER(argc);
  238. if( sqlite3_value_type(argv[0])==SQLITE_FLOAT
  239. || sqlite3_value_type(argv[0])==SQLITE_INTEGER
  240. ){
  241. double r = sqlite3_value_double(argv[0]);
  242. sqlite3_uint64 v;
  243. unsigned char a[sizeof(r)];
  244. unsigned int i;
  245. memcpy(&v, &r, sizeof(r));
  246. for(i=1; i<=sizeof(r); i++){
  247. a[sizeof(r)-i] = v&0xff;
  248. v >>= 8;
  249. }
  250. sqlite3_result_blob(context, a, sizeof(r), SQLITE_TRANSIENT);
  251. }
  252. }
  253. /*
  254. ** SQL Function: ieee754_inc(r,N)
  255. **
  256. ** Move the floating point value r by N quantums and return the new
  257. ** values.
  258. **
  259. ** Behind the scenes: this routine merely casts r into a 64-bit unsigned
  260. ** integer, adds N, then casts the value back into float.
  261. **
  262. ** Example: To find the smallest positive number:
  263. **
  264. ** SELECT ieee754_inc(0.0,+1);
  265. */
  266. static void ieee754inc(
  267. sqlite3_context *context,
  268. int argc,
  269. sqlite3_value **argv
  270. ){
  271. double r;
  272. sqlite3_int64 N;
  273. sqlite3_uint64 m1, m2;
  274. double r2;
  275. UNUSED_PARAMETER(argc);
  276. r = sqlite3_value_double(argv[0]);
  277. N = sqlite3_value_int64(argv[1]);
  278. memcpy(&m1, &r, 8);
  279. m2 = m1 + N;
  280. memcpy(&r2, &m2, 8);
  281. sqlite3_result_double(context, r2);
  282. }
  283. #ifdef _WIN32
  284. __declspec(dllexport)
  285. #endif
  286. int sqlite3_ieee_init(
  287. sqlite3 *db,
  288. char **pzErrMsg,
  289. const sqlite3_api_routines *pApi
  290. ){
  291. static const struct {
  292. char *zFName;
  293. int nArg;
  294. int iAux;
  295. void (*xFunc)(sqlite3_context*,int,sqlite3_value**);
  296. } aFunc[] = {
  297. { "ieee754", 1, 0, ieee754func },
  298. { "ieee754", 2, 0, ieee754func },
  299. { "ieee754_mantissa", 1, 1, ieee754func },
  300. { "ieee754_exponent", 1, 2, ieee754func },
  301. { "ieee754_to_blob", 1, 0, ieee754func_to_blob },
  302. { "ieee754_from_blob", 1, 0, ieee754func_from_blob },
  303. { "ieee754_inc", 2, 0, ieee754inc },
  304. };
  305. unsigned int i;
  306. int rc = SQLITE_OK;
  307. SQLITE_EXTENSION_INIT2(pApi);
  308. (void)pzErrMsg; /* Unused parameter */
  309. for(i=0; i<sizeof(aFunc)/sizeof(aFunc[0]) && rc==SQLITE_OK; i++){
  310. rc = sqlite3_create_function(db, aFunc[i].zFName, aFunc[i].nArg,
  311. SQLITE_UTF8|SQLITE_INNOCUOUS,
  312. (void*)&aFunc[i].iAux,
  313. aFunc[i].xFunc, 0, 0);
  314. }
  315. return rc;
  316. }