1 #include "lib9.h" 2 #include <float.h> 3 #include "mathi.h" 4 5 void 6 FPinit(void) 7 { 8 _controlfp(_EM_INEXACT,_MCW_EM); // abort on underflow, etc. 9 } 10 11 ulong 12 getFPstatus(void) 13 { 14 ulong fsr = 0, fsr32 = _statusfp(); 15 if(fsr32&_SW_INEXACT) fsr |= INEX; 16 if(fsr32&_SW_OVERFLOW) fsr |= OVFL; 17 if(fsr32&_SW_UNDERFLOW) fsr |= UNFL; 18 if(fsr32&_SW_ZERODIVIDE) fsr |= ZDIV; 19 if(fsr32&_SW_INVALID) fsr |= INVAL; 20 return fsr; 21 } 22 23 ulong 24 FPstatus(ulong fsr, ulong mask) 25 { 26 ulong old = getFPstatus(); 27 fsr = (fsr&mask) | (old&~mask); 28 if(fsr!=old){ 29 _clearfp(); 30 if(fsr){ 31 ulong fcr = _controlfp(0,0); 32 double x = 1., y = 1e200, z = 0.; 33 _controlfp(_MCW_EM,_MCW_EM); 34 if(fsr&INEX) z = x + y; 35 if(fsr&OVFL) z = y*y; 36 if(fsr&UNFL) z = (x/y)/y; 37 if(fsr&ZDIV) z = x/z; 38 if(fsr&INVAL) z = z/z; 39 _controlfp(fcr,_MCW_EM); 40 } 41 } 42 return(old&mask); 43 } 44 45 ulong 46 getFPcontrol(void) 47 { 48 ulong fcr, fcr32 = _controlfp(0,0); 49 switch(fcr32&_MCW_RC){ 50 case _RC_NEAR: fcr = RND_NR; break; 51 case _RC_DOWN: fcr = RND_NINF; break; 52 case _RC_UP: fcr = RND_PINF; break; 53 case _RC_CHOP: fcr = RND_Z; break; 54 } 55 if(!(fcr32&_EM_INEXACT)) fcr |= INEX; 56 if(!(fcr32&_EM_OVERFLOW)) fcr |= OVFL; 57 if(!(fcr32&_EM_UNDERFLOW)) fcr |= UNFL; 58 if(!(fcr32&_EM_ZERODIVIDE)) fcr |= ZDIV; 59 if(!(fcr32&_EM_INVALID)) fcr |= INVAL; 60 return fcr; 61 } 62 63 ulong 64 FPcontrol(ulong fcr, ulong mask) 65 { 66 ulong old = getFPcontrol(); 67 ulong fcr32 = _MCW_EM, mask32 = _MCW_RC|_MCW_EM; 68 fcr = (fcr&mask) | (old&~mask); 69 if(fcr&INEX) fcr32 ^= _EM_INEXACT; 70 if(fcr&OVFL) fcr32 ^= _EM_OVERFLOW; 71 if(fcr&UNFL) fcr32 ^= _EM_UNDERFLOW; 72 if(fcr&ZDIV) fcr32 ^= _EM_ZERODIVIDE; 73 if(fcr&INVAL) fcr32 ^= _EM_INVALID; 74 switch(fcr&RND_MASK){ 75 case RND_NR: fcr32 |= _RC_NEAR; break; 76 case RND_NINF: fcr32 |= _RC_DOWN; break; 77 case RND_PINF: fcr32 |= _RC_UP; break; 78 case RND_Z: fcr32 |= _RC_CHOP; break; 79 } 80 _controlfp(fcr32,mask32); 81 return(old&mask); 82 } 83