Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | Download | RSS feed

  1. /****************************************************************
  2.  *
  3.  * The author of this software is David M. Gay.
  4.  *
  5.  * Copyright (c) 1991 by AT&T.
  6.  *
  7.  * Permission to use, copy, modify, and distribute this software for any
  8.  * purpose without fee is hereby granted, provided that this entire notice
  9.  * is included in all copies of any software which is or includes a copy
  10.  * or modification of this software and in all copies of the supporting
  11.  * documentation for such software.
  12.  *
  13.  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
  14.  * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY
  15.  * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
  16.  * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
  17.  *
  18.  ***************************************************************/
  19.  
  20. /* Please send bug reports to
  21.         David M. Gay
  22.         AT&T Bell Laboratories, Room 2C-463
  23.         600 Mountain Avenue
  24.         Murray Hill, NJ 07974-2070
  25.         U.S.A.
  26.         dmg@research.att.com or research!dmg
  27.  */
  28.  
  29. #include <ieeefp.h>
  30. #include <math.h>
  31. #include <float.h>
  32. #include <errno.h>
  33. #include <sys/config.h>
  34. #include <sys/types.h>
  35.  
  36. #ifdef __IEEE_LITTLE_ENDIAN
  37. #define IEEE_8087
  38. #endif
  39.  
  40. #ifdef __IEEE_BIG_ENDIAN
  41. #define IEEE_MC68k
  42. #endif
  43.  
  44. #ifdef __Z8000__
  45. #define Just_16
  46. #endif
  47.  
  48. #ifdef DEBUG
  49. #include "stdio.h"
  50. #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
  51. #endif
  52.  
  53. #ifdef Unsigned_Shifts
  54. #define Sign_Extend(a,b) if (b < 0) a |= (__uint32_t)0xffff0000;
  55. #else
  56. #define Sign_Extend(a,b) /*no-op*/
  57. #endif
  58.  
  59. #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
  60. Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
  61. #endif
  62.  
  63. /* If we are going to examine or modify specific bits in a double using
  64.    the word0 and/or word1 macros, then we must wrap the double inside
  65.    a union.  This is necessary to avoid undefined behavior according to
  66.    the ANSI C spec.  */
  67. union double_union
  68. {
  69.   double d;
  70.   __uint32_t i[2];
  71. };
  72.  
  73. #ifdef IEEE_8087
  74. #define word0(x) (x.i[1])
  75. #define word1(x) (x.i[0])
  76. #else
  77. #define word0(x) (x.i[0])
  78. #define word1(x) (x.i[1])
  79. #endif
  80.  
  81. /* The following is taken from gdtoaimp.h for use with new strtod, but
  82.    adjusted to avoid invalid type-punning.  */
  83. typedef __int32_t Long;
  84.  
  85. /* Unfortunately, because __ULong might be a different type than
  86.    __uint32_t, we can't re-use union double_union as-is without
  87.    further edits in strtod.c.  */
  88. typedef union { double d; __ULong i[2]; } U;
  89.  
  90. #define dword0(x) word0(x)
  91. #define dword1(x) word1(x)
  92. #define dval(x) (x.d)
  93.  
  94. #undef SI
  95. #ifdef Sudden_Underflow
  96. #define SI 1
  97. #else
  98. #define SI 0
  99. #endif
  100.  
  101. #define Storeinc(a,b,c) (*(a)++ = (b) << 16 | (c) & 0xffff)
  102.  
  103. /* #define P DBL_MANT_DIG */
  104. /* Ten_pmax = floor(P*log(2)/log(5)) */
  105. /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
  106. /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
  107. /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
  108.  
  109. #if defined(IEEE_8087) + defined(IEEE_MC68k)
  110. #if defined (_DOUBLE_IS_32BITS)
  111. #define Exp_shift   23
  112. #define Exp_shift1  23
  113. #define Exp_msk1    ((__uint32_t)0x00800000L)
  114. #define Exp_msk11   ((__uint32_t)0x00800000L)
  115. #define Exp_mask    ((__uint32_t)0x7f800000L)
  116. #define P           24
  117. #define Bias        127
  118. #define NO_HEX_FP   /* not supported in this case */
  119. #define IEEE_Arith
  120. #define Emin        (-126)
  121. #define Exp_1       ((__uint32_t)0x3f800000L)
  122. #define Exp_11      ((__uint32_t)0x3f800000L)
  123. #define Ebits       8
  124. #define Frac_mask   ((__uint32_t)0x007fffffL)
  125. #define Frac_mask1  ((__uint32_t)0x007fffffL)
  126. #define Ten_pmax    10
  127. #define Sign_bit    ((__uint32_t)0x80000000L)
  128. #define Ten_pmax    10
  129. #define Bletch      2
  130. #define Bndry_mask  ((__uint32_t)0x007fffffL)
  131. #define Bndry_mask1 ((__uint32_t)0x007fffffL)
  132. #define LSB 1
  133. #define Sign_bit    ((__uint32_t)0x80000000L)
  134. #define Log2P       1
  135. #define Tiny0       0
  136. #define Tiny1       1
  137. #define Quick_max   5
  138. #define Int_max     6
  139. #define Infinite(x) (word0(x) == ((__uint32_t)0x7f800000L))
  140. #undef word0
  141. #undef word1
  142. #undef dword0
  143. #undef dword1
  144.  
  145. #define word0(x) (x.i[0])
  146. #define word1(x) 0
  147. #define dword0(x) word0(x)
  148. #define dword1(x) 0
  149. #else
  150.  
  151. #define Exp_shift  20
  152. #define Exp_shift1 20
  153. #define Exp_msk1    ((__uint32_t)0x100000L)
  154. #define Exp_msk11   ((__uint32_t)0x100000L)
  155. #define Exp_mask  ((__uint32_t)0x7ff00000L)
  156. #define P 53
  157. #define Bias 1023
  158. #define IEEE_Arith
  159. #define Emin (-1022)
  160. #define Exp_1  ((__uint32_t)0x3ff00000L)
  161. #define Exp_11 ((__uint32_t)0x3ff00000L)
  162. #define Ebits 11
  163. #define Frac_mask  ((__uint32_t)0xfffffL)
  164. #define Frac_mask1 ((__uint32_t)0xfffffL)
  165. #define Ten_pmax 22
  166. #define Bletch 0x10
  167. #define Bndry_mask  ((__uint32_t)0xfffffL)
  168. #define Bndry_mask1 ((__uint32_t)0xfffffL)
  169. #define LSB 1
  170. #define Sign_bit ((__uint32_t)0x80000000L)
  171. #define Log2P 1
  172. #define Tiny0 0
  173. #define Tiny1 1
  174. #define Quick_max 14
  175. #define Int_max 14
  176. #define Infinite(x) (word0(x) == ((__uint32_t)0x7ff00000L)) /* sufficient test for here */
  177.  
  178. #endif /* !_DOUBLE_IS_32BITS */
  179.  
  180. #ifndef Flt_Rounds
  181. #ifdef FLT_ROUNDS
  182. #define Flt_Rounds FLT_ROUNDS
  183. #else
  184. #define Flt_Rounds 1
  185. #endif
  186. #endif /*Flt_Rounds*/
  187.  
  188. #else /* !IEEE_8087 && !IEEE_MC68k */
  189. #undef  Sudden_Underflow
  190. #define Sudden_Underflow
  191. #ifdef IBM
  192. #define Flt_Rounds 0
  193. #define Exp_shift  24
  194. #define Exp_shift1 24
  195. #define Exp_msk1   ((__uint32_t)0x1000000L)
  196. #define Exp_msk11  ((__uint32_t)0x1000000L)
  197. #define Exp_mask  ((__uint32_t)0x7f000000L)
  198. #define P 14
  199. #define Bias 65
  200. #define Exp_1  ((__uint32_t)0x41000000L)
  201. #define Exp_11 ((__uint32_t)0x41000000L)
  202. #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
  203. #define Frac_mask  ((__uint32_t)0xffffffL)
  204. #define Frac_mask1 ((__uint32_t)0xffffffL)
  205. #define Bletch 4
  206. #define Ten_pmax 22
  207. #define Bndry_mask  ((__uint32_t)0xefffffL)
  208. #define Bndry_mask1 ((__uint32_t)0xffffffL)
  209. #define LSB 1
  210. #define Sign_bit ((__uint32_t)0x80000000L)
  211. #define Log2P 4
  212. #define Tiny0 ((__uint32_t)0x100000L)
  213. #define Tiny1 0
  214. #define Quick_max 14
  215. #define Int_max 15
  216. #else /* VAX */
  217. #define Flt_Rounds 1
  218. #define Exp_shift  23
  219. #define Exp_shift1 7
  220. #define Exp_msk1    0x80
  221. #define Exp_msk11   ((__uint32_t)0x800000L)
  222. #define Exp_mask  ((__uint32_t)0x7f80L)
  223. #define P 56
  224. #define Bias 129
  225. #define Exp_1  ((__uint32_t)0x40800000L)
  226. #define Exp_11 ((__uint32_t)0x4080L)
  227. #define Ebits 8
  228. #define Frac_mask  ((__uint32_t)0x7fffffL)
  229. #define Frac_mask1 ((__uint32_t)0xffff007fL)
  230. #define Ten_pmax 24
  231. #define Bletch 2
  232. #define Bndry_mask  ((__uint32_t)0xffff007fL)
  233. #define Bndry_mask1 ((__uint32_t)0xffff007fL)
  234. #define LSB ((__uint32_t)0x10000L)
  235. #define Sign_bit ((__uint32_t)0x8000L)
  236. #define Log2P 1
  237. #define Tiny0 0x80
  238. #define Tiny1 0
  239. #define Quick_max 15
  240. #define Int_max 15
  241. #endif
  242. #endif
  243.  
  244. #ifndef IEEE_Arith
  245. #define ROUND_BIASED
  246. #else
  247. #define Scale_Bit 0x10
  248. #if defined(_DOUBLE_IS_32BITS) && defined(__v800)
  249. #define n_bigtens 2
  250. #else
  251. #define n_bigtens 5
  252. #endif
  253. #endif
  254.  
  255. #ifdef IBM
  256. #define n_bigtens 3
  257. #endif
  258.  
  259. #ifdef VAX
  260. #define n_bigtens 2
  261. #endif
  262.  
  263. #ifndef __NO_INFNAN_CHECK
  264. #define INFNAN_CHECK
  265. #endif
  266.  
  267. /*
  268.  * NAN_WORD0 and NAN_WORD1 are only referenced in strtod.c.  Prior to
  269.  * 20050115, they used to be hard-wired here (to 0x7ff80000 and 0,
  270.  * respectively), but now are determined by compiling and running
  271.  * qnan.c to generate gd_qnan.h, which specifies d_QNAN0 and d_QNAN1.
  272.  * Formerly gdtoaimp.h recommended supplying suitable -DNAN_WORD0=...
  273.  * and -DNAN_WORD1=...  values if necessary.  This should still work.
  274.  * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
  275.  */
  276. #ifdef IEEE_Arith
  277. #ifdef IEEE_MC68k
  278. #define _0 0
  279. #define _1 1
  280. #ifndef NAN_WORD0
  281. #define NAN_WORD0 d_QNAN0
  282. #endif
  283. #ifndef NAN_WORD1
  284. #define NAN_WORD1 d_QNAN1
  285. #endif
  286. #else
  287. #define _0 1
  288. #define _1 0
  289. #ifndef NAN_WORD0
  290. #define NAN_WORD0 d_QNAN1
  291. #endif
  292. #ifndef NAN_WORD1
  293. #define NAN_WORD1 d_QNAN0
  294. #endif
  295. #endif
  296. #else
  297. #undef INFNAN_CHECK
  298. #endif
  299.  
  300. #ifdef RND_PRODQUOT
  301. #define rounded_product(a,b) a = rnd_prod(a, b)
  302. #define rounded_quotient(a,b) a = rnd_quot(a, b)
  303. #ifdef KR_headers
  304. extern double rnd_prod(), rnd_quot();
  305. #else
  306. extern double rnd_prod(double, double), rnd_quot(double, double);
  307. #endif
  308. #else
  309. #define rounded_product(a,b) a *= b
  310. #define rounded_quotient(a,b) a /= b
  311. #endif
  312.  
  313. #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
  314. #define Big1 ((__uint32_t)0xffffffffL)
  315.  
  316. #ifndef Just_16
  317. /* When Pack_32 is not defined, we store 16 bits per 32-bit long.
  318.  * This makes some inner loops simpler and sometimes saves work
  319.  * during multiplications, but it often seems to make things slightly
  320.  * slower.  Hence the default is now to store 32 bits per long.
  321.  */
  322.  
  323.  #ifndef Pack_32
  324.   #define Pack_32
  325.  #endif
  326. #else  /* Just_16 */
  327.  #ifndef Pack_16
  328.   #define Pack_16
  329.  #endif
  330. #endif /* Just_16 */
  331.  
  332. #ifdef Pack_32
  333. #define ULbits 32
  334. #define kshift 5
  335. #define kmask 31
  336. #define ALL_ON 0xffffffff
  337. #else
  338. #define ULbits 16
  339. #define kshift 4
  340. #define kmask 15
  341. #define ALL_ON 0xffff
  342. #endif
  343.  
  344. #ifdef __cplusplus
  345. extern "C" double strtod(const char *s00, char **se);
  346. extern "C" char *dtoa(double d, int mode, int ndigits,
  347.                         int *decpt, int *sign, char **rve);
  348. #endif
  349.  
  350.  
  351. typedef struct _Bigint _Bigint;
  352.  
  353. #define Balloc  _Balloc
  354. #define Bfree   _Bfree
  355. #define multadd __multadd
  356. #define s2b     __s2b
  357. #define lo0bits __lo0bits
  358. #define hi0bits __hi0bits
  359. #define i2b     __i2b
  360. #define mult    __multiply
  361. #define pow5mult        __pow5mult
  362. #define lshift  __lshift
  363. #define cmp     __mcmp
  364. #define diff    __mdiff
  365. #define ulp     __ulp
  366. #define b2d     __b2d
  367. #define d2b     __d2b
  368. #define ratio   __ratio
  369. #define any_on  __any_on
  370. #define gethex  __gethex
  371. #define copybits        __copybits
  372. #define hexnan  __hexnan
  373. #define hexdig_init     __hexdig_init
  374.  
  375. #define hexdig  __hexdig
  376.  
  377. #define tens __mprec_tens
  378. #define bigtens __mprec_bigtens
  379. #define tinytens __mprec_tinytens
  380.  
  381. struct _reent ;
  382. struct FPI;
  383. double          _EXFUN(ulp,(double x));
  384. double          _EXFUN(b2d,(_Bigint *a , int *e));
  385. _Bigint *       _EXFUN(Balloc,(struct _reent *p, int k));
  386. void            _EXFUN(Bfree,(struct _reent *p, _Bigint *v));
  387. _Bigint *       _EXFUN(multadd,(struct _reent *p, _Bigint *, int, int));
  388. _Bigint *       _EXFUN(s2b,(struct _reent *, const char*, int, int, __ULong));
  389. _Bigint *       _EXFUN(i2b,(struct _reent *,int));
  390. _Bigint *       _EXFUN(mult, (struct _reent *, _Bigint *, _Bigint *));
  391. _Bigint *       _EXFUN(pow5mult, (struct _reent *, _Bigint *, int k));
  392. int             _EXFUN(hi0bits,(__ULong));
  393. int             _EXFUN(lo0bits,(__ULong *));
  394. _Bigint *       _EXFUN(d2b,(struct _reent *p, double d, int *e, int *bits));
  395. _Bigint *       _EXFUN(lshift,(struct _reent *p, _Bigint *b, int k));
  396. _Bigint *       _EXFUN(diff,(struct _reent *p, _Bigint *a, _Bigint *b));
  397. int             _EXFUN(cmp,(_Bigint *a, _Bigint *b));
  398. int             _EXFUN(gethex,(struct _reent *p, _CONST char **sp, struct FPI *fpi, Long *exp, _Bigint **bp, int sign));    
  399. double          _EXFUN(ratio,(_Bigint *a, _Bigint *b));
  400. __ULong         _EXFUN(any_on,(_Bigint *b, int k));
  401. void            _EXFUN(copybits,(__ULong *c, int n, _Bigint *b));
  402. void            _EXFUN(hexdig_init,(void));
  403. #ifdef INFNAN_CHECK
  404. int             _EXFUN(hexnan,(_CONST char **sp, struct FPI *fpi, __ULong *x0));
  405. #endif
  406.  
  407. #define Bcopy(x,y) memcpy((char *)&x->_sign, (char *)&y->_sign, y->_wds*sizeof(__Long) + 2*sizeof(int))
  408.  
  409. extern _CONST double tinytens[];
  410. extern _CONST double bigtens[];
  411. extern _CONST double tens[];
  412. extern unsigned char hexdig[];
  413.  
  414.  
  415. double _EXFUN(_mprec_log10,(int));
  416.