Sigrid Haflínudóttirfirst dirty version of the portd7ad0686 years agoHistory
6262 lines · 182.2 KB · C
RawDownload
/*	FPMATHEM.h 	Copyright (C) 2007 John R. Hauser, Stanislav Shwartsman,		Ross Martin, Paul C. Pratt 	You can redistribute this file and/or modify it under the terms	of version 2 of the GNU General Public License as published by	the Free Software Foundation.  You should have received a copy	of the license along with this file; see the file COPYING. 	This file is distributed in the hope that it will be useful,	but WITHOUT ANY WARRANTY; without even the implied warranty of	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the	license for more details.*/ /*	Floating Point MATH implemented with software EMulation 	Based on the SoftFloat IEC/IEEE Floating-point Arithmetic	Package, Release 2b, written by John R. Hauser. 	Also uses extensions to SoftFloat, written for	Bochs (x86 achitecture simulator), by Stanislav Shwartsman.*/ /*	original SoftFloat Copyright notice: 	Written by John R. Hauser.  This work was made possible in part by the	International Computer Science Institute, located at Suite 600, 1947 Center	Street, Berkeley, California 94704.  Funding was partially provided by the	National Science Foundation under grant MIP-9311980.  The original version	of this code was written as part of a project to build a fixed-point vector	processor in collaboration with the University of California at Berkeley,	overseen by Profs. Nelson Morgan and John Wawrzynek.  More information	is available through the Web page `http://www.cs.berkeley.edu/~jhauser/	arithmetic/SoftFloat.html'. 	THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable effort has	been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT TIMES	RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS RESTRICTED TO PERSONS	AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ALL LOSSES,	COSTS, OR OTHER PROBLEMS THEY INCUR DUE TO THE SOFTWARE, AND WHO FURTHERMORE	EFFECTIVELY INDEMNIFY JOHN HAUSER AND THE INTERNATIONAL COMPUTER SCIENCE	INSTITUTE (possibly via similar legal warning) AGAINST ALL LOSSES, COSTS, OR	OTHER PROBLEMS INCURRED BY THEIR CUSTOMERS AND CLIENTS DUE TO THE SOFTWARE. 	Derivative works are acceptable, even for commercial purposes, so long as	(1) the source code for the derivative work includes prominent notice that	the work is derivative, and (2) the source code includes prominent notice with	these four paragraphs for those parts of this code that are retained.*/ /*	original Stanislav Shwartsman Copyright notice: 	This source file is an extension to the SoftFloat IEC/IEEE Floating-point	Arithmetic Package, Release 2b, written for Bochs (x86 achitecture simulator)	floating point emulation. 	THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE.  Although reasonable effort has	been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT TIMES	RESULT IN INCORRECT BEHAVIOR.  USE OF THIS SOFTWARE IS RESTRICTED TO PERSONS	AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ALL LOSSES,	COSTS, OR OTHER PROBLEMS THEY INCUR DUE TO THE SOFTWARE, AND WHO FURTHERMORE	EFFECTIVELY INDEMNIFY JOHN HAUSER AND THE INTERNATIONAL COMPUTER SCIENCE	INSTITUTE (possibly via similar legal warning) AGAINST ALL LOSSES, COSTS, OR	OTHER PROBLEMS INCURRED BY THEIR CUSTOMERS AND CLIENTS DUE TO THE SOFTWARE. 	Derivative works are acceptable, even for commercial purposes, so long as	(1) the source code for the derivative work includes prominent notice that	the work is derivative, and (2) the source code includes prominent notice with	these four paragraphs for those parts of this code that are retained. 	.*============================================================================	* Written for Bochs (x86 achitecture simulator) by	*            Stanislav Shwartsman [sshwarts at sourceforge net]	* ==========================================================================*.*/  /*	ReportAbnormalID unused 0x0204 - 0x02FF*/  /* soft float stuff */ /*	should avoid 64 bit literals.*/ typedef ui3r flag; /* 0/1 */ /*	To fix: softfloat representation of	denormalized extended precision numbers	is different than on 68881.*/ #define cIncludeFPUUnused cIncludeUnused /* ----- from original file "softfloat.h" ----- */ /*====================================================================== This C header file is part of the SoftFloat IEC/IEEE Floating-pointArithmetic Package, Release 2b. ["original SoftFloat Copyright notice" went here, included near top ofthis file.] ======================================================================*/ /*----------------------------------------------------------------------------| The macro `FLOATX80' must be defined to enable the extended double-precision| floating-point format `floatx80'.  If this macro is not defined, the| `floatx80' type will not be defined, and none of the functions that either| input or output the `floatx80' type will be defined.  The same applies to| the `FLOAT128' macro and the quadruple-precision format `float128'.*----------------------------------------------------------------------------*/#define FLOATX80#define FLOAT128 /*----------------------------------------------------------------------------| Software IEC/IEEE floating-point types.*----------------------------------------------------------------------------*/ typedef struct {	ui6b low;	unsigned short high;} floatx80; #ifdef FLOAT128typedef struct {	ui6b low, high;} float128;#endif  /* ----- end from original file "softfloat.h" ----- */  /*----------------------------------------------------------------------------| Software IEC/IEEE floating-point rounding mode.*----------------------------------------------------------------------------*/enum {	float_round_nearest_even = 0,	float_round_down         = 1,	float_round_up           = 2,	float_round_to_zero      = 3}; /*----------------------------------------------------------------------------| Floating-point rounding mode, extended double-precision rounding precision,| and exception flags.*----------------------------------------------------------------------------*/ LOCALVAR si3r float_rounding_mode = float_round_nearest_even;  /*----------------------------------------------------------------------------| Software IEC/IEEE floating-point exception flags.*----------------------------------------------------------------------------*/ enum {	float_flag_invalid   =  1,	float_flag_divbyzero =  4,	float_flag_overflow  =  8,	float_flag_underflow = 16,	float_flag_inexact   = 32};LOCALVAR si3r float_exception_flags = 0; /*----------------------------------------------------------------------------| Software IEC/IEEE extended double-precision rounding precision.  Valid| values are 32, 64, and 80.*----------------------------------------------------------------------------*/ LOCALVAR si3r floatx80_rounding_precision = 80; /*----------------------------------------------------------------------------| Primitive arithmetic functions, including multi-word arithmetic, and| division and square root approximations.  (Can be specialized to target if| desired.)*----------------------------------------------------------------------------*/ /* ----- from original file "softfloat-macros" ----- */ /*============================================================================ This C source fragment is part of the SoftFloat IEC/IEEE Floating-pointArithmetic Package, Release 2b. ["original SoftFloat Copyright notice" went here, included near top of this file.] =============================================================================*/ /*----------------------------------------------------------------------------| Shifts `a' right by the number of bits given in `count'.  If any nonzero| bits are shifted off, they are ``jammed'' into the least significant bit of| the result by setting the least significant bit to 1.  The value of `count'| can be arbitrarily large; in particular, if `count' is greater than 32, the| result will be either 0 or 1, depending on whether `a' is zero or nonzero.| The result is stored in the location pointed to by `zPtr'.*----------------------------------------------------------------------------*/ LOCALINLINEPROC shift32RightJamming( ui5b a, si4r count, ui5b *zPtr ){	ui5b z; 	if ( count == 0 ) {		z = a;	}	else if ( count < 32 ) {		z = ( a>>count ) | ( ( a<<( ( - count ) & 31 ) ) != 0 );	}	else {		z = ( a != 0 );	}	*zPtr = z; } /*----------------------------------------------------------------------------| Shifts `a' right by the number of bits given in `count'.  If any nonzero| bits are shifted off, they are ``jammed'' into the least significant bit of| the result by setting the least significant bit to 1.  The value of `count'| can be arbitrarily large; in particular, if `count' is greater than 64, the| result will be either 0 or 1, depending on whether `a' is zero or nonzero.| The result is stored in the location pointed to by `zPtr'.*----------------------------------------------------------------------------*/ LOCALINLINEPROC shift64RightJamming( ui6b a, si4r count, ui6b *zPtr ){	ui6b z; 	if ( count == 0 ) {		z = a;	}	else if ( count < 64 ) {		z = ( a>>count ) | ( ( a<<( ( - count ) & 63 ) ) != 0 );	}	else {		z = ( a != 0 );	}	*zPtr = z; } /*----------------------------------------------------------------------------| Shifts the 128-bit value formed by concatenating `a0' and `a1' right by 64| _plus_ the number of bits given in `count'.  The shifted result is at most| 64 nonzero bits; this is stored at the location pointed to by `z0Ptr'.  The| bits shifted off form a second 64-bit result as follows:  The _last_ bit| shifted off is the most-significant bit of the extra result, and the other| 63 bits of the extra result are all zero if and only if _all_but_the_last_| bits shifted off were all zero.  This extra result is stored in the location| pointed to by `z1Ptr'.  The value of `count' can be arbitrarily large.|     (This routine makes more sense if `a0' and `a1' are considered to form| a fixed-point value with binary point between `a0' and `a1'.  This fixed-| point value is shifted right by the number of bits given in `count', and| the integer part of the result is returned at the location pointed to by| `z0Ptr'.  The fractional part of the result may be slightly corrupted as| described above, and is returned at the location pointed to by `z1Ptr'.)*----------------------------------------------------------------------------*/ LOCALINLINEPROC shift64ExtraRightJamming(	ui6b a0, ui6b a1, si4r count, ui6b *z0Ptr, ui6b *z1Ptr ){	ui6b z0, z1;	si3r negCount = ( - count ) & 63; 	if ( count == 0 ) {		z1 = a1;		z0 = a0;	}	else if ( count < 64 ) {		z1 = ( a0<<negCount ) | ( a1 != 0 );		z0 = a0>>count;	}	else {		if ( count == 64 ) {			z1 = a0 | ( a1 != 0 );		}		else {			z1 = ( ( a0 | a1 ) != 0 );		}		z0 = 0;	}	*z1Ptr = z1;	*z0Ptr = z0; } /*----------------------------------------------------------------------------| Shifts the 128-bit value formed by concatenating `a0' and `a1' right by the| number of bits given in `count'.  Any bits shifted off are lost.  The value| of `count' can be arbitrarily large; in particular, if `count' is greater| than 128, the result will be 0.  The result is broken into two 64-bit pieces| which are stored at the locations pointed to by `z0Ptr' and `z1Ptr'.*----------------------------------------------------------------------------*/ LOCALINLINEPROC shift128Right(	ui6b a0, ui6b a1, si4r count, ui6b *z0Ptr, ui6b *z1Ptr ){	ui6b z0, z1;	si3r negCount = ( - count ) & 63; 	if ( count == 0 ) {		z1 = a1;		z0 = a0;	}	else if ( count < 64 ) {		z1 = ( a0<<negCount ) | ( a1>>count );		z0 = a0>>count;	}	else {		z1 = ( count < 64 ) ? ( a0>>( count & 63 ) ) : 0;		z0 = 0;	}	*z1Ptr = z1;	*z0Ptr = z0; } /*----------------------------------------------------------------------------| Shifts the 128-bit value formed by concatenating `a0' and `a1' right by the| number of bits given in `count'.  If any nonzero bits are shifted off, they| are ``jammed'' into the least significant bit of the result by setting the| least significant bit to 1.  The value of `count' can be arbitrarily large;| in particular, if `count' is greater than 128, the result will be either| 0 or 1, depending on whether the concatenation of `a0' and `a1' is zero or| nonzero.  The result is broken into two 64-bit pieces which are stored at| the locations pointed to by `z0Ptr' and `z1Ptr'.*----------------------------------------------------------------------------*/ LOCALINLINEPROC shift128RightJamming(	ui6b a0, ui6b a1, si4r count, ui6b *z0Ptr, ui6b *z1Ptr ){	ui6b z0, z1;	si3r negCount = ( - count ) & 63; 	if ( count == 0 ) {		z1 = a1;		z0 = a0;	}	else if ( count < 64 ) {		z1 = ( a0<<negCount ) | ( a1>>count ) | ( ( a1<<negCount ) != 0 );		z0 = a0>>count;	}	else {		if ( count == 64 ) {			z1 = a0 | ( a1 != 0 );		}		else if ( count < 128 ) {			z1 = ( a0>>( count & 63 ) ) | ( ( ( a0<<negCount ) | a1 ) != 0 );		}		else {			z1 = ( ( a0 | a1 ) != 0 );		}		z0 = 0;	}	*z1Ptr = z1;	*z0Ptr = z0; } /*----------------------------------------------------------------------------| Shifts the 192-bit value formed by concatenating `a0', `a1', and `a2' right| by 64 _plus_ the number of bits given in `count'.  The shifted result is| at most 128 nonzero bits; these are broken into two 64-bit pieces which are| stored at the locations pointed to by `z0Ptr' and `z1Ptr'.  The bits shifted| off form a third 64-bit result as follows:  The _last_ bit shifted off is| the most-significant bit of the extra result, and the other 63 bits of the| extra result are all zero if and only if _all_but_the_last_ bits shifted off| were all zero.  This extra result is stored in the location pointed to by| `z2Ptr'.  The value of `count' can be arbitrarily large.|     (This routine makes more sense if `a0', `a1', and `a2' are considered| to form a fixed-point value with binary point between `a1' and `a2'.  This| fixed-point value is shifted right by the number of bits given in `count',| and the integer part of the result is returned at the locations pointed to| by `z0Ptr' and `z1Ptr'.  The fractional part of the result may be slightly| corrupted as described above, and is returned at the location pointed to by| `z2Ptr'.)*----------------------------------------------------------------------------*/ LOCALINLINEPROC shift128ExtraRightJamming(	ui6b a0,	ui6b a1,	ui6b a2,	si4r count,	ui6b *z0Ptr,	ui6b *z1Ptr,	ui6b *z2Ptr){	ui6b z0, z1, z2;	si3r negCount = ( - count ) & 63; 	if ( count == 0 ) {		z2 = a2;		z1 = a1;		z0 = a0;	}	else {		if ( count < 64 ) {			z2 = a1<<negCount;			z1 = ( a0<<negCount ) | ( a1>>count );			z0 = a0>>count;		}		else {			if ( count == 64 ) {				z2 = a1;				z1 = a0;			}			else {				a2 |= a1;				if ( count < 128 ) {					z2 = a0<<negCount;					z1 = a0>>( count & 63 );				}				else {					z2 = ( count == 128 ) ? a0 : ( a0 != 0 );					z1 = 0;				}			}			z0 = 0;		}		z2 |= ( a2 != 0 );	}	*z2Ptr = z2;	*z1Ptr = z1;	*z0Ptr = z0; } /*----------------------------------------------------------------------------| Shifts the 128-bit value formed by concatenating `a0' and `a1' left by the| number of bits given in `count'.  Any bits shifted off are lost.  The value| of `count' must be less than 64.  The result is broken into two 64-bit| pieces which are stored at the locations pointed to by `z0Ptr' and `z1Ptr'.*----------------------------------------------------------------------------*/ LOCALINLINEPROC shortShift128Left(	ui6b a0, ui6b a1, si4r count, ui6b *z0Ptr, ui6b *z1Ptr ){ 	*z1Ptr = a1<<count;	*z0Ptr =		( count == 0 ) ? a0 : ( a0<<count ) | ( a1>>( ( - count ) & 63 ) ); } /*----------------------------------------------------------------------------| Adds the 128-bit value formed by concatenating `a0' and `a1' to the 128-bit| value formed by concatenating `b0' and `b1'.  Addition is modulo 2^128, so| any carry out is lost.  The result is broken into two 64-bit pieces which| are stored at the locations pointed to by `z0Ptr' and `z1Ptr'.*----------------------------------------------------------------------------*/ LOCALINLINEPROC add128(	ui6b a0, ui6b a1, ui6b b0, ui6b b1, ui6b *z0Ptr, ui6b *z1Ptr ){	ui6b z1; 	z1 = a1 + b1;	*z1Ptr = z1;	*z0Ptr = a0 + b0 + ( z1 < a1 );} /*----------------------------------------------------------------------------| Adds the 192-bit value formed by concatenating `a0', `a1', and `a2' to the| 192-bit value formed by concatenating `b0', `b1', and `b2'.  Addition is| modulo 2^192, so any carry out is lost.  The result is broken into three| 64-bit pieces which are stored at the locations pointed to by `z0Ptr',| `z1Ptr', and `z2Ptr'.*----------------------------------------------------------------------------*/ LOCALINLINEPROC add192(	ui6b a0,	ui6b a1,	ui6b a2,	ui6b b0,	ui6b b1,	ui6b b2,	ui6b *z0Ptr,	ui6b *z1Ptr,	ui6b *z2Ptr){	ui6b z0, z1, z2;	si3r carry0, carry1; 	z2 = a2 + b2;	carry1 = ( z2 < a2 );	z1 = a1 + b1;	carry0 = ( z1 < a1 );	z0 = a0 + b0;	z1 += carry1;

Lines 1–500 of 6262.

Keyboard shortcuts

General

Show keyboard shortcuts
?
Focus the current search
/

Issues and pull requests

Open filters
f

Issues, pull requests and Inbox

Select next item
j↓
Select previous item
k↑
Open selected item
Enter
Dismiss filters or clear selection
Esc

Repository code

Go to file
t