gwj / at91sam9263 · Commits

Commit 158142c2c2df728cfa3b5320c65534921a764f26

Authored by bellard
1 parent 4f716dc6

soft float support 


git-svn-id: svn://svn.savannah.nongnu.org/qemu/trunk@1332 c046a42c-6fe2-441c-8c8c-71466251a162
Inline Side-by-side
Showing 6 changed files with 1774 additions and 3 deletions

Too many changes to show.
Reload with full diff Plain diff Email patch

To preserve performance only 6 of 8 files are displayed.

Makefile.target
  View file @158142c
... ... @@ -229,7 +229,7 @@ ifeq ($(TARGET_ARCH), i386)
229 229 OBJS+= vm86.o
230 230 endif
231 231 ifeq ($(TARGET_ARCH), arm)
232   -OBJS+=nwfpe/softfloat.o nwfpe/fpa11.o nwfpe/fpa11_cpdo.o \
  232 +OBJS+=nwfpe/fpa11.o nwfpe/fpa11_cpdo.o \
233 233 nwfpe/fpa11_cpdt.o nwfpe/fpa11_cprt.o nwfpe/fpopcode.o nwfpe/single_cpdo.o \
234 234 nwfpe/double_cpdo.o nwfpe/extended_cpdo.o
235 235 endif
... ... @@ -237,8 +237,14 @@ SRCS:= $(OBJS:.o=.c)
237 237 OBJS+= libqemu.a
238 238  
239 239 # cpu emulator library
240   -LIBOBJS=exec.o kqemu.o translate-all.o cpu-exec.o\
  240 +LIBOBJS=exec.o kqemu.o translate-op.o translate-all.o cpu-exec.o\
241 241 translate.o op.o
  242 +ifdef CONFIG_SOFTFLOAT
  243 +LIBOBJS+=fpu/softfloat.o
  244 +else
  245 +LIBOBJS+=fpu/softfloat-native.o
  246 +endif
  247 +DEFINES+=-I$(SRC_PATH)/fpu
242 248  
243 249 ifeq ($(TARGET_ARCH), i386)
244 250 LIBOBJS+=helper.o helper2.o
... ... @@ -399,7 +405,9 @@ libqemu.a: $(LIBOBJS)
399 405  
400 406 translate.o: translate.c gen-op.h opc.h cpu.h
401 407  
402   -translate-all.o: translate-all.c op.h opc.h cpu.h
  408 +translate-all.o: translate-all.c opc.h cpu.h
  409 +
  410 +translate-op.o: translate-all.c op.h opc.h cpu.h
403 411  
404 412 op.h: op.o $(DYNGEN)
405 413 $(DYNGEN) -o $@ $<
... ...
configure
  View file @158142c
... ... @@ -593,6 +593,7 @@ fi
593 593 #echo "Creating $config_mak, $config_h and $target_dir/Makefile"
594 594  
595 595 mkdir -p $target_dir
  596 +mkdir -p $target_dir/fpu
596 597 if test "$target" = "arm-user" -o "$target" = "armeb-user" ; then
597 598 mkdir -p $target_dir/nwfpe
598 599 fi
... ... @@ -658,6 +659,10 @@ if test &quot;$target_user_only&quot; = &quot;yes&quot; ; then
658 659 echo "#define CONFIG_USER_ONLY 1" >> $config_h
659 660 fi
660 661  
  662 +if test "$target_cpu" = "arm" -o "$target_cpu" = "armeb" ; then
  663 + echo "CONFIG_SOFTFLOAT=yes" >> $config_mak
  664 + echo "#define CONFIG_SOFTFLOAT 1" >> $config_h
  665 +fi
661 666 # sdl defines
662 667  
663 668 if test "$target_user_only" = "no"; then
... ...
fpu/softfloat-macros.h 0 → 100644
  View file @158142c
  1 +
  2 +/*============================================================================
  3 +
  4 +This C source fragment is part of the SoftFloat IEC/IEEE Floating-point
  5 +Arithmetic Package, Release 2b.
  6 +
  7 +Written by John R. Hauser. This work was made possible in part by the
  8 +International Computer Science Institute, located at Suite 600, 1947 Center
  9 +Street, Berkeley, California 94704. Funding was partially provided by the
  10 +National Science Foundation under grant MIP-9311980. The original version
  11 +of this code was written as part of a project to build a fixed-point vector
  12 +processor in collaboration with the University of California at Berkeley,
  13 +overseen by Profs. Nelson Morgan and John Wawrzynek. More information
  14 +is available through the Web page `http://www.cs.berkeley.edu/~jhauser/
  15 +arithmetic/SoftFloat.html'.
  16 +
  17 +THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE. Although reasonable effort has
  18 +been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT TIMES
  19 +RESULT IN INCORRECT BEHAVIOR. USE OF THIS SOFTWARE IS RESTRICTED TO PERSONS
  20 +AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ALL LOSSES,
  21 +COSTS, OR OTHER PROBLEMS THEY INCUR DUE TO THE SOFTWARE, AND WHO FURTHERMORE
  22 +EFFECTIVELY INDEMNIFY JOHN HAUSER AND THE INTERNATIONAL COMPUTER SCIENCE
  23 +INSTITUTE (possibly via similar legal notice) AGAINST ALL LOSSES, COSTS, OR
  24 +OTHER PROBLEMS INCURRED BY THEIR CUSTOMERS AND CLIENTS DUE TO THE SOFTWARE.
  25 +
  26 +Derivative works are acceptable, even for commercial purposes, so long as
  27 +(1) the source code for the derivative work includes prominent notice that
  28 +the work is derivative, and (2) the source code includes prominent notice with
  29 +these four paragraphs for those parts of this code that are retained.
  30 +
  31 +=============================================================================*/
  32 +
  33 +/*----------------------------------------------------------------------------
  34 +| Shifts `a' right by the number of bits given in `count'. If any nonzero
  35 +| bits are shifted off, they are ``jammed'' into the least significant bit of
  36 +| the result by setting the least significant bit to 1. The value of `count'
  37 +| can be arbitrarily large; in particular, if `count' is greater than 32, the
  38 +| result will be either 0 or 1, depending on whether `a' is zero or nonzero.
  39 +| The result is stored in the location pointed to by `zPtr'.
  40 +*----------------------------------------------------------------------------*/
  41 +
  42 +INLINE void shift32RightJamming( bits32 a, int16 count, bits32 *zPtr )
  43 +{
  44 + bits32 z;
  45 +
  46 + if ( count == 0 ) {
  47 + z = a;
  48 + }
  49 + else if ( count < 32 ) {
  50 + z = ( a>>count ) | ( ( a<<( ( - count ) & 31 ) ) != 0 );
  51 + }
  52 + else {
  53 + z = ( a != 0 );
  54 + }
  55 + *zPtr = z;
  56 +
  57 +}
  58 +
  59 +/*----------------------------------------------------------------------------
  60 +| Shifts `a' right by the number of bits given in `count'. If any nonzero
  61 +| bits are shifted off, they are ``jammed'' into the least significant bit of
  62 +| the result by setting the least significant bit to 1. The value of `count'
  63 +| can be arbitrarily large; in particular, if `count' is greater than 64, the
  64 +| result will be either 0 or 1, depending on whether `a' is zero or nonzero.
  65 +| The result is stored in the location pointed to by `zPtr'.
  66 +*----------------------------------------------------------------------------*/
  67 +
  68 +INLINE void shift64RightJamming( bits64 a, int16 count, bits64 *zPtr )
  69 +{
  70 + bits64 z;
  71 +
  72 + if ( count == 0 ) {
  73 + z = a;
  74 + }
  75 + else if ( count < 64 ) {
  76 + z = ( a>>count ) | ( ( a<<( ( - count ) & 63 ) ) != 0 );
  77 + }
  78 + else {
  79 + z = ( a != 0 );
  80 + }
  81 + *zPtr = z;
  82 +
  83 +}
  84 +
  85 +/*----------------------------------------------------------------------------
  86 +| Shifts the 128-bit value formed by concatenating `a0' and `a1' right by 64
  87 +| _plus_ the number of bits given in `count'. The shifted result is at most
  88 +| 64 nonzero bits; this is stored at the location pointed to by `z0Ptr'. The
  89 +| bits shifted off form a second 64-bit result as follows: The _last_ bit
  90 +| shifted off is the most-significant bit of the extra result, and the other
  91 +| 63 bits of the extra result are all zero if and only if _all_but_the_last_
  92 +| bits shifted off were all zero. This extra result is stored in the location
  93 +| pointed to by `z1Ptr'. The value of `count' can be arbitrarily large.
  94 +| (This routine makes more sense if `a0' and `a1' are considered to form
  95 +| a fixed-point value with binary point between `a0' and `a1'. This fixed-
  96 +| point value is shifted right by the number of bits given in `count', and
  97 +| the integer part of the result is returned at the location pointed to by
  98 +| `z0Ptr'. The fractional part of the result may be slightly corrupted as
  99 +| described above, and is returned at the location pointed to by `z1Ptr'.)
  100 +*----------------------------------------------------------------------------*/
  101 +
  102 +INLINE void
  103 + shift64ExtraRightJamming(
  104 + bits64 a0, bits64 a1, int16 count, bits64 *z0Ptr, bits64 *z1Ptr )
  105 +{
  106 + bits64 z0, z1;
  107 + int8 negCount = ( - count ) & 63;
  108 +
  109 + if ( count == 0 ) {
  110 + z1 = a1;
  111 + z0 = a0;
  112 + }
  113 + else if ( count < 64 ) {
  114 + z1 = ( a0<<negCount ) | ( a1 != 0 );
  115 + z0 = a0>>count;
  116 + }
  117 + else {
  118 + if ( count == 64 ) {
  119 + z1 = a0 | ( a1 != 0 );
  120 + }
  121 + else {
  122 + z1 = ( ( a0 | a1 ) != 0 );
  123 + }
  124 + z0 = 0;
  125 + }
  126 + *z1Ptr = z1;
  127 + *z0Ptr = z0;
  128 +
  129 +}
  130 +
  131 +/*----------------------------------------------------------------------------
  132 +| Shifts the 128-bit value formed by concatenating `a0' and `a1' right by the
  133 +| number of bits given in `count'. Any bits shifted off are lost. The value
  134 +| of `count' can be arbitrarily large; in particular, if `count' is greater
  135 +| than 128, the result will be 0. The result is broken into two 64-bit pieces
  136 +| which are stored at the locations pointed to by `z0Ptr' and `z1Ptr'.
  137 +*----------------------------------------------------------------------------*/
  138 +
  139 +INLINE void
  140 + shift128Right(
  141 + bits64 a0, bits64 a1, int16 count, bits64 *z0Ptr, bits64 *z1Ptr )
  142 +{
  143 + bits64 z0, z1;
  144 + int8 negCount = ( - count ) & 63;
  145 +
  146 + if ( count == 0 ) {
  147 + z1 = a1;
  148 + z0 = a0;
  149 + }
  150 + else if ( count < 64 ) {
  151 + z1 = ( a0<<negCount ) | ( a1>>count );
  152 + z0 = a0>>count;
  153 + }
  154 + else {
  155 + z1 = ( count < 64 ) ? ( a0>>( count & 63 ) ) : 0;
  156 + z0 = 0;
  157 + }
  158 + *z1Ptr = z1;
  159 + *z0Ptr = z0;
  160 +
  161 +}
  162 +
  163 +/*----------------------------------------------------------------------------
  164 +| Shifts the 128-bit value formed by concatenating `a0' and `a1' right by the
  165 +| number of bits given in `count'. If any nonzero bits are shifted off, they
  166 +| are ``jammed'' into the least significant bit of the result by setting the
  167 +| least significant bit to 1. The value of `count' can be arbitrarily large;
  168 +| in particular, if `count' is greater than 128, the result will be either
  169 +| 0 or 1, depending on whether the concatenation of `a0' and `a1' is zero or
  170 +| nonzero. The result is broken into two 64-bit pieces which are stored at
  171 +| the locations pointed to by `z0Ptr' and `z1Ptr'.
  172 +*----------------------------------------------------------------------------*/
  173 +
  174 +INLINE void
  175 + shift128RightJamming(
  176 + bits64 a0, bits64 a1, int16 count, bits64 *z0Ptr, bits64 *z1Ptr )
  177 +{
  178 + bits64 z0, z1;
  179 + int8 negCount = ( - count ) & 63;
  180 +
  181 + if ( count == 0 ) {
  182 + z1 = a1;
  183 + z0 = a0;
  184 + }
  185 + else if ( count < 64 ) {
  186 + z1 = ( a0<<negCount ) | ( a1>>count ) | ( ( a1<<negCount ) != 0 );
  187 + z0 = a0>>count;
  188 + }
  189 + else {
  190 + if ( count == 64 ) {
  191 + z1 = a0 | ( a1 != 0 );
  192 + }
  193 + else if ( count < 128 ) {
  194 + z1 = ( a0>>( count & 63 ) ) | ( ( ( a0<<negCount ) | a1 ) != 0 );
  195 + }
  196 + else {
  197 + z1 = ( ( a0 | a1 ) != 0 );
  198 + }
  199 + z0 = 0;
  200 + }
  201 + *z1Ptr = z1;
  202 + *z0Ptr = z0;
  203 +
  204 +}
  205 +
  206 +/*----------------------------------------------------------------------------
  207 +| Shifts the 192-bit value formed by concatenating `a0', `a1', and `a2' right
  208 +| by 64 _plus_ the number of bits given in `count'. The shifted result is
  209 +| at most 128 nonzero bits; these are broken into two 64-bit pieces which are
  210 +| stored at the locations pointed to by `z0Ptr' and `z1Ptr'. The bits shifted
  211 +| off form a third 64-bit result as follows: The _last_ bit shifted off is
  212 +| the most-significant bit of the extra result, and the other 63 bits of the
  213 +| extra result are all zero if and only if _all_but_the_last_ bits shifted off
  214 +| were all zero. This extra result is stored in the location pointed to by
  215 +| `z2Ptr'. The value of `count' can be arbitrarily large.
  216 +| (This routine makes more sense if `a0', `a1', and `a2' are considered
  217 +| to form a fixed-point value with binary point between `a1' and `a2'. This
  218 +| fixed-point value is shifted right by the number of bits given in `count',
  219 +| and the integer part of the result is returned at the locations pointed to
  220 +| by `z0Ptr' and `z1Ptr'. The fractional part of the result may be slightly
  221 +| corrupted as described above, and is returned at the location pointed to by
  222 +| `z2Ptr'.)
  223 +*----------------------------------------------------------------------------*/
  224 +
  225 +INLINE void
  226 + shift128ExtraRightJamming(
  227 + bits64 a0,
  228 + bits64 a1,
  229 + bits64 a2,
  230 + int16 count,
  231 + bits64 *z0Ptr,
  232 + bits64 *z1Ptr,
  233 + bits64 *z2Ptr
  234 + )
  235 +{
  236 + bits64 z0, z1, z2;
  237 + int8 negCount = ( - count ) & 63;
  238 +
  239 + if ( count == 0 ) {
  240 + z2 = a2;
  241 + z1 = a1;
  242 + z0 = a0;
  243 + }
  244 + else {
  245 + if ( count < 64 ) {
  246 + z2 = a1<<negCount;
  247 + z1 = ( a0<<negCount ) | ( a1>>count );
  248 + z0 = a0>>count;
  249 + }
  250 + else {
  251 + if ( count == 64 ) {
  252 + z2 = a1;
  253 + z1 = a0;
  254 + }
  255 + else {
  256 + a2 |= a1;
  257 + if ( count < 128 ) {
  258 + z2 = a0<<negCount;
  259 + z1 = a0>>( count & 63 );
  260 + }
  261 + else {
  262 + z2 = ( count == 128 ) ? a0 : ( a0 != 0 );
  263 + z1 = 0;
  264 + }
  265 + }
  266 + z0 = 0;
  267 + }
  268 + z2 |= ( a2 != 0 );
  269 + }
  270 + *z2Ptr = z2;
  271 + *z1Ptr = z1;
  272 + *z0Ptr = z0;
  273 +
  274 +}
  275 +
  276 +/*----------------------------------------------------------------------------
  277 +| Shifts the 128-bit value formed by concatenating `a0' and `a1' left by the
  278 +| number of bits given in `count'. Any bits shifted off are lost. The value
  279 +| of `count' must be less than 64. The result is broken into two 64-bit
  280 +| pieces which are stored at the locations pointed to by `z0Ptr' and `z1Ptr'.
  281 +*----------------------------------------------------------------------------*/
  282 +
  283 +INLINE void
  284 + shortShift128Left(
  285 + bits64 a0, bits64 a1, int16 count, bits64 *z0Ptr, bits64 *z1Ptr )
  286 +{
  287 +
  288 + *z1Ptr = a1<<count;
  289 + *z0Ptr =
  290 + ( count == 0 ) ? a0 : ( a0<<count ) | ( a1>>( ( - count ) & 63 ) );
  291 +
  292 +}
  293 +
  294 +/*----------------------------------------------------------------------------
  295 +| Shifts the 192-bit value formed by concatenating `a0', `a1', and `a2' left
  296 +| by the number of bits given in `count'. Any bits shifted off are lost.
  297 +| The value of `count' must be less than 64. The result is broken into three
  298 +| 64-bit pieces which are stored at the locations pointed to by `z0Ptr',
  299 +| `z1Ptr', and `z2Ptr'.
  300 +*----------------------------------------------------------------------------*/
  301 +
  302 +INLINE void
  303 + shortShift192Left(
  304 + bits64 a0,
  305 + bits64 a1,
  306 + bits64 a2,
  307 + int16 count,
  308 + bits64 *z0Ptr,
  309 + bits64 *z1Ptr,
  310 + bits64 *z2Ptr
  311 + )
  312 +{
  313 + bits64 z0, z1, z2;
  314 + int8 negCount;
  315 +
  316 + z2 = a2<<count;
  317 + z1 = a1<<count;
  318 + z0 = a0<<count;
  319 + if ( 0 < count ) {
  320 + negCount = ( ( - count ) & 63 );
  321 + z1 |= a2>>negCount;
  322 + z0 |= a1>>negCount;
  323 + }
  324 + *z2Ptr = z2;
  325 + *z1Ptr = z1;
  326 + *z0Ptr = z0;
  327 +
  328 +}
  329 +
  330 +/*----------------------------------------------------------------------------
  331 +| Adds the 128-bit value formed by concatenating `a0' and `a1' to the 128-bit
  332 +| value formed by concatenating `b0' and `b1'. Addition is modulo 2^128, so
  333 +| any carry out is lost. The result is broken into two 64-bit pieces which
  334 +| are stored at the locations pointed to by `z0Ptr' and `z1Ptr'.
  335 +*----------------------------------------------------------------------------*/
  336 +
  337 +INLINE void
  338 + add128(
  339 + bits64 a0, bits64 a1, bits64 b0, bits64 b1, bits64 *z0Ptr, bits64 *z1Ptr )
  340 +{
  341 + bits64 z1;
  342 +
  343 + z1 = a1 + b1;
  344 + *z1Ptr = z1;
  345 + *z0Ptr = a0 + b0 + ( z1 < a1 );
  346 +
  347 +}
  348 +
  349 +/*----------------------------------------------------------------------------
  350 +| Adds the 192-bit value formed by concatenating `a0', `a1', and `a2' to the
  351 +| 192-bit value formed by concatenating `b0', `b1', and `b2'. Addition is
  352 +| modulo 2^192, so any carry out is lost. The result is broken into three
  353 +| 64-bit pieces which are stored at the locations pointed to by `z0Ptr',
  354 +| `z1Ptr', and `z2Ptr'.
  355 +*----------------------------------------------------------------------------*/
  356 +
  357 +INLINE void
  358 + add192(
  359 + bits64 a0,
  360 + bits64 a1,
  361 + bits64 a2,
  362 + bits64 b0,
  363 + bits64 b1,
  364 + bits64 b2,
  365 + bits64 *z0Ptr,
  366 + bits64 *z1Ptr,
  367 + bits64 *z2Ptr
  368 + )
  369 +{
  370 + bits64 z0, z1, z2;
  371 + int8 carry0, carry1;
  372 +
  373 + z2 = a2 + b2;
  374 + carry1 = ( z2 < a2 );
  375 + z1 = a1 + b1;
  376 + carry0 = ( z1 < a1 );
  377 + z0 = a0 + b0;
  378 + z1 += carry1;
  379 + z0 += ( z1 < carry1 );
  380 + z0 += carry0;
  381 + *z2Ptr = z2;
  382 + *z1Ptr = z1;
  383 + *z0Ptr = z0;
  384 +
  385 +}
  386 +
  387 +/*----------------------------------------------------------------------------
  388 +| Subtracts the 128-bit value formed by concatenating `b0' and `b1' from the
  389 +| 128-bit value formed by concatenating `a0' and `a1'. Subtraction is modulo
  390 +| 2^128, so any borrow out (carry out) is lost. The result is broken into two
  391 +| 64-bit pieces which are stored at the locations pointed to by `z0Ptr' and
  392 +| `z1Ptr'.
  393 +*----------------------------------------------------------------------------*/
  394 +
  395 +INLINE void
  396 + sub128(
  397 + bits64 a0, bits64 a1, bits64 b0, bits64 b1, bits64 *z0Ptr, bits64 *z1Ptr )
  398 +{
  399 +
  400 + *z1Ptr = a1 - b1;
  401 + *z0Ptr = a0 - b0 - ( a1 < b1 );
  402 +
  403 +}
  404 +
  405 +/*----------------------------------------------------------------------------
  406 +| Subtracts the 192-bit value formed by concatenating `b0', `b1', and `b2'
  407 +| from the 192-bit value formed by concatenating `a0', `a1', and `a2'.
  408 +| Subtraction is modulo 2^192, so any borrow out (carry out) is lost. The
  409 +| result is broken into three 64-bit pieces which are stored at the locations
  410 +| pointed to by `z0Ptr', `z1Ptr', and `z2Ptr'.
  411 +*----------------------------------------------------------------------------*/
  412 +
  413 +INLINE void
  414 + sub192(
  415 + bits64 a0,
  416 + bits64 a1,
  417 + bits64 a2,
  418 + bits64 b0,
  419 + bits64 b1,
  420 + bits64 b2,
  421 + bits64 *z0Ptr,
  422 + bits64 *z1Ptr,
  423 + bits64 *z2Ptr
  424 + )
  425 +{
  426 + bits64 z0, z1, z2;
  427 + int8 borrow0, borrow1;
  428 +
  429 + z2 = a2 - b2;
  430 + borrow1 = ( a2 < b2 );
  431 + z1 = a1 - b1;
  432 + borrow0 = ( a1 < b1 );
  433 + z0 = a0 - b0;
  434 + z0 -= ( z1 < borrow1 );
  435 + z1 -= borrow1;
  436 + z0 -= borrow0;
  437 + *z2Ptr = z2;
  438 + *z1Ptr = z1;
  439 + *z0Ptr = z0;
  440 +
  441 +}
  442 +
  443 +/*----------------------------------------------------------------------------
  444 +| Multiplies `a' by `b' to obtain a 128-bit product. The product is broken
  445 +| into two 64-bit pieces which are stored at the locations pointed to by
  446 +| `z0Ptr' and `z1Ptr'.
  447 +*----------------------------------------------------------------------------*/
  448 +
  449 +INLINE void mul64To128( bits64 a, bits64 b, bits64 *z0Ptr, bits64 *z1Ptr )
  450 +{
  451 + bits32 aHigh, aLow, bHigh, bLow;
  452 + bits64 z0, zMiddleA, zMiddleB, z1;
  453 +
  454 + aLow = a;
  455 + aHigh = a>>32;
  456 + bLow = b;
  457 + bHigh = b>>32;
  458 + z1 = ( (bits64) aLow ) * bLow;
  459 + zMiddleA = ( (bits64) aLow ) * bHigh;
  460 + zMiddleB = ( (bits64) aHigh ) * bLow;
  461 + z0 = ( (bits64) aHigh ) * bHigh;
  462 + zMiddleA += zMiddleB;
  463 + z0 += ( ( (bits64) ( zMiddleA < zMiddleB ) )<<32 ) + ( zMiddleA>>32 );
  464 + zMiddleA <<= 32;
  465 + z1 += zMiddleA;
  466 + z0 += ( z1 < zMiddleA );
  467 + *z1Ptr = z1;
  468 + *z0Ptr = z0;
  469 +
  470 +}
  471 +
  472 +/*----------------------------------------------------------------------------
  473 +| Multiplies the 128-bit value formed by concatenating `a0' and `a1' by
  474 +| `b' to obtain a 192-bit product. The product is broken into three 64-bit
  475 +| pieces which are stored at the locations pointed to by `z0Ptr', `z1Ptr', and
  476 +| `z2Ptr'.
  477 +*----------------------------------------------------------------------------*/
  478 +
  479 +INLINE void
  480 + mul128By64To192(
  481 + bits64 a0,
  482 + bits64 a1,
  483 + bits64 b,
  484 + bits64 *z0Ptr,
  485 + bits64 *z1Ptr,
  486 + bits64 *z2Ptr
  487 + )
  488 +{
  489 + bits64 z0, z1, z2, more1;
  490 +
  491 + mul64To128( a1, b, &z1, &z2 );
  492 + mul64To128( a0, b, &z0, &more1 );
  493 + add128( z0, more1, 0, z1, &z0, &z1 );
  494 + *z2Ptr = z2;
  495 + *z1Ptr = z1;
  496 + *z0Ptr = z0;
  497 +
  498 +}
  499 +
  500 +/*----------------------------------------------------------------------------
  501 +| Multiplies the 128-bit value formed by concatenating `a0' and `a1' to the
  502 +| 128-bit value formed by concatenating `b0' and `b1' to obtain a 256-bit
  503 +| product. The product is broken into four 64-bit pieces which are stored at
  504 +| the locations pointed to by `z0Ptr', `z1Ptr', `z2Ptr', and `z3Ptr'.
  505 +*----------------------------------------------------------------------------*/
  506 +
  507 +INLINE void
  508 + mul128To256(
  509 + bits64 a0,
  510 + bits64 a1,
  511 + bits64 b0,
  512 + bits64 b1,
  513 + bits64 *z0Ptr,
  514 + bits64 *z1Ptr,
  515 + bits64 *z2Ptr,
  516 + bits64 *z3Ptr
  517 + )
  518 +{
  519 + bits64 z0, z1, z2, z3;
  520 + bits64 more1, more2;
  521 +
  522 + mul64To128( a1, b1, &z2, &z3 );
  523 + mul64To128( a1, b0, &z1, &more2 );
  524 + add128( z1, more2, 0, z2, &z1, &z2 );
  525 + mul64To128( a0, b0, &z0, &more1 );
  526 + add128( z0, more1, 0, z1, &z0, &z1 );
  527 + mul64To128( a0, b1, &more1, &more2 );
  528 + add128( more1, more2, 0, z2, &more1, &z2 );
  529 + add128( z0, z1, 0, more1, &z0, &z1 );
  530 + *z3Ptr = z3;
  531 + *z2Ptr = z2;
  532 + *z1Ptr = z1;
  533 + *z0Ptr = z0;
  534 +
  535 +}
  536 +
  537 +/*----------------------------------------------------------------------------
  538 +| Returns an approximation to the 64-bit integer quotient obtained by dividing
  539 +| `b' into the 128-bit value formed by concatenating `a0' and `a1'. The
  540 +| divisor `b' must be at least 2^63. If q is the exact quotient truncated
  541 +| toward zero, the approximation returned lies between q and q + 2 inclusive.
  542 +| If the exact quotient q is larger than 64 bits, the maximum positive 64-bit
  543 +| unsigned integer is returned.
  544 +*----------------------------------------------------------------------------*/
  545 +
  546 +static bits64 estimateDiv128To64( bits64 a0, bits64 a1, bits64 b )
  547 +{
  548 + bits64 b0, b1;
  549 + bits64 rem0, rem1, term0, term1;
  550 + bits64 z;
  551 +
  552 + if ( b <= a0 ) return LIT64( 0xFFFFFFFFFFFFFFFF );
  553 + b0 = b>>32;
  554 + z = ( b0<<32 <= a0 ) ? LIT64( 0xFFFFFFFF00000000 ) : ( a0 / b0 )<<32;
  555 + mul64To128( b, z, &term0, &term1 );
  556 + sub128( a0, a1, term0, term1, &rem0, &rem1 );
  557 + while ( ( (sbits64) rem0 ) < 0 ) {
  558 + z -= LIT64( 0x100000000 );
  559 + b1 = b<<32;
  560 + add128( rem0, rem1, b0, b1, &rem0, &rem1 );
  561 + }
  562 + rem0 = ( rem0<<32 ) | ( rem1>>32 );
  563 + z |= ( b0<<32 <= rem0 ) ? 0xFFFFFFFF : rem0 / b0;
  564 + return z;
  565 +
  566 +}
  567 +
  568 +/*----------------------------------------------------------------------------
  569 +| Returns an approximation to the square root of the 32-bit significand given
  570 +| by `a'. Considered as an integer, `a' must be at least 2^31. If bit 0 of
  571 +| `aExp' (the least significant bit) is 1, the integer returned approximates
  572 +| 2^31*sqrt(`a'/2^31), where `a' is considered an integer. If bit 0 of `aExp'
  573 +| is 0, the integer returned approximates 2^31*sqrt(`a'/2^30). In either
  574 +| case, the approximation returned lies strictly within +/-2 of the exact
  575 +| value.
  576 +*----------------------------------------------------------------------------*/
  577 +
  578 +static bits32 estimateSqrt32( int16 aExp, bits32 a )
  579 +{
  580 + static const bits16 sqrtOddAdjustments[] = {
  581 + 0x0004, 0x0022, 0x005D, 0x00B1, 0x011D, 0x019F, 0x0236, 0x02E0,
  582 + 0x039C, 0x0468, 0x0545, 0x0631, 0x072B, 0x0832, 0x0946, 0x0A67
  583 + };
  584 + static const bits16 sqrtEvenAdjustments[] = {
  585 + 0x0A2D, 0x08AF, 0x075A, 0x0629, 0x051A, 0x0429, 0x0356, 0x029E,
  586 + 0x0200, 0x0179, 0x0109, 0x00AF, 0x0068, 0x0034, 0x0012, 0x0002
  587 + };
  588 + int8 index;
  589 + bits32 z;
  590 +
  591 + index = ( a>>27 ) & 15;
  592 + if ( aExp & 1 ) {
  593 + z = 0x4000 + ( a>>17 ) - sqrtOddAdjustments[ index ];
  594 + z = ( ( a / z )<<14 ) + ( z<<15 );
  595 + a >>= 1;
  596 + }
  597 + else {
  598 + z = 0x8000 + ( a>>17 ) - sqrtEvenAdjustments[ index ];
  599 + z = a / z + z;
  600 + z = ( 0x20000 <= z ) ? 0xFFFF8000 : ( z<<15 );
  601 + if ( z <= a ) return (bits32) ( ( (sbits32) a )>>1 );
  602 + }
  603 + return ( (bits32) ( ( ( (bits64) a )<<31 ) / z ) ) + ( z>>1 );
  604 +
  605 +}
  606 +
  607 +/*----------------------------------------------------------------------------
  608 +| Returns the number of leading 0 bits before the most-significant 1 bit of
  609 +| `a'. If `a' is zero, 32 is returned.
  610 +*----------------------------------------------------------------------------*/
  611 +
  612 +static int8 countLeadingZeros32( bits32 a )
  613 +{
  614 + static const int8 countLeadingZerosHigh[] = {
  615 + 8, 7, 6, 6, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4,
  616 + 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
  617 + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
  618 + 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
  619 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
  620 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
  621 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
  622 + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
  623 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  624 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  625 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  626 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  627 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  628 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  629 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  630 + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
  631 + };
  632 + int8 shiftCount;
  633 +
  634 + shiftCount = 0;
  635 + if ( a < 0x10000 ) {
  636 + shiftCount += 16;
  637 + a <<= 16;
  638 + }
  639 + if ( a < 0x1000000 ) {
  640 + shiftCount += 8;
  641 + a <<= 8;
  642 + }
  643 + shiftCount += countLeadingZerosHigh[ a>>24 ];
  644 + return shiftCount;
  645 +
  646 +}
  647 +
  648 +/*----------------------------------------------------------------------------
  649 +| Returns the number of leading 0 bits before the most-significant 1 bit of
  650 +| `a'. If `a' is zero, 64 is returned.
  651 +*----------------------------------------------------------------------------*/
  652 +
  653 +static int8 countLeadingZeros64( bits64 a )
  654 +{
  655 + int8 shiftCount;
  656 +
  657 + shiftCount = 0;
  658 + if ( a < ( (bits64) 1 )<<32 ) {
  659 + shiftCount += 32;
  660 + }
  661 + else {
  662 + a >>= 32;
  663 + }
  664 + shiftCount += countLeadingZeros32( a );
  665 + return shiftCount;
  666 +
  667 +}
  668 +
  669 +/*----------------------------------------------------------------------------
  670 +| Returns 1 if the 128-bit value formed by concatenating `a0' and `a1'
  671 +| is equal to the 128-bit value formed by concatenating `b0' and `b1'.
  672 +| Otherwise, returns 0.
  673 +*----------------------------------------------------------------------------*/
  674 +
  675 +INLINE flag eq128( bits64 a0, bits64 a1, bits64 b0, bits64 b1 )
  676 +{
  677 +
  678 + return ( a0 == b0 ) && ( a1 == b1 );
  679 +
  680 +}
  681 +
  682 +/*----------------------------------------------------------------------------
  683 +| Returns 1 if the 128-bit value formed by concatenating `a0' and `a1' is less
  684 +| than or equal to the 128-bit value formed by concatenating `b0' and `b1'.
  685 +| Otherwise, returns 0.
  686 +*----------------------------------------------------------------------------*/
  687 +
  688 +INLINE flag le128( bits64 a0, bits64 a1, bits64 b0, bits64 b1 )
  689 +{
  690 +
  691 + return ( a0 < b0 ) || ( ( a0 == b0 ) && ( a1 <= b1 ) );
  692 +
  693 +}
  694 +
  695 +/*----------------------------------------------------------------------------
  696 +| Returns 1 if the 128-bit value formed by concatenating `a0' and `a1' is less
  697 +| than the 128-bit value formed by concatenating `b0' and `b1'. Otherwise,
  698 +| returns 0.
  699 +*----------------------------------------------------------------------------*/
  700 +
  701 +INLINE flag lt128( bits64 a0, bits64 a1, bits64 b0, bits64 b1 )
  702 +{
  703 +
  704 + return ( a0 < b0 ) || ( ( a0 == b0 ) && ( a1 < b1 ) );
  705 +
  706 +}
  707 +
  708 +/*----------------------------------------------------------------------------
  709 +| Returns 1 if the 128-bit value formed by concatenating `a0' and `a1' is
  710 +| not equal to the 128-bit value formed by concatenating `b0' and `b1'.
  711 +| Otherwise, returns 0.
  712 +*----------------------------------------------------------------------------*/
  713 +
  714 +INLINE flag ne128( bits64 a0, bits64 a1, bits64 b0, bits64 b1 )
  715 +{
  716 +
  717 + return ( a0 != b0 ) || ( a1 != b1 );
  718 +
  719 +}
  720 +
... ...
fpu/softfloat-native.c 0 → 100644
  View file @158142c
  1 +/* Native implementation of soft float functions. Only a single status
  2 + context is supported */
  3 +#include "softfloat.h"
  4 +#include <math.h>
  5 +
  6 +void set_float_rounding_mode(int val STATUS_PARAM)
  7 +{
  8 + STATUS(float_rounding_mode) = val;
  9 +#if defined(_BSD) && !defined(__APPLE__)
  10 + fpsetround(val);
  11 +#elif defined(__arm__)
  12 + /* nothing to do */
  13 +#else
  14 + fesetround(val);
  15 +#endif
  16 +}
  17 +
  18 +#ifdef FLOATX80
  19 +void set_floatx80_rounding_precision(int val STATUS_PARAM)
  20 +{
  21 + STATUS(floatx80_rounding_precision) = val;
  22 +}
  23 +#endif
  24 +
  25 +#if defined(_BSD)
  26 +#define lrint(d) ((int32_t)rint(d))
  27 +#define llrint(d) ((int64_t)rint(d))
  28 +#endif
  29 +
  30 +#if defined(__powerpc__)
  31 +
  32 +/* correct (but slow) PowerPC rint() (glibc version is incorrect) */
  33 +double qemu_rint(double x)
  34 +{
  35 + double y = 4503599627370496.0;
  36 + if (fabs(x) >= y)
  37 + return x;
  38 + if (x < 0)
  39 + y = -y;
  40 + y = (x + y) - y;
  41 + if (y == 0.0)
  42 + y = copysign(y, x);
  43 + return y;
  44 +}
  45 +
  46 +#define rint qemu_rint
  47 +#endif
  48 +
  49 +/*----------------------------------------------------------------------------
  50 +| Software IEC/IEEE integer-to-floating-point conversion routines.
  51 +*----------------------------------------------------------------------------*/
  52 +float32 int32_to_float32(int v STATUS_PARAM)
  53 +{
  54 + return (float32)v;
  55 +}
  56 +
  57 +float64 int32_to_float64(int v STATUS_PARAM)
  58 +{
  59 + return (float64)v;
  60 +}
  61 +
  62 +#ifdef FLOATX80
  63 +floatx80 int32_to_floatx80(int v STATUS_PARAM)
  64 +{
  65 + return (floatx80)v;
  66 +}
  67 +#endif
  68 +float32 int64_to_float32( int64_t v STATUS_PARAM)
  69 +{
  70 + return (float32)v;
  71 +}
  72 +float64 int64_to_float64( int64_t v STATUS_PARAM)
  73 +{
  74 + return (float64)v;
  75 +}
  76 +#ifdef FLOATX80
  77 +floatx80 int64_to_floatx80( int64_t v STATUS_PARAM)
  78 +{
  79 + return (floatx80)v;
  80 +}
  81 +#endif
  82 +
  83 +/*----------------------------------------------------------------------------
  84 +| Software IEC/IEEE single-precision conversion routines.
  85 +*----------------------------------------------------------------------------*/
  86 +int float32_to_int32( float32 a STATUS_PARAM)
  87 +{
  88 + return lrintf(a);
  89 +}
  90 +int float32_to_int32_round_to_zero( float32 a STATUS_PARAM)
  91 +{
  92 + return (int)a;
  93 +}
  94 +int64_t float32_to_int64( float32 a STATUS_PARAM)
  95 +{
  96 + return llrintf(a);
  97 +}
  98 +
  99 +int64_t float32_to_int64_round_to_zero( float32 a STATUS_PARAM)
  100 +{
  101 + return (int64_t)a;
  102 +}
  103 +
  104 +float64 float32_to_float64( float32 a STATUS_PARAM)
  105 +{
  106 + return a;
  107 +}
  108 +#ifdef FLOATX80
  109 +floatx80 float32_to_floatx80( float32 a STATUS_PARAM)
  110 +{
  111 + return a;
  112 +}
  113 +#endif
  114 +
  115 +/*----------------------------------------------------------------------------
  116 +| Software IEC/IEEE single-precision operations.
  117 +*----------------------------------------------------------------------------*/
  118 +float32 float32_round_to_int( float32 a STATUS_PARAM)
  119 +{
  120 + return rintf(a);
  121 +}
  122 +
  123 +float32 float32_sqrt( float32 a STATUS_PARAM)
  124 +{
  125 + return sqrtf(a);
  126 +}
  127 +char float32_is_signaling_nan( float32 a1)
  128 +{
  129 + float32u u;
  130 + uint32_t a;
  131 + u.f = a1;
  132 + a = u.i;
  133 + return ( ( ( a>>22 ) & 0x1FF ) == 0x1FE ) && ( a & 0x003FFFFF );
  134 +}
  135 +
  136 +/*----------------------------------------------------------------------------
  137 +| Software IEC/IEEE double-precision conversion routines.
  138 +*----------------------------------------------------------------------------*/
  139 +int float64_to_int32( float64 a STATUS_PARAM)
  140 +{
  141 + return lrint(a);
  142 +}
  143 +int float64_to_int32_round_to_zero( float64 a STATUS_PARAM)
  144 +{
  145 + return (int)a;
  146 +}
  147 +int64_t float64_to_int64( float64 a STATUS_PARAM)
  148 +{
  149 + return llrint(a);
  150 +}
  151 +int64_t float64_to_int64_round_to_zero( float64 a STATUS_PARAM)
  152 +{
  153 + return (int64_t)a;
  154 +}
  155 +float32 float64_to_float32( float64 a STATUS_PARAM)
  156 +{
  157 + return a;
  158 +}
  159 +#ifdef FLOATX80
  160 +floatx80 float64_to_floatx80( float64 a STATUS_PARAM)
  161 +{
  162 + return a;
  163 +}
  164 +#endif
  165 +#ifdef FLOAT128
  166 +float128 float64_to_float128( float64 a STATUS_PARAM)
  167 +{
  168 + return a;
  169 +}
  170 +#endif
  171 +
  172 +/*----------------------------------------------------------------------------
  173 +| Software IEC/IEEE double-precision operations.
  174 +*----------------------------------------------------------------------------*/
  175 +float64 float64_round_to_int( float64 a STATUS_PARAM )
  176 +{
  177 +#if defined(__arm__)
  178 + switch(STATUS(float_rounding_mode)) {
  179 + default:
  180 + case float_round_nearest_even:
  181 + asm("rndd %0, %1" : "=f" (a) : "f"(a));
  182 + break;
  183 + case float_round_down:
  184 + asm("rnddm %0, %1" : "=f" (a) : "f"(a));
  185 + break;
  186 + case float_round_up:
  187 + asm("rnddp %0, %1" : "=f" (a) : "f"(a));
  188 + break;
  189 + case float_round_to_zero:
  190 + asm("rnddz %0, %1" : "=f" (a) : "f"(a));
  191 + break;
  192 + }
  193 +#else
  194 + return rint(a);
  195 +#endif
  196 +}
  197 +
  198 +float64 float64_sqrt( float64 a STATUS_PARAM)
  199 +{
  200 + return sqrt(a);
  201 +}
  202 +char float64_is_signaling_nan( float64 a1)
  203 +{
  204 + float64u u;
  205 + uint64_t a;
  206 + u.f = a1;
  207 + a = u.i;
  208 + return
  209 + ( ( ( a>>51 ) & 0xFFF ) == 0xFFE )
  210 + && ( a & LIT64( 0x0007FFFFFFFFFFFF ) );
  211 +
  212 +}
  213 +
  214 +#ifdef FLOATX80
  215 +
  216 +/*----------------------------------------------------------------------------
  217 +| Software IEC/IEEE extended double-precision conversion routines.
  218 +*----------------------------------------------------------------------------*/
  219 +int floatx80_to_int32( floatx80 a STATUS_PARAM)
  220 +{
  221 + return lrintl(a);
  222 +}
  223 +int floatx80_to_int32_round_to_zero( floatx80 a STATUS_PARAM)
  224 +{
  225 + return (int)a;
  226 +}
  227 +int64_t floatx80_to_int64( floatx80 a STATUS_PARAM)
  228 +{
  229 + return llrintl(a);
  230 +}
  231 +int64_t floatx80_to_int64_round_to_zero( floatx80 a STATUS_PARAM)
  232 +{
  233 + return (int64_t)a;
  234 +}
  235 +float32 floatx80_to_float32( floatx80 a STATUS_PARAM)
  236 +{
  237 + return a;
  238 +}
  239 +float64 floatx80_to_float64( floatx80 a STATUS_PARAM)
  240 +{
  241 + return a;
  242 +}
  243 +
  244 +/*----------------------------------------------------------------------------
  245 +| Software IEC/IEEE extended double-precision operations.
  246 +*----------------------------------------------------------------------------*/
  247 +floatx80 floatx80_round_to_int( floatx80 a STATUS_PARAM)
  248 +{
  249 + return rintl(a);
  250 +}
  251 +floatx80 floatx80_sqrt( floatx80 a STATUS_PARAM)
  252 +{
  253 + return sqrtl(a);
  254 +}
  255 +char floatx80_is_signaling_nan( floatx80 a1)
  256 +{
  257 + floatx80u u;
  258 + u.f = a1;
  259 + return ( ( u.i.high & 0x7FFF ) == 0x7FFF ) && (bits64) ( u.i.low<<1 );
  260 +}
  261 +
  262 +#endif
... ...
fpu/softfloat-native.h 0 → 100644
  View file @158142c
  1 +/* Native implementation of soft float functions */
  2 +#include <math.h>
  3 +#if defined(_BSD) && !defined(__APPLE__)
  4 +#include <ieeefp.h>
  5 +#else
  6 +#include <fenv.h>
  7 +#endif
  8 +
  9 +typedef float float32;
  10 +typedef double float64;
  11 +#ifdef FLOATX80
  12 +typedef long double floatx80;
  13 +#endif
  14 +
  15 +typedef union {
  16 + float32 f;
  17 + uint32_t i;
  18 +} float32u;
  19 +typedef union {
  20 + float64 f;
  21 + uint64_t i;
  22 +} float64u;
  23 +#ifdef FLOATX80
  24 +typedef union {
  25 + floatx80 f;
  26 + struct {
  27 + uint64_t low;
  28 + uint16_t high;
  29 + } i;
  30 +} floatx80u;
  31 +#endif
  32 +
  33 +/*----------------------------------------------------------------------------
  34 +| Software IEC/IEEE floating-point rounding mode.
  35 +*----------------------------------------------------------------------------*/
  36 +#if defined(_BSD) && !defined(__APPLE__)
  37 +enum {
  38 + float_round_nearest_even = FP_RN,
  39 + float_round_down = FE_RM,
  40 + float_round_up = FE_RP,
  41 + float_round_to_zero = FE_RZ
  42 +};
  43 +#elif defined(__arm__)
  44 +enum {
  45 + float_round_nearest_even = 0,
  46 + float_round_down = 1,
  47 + float_round_up = 2,
  48 + float_round_to_zero = 3
  49 +};
  50 +#else
  51 +enum {
  52 + float_round_nearest_even = FE_TONEAREST,
  53 + float_round_down = FE_DOWNWARD,
  54 + float_round_up = FE_UPWARD,
  55 + float_round_to_zero = FE_TOWARDZERO
  56 +};
  57 +#endif
  58 +
  59 +typedef struct float_status {
  60 + signed char float_rounding_mode;
  61 +#ifdef FLOATX80
  62 + signed char floatx80_rounding_precision;
  63 +#endif
  64 +} float_status;
  65 +
  66 +void set_float_rounding_mode(int val STATUS_PARAM);
  67 +#ifdef FLOATX80
  68 +void set_floatx80_rounding_precision(int val STATUS_PARAM);
  69 +#endif
  70 +
  71 +/*----------------------------------------------------------------------------
  72 +| Software IEC/IEEE integer-to-floating-point conversion routines.
  73 +*----------------------------------------------------------------------------*/
  74 +float32 int32_to_float32( int STATUS_PARAM);
  75 +float64 int32_to_float64( int STATUS_PARAM);
  76 +#ifdef FLOATX80
  77 +floatx80 int32_to_floatx80( int STATUS_PARAM);
  78 +#endif
  79 +#ifdef FLOAT128
  80 +float128 int32_to_float128( int STATUS_PARAM);
  81 +#endif
  82 +float32 int64_to_float32( int64_t STATUS_PARAM);
  83 +float64 int64_to_float64( int64_t STATUS_PARAM);
  84 +#ifdef FLOATX80
  85 +floatx80 int64_to_floatx80( int64_t STATUS_PARAM);
  86 +#endif
  87 +#ifdef FLOAT128
  88 +float128 int64_to_float128( int64_t STATUS_PARAM);
  89 +#endif
  90 +
  91 +/*----------------------------------------------------------------------------
  92 +| Software IEC/IEEE single-precision conversion routines.
  93 +*----------------------------------------------------------------------------*/
  94 +int float32_to_int32( float32 STATUS_PARAM);
  95 +int float32_to_int32_round_to_zero( float32 STATUS_PARAM);
  96 +int64_t float32_to_int64( float32 STATUS_PARAM);
  97 +int64_t float32_to_int64_round_to_zero( float32 STATUS_PARAM);
  98 +float64 float32_to_float64( float32 STATUS_PARAM);
  99 +#ifdef FLOATX80
  100 +floatx80 float32_to_floatx80( float32 STATUS_PARAM);
  101 +#endif
  102 +#ifdef FLOAT128
  103 +float128 float32_to_float128( float32 STATUS_PARAM);
  104 +#endif
  105 +
  106 +/*----------------------------------------------------------------------------
  107 +| Software IEC/IEEE single-precision operations.
  108 +*----------------------------------------------------------------------------*/
  109 +float32 float32_round_to_int( float32 STATUS_PARAM);
  110 +INLINE float32 float32_add( float32 a, float32 b STATUS_PARAM)
  111 +{
  112 + return a + b;
  113 +}
  114 +INLINE float32 float32_sub( float32 a, float32 b STATUS_PARAM)
  115 +{
  116 + return a - b;
  117 +}
  118 +INLINE float32 float32_mul( float32 a, float32 b STATUS_PARAM)
  119 +{
  120 + return a * b;
  121 +}
  122 +INLINE float32 float32_div( float32 a, float32 b STATUS_PARAM)
  123 +{
  124 + return a / b;
  125 +}
  126 +float32 float32_rem( float32, float32 STATUS_PARAM);
  127 +float32 float32_sqrt( float32 STATUS_PARAM);
  128 +INLINE char float32_eq( float32 a, float32 b STATUS_PARAM)
  129 +{
  130 + /* XXX: incorrect because it can raise an exception */
  131 + return a == b;
  132 +}
  133 +INLINE char float32_le( float32 a, float32 b STATUS_PARAM)
  134 +{
  135 + return a <= b;
  136 +}
  137 +INLINE char float32_lt( float32 a, float32 b STATUS_PARAM)
  138 +{
  139 + return a < b;
  140 +}
  141 +INLINE char float32_eq_signaling( float32 a, float32 b STATUS_PARAM)
  142 +{
  143 + return a == b;
  144 +}
  145 +INLINE char float32_le_quiet( float32 a, float32 b STATUS_PARAM)
  146 +{
  147 + return islessequal(a, b);
  148 +}
  149 +INLINE char float32_lt_quiet( float32 a, float32 b STATUS_PARAM)
  150 +{
  151 + return isless(a, b);
  152 +}
  153 +char float32_is_signaling_nan( float32 );
  154 +
  155 +INLINE float32 float32_abs(float32 a)
  156 +{
  157 + return fabsf(a);
  158 +}
  159 +
  160 +INLINE float32 float32_chs(float32 a)
  161 +{
  162 + return -a;
  163 +}
  164 +
  165 +/*----------------------------------------------------------------------------
  166 +| Software IEC/IEEE double-precision conversion routines.
  167 +*----------------------------------------------------------------------------*/
  168 +int float64_to_int32( float64 STATUS_PARAM );
  169 +int float64_to_int32_round_to_zero( float64 STATUS_PARAM );
  170 +int64_t float64_to_int64( float64 STATUS_PARAM );
  171 +int64_t float64_to_int64_round_to_zero( float64 STATUS_PARAM );
  172 +float32 float64_to_float32( float64 STATUS_PARAM );
  173 +#ifdef FLOATX80
  174 +floatx80 float64_to_floatx80( float64 STATUS_PARAM );
  175 +#endif
  176 +#ifdef FLOAT128
  177 +float128 float64_to_float128( float64 STATUS_PARAM );
  178 +#endif
  179 +
  180 +/*----------------------------------------------------------------------------
  181 +| Software IEC/IEEE double-precision operations.
  182 +*----------------------------------------------------------------------------*/
  183 +float64 float64_round_to_int( float64 STATUS_PARAM );
  184 +INLINE float64 float64_add( float64 a, float64 b STATUS_PARAM)
  185 +{
  186 + return a + b;
  187 +}
  188 +INLINE float64 float64_sub( float64 a, float64 b STATUS_PARAM)
  189 +{
  190 + return a - b;
  191 +}
  192 +INLINE float64 float64_mul( float64 a, float64 b STATUS_PARAM)
  193 +{
  194 + return a * b;
  195 +}
  196 +INLINE float64 float64_div( float64 a, float64 b STATUS_PARAM)
  197 +{
  198 + return a / b;
  199 +}
  200 +float64 float64_rem( float64, float64 STATUS_PARAM );
  201 +float64 float64_sqrt( float64 STATUS_PARAM );
  202 +INLINE char float64_eq( float64 a, float64 b STATUS_PARAM)
  203 +{
  204 + return a == b;
  205 +}
  206 +INLINE char float64_le( float64 a, float64 b STATUS_PARAM)
  207 +{
  208 + return a <= b;
  209 +}
  210 +INLINE char float64_lt( float64 a, float64 b STATUS_PARAM)
  211 +{
  212 + return a < b;
  213 +}
  214 +INLINE char float64_eq_signaling( float64 a, float64 b STATUS_PARAM)
  215 +{
  216 + return a == b;
  217 +}
  218 +INLINE char float64_le_quiet( float64 a, float64 b STATUS_PARAM)
  219 +{
  220 + return islessequal(a, b);
  221 +}
  222 +INLINE char float64_lt_quiet( float64 a, float64 b STATUS_PARAM)
  223 +{
  224 + return isless(a, b);
  225 +
  226 +}
  227 +char float64_is_signaling_nan( float64 );
  228 +
  229 +INLINE float64 float64_abs(float64 a)
  230 +{
  231 + return fabs(a);
  232 +}
  233 +
  234 +INLINE float64 float64_chs(float64 a)
  235 +{
  236 + return -a;
  237 +}
  238 +
  239 +#ifdef FLOATX80
  240 +
  241 +/*----------------------------------------------------------------------------
  242 +| Software IEC/IEEE extended double-precision conversion routines.
  243 +*----------------------------------------------------------------------------*/
  244 +int floatx80_to_int32( floatx80 STATUS_PARAM );
  245 +int floatx80_to_int32_round_to_zero( floatx80 STATUS_PARAM );
  246 +int64_t floatx80_to_int64( floatx80 STATUS_PARAM);
  247 +int64_t floatx80_to_int64_round_to_zero( floatx80 STATUS_PARAM);
  248 +float32 floatx80_to_float32( floatx80 STATUS_PARAM );
  249 +float64 floatx80_to_float64( floatx80 STATUS_PARAM );
  250 +#ifdef FLOAT128
  251 +float128 floatx80_to_float128( floatx80 STATUS_PARAM );
  252 +#endif
  253 +
  254 +/*----------------------------------------------------------------------------
  255 +| Software IEC/IEEE extended double-precision operations.
  256 +*----------------------------------------------------------------------------*/
  257 +floatx80 floatx80_round_to_int( floatx80 STATUS_PARAM );
  258 +INLINE floatx80 floatx80_add( floatx80 a, floatx80 b STATUS_PARAM)
  259 +{
  260 + return a + b;
  261 +}
  262 +INLINE floatx80 floatx80_sub( floatx80 a, floatx80 b STATUS_PARAM)
  263 +{
  264 + return a - b;
  265 +}
  266 +INLINE floatx80 floatx80_mul( floatx80 a, floatx80 b STATUS_PARAM)
  267 +{
  268 + return a * b;
  269 +}
  270 +INLINE floatx80 floatx80_div( floatx80 a, floatx80 b STATUS_PARAM)
  271 +{
  272 + return a / b;
  273 +}
  274 +floatx80 floatx80_rem( floatx80, floatx80 STATUS_PARAM );
  275 +floatx80 floatx80_sqrt( floatx80 STATUS_PARAM );
  276 +INLINE char floatx80_eq( floatx80 a, floatx80 b STATUS_PARAM)
  277 +{
  278 + return a == b;
  279 +}
  280 +INLINE char floatx80_le( floatx80 a, floatx80 b STATUS_PARAM)
  281 +{
  282 + return a <= b;
  283 +}
  284 +INLINE char floatx80_lt( floatx80 a, floatx80 b STATUS_PARAM)
  285 +{
  286 + return a < b;
  287 +}
  288 +INLINE char floatx80_eq_signaling( floatx80 a, floatx80 b STATUS_PARAM)
  289 +{
  290 + return a == b;
  291 +}
  292 +INLINE char floatx80_le_quiet( floatx80 a, floatx80 b STATUS_PARAM)
  293 +{
  294 + return islessequal(a, b);
  295 +}
  296 +INLINE char floatx80_lt_quiet( floatx80 a, floatx80 b STATUS_PARAM)
  297 +{
  298 + return isless(a, b);
  299 +
  300 +}
  301 +char floatx80_is_signaling_nan( floatx80 );
  302 +
  303 +INLINE floatx80 floatx80_abs(floatx80 a)
  304 +{
  305 + return fabsl(a);
  306 +}
  307 +
  308 +INLINE floatx80 floatx80_chs(floatx80 a)
  309 +{
  310 + return -a;
  311 +}
  312 +#endif
... ...
fpu/softfloat-specialize.h 0 → 100644
  View file @158142c
  1 +
  2 +/*============================================================================
  3 +
  4 +This C source fragment is part of the SoftFloat IEC/IEEE Floating-point
  5 +Arithmetic Package, Release 2b.
  6 +
  7 +Written by John R. Hauser. This work was made possible in part by the
  8 +International Computer Science Institute, located at Suite 600, 1947 Center
  9 +Street, Berkeley, California 94704. Funding was partially provided by the
  10 +National Science Foundation under grant MIP-9311980. The original version
  11 +of this code was written as part of a project to build a fixed-point vector
  12 +processor in collaboration with the University of California at Berkeley,
  13 +overseen by Profs. Nelson Morgan and John Wawrzynek. More information
  14 +is available through the Web page `http://www.cs.berkeley.edu/~jhauser/
  15 +arithmetic/SoftFloat.html'.
  16 +
  17 +THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE. Although reasonable effort has
  18 +been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT TIMES
  19 +RESULT IN INCORRECT BEHAVIOR. USE OF THIS SOFTWARE IS RESTRICTED TO PERSONS
  20 +AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ALL LOSSES,
  21 +COSTS, OR OTHER PROBLEMS THEY INCUR DUE TO THE SOFTWARE, AND WHO FURTHERMORE
  22 +EFFECTIVELY INDEMNIFY JOHN HAUSER AND THE INTERNATIONAL COMPUTER SCIENCE
  23 +INSTITUTE (possibly via similar legal warning) AGAINST ALL LOSSES, COSTS, OR
  24 +OTHER PROBLEMS INCURRED BY THEIR CUSTOMERS AND CLIENTS DUE TO THE SOFTWARE.
  25 +
  26 +Derivative works are acceptable, even for commercial purposes, so long as
  27 +(1) the source code for the derivative work includes prominent notice that
  28 +the work is derivative, and (2) the source code includes prominent notice with
  29 +these four paragraphs for those parts of this code that are retained.
  30 +
  31 +=============================================================================*/
  32 +
  33 +/*----------------------------------------------------------------------------
  34 +| Underflow tininess-detection mode, statically initialized to default value.
  35 +| (The declaration in `softfloat.h' must match the `int8' type here.)
  36 +*----------------------------------------------------------------------------*/
  37 +int8 float_detect_tininess = float_tininess_after_rounding;
  38 +
  39 +/*----------------------------------------------------------------------------
  40 +| Raises the exceptions specified by `flags'. Floating-point traps can be
  41 +| defined here if desired. It is currently not possible for such a trap
  42 +| to substitute a result value. If traps are not implemented, this routine
  43 +| should be simply `float_exception_flags |= flags;'.
  44 +*----------------------------------------------------------------------------*/
  45 +
  46 +void float_raise( int8 flags STATUS_PARAM )
  47 +{
  48 +
  49 + STATUS(float_exception_flags) |= flags;
  50 +
  51 +}
  52 +
  53 +/*----------------------------------------------------------------------------
  54 +| Internal canonical NaN format.
  55 +*----------------------------------------------------------------------------*/
  56 +typedef struct {
  57 + flag sign;
  58 + bits64 high, low;
  59 +} commonNaNT;
  60 +
  61 +/*----------------------------------------------------------------------------
  62 +| The pattern for a default generated single-precision NaN.
  63 +*----------------------------------------------------------------------------*/
  64 +#define float32_default_nan 0xFFC00000
  65 +
  66 +/*----------------------------------------------------------------------------
  67 +| Returns 1 if the single-precision floating-point value `a' is a NaN;
  68 +| otherwise returns 0.
  69 +*----------------------------------------------------------------------------*/
  70 +
  71 +flag float32_is_nan( float32 a )
  72 +{
  73 +
  74 + return ( 0xFF000000 < (bits32) ( a<<1 ) );
  75 +
  76 +}
  77 +
  78 +/*----------------------------------------------------------------------------
  79 +| Returns 1 if the single-precision floating-point value `a' is a signaling
  80 +| NaN; otherwise returns 0.
  81 +*----------------------------------------------------------------------------*/
  82 +
  83 +flag float32_is_signaling_nan( float32 a )
  84 +{
  85 +
  86 + return ( ( ( a>>22 ) & 0x1FF ) == 0x1FE ) && ( a & 0x003FFFFF );
  87 +
  88 +}
  89 +
  90 +/*----------------------------------------------------------------------------
  91 +| Returns the result of converting the single-precision floating-point NaN
  92 +| `a' to the canonical NaN format. If `a' is a signaling NaN, the invalid
  93 +| exception is raised.
  94 +*----------------------------------------------------------------------------*/
  95 +
  96 +static commonNaNT float32ToCommonNaN( float32 a STATUS_PARAM )
  97 +{
  98 + commonNaNT z;
  99 +
  100 + if ( float32_is_signaling_nan( a ) ) float_raise( float_flag_invalid STATUS_VAR );
  101 + z.sign = a>>31;
  102 + z.low = 0;
  103 + z.high = ( (bits64) a )<<41;
  104 + return z;
  105 +
  106 +}
  107 +
  108 +/*----------------------------------------------------------------------------
  109 +| Returns the result of converting the canonical NaN `a' to the single-
  110 +| precision floating-point format.
  111 +*----------------------------------------------------------------------------*/
  112 +
  113 +static float32 commonNaNToFloat32( commonNaNT a )
  114 +{
  115 +
  116 + return ( ( (bits32) a.sign )<<31 ) | 0x7FC00000 | ( a.high>>41 );
  117 +
  118 +}
  119 +
  120 +/*----------------------------------------------------------------------------
  121 +| Takes two single-precision floating-point values `a' and `b', one of which
  122 +| is a NaN, and returns the appropriate NaN result. If either `a' or `b' is a
  123 +| signaling NaN, the invalid exception is raised.
  124 +*----------------------------------------------------------------------------*/
  125 +
  126 +static float32 propagateFloat32NaN( float32 a, float32 b STATUS_PARAM)
  127 +{
  128 + flag aIsNaN, aIsSignalingNaN, bIsNaN, bIsSignalingNaN;
  129 +
  130 + aIsNaN = float32_is_nan( a );
  131 + aIsSignalingNaN = float32_is_signaling_nan( a );
  132 + bIsNaN = float32_is_nan( b );
  133 + bIsSignalingNaN = float32_is_signaling_nan( b );
  134 + a |= 0x00400000;
  135 + b |= 0x00400000;
  136 + if ( aIsSignalingNaN | bIsSignalingNaN ) float_raise( float_flag_invalid STATUS_VAR);
  137 + if ( aIsSignalingNaN ) {
  138 + if ( bIsSignalingNaN ) goto returnLargerSignificand;
  139 + return bIsNaN ? b : a;
  140 + }
  141 + else if ( aIsNaN ) {
  142 + if ( bIsSignalingNaN | ! bIsNaN ) return a;
  143 + returnLargerSignificand:
  144 + if ( (bits32) ( a<<1 ) < (bits32) ( b<<1 ) ) return b;
  145 + if ( (bits32) ( b<<1 ) < (bits32) ( a<<1 ) ) return a;
  146 + return ( a < b ) ? a : b;
  147 + }
  148 + else {
  149 + return b;
  150 + }
  151 +
  152 +}
  153 +
  154 +/*----------------------------------------------------------------------------
  155 +| The pattern for a default generated double-precision NaN.
  156 +*----------------------------------------------------------------------------*/
  157 +#define float64_default_nan LIT64( 0xFFF8000000000000 )
  158 +
  159 +/*----------------------------------------------------------------------------
  160 +| Returns 1 if the double-precision floating-point value `a' is a NaN;
  161 +| otherwise returns 0.
  162 +*----------------------------------------------------------------------------*/
  163 +
  164 +flag float64_is_nan( float64 a )
  165 +{
  166 +
  167 + return ( LIT64( 0xFFE0000000000000 ) < (bits64) ( a<<1 ) );
  168 +
  169 +}
  170 +
  171 +/*----------------------------------------------------------------------------
  172 +| Returns 1 if the double-precision floating-point value `a' is a signaling
  173 +| NaN; otherwise returns 0.
  174 +*----------------------------------------------------------------------------*/
  175 +
  176 +flag float64_is_signaling_nan( float64 a )
  177 +{
  178 +
  179 + return
  180 + ( ( ( a>>51 ) & 0xFFF ) == 0xFFE )
  181 + && ( a & LIT64( 0x0007FFFFFFFFFFFF ) );
  182 +
  183 +}
  184 +
  185 +/*----------------------------------------------------------------------------
  186 +| Returns the result of converting the double-precision floating-point NaN
  187 +| `a' to the canonical NaN format. If `a' is a signaling NaN, the invalid
  188 +| exception is raised.
  189 +*----------------------------------------------------------------------------*/
  190 +
  191 +static commonNaNT float64ToCommonNaN( float64 a STATUS_PARAM)
  192 +{
  193 + commonNaNT z;
  194 +
  195 + if ( float64_is_signaling_nan( a ) ) float_raise( float_flag_invalid STATUS_VAR);
  196 + z.sign = a>>63;
  197 + z.low = 0;
  198 + z.high = a<<12;
  199 + return z;
  200 +
  201 +}
  202 +
  203 +/*----------------------------------------------------------------------------
  204 +| Returns the result of converting the canonical NaN `a' to the double-
  205 +| precision floating-point format.
  206 +*----------------------------------------------------------------------------*/
  207 +
  208 +static float64 commonNaNToFloat64( commonNaNT a )
  209 +{
  210 +
  211 + return
  212 + ( ( (bits64) a.sign )<<63 )
  213 + | LIT64( 0x7FF8000000000000 )
  214 + | ( a.high>>12 );
  215 +
  216 +}
  217 +
  218 +/*----------------------------------------------------------------------------
  219 +| Takes two double-precision floating-point values `a' and `b', one of which
  220 +| is a NaN, and returns the appropriate NaN result. If either `a' or `b' is a
  221 +| signaling NaN, the invalid exception is raised.
  222 +*----------------------------------------------------------------------------*/
  223 +
  224 +static float64 propagateFloat64NaN( float64 a, float64 b STATUS_PARAM)
  225 +{
  226 + flag aIsNaN, aIsSignalingNaN, bIsNaN, bIsSignalingNaN;
  227 +
  228 + aIsNaN = float64_is_nan( a );
  229 + aIsSignalingNaN = float64_is_signaling_nan( a );
  230 + bIsNaN = float64_is_nan( b );
  231 + bIsSignalingNaN = float64_is_signaling_nan( b );
  232 + a |= LIT64( 0x0008000000000000 );
  233 + b |= LIT64( 0x0008000000000000 );
  234 + if ( aIsSignalingNaN | bIsSignalingNaN ) float_raise( float_flag_invalid STATUS_VAR);
  235 + if ( aIsSignalingNaN ) {
  236 + if ( bIsSignalingNaN ) goto returnLargerSignificand;
  237 + return bIsNaN ? b : a;
  238 + }
  239 + else if ( aIsNaN ) {
  240 + if ( bIsSignalingNaN | ! bIsNaN ) return a;
  241 + returnLargerSignificand:
  242 + if ( (bits64) ( a<<1 ) < (bits64) ( b<<1 ) ) return b;
  243 + if ( (bits64) ( b<<1 ) < (bits64) ( a<<1 ) ) return a;
  244 + return ( a < b ) ? a : b;
  245 + }
  246 + else {
  247 + return b;
  248 + }
  249 +
  250 +}
  251 +
  252 +#ifdef FLOATX80
  253 +
  254 +/*----------------------------------------------------------------------------
  255 +| The pattern for a default generated extended double-precision NaN. The
  256 +| `high' and `low' values hold the most- and least-significant bits,
  257 +| respectively.
  258 +*----------------------------------------------------------------------------*/
  259 +#define floatx80_default_nan_high 0xFFFF
  260 +#define floatx80_default_nan_low LIT64( 0xC000000000000000 )
  261 +
  262 +/*----------------------------------------------------------------------------
  263 +| Returns 1 if the extended double-precision floating-point value `a' is a
  264 +| NaN; otherwise returns 0.
  265 +*----------------------------------------------------------------------------*/
  266 +
  267 +flag floatx80_is_nan( floatx80 a )
  268 +{
  269 +
  270 + return ( ( a.high & 0x7FFF ) == 0x7FFF ) && (bits64) ( a.low<<1 );
  271 +
  272 +}
  273 +
  274 +/*----------------------------------------------------------------------------
  275 +| Returns 1 if the extended double-precision floating-point value `a' is a
  276 +| signaling NaN; otherwise returns 0.
  277 +*----------------------------------------------------------------------------*/
  278 +
  279 +flag floatx80_is_signaling_nan( floatx80 a )
  280 +{
  281 + bits64 aLow;
  282 +
  283 + aLow = a.low & ~ LIT64( 0x4000000000000000 );
  284 + return
  285 + ( ( a.high & 0x7FFF ) == 0x7FFF )
  286 + && (bits64) ( aLow<<1 )
  287 + && ( a.low == aLow );
  288 +
  289 +}
  290 +
  291 +/*----------------------------------------------------------------------------
  292 +| Returns the result of converting the extended double-precision floating-
  293 +| point NaN `a' to the canonical NaN format. If `a' is a signaling NaN, the
  294 +| invalid exception is raised.
  295 +*----------------------------------------------------------------------------*/
  296 +
  297 +static commonNaNT floatx80ToCommonNaN( floatx80 a STATUS_PARAM)
  298 +{
  299 + commonNaNT z;
  300 +
  301 + if ( floatx80_is_signaling_nan( a ) ) float_raise( float_flag_invalid STATUS_VAR);
  302 + z.sign = a.high>>15;
  303 + z.low = 0;
  304 + z.high = a.low<<1;
  305 + return z;
  306 +
  307 +}
  308 +
  309 +/*----------------------------------------------------------------------------
  310 +| Returns the result of converting the canonical NaN `a' to the extended
  311 +| double-precision floating-point format.
  312 +*----------------------------------------------------------------------------*/
  313 +
  314 +static floatx80 commonNaNToFloatx80( commonNaNT a )
  315 +{
  316 + floatx80 z;
  317 +
  318 + z.low = LIT64( 0xC000000000000000 ) | ( a.high>>1 );
  319 + z.high = ( ( (bits16) a.sign )<<15 ) | 0x7FFF;
  320 + return z;
  321 +
  322 +}
  323 +
  324 +/*----------------------------------------------------------------------------
  325 +| Takes two extended double-precision floating-point values `a' and `b', one
  326 +| of which is a NaN, and returns the appropriate NaN result. If either `a' or
  327 +| `b' is a signaling NaN, the invalid exception is raised.
  328 +*----------------------------------------------------------------------------*/
  329 +
  330 +static floatx80 propagateFloatx80NaN( floatx80 a, floatx80 b STATUS_PARAM)
  331 +{
  332 + flag aIsNaN, aIsSignalingNaN, bIsNaN, bIsSignalingNaN;
  333 +
  334 + aIsNaN = floatx80_is_nan( a );
  335 + aIsSignalingNaN = floatx80_is_signaling_nan( a );
  336 + bIsNaN = floatx80_is_nan( b );
  337 + bIsSignalingNaN = floatx80_is_signaling_nan( b );
  338 + a.low |= LIT64( 0xC000000000000000 );
  339 + b.low |= LIT64( 0xC000000000000000 );
  340 + if ( aIsSignalingNaN | bIsSignalingNaN ) float_raise( float_flag_invalid STATUS_VAR);
  341 + if ( aIsSignalingNaN ) {
  342 + if ( bIsSignalingNaN ) goto returnLargerSignificand;
  343 + return bIsNaN ? b : a;
  344 + }
  345 + else if ( aIsNaN ) {
  346 + if ( bIsSignalingNaN | ! bIsNaN ) return a;
  347 + returnLargerSignificand:
  348 + if ( a.low < b.low ) return b;
  349 + if ( b.low < a.low ) return a;
  350 + return ( a.high < b.high ) ? a : b;
  351 + }
  352 + else {
  353 + return b;
  354 + }
  355 +
  356 +}
  357 +
  358 +#endif
  359 +
  360 +#ifdef FLOAT128
  361 +
  362 +/*----------------------------------------------------------------------------
  363 +| The pattern for a default generated quadruple-precision NaN. The `high' and
  364 +| `low' values hold the most- and least-significant bits, respectively.
  365 +*----------------------------------------------------------------------------*/
  366 +#define float128_default_nan_high LIT64( 0xFFFF800000000000 )
  367 +#define float128_default_nan_low LIT64( 0x0000000000000000 )
  368 +
  369 +/*----------------------------------------------------------------------------
  370 +| Returns 1 if the quadruple-precision floating-point value `a' is a NaN;
  371 +| otherwise returns 0.
  372 +*----------------------------------------------------------------------------*/
  373 +
  374 +flag float128_is_nan( float128 a )
  375 +{
  376 +
  377 + return
  378 + ( LIT64( 0xFFFE000000000000 ) <= (bits64) ( a.high<<1 ) )
  379 + && ( a.low || ( a.high & LIT64( 0x0000FFFFFFFFFFFF ) ) );
  380 +
  381 +}
  382 +
  383 +/*----------------------------------------------------------------------------
  384 +| Returns 1 if the quadruple-precision floating-point value `a' is a
  385 +| signaling NaN; otherwise returns 0.
  386 +*----------------------------------------------------------------------------*/
  387 +
  388 +flag float128_is_signaling_nan( float128 a )
  389 +{
  390 +
  391 + return
  392 + ( ( ( a.high>>47 ) & 0xFFFF ) == 0xFFFE )
  393 + && ( a.low || ( a.high & LIT64( 0x00007FFFFFFFFFFF ) ) );
  394 +
  395 +}
  396 +
  397 +/*----------------------------------------------------------------------------
  398 +| Returns the result of converting the quadruple-precision floating-point NaN
  399 +| `a' to the canonical NaN format. If `a' is a signaling NaN, the invalid
  400 +| exception is raised.
  401 +*----------------------------------------------------------------------------*/
  402 +
  403 +static commonNaNT float128ToCommonNaN( float128 a STATUS_PARAM)
  404 +{
  405 + commonNaNT z;
  406 +
  407 + if ( float128_is_signaling_nan( a ) ) float_raise( float_flag_invalid STATUS_VAR);
  408 + z.sign = a.high>>63;
  409 + shortShift128Left( a.high, a.low, 16, &z.high, &z.low );
  410 + return z;
  411 +
  412 +}
  413 +
  414 +/*----------------------------------------------------------------------------
  415 +| Returns the result of converting the canonical NaN `a' to the quadruple-
  416 +| precision floating-point format.
  417 +*----------------------------------------------------------------------------*/
  418 +
  419 +static float128 commonNaNToFloat128( commonNaNT a )
  420 +{
  421 + float128 z;
  422 +
  423 + shift128Right( a.high, a.low, 16, &z.high, &z.low );
  424 + z.high |= ( ( (bits64) a.sign )<<63 ) | LIT64( 0x7FFF800000000000 );
  425 + return z;
  426 +
  427 +}
  428 +
  429 +/*----------------------------------------------------------------------------
  430 +| Takes two quadruple-precision floating-point values `a' and `b', one of
  431 +| which is a NaN, and returns the appropriate NaN result. If either `a' or
  432 +| `b' is a signaling NaN, the invalid exception is raised.
  433 +*----------------------------------------------------------------------------*/
  434 +
  435 +static float128 propagateFloat128NaN( float128 a, float128 b STATUS_PARAM)
  436 +{
  437 + flag aIsNaN, aIsSignalingNaN, bIsNaN, bIsSignalingNaN;
  438 +
  439 + aIsNaN = float128_is_nan( a );
  440 + aIsSignalingNaN = float128_is_signaling_nan( a );
  441 + bIsNaN = float128_is_nan( b );
  442 + bIsSignalingNaN = float128_is_signaling_nan( b );
  443 + a.high |= LIT64( 0x0000800000000000 );
  444 + b.high |= LIT64( 0x0000800000000000 );
  445 + if ( aIsSignalingNaN | bIsSignalingNaN ) float_raise( float_flag_invalid STATUS_VAR);
  446 + if ( aIsSignalingNaN ) {
  447 + if ( bIsSignalingNaN ) goto returnLargerSignificand;
  448 + return bIsNaN ? b : a;
  449 + }
  450 + else if ( aIsNaN ) {
  451 + if ( bIsSignalingNaN | ! bIsNaN ) return a;
  452 + returnLargerSignificand:
  453 + if ( lt128( a.high<<1, a.low, b.high<<1, b.low ) ) return b;
  454 + if ( lt128( b.high<<1, b.low, a.high<<1, a.low ) ) return a;
  455 + return ( a.high < b.high ) ? a : b;
  456 + }
  457 + else {
  458 + return b;
  459 + }
  460 +
  461 +}
  462 +
  463 +#endif
  464 +
... ...