ioldouble

Post your FreeBASIC tips and tricks here. Please don’t post your code without including an explanation.
srvaldez
Posts: 2525
Joined: Sep 25, 2005 21:54

ioldouble

Postby srvaldez » Jun 20, 2011 1:12

part 1 of 3

Code: Select all

 /' Extended precision arithmetic functions for long double I/O.
    This program has been placed in the public domain by Stephen L Moshier.
    http://www.moshier.net/#Rubid_pc
    file smldbl12.zip '/
   
' partial translation to FreeBASIC by srvaldez also known as jack on the
' thinBasic forum.
' I was going to incorporate these routines into the real10 lib but just don't
' have the time and decided to post the code, maybe it will be usefull to
' someone.

#define NE 10
#define NI (NE+3)
const EXPONENT = 1
const MANTISSA = 2
#define NBITS ((NI-4)*16)
#define NDEC (NBITS*8/27)
#define EXONE (&h3fff)
#define ANSIC 1
#define VOLATILE
#define DOMAIN      1   ' argument domain error
#define SING      2   ' argument singularity
#define OVERFLOW   3   ' overflow range error
#define UNDERFLOW   4   ' underflow range error
#define TLOSS      5   ' total loss of precision
#define PLOSS      6   ' partial loss of precision

#define EDOM      33
#define ERANGE      34
#define XPD 0
#define NANS
#define INFINITY
#define LONGBITS (8 * sizeof(long))
#define LLONGBITS (8 * sizeof(longint))

Type efloat
    ef(0 to NI-1) As ushort
End Type

'Extern ezero alias "ezero" as efloat
'Extern eone  alias "eone"  as efloat
'Extern epi  alias "epi"  as ushort ptr
dim shared ezero(NE-1) as ushort

dim shared eone(NE-1) as ushort
            eone(NE-1) = &h3fff
            eone(NE-2) = &h8000

dim shared etwo(NE-1) as ushort
            etwo(NE-1) = &h4000
            etwo(NE-2) = &h8000

dim shared ethree(NE-1) as ushort
           ethree(NE-1) = &h4000
           ethree(NE-2) = &hC000

dim shared eten(NE-1) as ushort
           eten(NE-1) = &h4002
           eten(NE-2) = &hA000
dim shared epi(NE-1) as ushort = {&h2902, &h1cd1, &h80dc, &h628b,_
  &hc4c6, &hc234, &h2168, &hdaa2, &hc90f, &h4000}
dim shared esqrt2(NE-1) as ushort = {&h1d6f, &hbe9f, &h754a, &h89b3,_
  &h597d, &h6484, &hf9de, &hf333, &hb504, &h3fff}
dim shared rndprc as integer = NBITS
dim shared nan113(8-1) as ushort = {0, 0, 0, 0, 0, 0, &hc000, &hffff}
dim shared nan64(6-1) as ushort = {0, 0, 0, &hc000, &hffff, 0}
dim shared nan53(4-1) as ushort = {0, 0, 0, &hfff8}
dim shared nan24(2-1) as ushort = {0, &hffc0}
dim shared equot(NI) as ushort = {0}
dim shared merror as integer = 0

dim shared rlast as integer = -1
dim shared rw as integer = 0
dim shared rmsk as ushort = 0
dim shared rmbit as ushort = 0
dim shared rebit as ushort = 0
dim shared re as integer = 0
dim shared rbit(NI) as ushort = {0,0,0,0,0,0,0,0}
dim shared subflg as integer = 0
dim shared lenldstr as integer

' Notice: the order of appearance of the following
' messages is bound to the error codes defined
' in mconf.h.

dim shared as string*25 ermsg(7-1) = {_
"unknown",_
"domain",_
"singularity",_
"overflow",_
"underflow",_
"total loss of precision",_
"partial loss of precision"}

type cmplx
   as double r
   as double i
end type


type longdouble
    ar(4) as ushort
end type


/'
; Clear out entire external format number.
;
; unsigned short x[];
; eclear( x );
'/
sub eclear(x as ushort ptr)
    dim as integer i

    for i=0 to NE-1
        *x = 0
        x += 1
    next
end sub

/' Move external format number from a to b.
 *
 * emov( a, b );
 '/

sub emov(a as ushort ptr, b as ushort ptr)
    dim as integer i

    for i=0 to NE-1
        *b = *a
        b += 1
        a += 1
    next
end sub

' Check if e-type number is not a number.

function eisnan(x as ushort ptr) as integer
    dim as integer i
    ' NaN has maximum exponent
    if (x[NE-1] and &h7fff) <> &h7fff then return 0
    ' ... and non-zero significand field.
    for i=0 to NE-2
        x += 1
        if *(x-1) <> 0 then
            x += 1
            return 1
        endif
    next
    return 0
end function

/'
;   Negate external format number
;
;   unsigned short x[NE];
;   eneg( x );
'/

sub eneg(x as ushort ptr)
    if eisnan(x) then return
    x[NE-1] xor= &h8000 ' Toggle the sign bit
end sub

/' Return 1 if external format number is negative,
 * else return zero.
 '/
function eisneg(x as ushort ptr) as integer
    if eisnan(x) then return 0
    if x[NE-1] and &h8000 then
        return 1
    else
        return 0
    endif
end function

/'
; Fill entire number, including exponent and significand, with
; largest possible number.  These programs implement a saturation
; value that is an ordinary, legal number.  A special value
; "infinity" may also be implemented; this would require tests
; for that value and implementation of special rules for arithmetic
; operations involving inifinity.
'/

sub einfin(x as ushort ptr)
    dim as integer i

    for i=0 to NE-2
        *x = 0
        x += 1
    next
    *x or= 32767
end sub

/' Return 1 if external format number has maximum possible exponent,
 * else return zero.
 '/
function eisinf(x as ushort ptr) as integer
    if (x[NE-1] and &h7fff) = &h7fff then
        if eisnan(x) then return 0
        return 1
    else
        return 0
    endif
end function

/' Move in external format number,
 * converting it to internal format.
 '/
sub emovi(a as ushort ptr, b as ushort ptr)
    dim  as ushort ptr p, q
    dim as integer i

    q = b
    p = a + (NE-1)   ' point to last word of external number
    ' get the sign bit
    q += 1
    if *p and &h8000 then
        *(q-1) = &hffff
    else
        *(q-1) = 0
    endif
    ' get the exponent
    *q = *p
    p -= 1
    *q and= &h7fff   ' delete the sign bit
    q += 1
    if (*(q-1) and &h7fff) = &h7fff then
        if eisnan(a) then
            *q = 0
            q += 1
            for i=3 to NI-1
                *q = *p
                q += 1
                p -= 1
            next
            return
        endif
        for i=2 to NI-1
            *q = 0
            q += 1
        next
        return
    endif

    ' clear high guard word
    *q = 0
    q += 1
    ' move in the significand
    for i=0 to NE-2
        *q = *p
        q += 1
        p -= 1
    next
    ' clear low guard word
    *q = 0
end sub

' Return nonzero if internal format number is a NaN.

function eiisnan (x as ushort ptr) as integer
    dim as integer i

    if (x[EXPONENT] and &h7fff) = &h7fff then
        for i=MANTISSA+1 to NI-1
            if x[i] <> 0 then return 1
        next
    endif
    return 0
end function

function mtherr(byval nam as string, byval code as integer) as integer

    ' Display string passed by calling program,
    ' which is supposed to be the name of the
    ' function in which the error occurred:

    print nam;

    ' Set global error message word
    merror = code

    ' Display error message defined
    ' by the code argument.

    if (code <= 0) or (code >= 6) then code = 0
    print " error "; ermsg(code)

    ' Return to calling
    ' program

    return 0
end function

' NaN bit patterns

sub enan (nan as ushort ptr, size as integer)
    dim as integer i, n
    dim as ushort ptr p

    select case size
        case 113
            n = 8
            p = @nan113(0)
            exit select

        case 64
            n = 6
            p = @nan64(0)
            exit select

        case 53
            n = 4
            p = @nan53(0)
            exit select

        case 24
            n = 2
            p = @nan24(0)
            exit select

        case NBITS
            for i=0 to NE-3
                *nan = 0
                nan += 1
            next
            *nan = &hc000
            nan += 1
            *nan = &h7fff
            nan += 1
            return

        case NI*16
            *nan = 0
            nan += 1
            *nan = &h7fff
            nan += 1
            *nan = 0
            nan += 1
            *nan = &hc000
            nan += 1
            for i=4 to NI-1
                *nan = 0
                nan += 1
            next
            return
        case else
            mtherr( "enan", DOMAIN )
            return
    end select
    for i=0 to n-1
        *nan = *p
        nan += 1
        p += 1
    next
end sub

/' Move internal format number out,
 * converting it to external format.
 '/
sub emovo(a as ushort ptr, b as ushort ptr)
    dim as ushort ptr p, q
    dim as ushort i

    p = a
    q = b + (NE-1) ' point to output exponent
    ' combine sign and exponent
    i = *p
    p += 1
    if i then
        *q = *p or &h8000
        q -= 1
        p += 1
    else
        *q = *p
        q -= 1
        p += 1
    endif
    if *(p-1) = &h7fff then
        if eiisnan(a) then
            enan( b, NBITS )
            return
        endif
        einfin(b)
        return
    endif
    ' skip over guard word
    p += 1
    ' move the significand
    for i=0 to NE-2
        *q = *p
        q -= 1
        p += 1
    next
end sub

' Clear out internal format number.

sub ecleaz(xi as ushort ptr)
    dim as integer i

    for i=0 to NI-1
        *xi = 0
        xi += 1
    next
end sub

' same, but don't touch the sign.

sub ecleazs(xi as ushort ptr)
    dim as integer i

    xi += 1
    for i=0 to NI-2
        *xi = 0
        xi += 1
    next
end sub

' Move internal format number from a to b.

sub emovz(a as ushort ptr, b as ushort ptr)
    dim as integer i

    for i=0 to NI-2
        *b = *a
        b += 1
        a += 1
    next
    ' clear low guard word
    *b = 0
end sub


/'
;   Compare significands of numbers in internal format.
;   Guard words are included in the comparison.
;
;   unsigned short a[NI], b[NI];
;   cmpm( a, b );
;
;   for the significands:
;   returns   +1 if a > b
;       0 if a == b
;      -1 if a < b
'/

function ecmpm(a as ushort ptr, b as ushort ptr) as integer
    dim as integer i

    a += MANTISSA ' skip up to significand area
    b += MANTISSA
    for i=MANTISSA to NI-1
        a += 1
        b += 1
        if *(a-1) <> *(b-1) then goto difrnt
    next
    return 0

    difrnt:
    a -= 1
    b -= 1
    if *a > *b then
        return 1
    else
        return -1
    endif
end function

'   Shift significand down by 1 bit

sub eshdn1(x as ushort ptr)
    dim as ushort bits
    dim as integer i

    x += MANTISSA   ' point to significand area

    bits = 0
    for i=MANTISSA to NI-1
        if *x and 1 then bits or= 1
        *x shr= 1
        if bits and 2 then *x or= &h8000
        bits shl= 1
        x += 1
    next   
end sub


'   Shift significand up by 1 bit

sub eshup1(x as ushort ptr)
    dim as ushort bits
    dim as integer i

    x += NI-1
    bits = 0

    for i=MANTISSA to NI-1
        if *x and &h8000 then bits or= 1
        *x shl= 1
        if bits and 2 then *x or= 1
        bits shl= 1
        x -= 1
    next
end sub

'   Shift significand down by 8 bits

sub eshdn8(x as ushort ptr)
    dim as ushort newbyt, oldbyt
    dim as integer i

    x += MANTISSA
    oldbyt = 0
    for i=MANTISSA to NI-1
        newbyt = *x shl 8
        *x shr= 8
        *x or= oldbyt
        oldbyt = newbyt
        x += 1
    next
end sub

'   Shift significand up by 8 bits

sub eshup8(x as ushort ptr)
    dim as integer i
    dim as ushort newbyt, oldbyt

    x += NI-1
    oldbyt = 0

    for i=MANTISSA to NI-1
        newbyt = *x shr 8
        *x shl= 8
        *x or= oldbyt
        oldbyt = newbyt
        x -= 1
    next
end sub

'   Shift significand up by 16 bits

sub eshup6(x as ushort ptr)
    dim as integer i
    dim as ushort ptr p

    p = x + MANTISSA
    x += MANTISSA + 1

    for i=MANTISSA to NI-2
        *p = *x
        p += 1
        x += 1
    next

    *p = 0
end sub

'   Shift significand down by 16 bits

sub eshdn6(x as ushort ptr)
    dim as integer i
    dim as ushort ptr p

    x += NI-1
    p = x + 1

    for i=MANTISSA to NI-2
        p -= 1
        x -= 1
        *p = *x
    next
    p -= 1
    *p = 0
end sub

/'
;   Add significands
;   x + y replaces y
'/

sub eaddm(x as ushort ptr, y as ushort ptr)
    dim as ulong a
    dim as integer i
    dim as uinteger carry

    x += NI-1
    y += NI-1
    carry = 0
    for i=MANTISSA to NI-1
        a = cast(ulong ,*x) + cast(ulong, *y) + carry
        if a and &h10000 then
            carry = 1
        else
            carry = 0
        endif
        *y = cast(ushort, a)
        x -= 1
        y -= 1
    next
end sub

/'
;   Subtract significands
;   y - x replaces y
'/

sub esubm(x as ushort ptr, y as ushort ptr)
    dim as ulong a
    dim as integer i
    dim as uinteger carry

    x += NI-1
    y += NI-1
    carry = 0
    for i=MANTISSA to NI-1
        a = cast(ulong, *y) - cast(ulong, *x) - carry
        if a and &h10000 then
            carry = 1
        else
            carry = 0
        endif
        *y = cast(ushort, a)
        x -= 1
        y -= 1
    next
end sub

' Divide significands

'static unsigned short equot[NI] = {0}; /* was static */


/' Multiply significand of e-type number b
by 16-bit quantity a, e-type result to c. '/

sub m16m(a as ushort, b as ushort ptr, c as ushort ptr)
    dim as ushort ptr pp
    dim as ulong carry
    dim as ushort ptr ps
    dim as ushort p(NI)
    dim as ulong aa, m
    dim as integer i

    aa = a
    pp = @p(NI-2)
    *pp = 0
    pp += 1
    *pp = 0
    ps = @b[NI-1]

    for i=MANTISSA+1 to NI-1
        if *ps = 0 then
            ps -= 1
            pp -= 1
            *(pp-1) = 0
        else
            m = cast(ulong, aa * *ps)
            ps -= 1
            carry = (m and &hffff) + *pp
            *pp = cast(ushort, carry)
            pp -= 1
            carry = (carry shr 16) + (m shr 16) + *pp
            *pp = cast(ushort, carry)
            *(pp-1) = carry shr 16
        endif
    next
    for i=MANTISSA to NI-1
        c[i] = p(i)
    next
end sub


/' Divide significands. Neither the numerator nor the denominator
is permitted to have its high guard word nonzero.  '/


function edivm(den as ushort ptr, num as ushort ptr) as integer
    dim as integer i
    dim as ushort ptr p
    dim as ulong tnum
    dim as ushort j, tdenm, tquot
    dim as ushort tprod(NI+1)

    p = @equot(0)
    *p = num[0]
    p += 1
    *p = num[1]
    p += 1

    for i=MANTISSA to NI-1
        *p = 0
        p += 1
    next
    eshdn1( num )
    tdenm = den[MANTISSA+1]
    for i=MANTISSA to NI-1
        ' Find trial quotient digit (the radix is 65536).
        tnum = ((cast(ulong, num[MANTISSA])) shl 16) + num[MANTISSA+1]
        ' Do not execute the divide instruction if it will overflow. */
        if (cast(ulong, tdenm) * &hffffL) < tnum then
            tquot = &hffff
        else
            tquot = tnum \ tdenm
        endif

'        ? "tquot = ------>";tquot
        ' Prove that the divide worked.
    /'
        tcheck = (unsigned long )tquot * tdenm;
        if( tnum - tcheck > tdenm )
            tquot = 0xffff;
    '/
        ' Multiply denominator by trial quotient digit.
        m16m( tquot, den, @tprod(0) )
        ' The quotient digit may have been overestimated.
        if ecmpm( @tprod(0), num ) > 0 then
            tquot -= 1
            esubm( den, @tprod(0) )
            if ecmpm( @tprod(0), num ) > 0 then
                tquot -= 1
                esubm( den, @tprod(0) )
            endif
        endif
    /'
        if( ecmpm( tprod, num ) > 0 )
            {
            eshow( "tprod", tprod );
            eshow( "num  ", num );
            printf( "tnum = %08lx, tden = %04x, tquot = %04x\n",
                 tnum, den[MANTISSA+1], tquot );
            }
    '/
        esubm( @tprod(0), num )
    /'
        if( ecmpm( num, den ) >= 0 )
            {
            eshow( "num  ", num );
            eshow( "den  ", den );
            printf( "tnum = %08lx, tden = %04x, tquot = %04x\n",
                 tnum, den[MANTISSA+1], tquot );
            }
    '/
        equot(i) = tquot
        eshup6(num)
    next
    ' test for nonzero remainder after roundoff bit
    p = @num[MANTISSA]
    j = 0
    for i=MANTISSA to NI-1
        j or= *p
        p += 1
    next
    if j then j = 1

    for i=0 to NI-1
        num[i] = equot(i)
    next

    return cast(integer, j )
end function


/' Divide significands. Neither the numerator nor the denominator
is permitted to have its high guard word nonzero.  '/

' Multiply significands
function emulm(a as ushort ptr, b as ushort ptr) as integer
    dim as ushort ptr p, q
    dim as ushort pprod(NI)
    dim as ushort j
    dim as integer i

    equot(0) = b[0]
    equot(1) = b[1]
    for i=MANTISSA to NI-1
        equot(i) = 0
    next

    j = 0
    p = @a[NI-1]
    q = @equot(NI-1)
    for i=MANTISSA+1 to NI-1
        if *p = 0 then
            p -= 1
        else
            m16m( *p, b, @pprod(0) )
            p -= 1
            eaddm(@pprod(0), @equot(0))
        endif
        j or= *q
        eshdn6(@equot(0))
    next

    for i=0 to NI-1
        b[i] = equot(i)
    next

    ' return flag for lost nonzero bits
    return cast(integer,j )
end function

/' Compare two e type numbers.
 *
 * unsigned short a[NE], b[NE];
 * ecmp( a, b );
 *
 *  returns +1 if a > b
 *           0 if a == b
 *          -1 if a < b
 *          -2 if either a or b is a NaN.
 '/
function ecmp(a as ushort ptr, b as ushort ptr) as integer
    dim as ushort ai(NI), bi(NI)
    dim as ushort ptr p, q
    dim as integer i, msign

    if (eisnan (a)  or eisnan (b)) then return -2
    emovi( a, @ai(0) )
    p = @ai(0)
    emovi( b, @bi(0) )
    q = @bi(0)

    if *p <> *q then
        ' the signs are different
        ' -0 equals + 0
        for i=1 to NI-2
            if ai(i) <> 0 then goto nzro
            if bi(i) <> 0 then goto nzro
        next
        return 0
    nzro:
        if *p = 0 then
            return 1
        else
            return -1
        endif
    endif
    ' both are the same sign
    if *p = 0 then
        msign = 1
    else
        msign = -1
    endif
    i = NI-1
    do
        p += 1
        q += 1
        if *(p-1) <> *(q-1) then goto diff
        i -= 1
    loop while( i > 0 )

    return 0   ' equality



    diff:
    p -= 1
    q -= 1
    if *p > *q then
        return msign      ' p is bigger
    else
        return -msign   ' p is littler
    endif
end function

/'
;   Shift significand
;
;   Shifts significand area up or down by the number of bits
;   given by the variable sc.
'/
function eshift(x as ushort ptr, sc as integer) as integer
    dim as ushort lost
    dim as ushort ptr p

    if sc = 0 then return 0

    lost = 0
    p = x + NI-1

    if sc < 0 then
        sc = -sc
        while sc >= 16
            lost or= *p   ' remember lost bits
            eshdn6(x)
            sc -= 16
        wend

        while sc >= 8
            lost or= *p and &hff
            eshdn8(x)
            sc -= 8
        wend

        while sc > 0
            lost or= *p and 1
            eshdn1(x)
            sc -= 1
        wend
    else
        while sc >= 16
            eshup6(x)
            sc -= 16
        wend

        while sc >= 8
            eshup8(x)
            sc -= 8
        wend

        while sc > 0
            eshup1(x)
            sc -= 1
        wend
    endif
    if lost then lost = 1
    return cast(integer, lost )
end function

/'
;   normalize
;
; Shift normalizes the significand area pointed to by argument
; shift count (up = positive) is returned.
'/
function enormlz(x as ushort ptr) as integer
    dim as ushort ptr p
    dim as integer sc

    sc = 0
    p = @x[MANTISSA]
    if *p <> 0 then   goto normdn
    p += 1
    if *p and &h8000 then   return 0   ' already normalized
    while *p = 0
        eshup6(x)
        sc += 16
    /' With guard word, there are NBITS+16 bits available.
     * return true if all are zero.
     '/
        if sc > NBITS then return sc
    wend
    ' see if high byte is zero
    while (*p and &hff00) = 0
        eshup8(x)
        sc += 8
    wend
    ' now shift 1 bit at a time
    while (*p  and &h8000) = 0
        eshup1(x)
        sc += 1
        if sc > (NBITS+16) then
            mtherr( "enormlz", UNDERFLOW )
            return sc
        endif
    wend
    return sc

    /' Normalize by shifting down out of the high guard word
       of the significand '/
    normdn:

    if *p and &hff00 then
        eshdn8(x)
        sc -= 8
    endif
    while *p <> 0
        eshdn1(x)
        sc -= 1

        if sc < -NBITS then
            mtherr( "enormlz", OVERFLOW )
            return( sc )
        endif
    wend
    return sc
end function

/'
 * Normalize and round off.
 *
 * The internal format number to be rounded is "s".
 * Input "lost" indicates whether the number is exact.
 * This is the so-called sticky bit.
 *
 * Input "subflg" indicates whether the number was obtained
 * by a subtraction operation.  In that case if lost is nonzero
 * then the number is slightly smaller than indicated.
 *
 * Input "exp" is the biased exponent, which may be negative.
 * the exponent field of "s" is ignored but is replaced by
 * "exp" as adjusted by normalization and rounding.
 *
 * Input "rcntrl" is the rounding control.
 '/
srvaldez
Posts: 2525
Joined: Sep 25, 2005 21:54

Postby srvaldez » Jun 20, 2011 1:14

part 2 of 3

Code: Select all


sub emdnorm(s as ushort ptr, lost as integer, isubflg as integer, exp1  as long, rcntrl as integer)
    dim as integer i, j
    dim as unsigned short r

    ' Normalize
    j = enormlz( s )

    ' a blank significand could mean either zero or infinity.

    exp1 -= j

    if (j > NBITS) and (exp1 < 32767L) then
        ecleazs( s )
        return
    endif

    if exp1 < 0L then
        if exp1 > cast(long ,(-NBITS-1) ) then
            j = cast(integer, exp1)
            i = eshift( s, j )
            if i then lost = 1
        else
            ecleazs( s )
            return
        endif
    endif
    ' Round off, unless told not to by rcntrl.
    if rcntrl = 0 then goto mdfin
    ' Set up rounding parameters if the control register changed.
    if rndprc <> rlast then
        ecleaz( @rbit(0) )
        select case rndprc
            default:
            case NBITS:
                rw = NI-1 ' low guard word
                rmsk = &hffff
                rmbit = &h8000
                rebit = 1
                re = rw - 1
                exit select
            case 113:
                rw = 10
                rmsk = &h7fff
                rmbit = &h4000
                rebit = &h8000
                re = rw
                exit select
            case 64:
                rw = 7
                rmsk = &hffff
                rmbit = &h8000
                rebit = 1
                re = rw-1
                exit select
    ' For DEC arithmetic
            case 56:
                rw = 6
                rmsk = &hff
                rmbit = &h80
                rebit = &h100
                re = rw
                exit select
            case 53:
                rw = 6
                rmsk = &h7ff
                rmbit = &h0400
                rebit = &h800
                re = rw
                exit select
            case 24:
                rw = 4
                rmsk = &hff
                rmbit = &h80
                rebit = &h100
                re = rw
                exit select
        end select
        rbit(re) = rebit
        rlast = rndprc
    endif

    /' Shift down 1 temporarily if the data structure has an implied
     * most significant bit and the number is denormal.
     * For rndprc = 64 or NBITS, there is no implied bit.
     '/
    if (exp1 <= 0) and (rndprc <> 64) and (rndprc <> NBITS) then
        lost or= s[NI-1] and 1
        eshdn1(s)
    endif
    /' Clear out all bits below the rounding bit,
     * remembering in r if any were nonzero.
     '/
    r = s[rw] and rmsk
    if rndprc < NBITS then
        i = rw + 1
        while i < NI
            if s[i] then r or= 1
            s[i] = 0
            i += 1
        wend
    endif
    s[rw] and= (not rmsk)
    if (r and rmbit) <> 0 then
        if r = rmbit then
            if lost = 0 then
                ' round to even
                if (s[re] and rebit) = 0 then goto mddone
            else
                if subflg <> 0 then   goto mddone
            endif
        endif
        eaddm( @rbit(0), s )
    endif
    mddone:
    if (exp1 <= 0) and (rndprc <> 64) and (rndprc <> NBITS) then
        eshup1(s)
    endif
    if s[2] <> 0 then
        ' overflow on roundoff
        eshdn1(s)
        exp1 += 1
    endif
    mdfin:
    s[NI-1] = 0
    if exp1 >= 32767L then
        s[1] = 32767
        for i=2 to NI-2
            s[i] = 0
        next

        return
    endif
    if exp1 < 0 then
        s[1] = 0
    else
        s[1] = cast(ushort, exp1)
    endif
end sub

sub eadd1(a as ushort ptr, b as ushort ptr, c as ushort ptr)
    dim as ushort ai(NI), bi(NI), ci(NI)
    dim as integer i, lost, j, k
    dim as long lt, lta, ltb

    if eisinf(a) then
        emov(a,c)
        if subflg then eneg(c)
        return
    endif
    if eisinf(b) then
        emov(b,c)
        return
    endif
    emovi( a, @ai(0) )
    emovi( b, @bi(0) )
    if subflg then ai(0) = not ai(0)

    'compare exponents
    lta = ai(EXPONENT)
    ltb = bi(EXPONENT)
    lt = lta - ltb
    if lt > 0L then
        ' put the larger number in bi
        emovz( @bi(0), @ci(0) )
        emovz( @ai(0), @bi(0) )
        emovz( @ci(0), @ai(0) )
        ltb = bi(EXPONENT)
        lt = -lt
    endif
    lost = 0
    if lt <> 0L then
        if lt < cast(long, (-NBITS-1) )   goto done ' answer same as larger addend
        k = cast(integer, lt)
        lost = eshift( @ai(0), k ) ' shift the smaller number down
    else
    ' exponents were the same, so must compare significands
        i = ecmpm( @ai(0), @bi(0) )
        if i = 0 then
            ' the numbers are identical in magnitude
            ' if different signs, result is zero
            if ai(0) <> bi(0) then
                eclear(c)
                return
            endif
            ' if same sign, result is double
            ' double denomalized tiny number
            if (bi(EXPONENT) = 0) and ((bi(3) and &h8000) = 0) then
                eshup1( @bi(0) )
                goto done
            endif
            ' add 1 to exponent unless both are zero!
            for j=1 to NI-2
                if bi(j) <> 0 then
    ' This could overflow, but let emovo take care of that.
                    ltb += 1
                    exit for
                endif
            next
            bi(EXPONENT) = cast(ushort, ltb)
            goto done
        endif
        if i > 0 then
            ' put the larger number in bi
            emovz( @bi(0), @ci(0) )
            emovz( @ai(0), @bi(0) )
            emovz( @ci(0), @ai(0) )
        endif
    endif
    if ai(0) = bi(0) then
        eaddm( @ai(0), @bi(0) )
        subflg = 0
    else
        esubm( @ai(0), @bi(0) )
        subflg = 1
    endif
    emdnorm( @bi(0), lost, subflg, ltb, 64 )

    done:
    emovo( @bi(0), c )
end sub

/'
;   Add.
;
;   unsigned short a[NE], b[NE], c[NE];
;   eadd( a, b, c );    c = b + a
'/

sub eadd(a as ushort ptr, b as ushort ptr, c as ushort ptr)
    ' NaN plus anything is a NaN.
    if eisnan(a) then
        emov(a,c)
        return
    endif
    if eisnan(b) then
        emov(b,c)
        return
    endif
    /' Infinity minus infinity is a NaN.
     * Test for adding infinities of opposite signs.
     '/
    if eisinf(a) and eisinf(b) and ((eisneg(a) xor eisneg(b)) <> 0) then
        mtherr( "eadd", DOMAIN )
        enan( c, NBITS )
        return
    endif
    subflg = 0
    eadd1( a, b, c )
end sub

/'
;   Subtract external format numbers.
;
;   unsigned short a[NE], b[NE], c[NE];
;   esub( a, b, c );    c = b - a
'/

'static int subflg = 0;

sub esub(a as ushort ptr, b as ushort ptr, c as ushort ptr)

    if eisnan(a) then
        emov (a, c)
        return
    endif
    if eisnan(b) then
        emov(b,c)
        return
    endif
    /' Infinity minus infinity is a NaN.
     * Test for subtracting infinities of the same sign.
     '/
    if (eisinf(a) and eisinf(b) and ((eisneg (a) xor eisneg (b)) = 0)) then
        mtherr( "esub", DOMAIN )
        enan( c, NBITS )
        return
    endif
    subflg = 1
    eadd1( a, b, c )
end sub

/'
;   Divide.
;
;   unsigned short a[NE], b[NE], c[NE];
;   ediv( a, b, c );   c = b / a
'/
sub ediv(a as ushort ptr, b as ushort ptr, c as ushort ptr)

    dim as ushort ai(NI), bi(NI)
    dim as integer i
    dim as long lt, lta, ltb

    ' Return any NaN input.
    if eisnan(a) then
        emov(a,c)
        return
    endif
    if eisnan(b) then
        emov(b,c)
        return
    endif
    ' Zero over zero, or infinity over infinity, is a NaN.
    if ((ecmp(a,@ezero(0)) = 0) and (ecmp(b,@ezero(0)) = 0)) or (eisinf (a) and eisinf (b)) then
        mtherr( "ediv", DOMAIN )
        enan( c, NBITS )
        return
    endif

    ' Infinity over anything else is infinity.
    if eisinf(b) then
        if (eisneg(a) xor eisneg(b)) then
            *(c+(NE-1)) = &h8000
        else
            *(c+(NE-1)) = 0
        endif
        einfin(c)
        return
    endif
    if eisinf(a) then
        eclear(c)
        return
    endif
    emovi( a, @ai(0) )
    emovi( b, @bi(0) )
    lta = ai(EXPONENT)
    ltb = bi(EXPONENT)
    if bi(EXPONENT) = 0 then
        ' See if numerator is zero.
        for i=1 to NI-2
            if bi(i) <> 0 then
                ltb -= enormlz( @bi(0) )
                goto dnzro1
            endif
        next
        eclear(c)
        return
    endif
    dnzro1:

    if ai(EXPONENT) = 0 then
        ' possible divide by zero
        for i=1 to NI-2
            if ai(i) <> 0 then
                lta -= enormlz( @ai(0) )
                goto dnzro2
            endif
        next
        if ai(0) = bi(0) then
            *(c+(NE-1)) = 0
        else
            *(c+(NE-1)) = &h8000
        endif
        einfin(c)
        mtherr( "ediv", SING )
        return
    endif
    dnzro2:

    i = edivm( @ai(0), @bi(0) )
    ' calculate exponent
    lt = ltb - lta + EXONE
    emdnorm( @bi(0), i, 0, lt, 64 )
    ' set the sign
    if ai(0) = bi(0) then
        bi(0) = 0
    else
        bi(0) = &hffff
    endif
    emovo( @bi(0), c )
end sub

/'
;   Multiply.
;
;   unsigned short a[NE], b[NE], c[NE];
;   emul( a, b, c );   c = b * a
'/
sub emul(a as ushort ptr, b as ushort ptr, c as ushort ptr)
    dim as ushort ai(NI), bi(NI)
    dim as integer i, j
    dim as long lt, lta, ltb

    ' NaN times anything is the same NaN.
    if eisnan(a) then
        emov(a,c)
        return
    endif
    if eisnan(b) then
        emov(b,c)
        return
    endif
    ' Zero times infinity is a NaN.
    if (eisinf(a) and (ecmp(b,@ezero(0)) = 0)) or (eisinf(b) and (ecmp(a,@ezero(0)) = 0)) then
        mtherr( "emul", DOMAIN )
        enan( c, NBITS )
        return
    endif
    'Infinity times anything else is infinity.
    if eisinf(a) or eisinf(b) then
        if eisneg(a) xor eisneg(b) then
            *(c+(NE-1)) = &h8000
        else
            *(c+(NE-1)) = 0
        endif
        einfin(c)
        return
    endif
    emovi( a, @ai(0) )
    emovi( b, @bi(0) )
    lta = ai(EXPONENT)
    ltb = bi(EXPONENT)
    if ai(EXPONENT) = 0 then
        for i=1 to NI-2
            if ai(i) <> 0 then
                lta -= enormlz( @ai(0) )
                goto mnzer1
            endif
        next
        eclear(c)
        return
    endif
    mnzer1:

    if bi(EXPONENT) = 0 then
        for i=1 to NI-2
            if bi(i) <> 0 then
                ltb -= enormlz( @bi(0) )
                goto mnzer2
            endif
        next
        eclear(c)
        return
    endif
    mnzer2:

    ' Multiply significands
    j = emulm( @ai(0), @bi(0) )
    ' calculate exponent
    lt = lta + ltb - (EXONE - 1)
    emdnorm( @bi(0), j, 0, lt, 64 )
    ' calculate sign of product
    if ai(0) = bi(0) then
        bi(0) = 0
    else
        bi(0) = &hffff
    endif
    emovo( @bi(0), c )
end sub


/'
; Convert IEEE double precision to e type
;   double d;
;   unsigned short x[N+2];
;   e53toe( &d, x );
'/

sub e64toe(pe as ushort ptr, y as ushort ptr)
    dim as ushort yy(NI)
    dim as ushort ptr p, q, e1
    dim as integer i

    e1 = pe
    p = @yy(0)
    for i=0 to NE-6
        *p = 0
        p += 1
    next

    for i=0 to 4
        *p = *e1
        p += 1
        e1 += 1
    next

    p = @yy(0)
    q = y

    if *p = &h7fff then
        for i=0 to 3
            if pe[i] <> 0 then
                enan( y, NBITS )
                return
            endif
        next

        eclear( y )
        einfin( y )
        if *p and &h8000 then eneg(y)
        return
    endif

    for i=0 to NE-1
        *q = *p
        q += 1
        p += 1
    next
end sub

' move out internal format to ieee long double
sub toe64(a as ushort ptr, b as ushort ptr)
    dim as ushort ptr p, q
    dim as ushort i


    if eiisnan(a) then
        enan( b, 64 )
        return
    endif

    p = a
    q = b + 4 ' point to output exponent

    ' NOTE: if data type is 96 bits wide, clear the last word here.
    *(q+1)= 0

    ' combine sign and exponent
    i = *p
    p += 1

    if i then
        *q = *p or &h8000
        q -= 1
        p += 1
    else
        *q = *p
        q -= 1
        p += 1
    endif

    ' skip over guard word
    p += 1
    ' move the significand

    for i=0 to 3
        *q = *p
        q -= 1
        p += 1
    next
end sub

sub etoe64(x as ushort ptr, ex as ushort ptr)
    dim as ushort xi(NI-1)
    dim as long expn
    dim as integer rndsav

    if eisnan(x) then
        enan( ex, 64 )
        return
    endif
    emovi( x, @xi(0) )
    expn = cast(long, xi(EXPONENT)) /' adjust exponent for offset '/
    if eisinf(x) then
        goto nonorm
    endif
    /' round off to nearest or even '/
    rndsav = rndprc
    rndprc = 64
    emdnorm( @xi(0), 0, 0, expn, 64 )
    rndprc = rndsav
    nonorm:
    toe64( @xi(0), ex )
end sub

/'
; Convert IEEE double precision to e type
;   double d;
;   unsigned short x[N+2];
;   e53toe( &d, x );
'/
sub e53toe(pe as ushort ptr, y as ushort ptr)
    dim as ushort r
    dim as ushort ptr p, ex
    dim as ushort yy(NI-1)
    dim as integer denorm, k

    ex = pe
    denorm = 0   /' flag if denormalized number '/
    ecleaz(@yy(0))
    ex += 3
    r = *ex
    yy(0) = 0
    if (r and &h8000) then
        yy(0) = &hffff
    end if
    yy(MANTISSA) = (r and &h0f) or &h10
    r and= not(&h800f)   /' strip sign and 4 significand bits '/
    if r = &h7ff0 then
        if ((pe[3] and &hf) <> 0) or (pe[2] <> 0) or (pe[1] <> 0) or (pe[0] <> 0) then
            enan( y, NBITS )
            return
        end if
        eclear( y )
        einfin( y )
        if yy(0) then
            eneg(y)
        end if
        return
    end if
    r shr= 4
    /' If zero exponent, then the significand is denormalized.
     * So, take back the understood high significand bit. '/
    if r = 0 then
        denorm = 1
        yy(MANTISSA) and= not(&h10)
    end if
    r += EXONE - &h3ff
    yy(EXPONENT) = r
    p = @yy(MANTISSA+1)
    ex -= 1
    *p = *ex
    p += 1
    ex -= 1
    *p = *ex
    p += 1
    ex -= 1
    *p = *ex
    p += 1

    eshift( @yy(0), -5 )
    if denorm then
        /' if zero exponent, then normalize the significand '/
        if (k = enormlz(@yy(0))) > NBITS then
            ecleazs(@yy(0))
        else
            yy(EXPONENT) -= cast(ushort, (k-1))
        end if
    end if
    emovo( @yy(0), y )
end sub


/'
; e type to IEEE double precision
;   double d;
;   unsigned short x[NE];
;   etoe53( x, &d );
'/


sub toe53(x as ushort ptr, y as ushort ptr)
    dim as ushort i
    dim as ushort ptr p


    if eiisnan(x) then
        enan( y, 53 )
        return
    endif
    p = @x[0]
    y += 3
    *y = 0   ' output high order
    p += 1
    if *(p-1) then *y = &h8000   ' output sign bit

    i = *p
    p += 1
    if i >= cast(uinteger, 2047) then
        ' Saturate at largest number less than infinity.

        *y or= &h7ff0
        y -= 1
        *y = 0
        y -= 1
        *y = 0
        y -= 1
        *y = 0

        return
    endif
    if i = 0 then
        eshift( x, 4 )
    else
        i shl= 4
        eshift( x, 5 )
    endif
    i or= *p and cast(ushort, &h0f)   ' *p = xi[MANTISSA]
    p += 1
    *y or= cast(ushort, i) ' high order output already has sign bit set
    y -= 1
    *y = *p
    p += 1
    y -= 1
    *y = *p
    p += 1
    y -= 1
    *y = *p
end sub

sub etoe53(x as ushort ptr, ex as ushort ptr)
    dim as ushort xi(NI-1)
    dim as long expn
    dim as integer rndsav

    if eisnan(x) then
        enan( ex, 53 )
        return
    endif
    emovi( x, @xi(0) )
    expn = cast(long, xi(EXPONENT)) - (EXONE - &h3ff) /' adjust exponent for offsets '/
    if eisinf(x) then
        goto nonorm
    endif
    /' round off to nearest or even '/
    rndsav = rndprc
    rndprc = 53
    emdnorm( @xi(0), 0, 0, expn, 64 )
    rndprc = rndsav
    nonorm:
    toe53( @xi(0), ex )
end sub

/'
 Convert IEEE single precision to e type
   float d
   unsigned short x[N+2]
   dtox( &d, x )
'/

sub e24toe(pe as ushort ptr, y as ushort ptr)
dim as ushort r
dim as ushort ptr p, e
dim as ushort yy(NI-1)
dim as integer denorm, k

e = pe
denorm = 0   /' flag if denormalized number '/
ecleaz(@yy(0))

e += 1

r = *e
yy(0) = 0
if r and &h8000 then
   yy(0) = &hffff
end if
yy(MANTISSA) = (r and &h7f) or &h80
r and= not(&h807f)   /' strip sign and 7 significand bits '/
if r = &h7f80 then
   if ((pe[1] and &h7f) <> 0) or (pe[0] <> 0) then
      enan( y, NBITS )
      return
    end if
   eclear( y )
   einfin( y )
   if yy(0) then
      eneg(y)
    end if
   return
end if
r shr= 7
/' If zero exponent, then the significand is denormalized.
 * So, take back the understood high significand bit. '/
if r = 0 then
   denorm = 1
   yy(MANTISSA) and= not(&h80)
end if
r += EXONE - &h7f
yy(EXPONENT) = r
p = @yy(MANTISSA+1)
e -= 1
*p = *e
p += 1
eshift( @yy(0), -8 )
if denorm then
   /' if zero exponent, then normalize the significand '/
    k = enormlz(@yy(0))
   if k > NBITS then
      ecleazs(@yy(0))
   else
      yy(EXPONENT) -= cast(ushort, (k-1))
    end if
end if
emovo( @yy(0), y )
end sub

sub toe24(x as ushort ptr, y as ushort ptr)
dim as ushort i
dim as ushort ptr p

if eiisnan(x) then
   enan( y, 24 )
   return
end if
p = @x[0]
y += 1
*y = 0   /' output high order '/
p += 1
if *(p-1) then
   *y = &h8000   /' output sign bit '/
end if
i = *p
p += 1
if i >= 255 then
   /' Saturate at largest number less than infinity. '/
   *y or= cast(ushort, &h7f80)
   *(y-1) = 0
   *y = 0
   return
end if
if i = 0 then
   eshift( x, 7 )
else
   i shl= 7
   eshift( x, 8 )
end if
i or= *p and cast(ushort, &h7f)   /' *p = xi[M] '/
p += 1
*y or= i   /' high order output already has sign bit set '/
*(y-1) = *p
end sub

/'
 e type to IEEE single precision
   float d;
   unsigned short x[N+2];
   xtod( x, &d );
'/
sub etoe24(x as ushort ptr, e as ushort ptr)
dim as long exp1
dim as ushort xi(NI-1)
dim as integer rndsav

if eisnan(x) then
   enan( e, 24 )
   return
end if

emovi( x, @xi(0) )
exp1 = cast(long, xi(EXPONENT)) - (EXONE - &h7f) /' adjust exponent for offsets '/

if eisinf(x) then
   goto nonorm
end if
/' round off to nearest or even '/
rndsav = rndprc
rndprc = 24
emdnorm( @xi(0), 0, 0, exp1, 64 )
rndprc = rndsav
nonorm:
toe24( @xi(0), e )
end sub

/'
 convert long (32-bit) integer to e type

   long l
   unsigned short x[NE]
   ltoe( &l, x )
 note &l is the memory address of l
'/

sub ltoe(lp as long ptr, y as ushort ptr)
    dim as ushort yi(NI-1)
    dim as ulong ll
    dim as integer k

    ecleaz( @yi(0) )
    if *lp < 0 then
        ll =  cast(ulong, (-(*lp)) ) /' make it positive '/
        yi(0) = &hffff /' put correct sign in the e type number '/
    else
        ll = cast(ulong, *lp )
    end if
    /' move the long integer to yi significand area '/
    if sizeof(long) = 8 then
        yi(MANTISSA) = cast(ushort, (ll shr (LONGBITS - 16)))
        yi(MANTISSA + 1) = cast(ushort, (ll shr (LONGBITS - 32)))
        yi(MANTISSA + 2) = cast(ushort, (ll shr 16))
        yi(MANTISSA + 3) = cast(ushort, ll)
        yi(EXPONENT) = EXONE + 47 /' exponent if normalize shift count were 0 '/
    else
        yi(MANTISSA) = cast(ushort, (ll shr 16))
        yi(MANTISSA+1) = cast(ushort, ll)
        yi(EXPONENT) = EXONE + 15 /' exponent if normalize shift count were 0 '/
    end if
    k = enormlz( @yi(0) )
    if k > NBITS then /' normalize the significand '/
        ecleaz( @yi(0) )   /' it was zero '/
    else
        yi(EXPONENT) -= cast(ushort, k) /' subtract shift count from exponent '/
    end if   
    emovo( @yi(0), y )   /' output the answer '/
end sub

/'
 convert longint (64-bit) integer to e type

   longint l
   unsigned short x[NE]
   ltoe( &l, x )
 note &l is the memory address of l
'/

sub lltoe(lp as longint ptr, y as ushort ptr)
    dim as ushort yi(NI-1)
    dim as ulongint ll
    dim as integer k

    ecleaz( @yi(0) )
    if *lp < 0 then
        ll =  cast(ulongint, (-(*lp)) ) /' make it positive '/
        yi(0) = &hffff /' put correct sign in the e type number '/
    else
        ll = cast(ulongint, *lp )
    end if
    /' move the long integer to yi significand area '/
    if sizeof(longint) = 8 then
        yi(MANTISSA) = cast(ushort, (ll shr (LLONGBITS - 16)))
        yi(MANTISSA + 1) = cast(ushort, (ll shr (LLONGBITS - 32)))
        yi(MANTISSA + 2) = cast(ushort, (ll shr 16))
        yi(MANTISSA + 3) = cast(ushort, ll)
        yi(EXPONENT) = EXONE + 47 /' exponent if normalize shift count were 0 '/
    else
        yi(MANTISSA) = cast(ushort, (ll shr 16))
        yi(MANTISSA+1) = cast(ushort, ll)
        yi(EXPONENT) = EXONE + 15 /' exponent if normalize shift count were 0 '/
    end if
    k = enormlz( @yi(0) )
    if k > NBITS then /' normalize the significand '/
        ecleaz( @yi(0) )   /' it was zero '/
    else
        yi(EXPONENT) -= cast(ushort, k) /' subtract shift count from exponent '/
    end if   
    emovo( @yi(0), y )   /' output the answer '/
end sub

srvaldez
Posts: 2525
Joined: Sep 25, 2005 21:54

Postby srvaldez » Jun 20, 2011 1:14

part 3 of 3

Code: Select all


/'
 convert unsigned long (32-bit) integer to e type

   unsigned long l
   unsigned short x(NE)
   ltox( &l, x )
 note &l is the memory address of l
'/

sub ultoe(lp as ulong ptr, y as ushort ptr)
                   /' lp is the memory address of a long integer '/
                     /' y is the address of a short '/
    dim as ushort yi(NI-1)
    dim as ulong ll
    dim as integer k

    ecleaz( @yi(0) )
    ll = *lp

    /' move the long integer to ayi significand area '/
    if sizeof(long) = 8 then
        yi(MANTISSA) = cast(ushort, (ll shr (LONGBITS - 16)))
        yi(MANTISSA + 1) = cast(ushort, (ll shr (LONGBITS - 32)))
        yi(MANTISSA + 2) = cast(ushort, (ll shr 16))
        yi(MANTISSA + 3) = cast(ushort, ll)
        yi(EXPONENT) = EXONE + 47 /' exponent if normalize shift count were 0 '/
    else
        yi(MANTISSA) = cast(ushort, (ll shr 16))
        yi(MANTISSA+1) = cast(ushort, ll)
        yi(EXPONENT) = EXONE + 15 /' exponent if normalize shift count were 0 '/
    end if
    k = enormlz( @yi(0) )
    if k > NBITS then /' normalize the significand '/
        ecleaz( @yi(0) )   /' it was zero '/
    else
        yi(EXPONENT) -= cast(ushort, k) /' subtract shift count from exponent '/
    end if
    emovo( @yi(0), y )   /' output the answer '/
end sub

/'
 convert unsigned longint (64-bit) integer to e type

   unsigned longint l
   unsigned short x(NE)
   ltox( &l, x )
 note &l is the memory address of l
'/

sub ulltoe(lp as ulongint ptr, y as ushort ptr)
                   /' lp is the memory address of a long integer '/
                     /' y is the address of a short '/
    dim as ushort yi(NI-1)
    dim as ulongint ll
    dim as integer k

    ecleaz( @yi(0) )
    ll = *lp
    /' move the long integer to ayi significand area '/
    if sizeof(longint) = 8 then
        yi(MANTISSA) = cast(ushort, (ll shr (LLONGBITS - 16)))
        yi(MANTISSA + 1) = cast(ushort, (ll shr (LLONGBITS - 32)))
        yi(MANTISSA + 2) = cast(ushort, (ll shr 16))
        yi(MANTISSA + 3) = cast(ushort, ll)
        yi(EXPONENT) = EXONE + 47 /' exponent if normalize shift count were 0 '/
    else
        yi(MANTISSA) = cast(ushort, (ll shr 16))
        yi(MANTISSA+1) = cast(ushort, ll)
        yi(EXPONENT) = EXONE + 15 /' exponent if normalize shift count were 0 '/
    end if
    k = enormlz( @yi(0) )
    if k > NBITS then /' normalize the significand '/
        ecleaz( @yi(0) )   /' it was zero '/
    else
        yi(EXPONENT) -= cast(ushort, k) /' subtract shift count from exponent '/
    end if
    emovo( @yi(0), y )   /' output the answer '/
end sub

sub eiremain(den as ushort ptr, num as ushort ptr)
    dim as long ld, ln
    dim as ushort j

    ld = den[EXPONENT]
    ld -= enormlz( den )
    ln = num[EXPONENT]
    ln -= enormlz( num )
    ecleaz( @equot(0) )
    while( ln >= ld )
        if ( ecmpm(den,num) <= 0 ) then
            esubm(den, num)
            j = 1
        else
            j = 0
        end if
        eshup1(@equot(0))
        equot(NI-1) or= j
        eshup1(num)
        ln -= 1
    wend
    emdnorm( num, 0, 0, ln, 0 )
end sub

sub efloor(x as ushort ptr, y as ushort ptr)

    dim as ushort ptr p
    dim as integer e2, expon, i
    dim as ushort f(NE-1)

    static bmask(16) as ushort = {&hffff, &hfffe, &hfffc, &hfff8,_
              &hfff0, &hffe0, &hffc0, &hff80, &hff00, &hfe00, &hfc00, _
              &hf800, &hf000, &he000, &hc000, &h8000, &h0000}
             
    emov( x, @f(0) ) /' leave in external format '/
    expon = cast(integer, f(NE-1))
    e2 = (expon and &h7fff) - (EXONE - 1)
    if( e2 <= 0 ) then
        eclear(y)
        goto isitneg
    endif
    /' number of bits to clear out '/
    e2 = NBITS - e2
    emov( @f(0), y )
    if ( e2 <= 0 ) then
        return
    endif
    p = y
    while( e2 >= 16 )
        *p = 0
        p += 1
        e2 -= 16
    wend

    /' clear the remaining bits '/
    *p and= bmask(e2)
    /' truncate negatives toward minus infinity '/
    isitneg:

    if( cast(ushort,expon) and cast(ushort, &h8000 )) then
        p = y
        for i=0 to NE-2
            if( f(i) <> *p ) then
                esub( @eone(0), y, y )
                exit sub
            endif
            p += 1
        next
    endif
end sub



/' Convert e type number to decimal format ASCII string.
 * The constants are for 64 bit precision.
 '/

#define NTEN 12
#define MAXP 4096


dim shared etens(NTEN + 1, NE-1) as ushort = _
{ _
  {&h6576, &h4a92, &h804a, &h153f, _
   &hc94c, &h979a, &h8a20, &h5202, &hc460, &h7525}, _   /' 10**4096 '/
  {&h6a32, &hce52, &h329a, &h28ce, _
   &ha74d, &h5de4, &hc53d, &h3b5d, &h9e8b, &h5a92}, _   /' 10**2048 '/
  {&h526c, &h50ce, &hf18b, &h3d28, _
   &h650d, &h0c17, &h8175, &h7586, &hc976, &h4d48}, _
  {&h9c66, &h58f8, &hbc50, &h5c54, _
   &hcc65, &h91c6, &ha60e, &ha0ae, &he319, &h46a3}, _
  {&h851e, &heab7, &h98fe, &h901b, _
   &hddbb, &hde8d, &h9df9, &hebfb, &haa7e, &h4351}, _
  {&h0235, &h0137, &h36b1, &h336c, _
   &hc66f, &h8cdf, &h80e9, &h47c9, &h93ba, &h41a8}, _
  {&h50f8, &h25fb, &hc76b, &h6b71, _
   &h3cbf, &ha6d5, &hffcf, &h1f49, &hc278, &h40d3}, _
  {&h0000, &h0000, &h0000, &h0000, _
   &hf020, &hb59d, &h2b70, &hada8, &h9dc5, &h4069}, _
  {&h0000, &h0000, &h0000, &h0000, _
   &h0000, &h0000, &h0400, &hc9bf, &h8e1b, &h4034}, _
  {&h0000, &h0000, &h0000, &h0000, _
   &h0000, &h0000, &h0000, &h2000, &hbebc, &h4019}, _
  {&h0000, &h0000, &h0000, &h0000, _
   &h0000, &h0000, &h0000, &h0000, &h9c40, &h400c}, _
  {&h0000, &h0000, &h0000, &h0000, _
   &h0000, &h0000, &h0000, &h0000, &hc800, &h4005}, _
  {&h0000, &h0000, &h0000, &h0000, _
   &h0000, &h0000, &h0000, &h0000, &ha000, &h4002} _   /' 10**1 '/
}

dim shared emtens(NTEN, NE-1) as ushort = _
{ _
  {&h2030, &hcffc, &ha1c3, &h8123, _
   &h2de3, &h9fde, &hd2ce, &h04c8, &ha6dd, &h0ad8}, _   /' 10**-4096 '/
  {&h8264, &hd2cb, &hf2ea, &h12d4, _
   &h4925, &h2de4, &h3436, &h534f, &hceae, &h256b}, _   /' 10**-2048 '/
  {&hf53f, &hf698, &h6bd3, &h0158, _
   &h87a6, &hc0bd, &hda57, &h82a5, &ha2a6, &h32b5}, _
  {&he731, &h04d4, &he3f2, &hd332, _
   &h7132, &hd21c, &hdb23, &hee32, &h9049, &h395a}, _
  {&ha23e, &h5308, &hfefb, &h1155, _
   &hfa91, &h1939, &h637a, &h4325, &hc031, &h3cac}, _
  {&he26d, &hdbde, &hd05d, &hb3f6, _
   &hac7c, &he4a0, &h64bc, &h467c, &hddd0, &h3e55}, _
  {&h2a20, &h6224, &h47b3, &h98d7, _
   &h3f23, &he9a5, &ha539, &hea27, &ha87f, &h3f2a}, _
  {&h0b5b, &h4af2, &ha581, &h18ed, _
   &h67de, &h94ba, &h4539, &h1ead, &hcfb1, &h3f94}, _
  {&hbf71, &ha9b3, &h7989, &hbe68, _
   &h4c2e, &he15b, &hc44d, &h94be, &he695, &h3fc9}, _
  {&h3d4d, &h7c3d, &h36ba, &h0d2b, _
   &hfdc2, &hcefc, &h8461, &h7711, &habcc, &h3fe4}, _
  {&hc155, &ha4a8, &h404e, &h6113, _
   &hd3c3, &h652b, &he219, &h1758, &hd1b7, &h3ff1}, _
  {&hd70a, &h70a3, &h0a3d, &ha3d7, _
   &h3d70, &hd70a, &h70a3, &h0a3d, &ha3d7, &h3ff8}, _
  {&hcccd, &hcccc, &hcccc, &hcccc, _
   &hcccc, &hcccc, &hcccc, &hcccc, &hcccc, &h3ffb} _   /' 10**-1 '/
}

dim shared outformat as integer = -1

function etoasc(byval x as ushort ptr, byval ndigits as integer) as string

dim digit as long
dim as ushort y(NI-1), t(NI-1), u(NI-1), w(NI-1)
dim as ushort ptr p, r, ten
dim as ushort sign
dim as integer i, j, k, expon, rndsav, ndigs, outexpon
dim as integer s, ss
dim as string*256 strin
dim as ushort m
dim as string oute

ndigs = ndigits
rndsav = rndprc

if eisnan(x) then
   strin = " NaN "
   expon = 9999
   goto bxit
end if

rndprc = NBITS      /' set to full precision '/
emov( x, @y(0) ) /' retain external format '/

if y(NE-1) and &h8000 then
   sign = &hffff
   y(NE-1) and= &h7fff
else
   sign = 0
end if
expon = 0
ten = @etens(NTEN,0)
emov( @eone(0), @t(0) )
/' Test for zero exponent '/
if y(NE-1) = 0 then
   for k=0 to NE-2
      if y(k) <> 0 then goto tnzro /' denormalized number '/
   next
   goto isone /' legal all zeros '/
end if
tnzro:

/' Test for infinity.
 '/
if y(NE-1) = &h7fff then
   if sign then
      strin = " -Infinity "
   else
      strin = " Infinity "
    end if
   expon = 9999
   goto bxit
end if

/' Test for exponent nonzero but significand denormalized.
 * This is an error condition.
 '/
if (y(NE-1) <> 0) and ((y(NE-2) and &h8000) = 0) then
   mtherr( "etoasc", DOMAIN )
   strin = "NaN"
   expon = 9999
   goto bxit
end if

/' Compare to 1.0 '/
i = ecmp( @eone(0), @y(0) )

if i = 0 then goto isone

if i < 0 then
/' Number is greater than 1 '/
/' Convert significand to an integer and strip trailing decimal zeros. '/
   emov( @y(0), @u(0) )
   u(NE-1) = EXONE + NBITS - 1
   p = @etens(NTEN-4,0)
   m = 16
do
   ediv( p, @u(0), @t(0) )
   efloor( @t(0), @w(0) )
   for j=0 to NE-2
      if t(j) <> w(j) then goto noint
   next
   emov( @t(0), @u(0) )
   expon += cast(integer, m)
noint:
   p += NE
   m shr= 1
loop while( m <> 0 )
/' Rescale from integer significand '/
   u(NE-1) += y(NE-1) - cast(uinteger ,EXONE + NBITS - 1)
   emov( @u(0), @y(0) )
/' Find power of 10 '/
   emov( @eone(0), @t(0) )
   m = MAXP
   p = @etens(0,0)
   while( ecmp( ten, @u(0) ) <= 0 )
      if ecmp( p, @u(0) ) <= 0 then
         ediv( p, @u(0), @u(0) )
         emul( p, @t(0), @t(0) )
         expon += cast(integer, m)
      end if
      m shr= 1
      if m = 0 then exit while
      p += NE
   wend
else
/' Number is less than 1.0 '/
/' Pad significand with trailing decimal zeros. '/
   if y(NE-1) = 0 then
      while( (y(NE-2) and &h8000) = 0 )
         emul( ten, @y(0), @y(0) )
         expon -= 1
      wend
   else
      emovi( @y(0), @w(0) )
      for i=0 to NDEC
         if (w(NI-1) and &h7) <> 0 then exit for
/' multiply by 10 '/
         emovz( @w(0), @u(0) )
         eshdn1( @u(0) )
         eshdn1( @u(0) )
         eaddm( @w(0), @u(0) )
         u(1) += 3
         while( u(2) <> 0 )
            eshdn1(@u(0))
            u(1) += 1
         wend
         if u(NI-1) <> 0 then exit for
         if eone(NE-1) <= u(1) then exit for
         emovz( @u(0), @w(0) )
         expon -= 1
      next
      emovo( @w(0), @y(0) )
   end if
   k = -MAXP
   p = @emtens(0,0)
   r = @etens(0,0)
   emov( @y(0), @w(0) )
   emov( @eone(0), @t(0) )
   while( ecmp( @eone(0), @w(0) ) > 0 )
      if ecmp( p, @w(0) ) >= 0 then
         emul( r, @w(0), @w(0) )
         emul( r, @t(0), @t(0) )
         expon += k
      end if
      k \= 2
      if k = 0 then exit while
      p += NE
      r += NE
   wend
   ediv( @t(0), @eone(0), @t(0) )
end if

isone:
/' Find the first (leading) digit. '/
emovi( @t(0), @w(0) )
emovz( @w(0), @t(0) )
emovi( @y(0), @w(0) )
emovz( @w(0), @y(0) )
eiremain( @t(0), @y(0) )
digit = equot(NI-1)

while( (digit = 0) and (ecmp(@y(0),@ezero(0)) <> 0) )
   eshup1( @y(0) )
   emovz( @y(0), @u(0) )
   eshup1( @u(0) )
   eshup1( @u(0) )
   eaddm( @u(0), @y(0) )
   eiremain( @t(0), @y(0) )
   digit = equot(NI-1)
   expon -= 1
wend

s = 0's = strin
if sign then
   strin[s] = asc("-")
    s += 1
else
   strin[s] = asc(" ")
    s += 1
end if
/' Examine number of digits requested by caller. '/
if outformat = 3 then ndigs += expon
/'
else if ndigs < 0 )
        ndigs = 0
'/
if ndigs > NDEC then ndigs = NDEC
if digit = 10 then
   strin[s] = asc("1")
    s += 1
   strin[s] = asc(".")
    s += 1
   if ndigs > 0 then
      strin[s] = asc("0")
        s += 1
      ndigs -= 1
   end if
   expon += 1
   if ndigs < 0 then
           ss = s
           goto doexp
   end if
else
   strin[s] = digit+48 'asc(48) = '0'
    s += 1
   strin[s] = asc(".")
    s += 1
end if
/' Generate digits after the decimal point. '/
for k=0 to ndigs
/' multiply current number by 10, without normalizing '/
   eshup1( @y(0) )
   emovz( @y(0), @u(0) )
   eshup1( @u(0) )
   eshup1( @u(0) )
   eaddm( @u(0), @y(0) )
   eiremain( @t(0), @y(0) )
   strin[s] = equot(NI-1)+48 'asc(48) = '0'
    s += 1
next

digit = equot(NI-1)
s -= 1
ss=s
/' round off the ASCII string '/
if digit > 4 then
/' Test for critical rounding case in ASCII output. '/
   if digit = 5 then
      emovo( @y(0), @t(0) )
      if ecmp(@t(0),@ezero(0)) <> 0 then   goto roun   /' round to nearest '/
      if (strin[s-1] and 1) = 0 then goto doexp   /' round to even '/
   end if
/' Round up and propagate carry-outs '/
roun:
   s -= 1
   k = strin[s] and &h7f
/' Carry out to most significant digit? '/
   if ndigs < 0 then
           /' This will print like "1E-6". '/
      strin[s] = asc("1")
      expon += 1
      goto doexp
   elseif k = asc("." ) then
      s -= 1
      k = strin[s]
      k += 1
      strin[s] = k
/' Most significant digit carries to 10? '/
      if k > asc("9") then
         expon += 1
         strin[s] = asc("1")
      end if
      goto doexp
   end if
/' Round up and carry out from less significant digits '/
   k += 1
   strin[s] = k
   if k > asc("9" ) then
      strin[s] = asc("0")
      goto roun
   end if
end if
doexp:
/'
if expon >= 0 )
   sprintf( ss, "e+%d", expon )
else
   sprintf( ss, "e%d", expon )
'/
'   sprintf( ss, "E%d", expon )
bxit:

oute = left(strin,ndigs+3)+"e"+str(expon)
rndprc = rndsav
outexpon =  expon

return oute

end function

sub asctoeg(byref ss1 as string, y as ushort ptr, oprec as integer)
dim as ushort yy(NI-1), xt(NI-1), tt(NI-1)
dim as integer esign, decflg, sgnflg, nexp, exp1, prec, lost
dim as integer k, trail, c, rndsav
dim as long lexp
dim as ushort nsign
dim as ushort ptr p
dim as string*256 lstr, ss
dim as integer sp, s

ss = ss1
c = len (ss1) + 2

/' Copy the input string. '/

'lstr = (char *) alloca (c)
s = 0
lenldstr = 0

while ss[s] = asc(" ") /' skip leading spaces '/
    s += 1
    lenldstr += 1
wend
sp = 0
for k=0 to c-1
    lstr[sp] = ss[s]
    sp += 1
    s  += 1
   if lstr[sp-1] = 0 then exit for
next

lstr[sp] = 0
s = 0

rndsav = rndprc
rndprc = NBITS /' Set to full precision '/
lost = 0
nsign = 0
decflg = 0
sgnflg = 0
nexp = 0
exp1 = 0
prec = 0
ecleaz( @yy(0) )
trail = 0

nxtcom:
k = lstr[s] - asc("0")

if (k >= 0) and (k <= 9) then
/' Ignore leading zeros '/
   if (prec = 0) and (decflg = 0) and (k = 0) then   goto donchr
/' Identify and strip trailing zeros after the decimal point. '/
   if (trail = 0) and (decflg <> 0) then
      sp = s
      while (lstr[sp] >= asc("0")) and (lstr[sp] <= asc("9"))
         sp += 1
        wend
/' Check for syntax error '/
      c = lstr[sp] and &h7f
      if (c <> asc("e")) and (c <> asc("E")) and (c <> 0)_
         and (c <> 10) and (c <> 13) and (c <> asc(" "))_
         and (c <> asc(","))_
         then goto errors
      sp -= 1
      while lstr[sp] = asc("0")
            lstr[sp] = asc("z")
         sp -= 1
        wend
      trail = 1
      if lstr[s] = asc("z") then goto donchr
   end if
/' If enough digits were given to more than fill up the yy register,
 * continuing until overflow into the high guard word yy[2]
 * guarantees that there will be a roundoff bit at the top
 * of the low guard word after normalization.
 '/
   if yy(2) = 0 then
      if decflg then nexp += 1 /' count digits after decimal point '/
      eshup1( @yy(0) )   /' multiply current number by 10 '/
      emovz( @yy(0), @xt(0) )
      eshup1( @xt(0) )
      eshup1( @xt(0) )
      eaddm( @xt(0), @yy(0) )
      ecleaz( @xt(0) )
      xt(NI-2) = cast(ushort, k)
      eaddm( @xt(0), @yy(0) )
   else
      /' Mark any lost non-zero digit.  '/
      lost or= k
      /' Count lost digits before the decimal point.  '/
      if (decflg = 0) then nexp -= 1
   end if
   prec += 1
   goto donchr
end if

select case as const lstr[s]
   case asc("z")
      exit select
   case asc("e")
        goto expnt
    case asc("E")
      goto expnt
   case asc(".")   /' decimal point '/
      if decflg then goto errors
      decflg += 1
      exit select
   case asc("-")
      nsign = &hffff
      if sgnflg then goto errors
      sgnflg += 1
      exit select
   case asc("+")
      if sgnflg then goto errors
      sgnflg += 1
      exit select
   case asc(",")
        goto daldone
   case asc(" ")
        goto daldone
   case 0
        goto daldone
   case 10
        goto daldone
   case 13
      goto daldone
   case asc("i")
        goto infinite
   case asc("I")
      goto infinite
   case else
   errors:
    enan( @yy(0), NI*16 )
    goto aexit
end select
donchr:
s += 1
goto nxtcom

/' Exponent interpretation '/
expnt:

esign = 1
exp1 = 0
s += 1
/' check for + or - '/
if lstr[s] = asc("-") then
   esign = -1
   s += 1
end if
if lstr[s] = asc("+") then s += 1
while (lstr[s] >= asc("0")) and (lstr[s] <= asc("9"))
   exp1 *= 10
   exp1 += lstr[s] - asc("0")
    s += 1
   if (exp1 > 4977) then
      if (esign < 0) then
         goto zero
      else
         goto infinite
        end if
   end if
wend
if esign < 0 then exp1 = -exp1
if exp1 > 4932 then
infinite:
   ecleaz(@yy(0))
   yy(EXPONENT) = &h7fff  /' infinity '/
   goto aexit
end if
if exp1 < -4977 then
zero:
   ecleaz(@yy(0))
   goto aexit
end if

daldone:
nexp = exp1 - nexp
/' Pad trailing zeros to minimize power of 10, per IEEE spec. '/
while (nexp > 0) and (yy(2) = 0)
   emovz( @yy(0), @xt(0) )
   eshup1( @xt(0) )
   eshup1( @xt(0) )
   eaddm( @yy(0), @xt(0) )
   eshup1( @xt(0) )
   if xt(2) <> 0 then exit while
   nexp -= 1
   emovz( @xt(0), @yy(0) )
wend
k = enormlz(@yy(0))
if k > NBITS then
   ecleaz(@yy(0))
   goto aexit
end if
lexp = (EXONE - 1 + NBITS) - k
emdnorm( @yy(0), lost, 0, lexp, 64 )
/' convert to external format '/


/' Multiply by 10**nexp.  If precision is 64 bits,
 * the maximum relative error incurred in forming 10**n
 * for 0 <= n <= 324 is 8.2e-20, at 10**180.
 * For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
 * For 0 >= n >= -999, it is -1.55e-19 at 10**-435.
 '/
lexp = yy(EXPONENT)
if nexp = 0 then
   k = 0
   goto expdon
end if
esign = 1
if nexp < 0 then
   nexp = -nexp
   esign = -1
   if nexp > 4096 then
      /' Punt.  Can't handle this without 2 divides. '/
      emovi( @etens(0,0), @tt(0) )
      lexp -= tt(EXPONENT)
      k = edivm( @tt(0), @yy(0) )
      lexp += EXONE
      nexp -= 4096
   end if
end if
p = @etens(NTEN,0)
emov( @eone(0), @xt(0) )
exp1 = 1
do
   if exp1 and nexp then
        emul( p, @xt(0), @xt(0) )
    end if
   p -= NE
   exp1 = exp1 + exp1
loop while( exp1 <= MAXP )

emovi( @xt(0), @tt(0) )
if esign < 0 then
   lexp -= tt(EXPONENT)
   k = edivm( @tt(0), @yy(0) )
   lexp += EXONE
else
   lexp += tt(EXPONENT)
   k = emulm( @tt(0), @yy(0) )
   lexp -= EXONE - 1
end if

expdon:

/' Round and convert directly to the destination type '/
if oprec = 53 then
   lexp -= EXONE - &h3ff
elseif oprec = 24 then
   lexp -= EXONE - &h7f
end if
rndprc = oprec
emdnorm( @yy(0), k, 0, lexp, 64 )

aexit:

rndprc = rndsav
yy(0) = nsign
select case oprec
   case 53
      toe53( @yy(0), y )
      exit select
'   case 24
'      toe24( @yy(0), y )
'      exit select
   case 64
      toe64( @yy(0), y )
      exit select
'   case 113
'      toe113( @yy(0), y )
'      exit select
   case NBITS
      emovo( @yy(0), y )
      'exit select
end select
'lenldstr += s - lstr

end sub

'###################################

dim as ushort w, ip(NI-1)
dim as string s = "2.7182818284590452353602874713526624977572470937000"
dim as integer i
dim as longdouble pi
dim as ushort carry(NI-1)
dim as double dd = 1234567
dim as single ss = 3.141592654
dim as ulong ul = 4294967295
dim as long  sl = -2147483648
dim as longint sll = -9223372036854775808
dim as ulongint ull = 18446744073709551615

'Const Wd = 800, Ht = 600
'ScreenRes Wd, Ht

'Width Wd\8, Ht\16 '' Use 8*16 font

'store the 80-bit value of Pi into the variable pi
asm
    fldpi
    fstp tbyte ptr [pi]
end asm
'convert the 80-bit value to extended format
e64toe(cast(ushort ptr, @pi), @ip(0))
'now convert and print the extended value
? "Pi converted from 80-bit to extended precision and then converted to string "
? etoasc(@ip(0), 18)
?
'convert the extended value of Pi to double
etoe53(@epi(0), cast(ushort ptr, @dd))
? "Pi converted from extended precision to double "
? dd
?
'convert double to extended precision
e53toe(cast(ushort ptr, @dd), @ip(0))
? "Pi converted from double to extended precision and then converted to string "
? etoasc(@ip(0), 20)
?
? "convert the signed long value -2147483648 to extended precision and print"
ltoe(@sl, @ip(0))
? etoasc(@ip(0), 20)
?
? "convert the unsigned long value 4294967295 to extended precision and print"
ultoe(@ul, @ip(0))
? etoasc(@ip(0), 20)
?
? "convert the signed longint value -9223372036854775808 to extended precision and print"
lltoe(@sll, @ip(0))
? etoasc(@ip(0), 20)
?
? "convert the unsigned longint value 18446744073709551615 to extended precision and print"
ulltoe(@ull, @ip(0))
? etoasc(@ip(0), 20)
?
? "convert the string 2.7182818 ... to extended precision and back to string"
asctoeg(s, @ip(0), NBITS)
? etoasc(@ip(0), 45)
?
? "convert the string 2.7182818 ... to double precision and print"
asctoeg(s, cast(ushort ptr, @dd), 53)
? dd
?
? "Pi converted from single to extended precision and then converted to string "
e24toe(cast(ushort ptr, @ss), @ip(0))
? etoasc(@ip(0), 45)
?
? "Pi converted from extended precision to single and then printed"
etoe24(@epi(0), cast(ushort ptr, @ss))
? ss
?
Print "Press RETURN to end ";

sleep


Return to “Tips and Tricks”

Who is online

Users browsing this forum: Bing [Bot] and 2 guests