General FreeBASIC programming questions.
srvaldez
Posts: 2850
Joined: Sep 25, 2005 21:54

I am still working on my big-float package and while examining the multiplication routine by Stephan Brunker https://github.com/stephanbrunker/big_integer I noticed that he was assigning a ulongint to a ulong, it turns out that the least significant portion of the ulongint is stored in the ulong so he saves from having to do a mod operation, he uses base 2^32 I think
he's multiplication routine is over 2 times faster than mine partially because of the larger base and the trick I mentioned, I plan on rewriting my bigfloat package using base 2^32 to take advantage of this trick, I also plan on implementing the toom 3 algorithm
dodicat
Posts: 7146
Joined: Jan 10, 2006 20:30
Location: Scotland

Hi srvaldez.
Many years ago in topic circles, Richard, Albert and myself (to a lesser degree) messed about with bignums.
Richard finally used base 1e9 for his method.
Albert used various bases (as those mults I posted a while back).
I stuck to high school methods using string indexing base 10.
My method could never get up to speed, but I stuck to it anyway, I still post bits and pieces now and then.
Maybe you can test your right/left shifts (> fb ulongint) with this, It looks and tests out OK.
I have kept the test numbers smallish to suit 32 bit gas.

Code: Select all

`Function plus(_num1 As String,_num2 As String) As String      Static As integer ADDQmod(19)={48,49,50,51,52,53,54,55,56,57,48,49,50,51,52,53,54,55,56,57}      Static As integer ADDbool(19)={0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1}    Var _flag=0,n_=0    Dim As integer addup=Any,addcarry=Any    #macro finish()    answer=Ltrim(answer,"0")    If _flag=1 Then Swap _num2,_num1    Return answer    #endmacro    If Len(_num2)>Len(_num1) Then         Swap _num2,_num1        _flag=1        Endif        Var diff=Len(_num1)-Len(_num2)        Var answer="0"+_num1        addcarry=0        For n_=Len(_num1)-1 To diff Step -1             addup=_num2[n_-diff]+_num1[n_]-96            answer[n_+1]=ADDQmod(addup+addcarry)            addcarry=ADDbool(addup+addcarry)        Next n_         If addcarry=0 Then             finish()        End If        If n_=-1 Then             answer[0]=addcarry+48            finish()            Endif            For n_=n_ To 0 Step -1                 addup=_num1[n_]-48                answer[n_+1]=ADDQmod(addup+addcarry)                addcarry=ADDbool(addup+addcarry)                If addcarry=0 Then Exit For            Next n_            answer[0]=addcarry+48            finish()        End Function                ' ============== BINARY TO DECIMAL =================================             Function Base10(BinaryNumber As String) As String            Dim As String sum            sum=Left(BinaryNumber,1)            For n As Integer=2 To Len(BinaryNumber)                sum=plus(plus(sum,sum),Chr(BinaryNumber[n-1]))             Next n            Return sum        End Function                ' ==================== DECIMAL TO BINARY ==================================                  Function base2(DecimalNumber As String) As String            Dim As String starter=DecimalNumber,ans,m,b            Dim As Ubyte main,carry,temp',c            Dim As Integer c            #macro reverse(s)                Scope                Var lens=Len(s)                For n As Integer=0 To Int((lens-1)/2):Swap s[n],s[lens-1-n]:Next                End Scope                #endmacro                #macro div2(s,m,c)                carry=0:ans=s                For z As Integer=0 To Len(s)-1                    temp=(s[z]-48+carry)                    main=temp Shr 1                    carry=(temp And 1) Shl 3 +(temp And 1) Shl 1                    ans[z]=main+48                Next z                c= carry\10                m= Ltrim(ans,"0")                #endmacro                Do                    c=0                    div2(starter,m,c)                    b=b+Str(c)                    starter=m                Loop Until m=""                reverse(b)                b=Str(m)+b                Return b            End Function                   Function shiftr(sd As String,k As Long ) As String                  Var b2=base2(sd)                  Var l=Len(b2)                  Dim As String s=String(k,"0")+b2                  var ans=base10(Left(s,l))                  Return iif(ans="","0",ans)       End Function       Function shiftl(sd As String,k As Long ) As String                Var b2=base2(sd)                Dim As String s=b2+String(k,"0")                var ans=base10(s)                Return iif(ans="","0",ans)       End Function              dim as integer a,b       print "loop test . . ."       for n as long=1 to 10000             a=20+rnd*10000             b=1+rnd*15             if a shl b<>val(shiftl(str(a),b)) then print n             if a shr b<>val(shiftr(str(a),b)) then print n             nextPrint Print shiftl("2021",2021)print "done"sleep `
srvaldez
Posts: 2850
Joined: Sep 25, 2005 21:54

Hi dodicat :-)
thank you very much, it's just what I needed
I got started on the 32-bit binary math routines and found out that the bit shifter with the asm code is way too slow so I used fxm suggestion and it's working ok so far, your conversion routines are just what I needed :-)
srvaldez
Posts: 2850
Joined: Sep 25, 2005 21:54

so far I have implemented add, subtract, multiply and divide_by_long, a big hurdle to overcome is the conversion from string to bigfloat and back
I tested the routines by calculating sin(0.5) to 1030 decimal digits verified with the pari/gp calculator
the binary routines were about 38% faster than my decimal routines, I am somewhat disappointed
srvaldez
Posts: 2850
Joined: Sep 25, 2005 21:54

just for fun, here's my try at multi-precision binary floating-point math
the input/output is crude but it's better than nothing, was having trouble figuring what was wrong with the conversion to string until it dawned on me that since I use double to extract the digits the value could easily be out of range

Code: Select all

`Extern "c"   Declare Function __builtin_clz (Byval x As Ulong) As Long   Declare Function __builtin_smull_overflow (Byval As Long, Byval As Long, Byref As Long) As Boolean   Declare Function __builtin_smulll_overflow (Byval As Longint, Byval As Longint, Byref As Longint) As Boolean   Declare Function __builtin_umull_overflow (Byval As Ulong, Byval As Ulong, Byref As Ulong) As Boolean   Declare Function __builtin_umulll_overflow (Byval As Ulongint, Byval As Ulongint, Byref As Ulongint) As Boolean   Declare Function __builtin_uaddll_overflow (Byval As Ulongint, Byval As Ulongint, Byref As Ulongint) As BooleanEnd ExternConst NUMBER_OF_DIGITS = 96Const NUMBER_OF_BITS   =  NUMBER_OF_DIGITS*3.321928094887362#If (NUMBER_OF_BITS Mod 32)<>0   Const NUM_BITS = (NUMBER_OF_BITS\32+1)*32#Else   Const NUM_BITS = NUMBER_OF_BITS#EndifConst NUM_DWORDS    =   NUM_BITS\32Const BIAS          =   1073741824 '2 ^ 30' Error definitions#Define DIVZ_ERR        1                    'Divide by zero#Define EXPO_ERR        2                    'Exponent overflow error#Define EXPU_ERR        3                    'Exponent underflow errorType bigfloat_struct   Declare Constructor ( )   Declare Constructor ( Byref rhs As bigfloat_struct )   Declare Destructor ( )   Declare Operator Let ( Byref rhs As bigfloat_struct )    As Short sign    As Ulong  exponent   As Ulong  Ptr mantissaEnd TypeConstructor bigfloat_struct ( )   If this.mantissa<>0 Then      Print "pointer is not 0"      End(1)   End If   this.mantissa=Callocate((NUM_DWORDS +1), 4)   If this.mantissa=0 Then      Print "unable to allocate memory"      End(4)   End If   this.sign=0   this.exponent=0End ConstructorConstructor bigfloat_struct ( Byref rhs As bigfloat_struct )   If this.mantissa<>0 Then      Print "pointer is not 0"      End(1)   End If   this.mantissa=Callocate((NUM_DWORDS +1), 4)   If this.mantissa=0 Then      Print "unable to allocate memory"      End(4)   End If   this.sign=rhs.sign   this.exponent=rhs.exponent   For i As Long=0 To NUM_DWORDS      this.mantissa[i]=rhs.mantissa[i]   NextEnd ConstructorOperator bigfloat_struct.let ( Byref rhs As bigfloat_struct )   this.sign=rhs.sign   this.exponent=rhs.exponent   For i As Long=0 To NUM_DWORDS      this.mantissa[i]=rhs.mantissa[i]   NextEnd OperatorDestructor bigfloat_struct ( )   If this.mantissa=0 Then      Print "unable to de-allocate memory"      End(-4)   Else      Deallocate(this.mantissa)   End If   this.mantissa=0End DestructorDim Shared As bigfloat_struct Log102Log102.sign=0Log102.exponent=(-2)+BIAS+1If NUM_DWORDS>=0 Then Log102.mantissa[0]=&h13441350If NUM_DWORDS>=1 Then Log102.mantissa[1]=&h9F79FEF3If NUM_DWORDS>=2 Then Log102.mantissa[2]=&h11F12B35If NUM_DWORDS>=3 Then Log102.mantissa[3]=&h816F922FPrivate Function log2_32(Byval value As Ulong) As Long'https://stackoverflow.com/questions/11376288/fast-computing-of-log2-for-64-bit-integers   Static tab32(0 To 31) As Const Long = {0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30, 8, 12, 20, 28, 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31}   value Or= Culng(value Shr 1)   value Or= Culng(value Shr 2)   value Or= Culng(value Shr 4)   value Or= Culng(value Shr 8)   value Or= Culng(value Shr 16)   Return tab32(Culng(Culng(value * &h07C4ACDD) Shr 27))End FunctionPrivate Function log2_64(Byval value As Ulongint) As Long'https://stackoverflow.com/questions/11376288/fast-computing-of-log2-for-64-bit-integers   Static tab64(0 To 63) As Const Long ={63, 0, 58, 1, 59, 47, 53, 2, 60, 39, 48, 27, 54, 33, 42, 3, 61, 51, 37, 40, 49, 18, 28, 20, 55, 30, 34, 11, 43, 14, 22, 4, 62, 57, 46, 52, 38, 26, 32, 41, 50, 36, 17, 19, 29, 10, 13, 21, 56, 45, 25, 31, 35, 16, 9, 12, 44, 24, 15, 8, 23, 7, 6, 5}   value Or= Culngint(value Shr 1)   value Or= Culngint(value Shr 2)   value Or= Culngint(value Shr 4)   value Or= Culngint(value Shr 8)   value Or= Culngint(value Shr 16)   value Or= Culngint(value Shr 32)   Return tab64((Culngint((value - (value Shr 1)) * &h07EDD5E59A4E28C2) Shr 58))End FunctionFunction ipower (Byval x As Ulongint, Byval e As Ulongint) As Ulongint    'take x to an integer power    Dim As Ulongint z, y, n   y = x    n = e    z = 1    While n > 0        While (n And 1) = 0            n = n \ 2            y = y * y        Wend        n = n - 1        z = y * z    Wend    Return zEnd FunctionFunction shl32(Byval n As Ulong, Byval k As Ubyte, Byref c As Ulong) As Ulong   If k>0 And k<32 Then      Dim As Ulong carry=n      Dim As Ubyte k32=32-k      Asm         mov cl, [k32]         Shr dword Ptr [carry], cl         mov cl, [k]         Shl dword Ptr [n], cl      End Asm      c=carry   End If   Return nEnd FunctionFunction shr32(Byval n As Ulong, Byval k As Ubyte, Byref c As Ulong) As Ulong   If k>0 And k<32 Then      Dim As Ulong carry=n      Dim As Ubyte k32=32-k      Asm         mov cl, [k32]         Shl dword Ptr [carry], cl         mov cl, [k]         Shr dword Ptr [n], cl      End Asm      c=carry   End If   Return nEnd FunctionSub shiftl(Byref fac1 As bigfloat_struct, Byval k As Long, Byval dwords As Long=NUM_DWORDS)   If k>0 And k<32 Then      Dim As Long i      Dim As Ulong n, carry, c=0      Dim As Long k32=32-k      For i=dwords To 0 Step -1         n=fac1.mantissa[i]         carry=n         carry=Cast(Ulong, (carry Shr k32))         n=Cast(Ulong, (n Shl k))         fac1.mantissa[i]=n+c         c=carry      Next   Elseif k=32 Then      Dim As Long i      For i=0 To dwords-1         fac1.mantissa[i]=fac1.mantissa[i+1]      Next      fac1.mantissa[dwords]=0   End IfEnd SubSub shiftr(Byref fac1 As bigfloat_struct, Byval k As Long, Byval dwords As Long=NUM_DWORDS)   If k>0 And k<32 Then      Dim As Long i      Dim As Ulong n, carry, c=0      Dim As Long k32=32-k      For i=0 To dwords         n=fac1.mantissa[i]         carry=n         carry=Cast(Ulong, (carry Shl k32))         n=Cast(Ulong, (n Shr k))            fac1.mantissa[i]=c+n         c=carry      Next   Elseif k=32 Then      Dim As Long i      For i=dwords To 1 Step -1         fac1.mantissa[i]=fac1.mantissa[i-1]      Next      fac1.mantissa[0]=0   End IfEnd SubFunction fpcmp(Byref x As bigfloat_struct, Byref y As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As Long   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   Dim As Long c, i    If x.sign < y.sign Then fpcmp = -1: Exit Function    If x.sign > y.sign Then fpcmp = 1: Exit Function   If x.sign=y.sign Then      If x.exponent=y.exponent Then         For i=0 To dwords            c=x.mantissa[i]-y.mantissa[i]            If c<>0 Then Exit For         Next         If c<0 Then Return -1         If c=0 Then Return 0         If c>0 Then Return 1      End If      If x.exponent<y.exponent Then Return -1      If x.exponent>y.exponent Then Return 1   End IfEnd FunctionFunction NORM_FAC1(Byref fac1 As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As Integer   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS    ' normalize the number in fac1    ' all routines exit through this one.    'see if the mantissa is all zeros.    'if so, set the exponent and sign equal to 0.    Dim As Integer i,er=0,f=0    For i=0 To dwords        If fac1.mantissa[i]>0 Then f=1    Next    If f=0 Then        fac1.exponent=0        fac1.sign=0        Exit Function   End If   While fac1.mantissa[0]=0      shiftl(fac1, 32, dwords)      fac1.exponent-=32      If fac1.exponent=0 Then         NORM_FAC1=EXPU_ERR         Exit Function      End If   Wend   i=__builtin_clz(fac1.mantissa[0])-3   If i>0 Then      shiftl(fac1, i, dwords)      fac1.exponent-=i   End If   'if the highmost Bit in fac1_man is nonzero,   'shift the mantissa right 1 Bit and   'increment the exponent          If fac1.mantissa[0]>(&h1ffffffful)Then      While fac1.mantissa[0]>&h1ffffffful         shiftr(fac1, 1, dwords)         fac1.exponent+=1      Wend   Elseif fac1.mantissa[0]<&h10000000ul Then   /' the following will probably never be executed '/      'now shift fac1_man 1 to the left until a      'nonzero Bit appears in the next-to-highest      'Bit of fac1_man.  decrement exponent for      'each shift.      While fac1.mantissa[0]<&h10000000ul         shiftl(fac1, 1, dwords)         fac1.exponent-=1      Wend    End If    'check for overflow/underflow    If fac1.exponent<0 Then        NORM_FAC1=EXPO_ERR    End IfEnd FunctionFunction si2fp(Byval m As Long, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   Dim As bigfloat_struct fac1   Dim As Long n=Abs(m)      If m=0 Then      Return fac1   End If   fac1.mantissa[0]=n   NORM_FAC1(fac1, NUM_DWORDS)   fac1.exponent=log2_32(n)+BIAS+1      If m<0 Then      fac1.sign=&H8000   Else      fac1.sign=0   End If   Return fac1End FunctionFunction ui2fp(Byval m As Ulong, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   Dim As bigfloat_struct fac1      If m=0 Then      Return fac1   End If   fac1.mantissa[0]=m   NORM_FAC1(fac1, NUM_DWORDS)   fac1.exponent=log2_32(m)+BIAS+1      Return fac1End FunctionSub fpadd_aux(Byref fac1 As bigfloat_struct, Byref fac2 As bigfloat_struct, Byval dwords As Long=NUM_DWORDS)   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   Dim As Long c, i   Dim As Longint v   c=0   For i=dwords To 1 Step -1      v=Clngint(fac2.mantissa[i])+Clngint(fac1.mantissa[i])+c      If v>(&h100000000ull-1) Then         v=v-&h100000000ull         c=1      Else         c=0      End If      fac1.mantissa[i]=v   Next   v=Clngint(fac1.mantissa[0])+Clngint(fac2.mantissa[0])+c   fac1.mantissa[0]=v    NORM_FAC1(fac1, dwords)End SubSub fpsub_aux(Byref fac1 As bigfloat_struct, Byref fac2 As bigfloat_struct, Byval dwords As Long=NUM_DWORDS)   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   Dim As Long c, i   Dim As Longint v   c=0   For i=dwords To 1 Step -1      v=Clngint(fac1.mantissa[i])-clngint(fac2.mantissa[i])-c      If v<0 Then         v=v+&h100000000ull         c=1      Else         c=0      End If      fac1.mantissa[i]=v   Next   v=Clngint(fac1.mantissa[0])-clngint(fac2.mantissa[0])-c   fac1.mantissa[0]=v    NORM_FAC1(fac1, dwords)End SubFunction fpadd(Byref x As bigfloat_struct, Byref y As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS       Dim As bigfloat_struct fac1,fac2    Dim As Long i, t, c, xsign, ysign       xsign=x.sign:x.sign=0    ysign=y.sign:y.sign=0    c=fpcmp(x, y, dwords)    x.sign=xsign    y.sign=ysign    If c<0 Then        fac1=y        fac2=x    Else        fac1=x        fac2=y    End If    t=fac1.exponent-fac2.exponent    If t<(NUM_BITS+32) Then        'The difference between the two        'exponents indicate how many times        'we have to multiply the mantissa        'of FAC2 by 10 (i.e., shift it right 1 place).        'If we have to shift more times than        'we have dwords, the result is already in FAC1.               If t>0 And t<(NUM_BITS+32) Then 'shift                i=t\32         While i>0            shiftr(fac2, 32, dwords)            t-=32            i-=1         Wend         'While t>0                  If t>0 Then   shiftr(fac2, t, dwords)         '   t-=1         'Wend        End If        'See if the signs of the two numbers        'are the same.  If so, add; if not, subtract.        If fac1.sign=fac2.sign Then 'add            fpadd_aux(fac1,fac2, dwords)        Else            fpsub_aux(fac1,fac2, dwords)        End If    Endif    Return fac1End FunctionFunction fpsub(Byref x As bigfloat_struct, Byref y As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS    Dim As bigfloat_struct fac1,fac2    Dim As Long s    fac1=x    fac2=y    fac2.sign=fac2.sign Xor &h8000    fac1=fpadd(fac1,fac2, dwords)    Return fac1End FunctionFunction fpmul(Byref x As bigfloat_struct, Byref y As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS    Dim As bigfloat_struct fac1,fac2    Dim As Integer i, j, ex, er, den, num    Dim As Ulongint digit, carry, tmp, prod    Dim As Ulong  fac3(0 To 2*dwords+1)       fac1=x    fac2=y    'check exponents.  if either is zero,    'the result is zero    If fac1.exponent=0 Or fac2.exponent=0 Then 'result is zero...clear fac1.        fac1.sign=0        fac1.exponent=0        For i=0 To dwords            fac1.mantissa[i]=0        Next        'NORM_FAC1(fac1)        Return fac1    Else               If ex<0 Then            er=EXPO_ERR            Return fac1 'Exit Function        End If        'clear fac3 mantissa        For i=0 To dwords            fac3(i)=0        Next      den=dwords      While fac2.mantissa[den]=0         den-=1      Wend      num=dwords      While fac1.mantissa[num]=0         num-=1      Wend            If num<den Then         Swap fac1, fac2         'fac1=y         'fac2=x                  Swap den, num      End If      For j=den To 0 Step -1         carry=0         digit=fac2.mantissa[j]         For i=num To 0 Step -1         /'   prod=fac3(i+j+1)+digit*fac1.mantissa[i]+carry '/         tmp=fac1.mantissa[i]            If __builtin_umulll_overflow (digit, tmp, prod) Then Print "mul overflow"         tmp=fac3(i+j+1)            If __builtin_uaddll_overflow (tmp, prod, prod) Then Print "add overflow"            If __builtin_uaddll_overflow (carry, prod, prod) Then Print "add overflow"            carry=prod   Shr 32'\&H100000000            fac3(i+j+1)=prod '(prod mod &H100000000)         Next         fac3(j)=carry      Next      For i=0 To dwords         fac1.mantissa[i]=fac3(i)      Next         End If   'now determine exponent of result.   'as you do...watch for overflow.   'ex=fac2.exponent-BIAS+fac1.exponent   'fac1.exponent=ex      ex=(fac2.exponent And &h7FFFFFFFul)-BIAS+1   ex=ex+(fac1.exponent And &h7FFFFFFFul)-BIAS+1   fac1.exponent=ex+BIAS+1      'determine the sign of the product   fac1.sign=fac1.sign Xor fac2.sign    NORM_FAC1(fac1, dwords)    Return fac1End FunctionFunction fpmul_si(Byref x As bigfloat_struct, Byval y As Long, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS    Dim As bigfloat_struct fac1,fac2    fac1=x   fac2=ui2fp(y, dwords)   fac1=fpmul(fac1, fac2, dwords)   Return fac1End FunctionFunction fpdiv_si(Byref num As bigfloat_struct, Byval den As Long, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   Dim As bigfloat_struct fac1    Dim As Ulongint carry, remder=0    Dim As Longint i, divisor    Dim As Longint quotient   Dim As Long j      divisor=Abs(den)   fac1=num   If divisor=1 Then      Return fac1   End If    If divisor = 0 Then        Print "error: divisor = 0"        Return fac1 'Exit function    End If   If divisor>2147483647 Then        Print "error: divisor too large"        Return fac1 'Exit function    End If       For i = 0 To dwords        quotient = fac1.mantissa[i] + remder * &h100000000ull        remder = quotient Mod divisor        fac1.mantissa[i]=quotient \ divisor    Next    quotient = remder * &h100000000ull    quotient=quotient \ divisor    carry=fac1.mantissa[0]   j=__builtin_clz(fac1.mantissa[0])-3   If j>0 Then      shiftl(fac1, j, dwords)      fac1.exponent-=j   End If    While fac1.mantissa[0]>(&h1ffffffful)      shiftr(fac1, 1, dwords)        fac1.exponent+=1    Wend    NORM_FAC1(fac1)    If den<0 Then      fac1.sign=fac1.sign Xor &h8000   End If   Return fac1End FunctionFunction fpdiv(Byref x As bigfloat_struct, Byref y As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS    Dim As bigfloat_struct fac1, fac2, one    Dim As Integer i, er, is_power_of_two   Dim As Short sign   fac1=x   fac2=y      one.exponent=(BIAS+1)   one.mantissa[0]=&h10000000ul   one=ui2fp(1)   sign=fac2.sign   fac2.sign=0   If fpcmp(fac2, one)=0 Then      Return fac1   Else      fac2.sign=sign   End If    If fac2.exponent=0 Then ' if fac2 = 0, return        ' a divide-by-zero error and        ' bail out.        fac1.mantissa[i]=&H1FFFFFFFul        For i=1 To dwords         fac1.mantissa[i]=&HFFFFFFFFul      Next      fac1.exponent=&HFFFFF+BIAS+1        er=DIVZ_ERR        Return fac1    Elseif fac1.exponent=0 Then 'fact1=0, just return        er=0        Return fac1    Else      'check to see if fac2 is a power of ten      is_power_of_two=0      If fac2.mantissa[0]=&H10000000ul Then         is_power_of_two=1         For i=1 To dwords            If fac2.mantissa[i]<>0 Then               is_power_of_two=0               Exit For            End If         Next      End If      'if fac2 is a power of ten then all we need to do is to adjust the sign and exponent and we are finished      If is_power_of_two=1 Then         fac1.sign=fac1.sign Xor fac2.sign         fac1.exponent=fac1.exponent-fac2.exponent+BIAS+1         Return fac1      End If      #Define min(a,b) Iif(((a)<(b)),(a),(b))      #Macro realw(w, j)         ((w(j - 1)*b + w(j))*b + w(j + 1))*b +Iif(Ubound(w)>=j+2,w(j+2),0)      #Endmacro      #Macro subtract(w, q, d, ka, kb)         For j=ka To kb            w(j) = w(j) - q*d(j - ka + 2)         Next      #Endmacro      #Macro normalize(w, ka, q)         w(ka) = w(ka) + w(ka - 1)*b         w(ka - 1) = q      #Endmacro      #Macro finalnorm(w, kb)         For j=kb To 3 Step -1            carry=Iif(w(j)<0, ((-w(j) - 1)\b) + 1, Iif(w(j) >= b, -(w(j)\b), 0))            w(j) = w(j) + carry*b            w(j - 1) = w(j - 1) - carry         Next      #Endmacro      Dim As Double result(1 To 2*dwords+3), n(1 To 2*dwords+3), d(1 To 2*dwords+3)      Const b=&H10000         Dim As Long j, last, laststep, q, t      Dim As Long stp, carry      Dim As Double xd, xn, rund      Dim As Double w(1 To Ubound(n)+4)      For j=0 To dwords         n(2*j+2)=fac1.mantissa[j]\&H10000         n(2*j+3)=fac1.mantissa[j] Mod &H10000         d(2*j+2)=fac2.mantissa[j]\&H10000         d(2*j+3)=fac2.mantissa[j] Mod &H10000      Next      n(1)=(fac1.exponent And &h7FFFFFFF)-BIAS-1      d(1)=(fac2.exponent And &h7FFFFFFF)-BIAS-1      For j=Ubound(n) To Ubound(w)         w(j)=0      Next      t=Ubound(n)-1      w(1)=n(1)-d(1)+1      w(2)=0      For j=2 To Ubound(n)         w(j+1)=n(j)      Next         xd = (d(2)*b + d(3))*b + d(4) + d(5)/b      laststep = t + 2      For stp=1 To laststep         xn=RealW(w, (stp + 2))         q=Int(xn/xd)         last = Min(stp + t + 1, Ubound(W))         subtract(w, q, d, (stp + 2), last)         normalize(w, (stp + 2), q)      Next      FinalNorm(w, (laststep + 1))      laststep = Iif(w(2) = 0, laststep, laststep - 1)      rund = w(laststep + 1)/b      w(laststep) = w(laststep) + Iif(rund >= 0.5, 1, 0)      If w(2)=0 Then         For j=1 To t+1            result(j)=w(j+1)         Next      Else         For j=1 To t+1            result(j)=w(j)         Next      End If      result(1) = Iif(w(2) = 0, w(1) - 1, w(1))            For j=0 To dwords         fac1.mantissa[j]=result(2*j+2)*&H10000+result(2*j+3)      Next      NORM_FAC1(fac1, dwords)      fac1.exponent=(result(1)+BIAS)   End If   fac1.sign=fac1.sign Xor fac2.sign    Return fac1End FunctionFunction fpipow(Byref x As bigfloat_struct, Byval e As Longint, Byval dwords As Long=NUM_DWORDS) As bigfloat_struct   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   'take x to an Long power   Dim As bigfloat_struct y=x   Dim As bigfloat_struct z, one   Dim As Longint n, c=0   Dim As Integer i      n = Abs(e)   one.exponent=(BIAS+1)   one.mantissa[0]=&h10000000ul   z=one   While n>0      While (n And 1)=0         n\=2         y=fpmul(y, y, dwords)         c+=1      Wend      n-=1      z=fpmul(y, z, dwords)      c+=1   Wend      If e<0 Then      z=fpdiv(one, z, dwords)   End If   Return zEnd FunctionFunction dbl2fp(Byval x As Double) As bigfloat_struct   Dim As bigfloat_struct bf   Dim As Ulongint n, e      n=Peek ( Ulongint, @x )   e=n   n=Cast(Ulongint, (n Shl 12))   n=Cast(Ulongint, (n Shr 4))   n+=&h1000000000000000ull   bf.mantissa[0]=Cast(Ulongint, (n Shr 32))   n=e   n=Cast(Ulongint, (n Shl 40))   n=Cast(Ulongint, (n Shr 32))   bf.mantissa[1]=n   n=n+&h1000000000000000ull   e=Cast(Ulongint, (e Shr 52))   e=e-&h3FFul   bf.exponent=e+BIAS+1   Return bfEnd FunctionFunction fp2dbl(Byref x As bigfloat_struct) As Double   Dim As Double dbl   Dim As Ulongint n, e   Dim As Ushort ex   Dim As Long sign   If x.sign<>0 Then      sign=-1   Else      sign=1   End If    If x.exponent>0 Then        ex=(x.exponent And &h7FFFFFFF)-BIAS-1    Else        ex=0    End If   n=x.mantissa[0]   n-=&h10000000ul   n=Cast(Ulongint, (n Shl 24))   e=&h3FF+ex   e=e*&h10000000000000ull   n=n+e   e=x.mantissa[1]\&h100   n=n+e   Poke Ulongint,@dbl,n   Return dbl*signEnd FunctionFunction fp2dbl2(Byref x As bigfloat_struct) As Double   Dim As Double dbl   Dim As Ulongint n, e   n=x.mantissa[0]   n-=&h10000000ul   n=Cast(Ulongint, (n Shl 24))   e=&h3FF+1   e=e*&h10000000000000ull   n=n+e   e=x.mantissa[1]\&h100   n=n+e   Poke Ulongint,@dbl,n   Return dblEnd FunctionFunction factorial(Byval n As Ulong) As bigfloat_struct   Dim As bigfloat_struct one, f, indx   one.mantissa[0]=&h10000000ul   one.exponent=0+BIAS+1   f=one   indx=one   For i As Long=1 To n      f=fpmul(f, indx)      indx=fpadd(one, indx)   Next   Return fEnd FunctionFunction factorial_inv(Byval n As Ulong) As bigfloat_struct   Dim As bigfloat_struct f   f.mantissa[0]=&h10000000ul   f.exponent=0+BIAS+1   For i As Long=1 To n      f=fpdiv_si(f, i)   Next   Return fEnd FunctionFunction fp2ui(Byref x As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As Ulong   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   Dim As bigfloat_struct fac1   Dim As Double n   fac1.exponent=0   fac1.mantissa[0]=x.mantissa[0]   n=fp2dbl(x)   Return Fix(n)End FunctionFunction fp2si(Byref x As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As Long   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   Dim As bigfloat_struct fac1   Dim As Long sign   Dim As Double n      fac1.exponent=0   fac1.mantissa[0]=x.mantissa[0]   n=fp2dbl(x)   Return Fix(n)End FunctionFunction str2fp_aux(Byval value As String, Byref s As Long, Byref ex As Long ) As String   Dim As Integer j,d,e,ep,es,i,f,fp,fln   Dim As String c,f1,f2,f3, ts   Dim As Ulong ulng   j=1   s=1   d=0   e=0   ep=0   ex=0   es=1   i=0   f=0   fp=0   f1=""   f2=""   f3=""   value=Ucase(value)   fln=Len(value)   While j<=fln      c=Mid(value,j,1)      If ep=1 Then         If c=" " Then            j=j+1            Continue While         Endif         If c="-" Then            es=-es            c=""         Endif         If c="+" Then            j=j+1            Continue While         Endif         If (c="0") And (f3="") Then            j=j+1            Continue While         Endif         If (c>"/") And (c<":") Then 'c is digit between 0 and 9            f3=f3+c            ex=10*ex+(Asc(c)-48)            j=j+1            Continue While         Endif      Endif      If c=" " Then         j=j+1         Continue While      Endif      If c="-" Then         s=-s         j=j+1         Continue While      Endif      If c="+" Then         j=j+1         Continue While      Endif      If c="." Then         If d=1 Then            j=j+1            Continue While         Endif         d=1      Endif      If (c>"/") And (c<":") Then 'c is digit between 0 and 9         If ((c="0") And (i=0)) Then            If d=0 Then               j=j+1               Continue While            Endif            If (d=1) And (f=0) Then               e=e-1               j=j+1               Continue While            Endif         Endif         If d=0 Then            f1=f1+c            i=i+1         Else            If (c>"0") Then               fp=1            Endif            f2=f2+c            f=f+1         Endif      Endif      If c="E" Or c="D" Then         ep=1      Endif      j=j+1   Wend   If fp=0 Then      f=0      f2=""   Endif   If s=-1 Then s=&h8000 Else s=0      f1=f1+f2   f1=Mid(f1,1,1)+Right(f1,Len(f1)-1)   fln=Len(f1)   If fp=0 Andalso ex=0 Then      ex=es*fln-1   Else      ex=es*ex-1+i+e   End If   If fln>((NUM_BITS+32)*0.3010299956639811) Then      f1=Mid(f1,1,((NUM_BITS+32)*0.3010299956639811)-1)   Endif   Return f1End FunctionFunction str2fp(Byref x As String) As bigfloat_struct   Dim As Long strlen, i, n, ex, sign   Dim As String s2, strin=x   Dim As bigfloat_struct y, z, pw, ten         strin=str2fp_aux(strin, sign, ex)      strlen=Len(strin)      n=1      ten=ui2fp(Valuint("1"+String(n,"0")))      If strlen>n Then         s2=Left(strin, n)         strin=Mid(strin, n+1)      Else         z=ui2fp(Valuint(strin))         If ex<>0 Then            pw=ui2fp(10)            pw=fpipow(pw, strlen-ex-1)            z=fpdiv(z, pw)         End If         Return z      End If      z=ui2fp(Valuint(s2))      While Len(strin)>=n         s2=Left(strin, n)         strin=Mid(strin, n+1)         y=ui2fp(Valuint(s2))         z=fpmul(z, ten)         z=fpadd(z,y)      Wend      pw=ui2fp(10)      pw=fpipow(pw,strlen-ex-1)      z=fpdiv(z, pw)      Return zEnd FunctionFunction fp2str(Byref zz As bigfloat_struct, Byval dwords As Long=NUM_DWORDS) As String   Dim As Long ex, i, n, powten, sign, tn   Dim As Long m   Dim As bigfloat_struct y, ten, z, ep   Dim As String s, s2      If zz.exponent=0 Then         Return "0"      End If      z=zz      z.sign=0      sign=zz.sign      If z.exponent<>0 Then         ex=(z.exponent And &h7FFFFFFFul)-BIAS-1      Else         ex=0      End If      n=NUM_BITS-ex      ep.exponent=-(n+15)+BIAS+1      ep.mantissa[0]=&h10000000ul      z=fpadd(z, ep)      If ex>0 Then         While ex>=3            z=fpdiv_si(z, 10)            powten+=1            ex=(z.exponent And &h7FFFFFFFul)-BIAS-1         Wend      Else         While ex<0            z=fpmul_si(z, 10)            powten-=1            ex=(z.exponent And &h7FFFFFFFul)-BIAS-1         Wend      End If            n=1      m=fp2ui(z)      If m=0 Then         z=fpmul_si(z, 10)         powten-=1         ex=(z.exponent And &h7FFFFFFFul)-BIAS-1         m=fp2ui(z)      End If      s=Trim(Str(m))+"."      If sign<>0 Then         s="-"+s      Else         s=" "+s      End If      ten=ui2fp(Clng("1"+String(n,"0")))      tn=Clng("1"+String(n,"0"))      For i=1 To (dwords*9.63)\n+1         y=ui2fp(m)         z=fpsub(z,y)         z=fpmul(z,ten)         m=fp2ui(z)         If m<0 Or m>tn Then            m=0         End If         s2=Trim(Str(m))         s2=String(n-len(s2),"0")+s2         s=s+s2      Next      s2=Trim(Str(Abs(powten)))      s2=String(5-len(s2),"0")+s2      If powten<0 Then         s2="-"+s2      Else         s2="+"+s2      End If      s=s+"e"+s2      Return sEnd Function'integer part of numFunction fpfix( Byref num As bigfloat_struct ) As bigfloat_struct   Dim As bigfloat_struct ip   Dim As Long ex, ex2, j, k, c      ex=(num.exponent And &h7FFFFFFFul)-BIAS   If ex<1 Then      Return ip   End If   If ex>=(NUM_BITS) Then      Return num   End If   ex2=ex\32   k=ex2   j=ex Mod 32   While ex2>0      ex2-=1      ip.mantissa[ex2]=num.mantissa[ex2]   Wend   ip.mantissa[k]=num.mantissa[k]   ip.mantissa[k]=shr32(ip.mantissa[k], (32-j)-3, c)   ip.mantissa[k]=shl32(ip.mantissa[k], (32-j)-3, c)   ip.exponent=ex+BIAS   ip.sign=num.sign   Return ipEnd Function'fractional part of numFunction fpfrac( Byref num As bigfloat_struct ) As bigfloat_struct   Dim As bigfloat_struct ip, fp   ip=fpfix(num)   fp=fpsub(num, ip)   Return fpEnd FunctionDim As bigfloat_struct x, x2, x3, y, z, term, fac, pi, epsDim As Long i, n, signDim As String s, spiDim As Double teps.exponent=-(NUM_BITS)+BIAS+1eps.mantissa[0]=&h10000000ulspi="3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821"pi=str2fp(spi)Print "Pi+Pi"y=fpadd(pi, pi)Print fp2str(y)Print "Pi/2"y=fpdiv_si(pi, 2)Print fp2str(y)Print "integer part of Pi/2"z=fpfix(y)Print fp2str(z)Print "fractional part of Pi/2"z=fpfrac(y)Print fp2str(z)Print "Pi*Pi"z=fpmul(pi, pi)Print fp2str(z)Print "factorial(10000)"z=factorial(10000)Print fp2str(z)Print "Pi^100"z=fpipow(pi, 100)Print fp2str(z)t=Timer'' sin(.5)x=str2fp(".5")y=xx3=xx2=fpmul(x, x)term=ui2fp(1)fac=termsign=-1i=1While fpcmp(term, eps)>0 ' fpcmp returns -1 if term<eps, 0 if equal and 1 if greater   x3=fpmul(x3, x2)   fac=fpdiv_si(fac, i+1)   fac=fpdiv_si(fac, i+2)   term=fpmul(x3, fac)   If sign<0 Then      y=fpsub(y, term)   Else      y=fpadd(y, term)   End If   i+=2   sign=-signWendt=timer-tPrint "sin(.5)"Print fp2str(y)Print "time to comput sine function ";t;" seconds"PrintPrint "you can specify the number of dwords to the fp2str function"Print fp2str(pi,2)Sleep`
srvaldez
Posts: 2850
Joined: Sep 25, 2005 21:54

Code: Select all

`Extern "c"   Declare Function __builtin_clz (Byval x As Ulong) As Long   Declare Function __builtin_smull_overflow (Byval As Long, Byval As Long, Byref As Long) As Boolean   Declare Function __builtin_smulll_overflow (Byval As Longint, Byval As Longint, Byref As Longint) As Boolean   Declare Function __builtin_umull_overflow (Byval As Ulong, Byval As Ulong, Byref As Ulong) As Boolean   Declare Function __builtin_umulll_overflow (Byval As Ulongint, Byval As Ulongint, Byref As Ulongint) As Boolean   Declare Function __builtin_uaddll_overflow (Byval As Ulongint, Byval As Ulongint, Byref As Ulongint) As BooleanEnd ExternConst NUMBER_OF_DIGITS = 96Const NUMBER_OF_BITS   =  NUMBER_OF_DIGITS*3.321928094887362#If (NUMBER_OF_BITS Mod 32)<>0   Const NUM_BITS = (NUMBER_OF_BITS\32+1)*32#Else   Const NUM_BITS = NUMBER_OF_BITS#EndifConst NUM_DWORDS    =   NUM_BITS\32Const BIAS          =   1073741824 '2 ^ 30' Error definitions#Define DIVZ_ERR        1                    'Divide by zero#Define EXPO_ERR        2                    'Exponent overflow error#Define EXPU_ERR        3                    'Exponent underflow errorType BigFloat_struct   Declare Constructor ( )   Declare Constructor ( Byref rhs As BigFloat_struct )   Declare Destructor ( )   Declare Operator Let ( Byref rhs As BigFloat_struct )   As Short sign   As Ulong  exponent   As Ulong  Ptr mantissaEnd TypeConstructor BigFloat_struct ( )   If this.mantissa<>0 Then      Print "pointer is not 0"      End(1)   End If   this.mantissa=Callocate((NUM_DWORDS +1), 4)   If this.mantissa=0 Then      Print "unable to allocate memory"      End(4)   End If   this.sign=0   this.exponent=0End ConstructorConstructor BigFloat_struct ( Byref rhs As BigFloat_struct )   If this.mantissa<>0 Then      Print "pointer is not 0"      End(1)   End If   this.mantissa=Callocate((NUM_DWORDS +1), 4)   If this.mantissa=0 Then      Print "unable to allocate memory"      End(4)   End If   this.sign=rhs.sign   this.exponent=rhs.exponent   For i As Long=0 To NUM_DWORDS      this.mantissa[i]=rhs.mantissa[i]   NextEnd ConstructorOperator BigFloat_struct.let ( Byref rhs As BigFloat_struct )   this.sign=rhs.sign   this.exponent=rhs.exponent   For i As Long=0 To NUM_DWORDS      this.mantissa[i]=rhs.mantissa[i]   NextEnd OperatorDestructor BigFloat_struct ( )   If this.mantissa=0 Then      Print "unable to de-allocate memory"      End(-4)   Else      Deallocate(this.mantissa)   End If   this.mantissa=0End DestructorType BigFloat   Declare Constructor ( )   Declare Constructor ( Byval rhs As Long )   Declare Constructor ( Byval rhs As Double )   Declare Constructor ( Byref rhs As String )   Declare Constructor ( Byref rhs As BigFloat )   Declare Destructor ( )      Declare Operator Let ( Byval rhs As Long )   Declare Operator Let ( Byval rhs As Double )   Declare Operator Let ( Byref rhs As String )   Declare Operator Let ( Byref rhs As BigFloat )   Declare Operator Cast ( ) As String   Declare Operator Cast ( ) As Long   Declare Operator Cast ( ) As Double      '----------------------------------------------   Declare Operator += (Byref rhs As BigFloat)   Declare Operator += (Byval rhs As Long)   Declare Operator += (Byval rhs As Double)   Declare Operator += (Byref rhs As String)   Declare Operator -= (Byref rhs As BigFloat)   Declare Operator -= (Byval rhs As Long)   Declare Operator -= (Byval rhs As Double)   Declare Operator -= (Byref rhs As String)   Declare Operator *= (Byref rhs As BigFloat)   Declare Operator *= (Byval rhs As Long)   Declare Operator *= (Byval rhs As Double)   Declare Operator *= (Byref rhs As String)   Declare Operator /= (Byref rhs As BigFloat)   Declare Operator /= (Byval rhs As Long)   Declare Operator /= (Byval rhs As Double)   Declare Operator /= (Byval rhs As Single)   Declare Operator /= (Byref rhs As String)   ' For Next Implicit step = +1   Declare Operator For ( )   Declare Operator Step( )   Declare Operator Next( Byref end_cond As BigFloat ) As Integer   ' For Next Exlicit step   Declare Operator For ( Byref stp As BigFloat )   Declare Operator Step( Byref stp As BigFloat )   Declare Operator Next( Byref end_cond As BigFloat, Byref step_var As BigFloat ) As Integer      Declare Function toString( byval places as long=NUM_DWORDS*9.63 ) As String   Declare Function toLong ( ) As Long   Declare Function toDouble ( ) As Double   BigNum As BigFloat_structEnd TypeDeclare Function fpcmp(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As LongDeclare Function fpadd(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_structDeclare Function fpsub(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_structDeclare Function fpmul(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_structDeclare Function fpmul_si(Byref x As BigFloat_struct, Byval y As Long, Byval dwords As Long=NUM_DWORDS) As BigFloat_structDeclare Function fpdiv_si(Byref num As BigFloat_struct, Byval den As Long, Byval dwords As Long=NUM_DWORDS) As BigFloat_structDeclare Function fpdiv(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_structDeclare Function fpipow(Byref x As BigFloat_struct, Byval e As Longint, Byval dwords As Long=NUM_DWORDS) As BigFloat_structDeclare Function dbl2fp(Byval x As Double) As BigFloat_structDeclare Function fp2dbl(Byref x As BigFloat_struct) As DoubleDeclare Function ui2fp(Byval m As Ulong, Byval dwords As Long=NUM_DWORDS) As BigFloat_structDeclare Function fp2ui(Byref x As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As UlongDeclare Function si2fp(Byval m As Long, Byval dwords As Long=NUM_DWORDS) As BigFloat_structDeclare Function fp2si(Byref x As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As LongDeclare Function str2fp(Byref x As String) As BigFloat_structDeclare Function fp2str(Byref zz As BigFloat_struct, Byval digits As Long=NUM_DWORDS*9.63) As StringDeclare Function fpfix( Byref num As BigFloat_struct ) As BigFloat_structDeclare Function fpfrac( Byref num As BigFloat_struct ) As BigFloat_structFunction Bigfloat.toString( byval places as long=NUM_DWORDS*9.63 ) As String   Function = fp2str( this.BigNum, places )End FunctionFunction Bigfloat.toLong ( ) As Long   dim as double x   x=fp2dbl(this.BigNum)   return clng(x)End FunctionFunction Bigfloat.toDouble ( ) As Double   Function = fp2dbl(this.BigNum)End FunctionConstructor Bigfloat ( )   this.BigNum=si2fp(0, NUM_DWORDS)End ConstructorConstructor Bigfloat ( Byval rhs As Long )   this.BigNum=si2fp( rhs, NUM_DWORDS )End ConstructorConstructor Bigfloat ( Byref rhs As String )   this.BigNum=str2fp( rhs )End ConstructorConstructor Bigfloat ( Byref rhs As Bigfloat)   this.BigNum.sign=rhs.BigNum.sign   this.BigNum.exponent=rhs.BigNum.exponent   for i as long=0 to NUM_DWORDS      this.BigNum.mantissa[i]=rhs.BigNum.mantissa[i]   nextEnd ConstructorConstructor Bigfloat ( Byval rhs As Double )   this.BigNum=str2fp( str( rhs ) )End ConstructorOperator Bigfloat.let ( Byref rhs As Bigfloat )   this.BigNum.sign=rhs.BigNum.sign   this.BigNum.exponent=rhs.BigNum.exponent   for i as long=0 to NUM_DWORDS      this.BigNum.mantissa[i]=rhs.BigNum.mantissa[i]   nextEnd OperatorOperator Bigfloat.let ( Byval rhs As Long )   this.BigNum=si2fp( rhs, NUM_DWORDS )End OperatorOperator Bigfloat.let ( Byref rhs As String )   this.BigNum=str2fp( rhs )End OperatorOperator Bigfloat.Let ( Byval rhs As Double )   this.BigNum=str2fp( str(rhs) )End OperatorOperator Bigfloat.cast ( ) As String   Operator = fp2str(this.BigNum)End OperatorOperator Bigfloat.cast ( ) As Long   dim as double x   x=fp2dbl(this.BigNum)   Operator = clng(x)End OperatorOperator Bigfloat.cast ( ) As Double   Operator = fp2dbl(this.BigNum)End OperatorDestructor Bigfloat ( )End Destructor'============================================================================Operator BigFloat.for ( )End Operator Operator BigFloat.step ( )   this.BigNum = fpadd(this.BigNum,si2fp(1, NUM_DWORDS), NUM_DWORDS) End Operator  Operator BigFloat.next ( Byref end_cond As BigFloat ) As Integer   Return fpcmp(This.BigNum, end_cond.BigNum, NUM_DWORDS)<=0End Operator  '' explicit step versions'' Operator BigFloat.for ( Byref step_var As BigFloat )End Operator Operator BigFloat.step ( Byref step_var As BigFloat )   this.BigNum = fpadd(this.BigNum, step_var.BigNum, NUM_DWORDS)   End Operator  Operator BigFloat.next ( Byref end_cond As BigFloat, Byref step_var As BigFloat ) As Integer   If fpcmp(step_var.BigNum, si2fp(0, NUM_DWORDS), NUM_DWORDS) < 0 Then      Return fpcmp(This.BigNum, end_cond.BigNum, NUM_DWORDS) >= 0   Else      Return fpcmp(This.BigNum, end_cond.BigNum, NUM_DWORDS) <= 0   End IfEnd Operator'============================================================================Operator + ( Byref lhs As Bigfloat, Byval rhs As Bigfloat ) As Bigfloat   Dim As Bigfloat result   result.BigNum= fpadd(lhs.BigNum, rhs.BigNum, NUM_DWORDS)   Operator = resultEnd OperatorOperator + ( Byref lhs As Bigfloat, Byval rhs As Double ) As Bigfloat   Dim As Bigfloat result   result.BigNum=str2fp( str(rhs) )   result.BigNum= fpadd(lhs.BigNum, result.BigNum, NUM_DWORDS)   Operator = resultEnd OperatorOperator + ( Byval lhs As Double, Byref rhs As Bigfloat ) As Bigfloat   Dim As Bigfloat result   result.BigNum=str2fp( str(lhs) )   result.BigNum= fpadd( result.BigNum, rhs.BigNum, NUM_DWORDS)   Operator = resultEnd OperatorOperator + (Byref lhs As Bigfloat, byref rhs as String) As Bigfloat   dim as Bigfloat result   result.BigNum=str2fp(rhs)   result.BigNum= fpadd(lhs.BigNum, result.BigNum, NUM_DWORDS)   Operator = resultEnd OperatorOperator - ( Byref rhs As Bigfloat ) As Bigfloat   Dim As Bigfloat result = rhs   result.BigNum.sign = result.BigNum.sign xor &h8000   Operator = resultEnd OperatorOperator - ( Byref lhs As Bigfloat, Byval rhs As Bigfloat ) As Bigfloat   Dim As Bigfloat result   result.BigNum= fpsub(lhs.BigNum, rhs.BigNum, NUM_DWORDS)   Operator = resultEnd OperatorOperator - ( Byref lhs As Bigfloat, Byval rhs As Double ) As Bigfloat   Dim As Bigfloat result   result.BigNum=str2fp( str(rhs) )   result.BigNum= fpsub(lhs.BigNum, result.BigNum, NUM_DWORDS)   Operator = resultEnd OperatorOperator - ( Byval lhs As Double, Byref rhs As Bigfloat ) As Bigfloat   Dim As Bigfloat result   result.BigNum=str2fp( str(lhs) )   result.BigNum= fpsub( result.BigNum, rhs.BigNum, NUM_DWORDS)   Operator = resultEnd OperatorOperator - ( Byref lhs As Bigfloat, Byval rhs As Long ) As Bigfloat   Dim As Bigfloat result   result.BigNum=si2fp(rhs, NUM_DWORDS)   result.BigNum= fpsub(lhs.BigNum, result.BigNum, NUM_DWORDS)   Operator = resultEnd OperatorOperator - ( Byval lhs As Long, byref rhs As Bigfloat  ) As Bigfloat   Dim As Bigfloat result   result.BigNum=si2fp(lhs, NUM_DWORDS)   result.BigNum= fpsub(result.BigNum, rhs.BigNum, NUM_DWORDS)   Operator = resultEnd OperatorOperator - (Byref lhs As Bigfloat, byref rhs as String) As Bigfloat   dim as Bigfloat result   result.BigNum=str2fp(rhs)   result.BigNum= fpsub(lhs.BigNum, result.BigNum, NUM_DWORDS)   Operator = resultEnd OperatorOperator * ( Byref lhs As Bigfloat, Byval rhs As Bigfloat ) As Bigfloat   Dim As Bigfloat result   result.BigNum= fpmul(lhs.BigNum, rhs.BigNum, NUM_DWORDS)   Operator = resultEnd OperatorOperator * ( Byref lhs As Bigfloat, Byval rhs As Double ) As Bigfloat   Dim As Bigfloat result   result.BigNum=str2fp( str(rhs) )   result.BigNum= fpmul(lhs.BigNum, result.BigNum, NUM_DWORDS)   Operator = resultEnd OperatorOperator * ( Byval lhs As Double, Byref rhs As Bigfloat ) As Bigfloat   Dim As Bigfloat result   result.BigNum=str2fp( str(lhs) )   result.BigNum= fpmul( result.BigNum, rhs.BigNum, NUM_DWORDS)   Operator = resultEnd OperatorOperator * ( Byref lhs As Bigfloat, Byval rhs As Long ) As Bigfloat   Dim As Bigfloat result   if abs(rhs)<&h100000000ull then      result.BigNum=fpmul_si(lhs.BigNum, rhs)   else      result.BigNum=si2fp(rhs, NUM_DWORDS)      result.BigNum= fpmul(lhs.BigNum, result.BigNum, NUM_DWORDS)   end if   Operator = resultEnd OperatorOperator * ( Byval lhs As Long, byref rhs As Bigfloat  ) As Bigfloat   Dim As Bigfloat result   result.BigNum=si2fp(lhs, NUM_DWORDS)   result.BigNum= fpmul(result.BigNum, rhs.BigNum, NUM_DWORDS)   Operator = resultEnd OperatorOperator * (Byref lhs As Bigfloat, byref rhs as String) As Bigfloat   dim as Bigfloat result   result.BigNum=str2fp(rhs)   result.BigNum= fpmul(lhs.BigNum, result.BigNum, NUM_DWORDS)   Operator = resultEnd OperatorOperator / ( Byref lhs As Bigfloat, Byval rhs As Bigfloat ) As Bigfloat   Dim As Bigfloat result   result.BigNum= fpdiv(lhs.BigNum, rhs.BigNum, NUM_DWORDS)   Operator = resultEnd OperatorOperator / ( Byref lhs As Bigfloat, Byval rhs As Double ) As Bigfloat   Dim As Bigfloat result   result.BigNum=str2fp( str(rhs) )   result.BigNum= fpdiv(lhs.BigNum, result.BigNum, NUM_DWORDS)   Operator = resultEnd OperatorOperator / ( Byval lhs As Double, Byref rhs As Bigfloat ) As Bigfloat   Dim As Bigfloat result   result.BigNum=str2fp( str(lhs) )   result.BigNum= fpdiv( result.BigNum, rhs.BigNum, NUM_DWORDS)   Operator = resultEnd OperatorOperator / ( Byref lhs As Bigfloat, Byval rhs As Long ) As Bigfloat   Dim As Bigfloat result   result.BigNum=si2fp(rhs, NUM_DWORDS)   result.BigNum= fpdiv(lhs.BigNum, result.BigNum, NUM_DWORDS)   Operator = resultEnd OperatorOperator / ( Byval lhs As Long, byref rhs As Bigfloat  ) As Bigfloat   Dim As Bigfloat result   result.BigNum=si2fp(lhs, NUM_DWORDS)   result.BigNum= fpdiv(result.BigNum, rhs.BigNum, NUM_DWORDS)   Operator = resultEnd OperatorOperator / (Byref lhs As Bigfloat, byref rhs as String) As Bigfloat   dim as Bigfloat result   result.BigNum=str2fp(rhs)   result.BigNum= fpdiv(lhs.BigNum, result.BigNum, NUM_DWORDS)   Operator = resultEnd OperatorOperator ^ ( Byref lhs As Bigfloat, Byval rhs As Longint ) As Bigfloat   Dim As Bigfloat lhs2   lhs2.BigNum=fpipow(lhs.BigNum, rhs, NUM_DWORDS)   Operator = lhs2end operatorOperator BigFloat.+= (Byref rhs As BigFloat)   this.BigNum = fpadd(this.BigNum, rhs.BigNum, NUM_DWORDS)End OperatorOperator BigFloat.+= (Byval rhs As Long)   Dim As BigFloat result   result.BigNum=si2fp(rhs, NUM_DWORDS)   this.BigNum = fpadd(this.BigNum, result.BigNum, NUM_DWORDS)End OperatorOperator BigFloat.+= (Byval rhs As Double)   Dim As BigFloat result   result.BigNum=str2fp( str(rhs) )   this.BigNum = fpadd(this.BigNum, result.BigNum, NUM_DWORDS)End OperatorOperator BigFloat.+= (Byref rhs As String)   Dim As BigFloat result   result.BigNum=str2fp(rhs)   this.BigNum = fpadd(this.BigNum, result.BigNum, NUM_DWORDS)End OperatorOperator BigFloat.-= (Byref rhs As BigFloat)   this.BigNum = fpsub(this.BigNum, rhs.BigNum, NUM_DWORDS)End OperatorOperator BigFloat.-= (Byval rhs As Long)   Dim As BigFloat result   result.BigNum=si2fp(rhs, NUM_DWORDS)   this.BigNum = fpsub(this.BigNum, result.BigNum, NUM_DWORDS)End OperatorOperator BigFloat.-= (Byval rhs As Double)   Dim As BigFloat result   result.BigNum=str2fp( str(rhs) )   this.BigNum = fpsub(this.BigNum, result.BigNum, NUM_DWORDS)End OperatorOperator BigFloat.-= (Byref rhs As String)   Dim As BigFloat result   result.BigNum=str2fp(rhs)   this.BigNum = fpsub(this.BigNum, result.BigNum, NUM_DWORDS)End OperatorOperator BigFloat.*= (Byref rhs As BigFloat)   this.BigNum = fpmul(this.BigNum, rhs.BigNum, NUM_DWORDS)End OperatorOperator BigFloat.*= (Byval rhs As Long)   Dim As BigFloat result   result.BigNum=si2fp(rhs, NUM_DWORDS)   this.BigNum = fpmul(this.BigNum, result.BigNum, NUM_DWORDS)End OperatorOperator BigFloat.*= (Byval rhs As Double)   Dim As BigFloat result   result.BigNum=str2fp( str(rhs) )   this.BigNum = fpmul(this.BigNum, result.BigNum, NUM_DWORDS)End OperatorOperator BigFloat.*= (Byref rhs As String)   Dim As BigFloat result   result.BigNum=str2fp(rhs)   this.BigNum = fpmul(this.BigNum, result.BigNum, NUM_DWORDS)End OperatorOperator BigFloat./= (Byref rhs As BigFloat)   this.BigNum = fpdiv(this.BigNum, rhs.BigNum, NUM_DWORDS)End OperatorOperator BigFloat./= (Byval rhs As Long)   Dim As BigFloat result   result.BigNum=si2fp(rhs, NUM_DWORDS)   this.BigNum = fpdiv(this.BigNum, result.BigNum, NUM_DWORDS)End OperatorOperator BigFloat./= (Byval rhs As Double)   Dim As BigFloat result   result.BigNum=str2fp( str(rhs) )   this.BigNum = fpdiv(this.BigNum, result.BigNum, NUM_DWORDS)End OperatorOperator BigFloat./= (Byref rhs As String)   Dim As BigFloat result   result.BigNum=str2fp(rhs)   this.BigNum = fpdiv(this.BigNum, result.BigNum, NUM_DWORDS)End OperatorOperator = ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long   Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)=0End OperatorOperator < ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long   Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)<0End OperatorOperator > ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long   Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)>0End OperatorOperator <= ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long   Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)<=0End OperatorOperator >= ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long   Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)>=0End OperatorOperator <> ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long   Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)<>0End OperatorOperator Abs(Byref rhs As BigFloat) As BigFloat   Dim As BigFloat result   result.BigNum=rhs.BigNum   result.BigNum.sign = 0   Operator = resultEnd OperatorOperator  Fix (Byref x As BigFloat) As BigFloat   Dim As BigFloat result   result.BigNum =fpfix(x.BigNum)   Operator = resultEnd OperatorOperator Frac(Byref x As BigFloat) As BigFloat   Dim As BigFloat result   result.BigNum=fpfrac(x.BigNum)   Operator = resultEnd OperatorDim Shared As BigFloat_struct Log102Log102.sign=0Log102.exponent=(-2)+BIAS+1If NUM_DWORDS>=0 Then Log102.mantissa[0]=&h13441350If NUM_DWORDS>=1 Then Log102.mantissa[1]=&h9F79FEF3If NUM_DWORDS>=2 Then Log102.mantissa[2]=&h11F12B35If NUM_DWORDS>=3 Then Log102.mantissa[3]=&h816F922FPrivate Function log2_32(Byval value As Ulong) As Long'https://stackoverflow.com/questions/11376288/fast-computing-of-log2-for-64-bit-integers   Static tab32(0 To 31) As Const Long = {0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30, 8, 12, 20, 28, 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31}   value Or= Culng(value Shr 1)   value Or= Culng(value Shr 2)   value Or= Culng(value Shr 4)   value Or= Culng(value Shr 8)   value Or= Culng(value Shr 16)   Return tab32(Culng(Culng(value * &h07C4ACDD) Shr 27))End FunctionPrivate Function log2_64(Byval value As Ulongint) As Long'https://stackoverflow.com/questions/11376288/fast-computing-of-log2-for-64-bit-integers   Static tab64(0 To 63) As Const Long ={63, 0, 58, 1, 59, 47, 53, 2, 60, 39, 48, 27, 54, 33, 42, 3, 61, 51, 37, 40, 49, 18, 28, 20, 55, 30, 34, 11, 43, 14, 22, 4, 62, 57, 46, 52, 38, 26, 32, 41, 50, 36, 17, 19, 29, 10, 13, 21, 56, 45, 25, 31, 35, 16, 9, 12, 44, 24, 15, 8, 23, 7, 6, 5}   value Or= Culngint(value Shr 1)   value Or= Culngint(value Shr 2)   value Or= Culngint(value Shr 4)   value Or= Culngint(value Shr 8)   value Or= Culngint(value Shr 16)   value Or= Culngint(value Shr 32)   Return tab64((Culngint((value - (value Shr 1)) * &h07EDD5E59A4E28C2) Shr 58))End FunctionFunction ipower (Byval x As Ulongint, Byval e As Ulongint) As Ulongint   'take x to an integer power   Dim As Ulongint z, y, n   y = x   n = e   z = 1   While n > 0      While (n And 1) = 0         n = n \ 2         y = y * y      Wend         n = n - 1         z = y * z   Wend   Return zEnd FunctionFunction shl32(Byval n As Ulong, Byval k As Ubyte, Byref c As Ulong) As Ulong   If k>0 And k<32 Then      Dim As Ulong carry=n      Dim As Ubyte k32=32-k      Asm         mov cl, [k32]         Shr dword Ptr [carry], cl         mov cl, [k]         Shl dword Ptr [n], cl      End Asm      c=carry   End If   Return nEnd FunctionFunction shr32(Byval n As Ulong, Byval k As Ubyte, Byref c As Ulong) As Ulong   If k>0 And k<32 Then      Dim As Ulong carry=n      Dim As Ubyte k32=32-k      Asm         mov cl, [k32]         Shl dword Ptr [carry], cl         mov cl, [k]         Shr dword Ptr [n], cl      End Asm      c=carry   End If   Return nEnd FunctionSub shiftl(Byref fac1 As BigFloat_struct, Byval k As Long, Byval dwords As Long=NUM_DWORDS)   If k>0 And k<32 Then      Dim As Long i      Dim As Ulong n, carry, c=0      Dim As Long k32=32-k      For i=dwords To 0 Step -1         n=fac1.mantissa[i]         carry=n         carry=Cast(Ulong, (carry Shr k32))         n=Cast(Ulong, (n Shl k))         fac1.mantissa[i]=n+c         c=carry      Next   Elseif k=32 Then      Dim As Long i      For i=0 To dwords-1         fac1.mantissa[i]=fac1.mantissa[i+1]      Next      fac1.mantissa[dwords]=0   End IfEnd SubSub shiftr(Byref fac1 As BigFloat_struct, Byval k As Long, Byval dwords As Long=NUM_DWORDS)   If k>0 And k<32 Then      Dim As Long i      Dim As Ulong n, carry, c=0      Dim As Long k32=32-k      For i=0 To dwords         n=fac1.mantissa[i]         carry=n         carry=Cast(Ulong, (carry Shl k32))         n=Cast(Ulong, (n Shr k))            fac1.mantissa[i]=c+n         c=carry      Next   Elseif k=32 Then      Dim As Long i      For i=dwords To 1 Step -1         fac1.mantissa[i]=fac1.mantissa[i-1]      Next      fac1.mantissa[0]=0   End IfEnd SubFunction fpcmp(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As Long   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS      Dim As Long c, i      If x.sign < y.sign Then fpcmp = -1: Exit Function      If x.sign > y.sign Then fpcmp = 1: Exit Function      If x.sign=y.sign Then         If x.exponent=y.exponent Then            For i=0 To dwords               c=x.mantissa[i]-y.mantissa[i]               If c<>0 Then Exit For            Next            If c<0 Then Return -1            If c=0 Then Return 0            If c>0 Then Return 1         End If         If x.exponent<y.exponent Then Return -1         If x.exponent>y.exponent Then Return 1      End IfEnd FunctionFunction NORM_FAC1(Byref fac1 As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As Integer   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   ' normalize the number in fac1   ' all routines exit through this one.   'see if the mantissa is all zeros.   'if so, set the exponent and sign equal to 0.   Dim As Integer i,er=0,f=0   For i=0 To dwords      If fac1.mantissa[i]>0 Then f=1   Next   If f=0 Then      fac1.exponent=0      fac1.sign=0      Exit Function   End If   While fac1.mantissa[0]=0      shiftl(fac1, 32, dwords)      fac1.exponent-=32      If fac1.exponent=0 Then         NORM_FAC1=EXPU_ERR         Exit Function      End If   Wend   i=__builtin_clz(fac1.mantissa[0])-3   If i>0 Then      shiftl(fac1, i, dwords)      fac1.exponent-=i   End If   'if the highmost Bit in fac1_man is nonzero,   'shift the mantissa right 1 Bit and   'increment the exponent         If fac1.mantissa[0]>(&h1ffffffful)Then      While fac1.mantissa[0]>&h1ffffffful      shiftr(fac1, 1, dwords)      fac1.exponent+=1      Wend   Elseif fac1.mantissa[0]<&h10000000ul Then      /' the following will probably never be executed '/      'now shift fac1_man 1 to the left until a      'nonzero Bit appears in the next-to-highest      'Bit of fac1_man.  decrement exponent for      'each shift.      While fac1.mantissa[0]<&h10000000ul         shiftl(fac1, 1, dwords)         fac1.exponent-=1      Wend   End If   'check for overflow/underflow   If fac1.exponent<0 Then   NORM_FAC1=EXPO_ERR   End IfEnd FunctionFunction si2fp(Byval m As Long, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   Dim As BigFloat_struct fac1   Dim As Long n=Abs(m)   If m=0 Then      Return fac1   End If   fac1.mantissa[0]=n   NORM_FAC1(fac1, NUM_DWORDS)   fac1.exponent=log2_32(n)+BIAS+1      If m<0 Then      fac1.sign=&H8000   Else      fac1.sign=0   End If   Return fac1End FunctionFunction ui2fp(Byval m As Ulong, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   Dim As BigFloat_struct fac1   If m=0 Then      Return fac1   End If   fac1.mantissa[0]=m   NORM_FAC1(fac1, NUM_DWORDS)   fac1.exponent=log2_32(m)+BIAS+1      Return fac1End FunctionSub fpadd_aux(Byref fac1 As BigFloat_struct, Byref fac2 As BigFloat_struct, Byval dwords As Long=NUM_DWORDS)   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   Dim As Long c, i   Dim As Longint v   c=0   For i=dwords To 1 Step -1      v=Clngint(fac2.mantissa[i])+Clngint(fac1.mantissa[i])+c      If v>(&h100000000ull-1) Then         v=v-&h100000000ull         c=1      Else         c=0      End If      fac1.mantissa[i]=v   Next   v=Clngint(fac1.mantissa[0])+Clngint(fac2.mantissa[0])+c   fac1.mantissa[0]=v   NORM_FAC1(fac1, dwords)End SubSub fpsub_aux(Byref fac1 As BigFloat_struct, Byref fac2 As BigFloat_struct, Byval dwords As Long=NUM_DWORDS)   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   Dim As Long c, i   Dim As Longint v   c=0   For i=dwords To 1 Step -1      v=Clngint(fac1.mantissa[i])-clngint(fac2.mantissa[i])-c      If v<0 Then         v=v+&h100000000ull         c=1      Else         c=0      End If      fac1.mantissa[i]=v   Next   v=Clngint(fac1.mantissa[0])-clngint(fac2.mantissa[0])-c   fac1.mantissa[0]=v   NORM_FAC1(fac1, dwords)End SubFunction fpadd(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   Dim As BigFloat_struct fac1,fac2   Dim As Long i, t, c, xsign, ysign   xsign=x.sign:x.sign=0   ysign=y.sign:y.sign=0   c=fpcmp(x, y, dwords)   x.sign=xsign   y.sign=ysign   If c<0 Then      fac1=y      fac2=x   Else      fac1=x      fac2=y   End If   t=fac1.exponent-fac2.exponent   If t<(NUM_BITS+32) Then      'The difference between the two      'exponents indicate how many times      'we have to multiply the mantissa      'of FAC2 by 10 (i.e., shift it right 1 place).      'If we have to shift more times than      'we have dwords, the result is already in FAC1.      If t>0 And t<(NUM_BITS+32) Then 'shift         i=t\32         While i>0            shiftr(fac2, 32, dwords)            t-=32            i-=1         Wend         'While t>0         If t>0 Then   shiftr(fac2, t, dwords)            't-=1            'Wend         End If         'See if the signs of the two numbers         'are the same.  If so, add; if not, subtract.         If fac1.sign=fac2.sign Then 'add            fpadd_aux(fac1,fac2, dwords)         Else            fpsub_aux(fac1,fac2, dwords)         End If   Endif   Return fac1End FunctionFunction fpsub(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   Dim As BigFloat_struct fac1,fac2   Dim As Long s   fac1=x   fac2=y   fac2.sign=fac2.sign Xor &h8000   fac1=fpadd(fac1,fac2, dwords)   Return fac1End FunctionFunction fpmul(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   Dim As BigFloat_struct fac1,fac2   Dim As Integer i, j, ex, er, den, num   Dim As Ulongint digit, carry, tmp, prod   Dim As Ulong  fac3(0 To 2*dwords+1)   fac1=x   fac2=y   'check exponents.  if either is zero,   'the result is zero   If fac1.exponent=0 Or fac2.exponent=0 Then 'result is zero...clear fac1.      fac1.sign=0      fac1.exponent=0      For i=0 To dwords         fac1.mantissa[i]=0      Next      'NORM_FAC1(fac1)      Return fac1   Else      If ex<0 Then         er=EXPO_ERR         Return fac1 'Exit Function      End If      'clear fac3 mantissa      For i=0 To dwords         fac3(i)=0      Next      den=dwords      While fac2.mantissa[den]=0         den-=1      Wend      num=dwords      While fac1.mantissa[num]=0         num-=1      Wend      If num<den Then         Swap fac1, fac2         'fac1=y         'fac2=x                  Swap den, num      End If      For j=den To 0 Step -1         carry=0         digit=fac2.mantissa[j]         For i=num To 0 Step -1            /'   prod=fac3(i+j+1)+digit*fac1.mantissa[i]+carry '/            tmp=fac1.mantissa[i]            If __builtin_umulll_overflow (digit, tmp, prod) Then Print "mul overflow"            tmp=fac3(i+j+1)            If __builtin_uaddll_overflow (tmp, prod, prod) Then Print "add overflow"            If __builtin_uaddll_overflow (carry, prod, prod) Then Print "add overflow"            carry=prod   Shr 32'\&H100000000            fac3(i+j+1)=prod '(prod mod &H100000000)         Next         fac3(j)=carry      Next      For i=0 To dwords         fac1.mantissa[i]=fac3(i)      Next         End If   'now determine exponent of result.   'as you do...watch for overflow.   'ex=fac2.exponent-BIAS+fac1.exponent   'fac1.exponent=ex   ex=(fac2.exponent And &h7FFFFFFFul)-BIAS+1   ex=ex+(fac1.exponent And &h7FFFFFFFul)-BIAS+1   fac1.exponent=ex+BIAS+1   'determine the sign of the product   fac1.sign=fac1.sign Xor fac2.sign   NORM_FAC1(fac1, dwords)   Return fac1End FunctionFunction fpmul_si(Byref x As BigFloat_struct, Byval y As Long, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   Dim As BigFloat_struct fac1,fac2   fac1=x   fac2=ui2fp(y, dwords)   fac1=fpmul(fac1, fac2, dwords)   Return fac1End FunctionFunction fpdiv_si(Byref num As BigFloat_struct, Byval den As Long, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   Dim As BigFloat_struct fac1   Dim As Ulongint carry, remder=0   Dim As Longint i, divisor   Dim As Longint quotient   Dim As Long j   divisor=Abs(den)   fac1=num   If divisor=1 Then      Return fac1   End If   If divisor = 0 Then      Print "error: divisor = 0"      Return fac1 'Exit function   End If   If divisor>2147483647 Then      Print "error: divisor too large"      Return fac1 'Exit function   End If   For i = 0 To dwords      quotient = fac1.mantissa[i] + remder * &h100000000ull      remder = quotient Mod divisor      fac1.mantissa[i]=quotient \ divisor   Next   quotient = remder * &h100000000ull   quotient=quotient \ divisor   carry=fac1.mantissa[0]   j=__builtin_clz(fac1.mantissa[0])-3   If j>0 Then      shiftl(fac1, j, dwords)      fac1.exponent-=j   End If   While fac1.mantissa[0]>(&h1ffffffful)      shiftr(fac1, 1, dwords)      fac1.exponent+=1   Wend   NORM_FAC1(fac1)   If den<0 Then   fac1.sign=fac1.sign Xor &h8000   End If   Return fac1End FunctionFunction fpdiv(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   Dim As BigFloat_struct fac1, fac2, one   Dim As Integer i, er, is_power_of_two   Dim As Short sign   fac1=x   fac2=y   one.exponent=(BIAS+1)   one.mantissa[0]=&h10000000ul   one=ui2fp(1)   sign=fac2.sign   fac2.sign=0   If fpcmp(fac2, one)=0 Then      Return fac1   Else      fac2.sign=sign   End If   If fac2.exponent=0 Then ' if fac2 = 0, return      ' a divide-by-zero error and      ' bail out.      fac1.mantissa[i]=&H1FFFFFFFul      For i=1 To dwords         fac1.mantissa[i]=&HFFFFFFFFul      Next      fac1.exponent=&HFFFFF+BIAS+1      er=DIVZ_ERR      Return fac1   Elseif fac1.exponent=0 Then 'fact1=0, just return      er=0      Return fac1   Else      'check to see if fac2 is a power of ten      is_power_of_two=0      If fac2.mantissa[0]=&H10000000ul Then         is_power_of_two=1         For i=1 To dwords            If fac2.mantissa[i]<>0 Then               is_power_of_two=0               Exit For            End If         Next      End If      'if fac2 is a power of ten then all we need to do is to adjust the sign and exponent and we are finished      If is_power_of_two=1 Then         fac1.sign=fac1.sign Xor fac2.sign         fac1.exponent=fac1.exponent-fac2.exponent+BIAS+1         Return fac1      End If      #Define min(a,b) Iif(((a)<(b)),(a),(b))      #Macro realw(w, j)         ((w(j - 1)*b + w(j))*b + w(j + 1))*b +Iif(Ubound(w)>=j+2,w(j+2),0)      #Endmacro      #Macro subtract(w, q, d, ka, kb)         For j=ka To kb            w(j) = w(j) - q*d(j - ka + 2)         Next      #Endmacro      #Macro normalize(w, ka, q)         w(ka) = w(ka) + w(ka - 1)*b         w(ka - 1) = q      #Endmacro      #Macro finalnorm(w, kb)         For j=kb To 3 Step -1            carry=Iif(w(j)<0, ((-w(j) - 1)\b) + 1, Iif(w(j) >= b, -(w(j)\b), 0))            w(j) = w(j) + carry*b            w(j - 1) = w(j - 1) - carry         Next      #Endmacro      Dim As Double result(1 To 2*dwords+3), n(1 To 2*dwords+3), d(1 To 2*dwords+3)      Const b=&H10000         Dim As Long j, last, laststep, q, t      Dim As Long stp, carry      Dim As Double xd, xn, rund      Dim As Double w(1 To Ubound(n)+4)      For j=0 To dwords         n(2*j+2)=fac1.mantissa[j]\&H10000         n(2*j+3)=fac1.mantissa[j] Mod &H10000         d(2*j+2)=fac2.mantissa[j]\&H10000         d(2*j+3)=fac2.mantissa[j] Mod &H10000      Next      n(1)=(fac1.exponent And &h7FFFFFFF)-BIAS-1      d(1)=(fac2.exponent And &h7FFFFFFF)-BIAS-1      For j=Ubound(n) To Ubound(w)         w(j)=0      Next      t=Ubound(n)-1      w(1)=n(1)-d(1)+1      w(2)=0      For j=2 To Ubound(n)         w(j+1)=n(j)      Next         xd = (d(2)*b + d(3))*b + d(4) + d(5)/b      laststep = t + 2      For stp=1 To laststep         xn=RealW(w, (stp + 2))         q=Int(xn/xd)         last = Min(stp + t + 1, Ubound(W))         subtract(w, q, d, (stp + 2), last)         normalize(w, (stp + 2), q)      Next      FinalNorm(w, (laststep + 1))      laststep = Iif(w(2) = 0, laststep, laststep - 1)      rund = w(laststep + 1)/b      w(laststep) = w(laststep) + Iif(rund >= 0.5, 1, 0)      If w(2)=0 Then         For j=1 To t+1            result(j)=w(j+1)         Next      Else         For j=1 To t+1            result(j)=w(j)         Next      End If      result(1) = Iif(w(2) = 0, w(1) - 1, w(1))      For j=0 To dwords         fac1.mantissa[j]=result(2*j+2)*&H10000+result(2*j+3)      Next      NORM_FAC1(fac1, dwords)      fac1.exponent=(result(1)+BIAS)   End If   fac1.sign=fac1.sign Xor fac2.sign   Return fac1End FunctionFunction fpipow(Byref x As BigFloat_struct, Byval e As Longint, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   'take x to an Long power   Dim As BigFloat_struct y=x   Dim As BigFloat_struct z, one   Dim As Longint n, c=0   Dim As Integer i   n = Abs(e)   one.exponent=(BIAS+1)   one.mantissa[0]=&h10000000ul   z=one   While n>0      While (n And 1)=0         n\=2         y=fpmul(y, y, dwords)         c+=1      Wend      n-=1      z=fpmul(y, z, dwords)      c+=1   Wend   If e<0 Then      z=fpdiv(one, z, dwords)   End If   Return zEnd FunctionFunction dbl2fp(Byval x As Double) As BigFloat_struct   Dim As BigFloat_struct bf   Dim As Ulongint n, e   n=Peek ( Ulongint, @x )   e=n   n=Cast(Ulongint, (n Shl 12))   n=Cast(Ulongint, (n Shr 4))   n+=&h1000000000000000ull   bf.mantissa[0]=Cast(Ulongint, (n Shr 32))   n=e   n=Cast(Ulongint, (n Shl 40))   n=Cast(Ulongint, (n Shr 32))   bf.mantissa[1]=n   n=n+&h1000000000000000ull   e=Cast(Ulongint, (e Shr 52))   e=e-&h3FFul   bf.exponent=e+BIAS+1   Return bfEnd FunctionFunction fp2dbl(Byref x As BigFloat_struct) As Double   Dim As Double dbl   Dim As Ulongint n, e   Dim As Ushort ex   Dim As Long sign   If x.sign<>0 Then      sign=-1   Else      sign=1   End If   If x.exponent>0 Then      ex=(x.exponent And &h7FFFFFFF)-BIAS-1   Else      ex=0   End If   n=x.mantissa[0]   n-=&h10000000ul   n=Cast(Ulongint, (n Shl 24))   e=&h3FF+ex   e=e*&h10000000000000ull   n=n+e   e=x.mantissa[1]\&h100   n=n+e   Poke Ulongint,@dbl,n   Return dbl*signEnd FunctionFunction factorial(Byval n As Ulong) As BigFloat_struct   Dim As BigFloat_struct one, f, indx   one.mantissa[0]=&h10000000ul   one.exponent=0+BIAS+1   f=one   indx=one   For i As Long=1 To n      f=fpmul(f, indx)      indx=fpadd(one, indx)   Next   Return fEnd FunctionFunction factorial_inv(Byval n As Ulong) As BigFloat_struct   Dim As BigFloat_struct f   f.mantissa[0]=&h10000000ul   f.exponent=0+BIAS+1   For i As Long=1 To n      f=fpdiv_si(f, i)   Next   Return fEnd FunctionFunction fp2ui(Byref x As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As Ulong   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   Dim As BigFloat_struct fac1   Dim As Double n   fac1.exponent=0   fac1.mantissa[0]=x.mantissa[0]   n=fp2dbl(x)   Return Fix(n)End FunctionFunction fp2si(Byref x As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As Long   If dwords>NUM_DWORDS Then dwords=NUM_DWORDS   Dim As BigFloat_struct fac1   Dim As Long sign   Dim As Double n   fac1.exponent=0   fac1.mantissa[0]=x.mantissa[0]   n=fp2dbl(x)   Return Fix(n)End FunctionFunction str2fp_aux(Byval value As String, Byref s As Long, Byref ex As Long ) As String   Dim As Integer j,d,e,ep,es,i,f,fp,fln   Dim As String c,f1,f2,f3, ts   Dim As Ulong ulng   j=1   s=1   d=0   e=0   ep=0   ex=0   es=1   i=0   f=0   fp=0   f1=""   f2=""   f3=""   value=Ucase(value)   fln=Len(value)   While j<=fln      c=Mid(value,j,1)      If ep=1 Then         If c=" " Then            j=j+1            Continue While         Endif         If c="-" Then            es=-es            c=""         Endif         If c="+" Then            j=j+1            Continue While         Endif         If (c="0") And (f3="") Then            j=j+1            Continue While         Endif         If (c>"/") And (c<":") Then 'c is digit between 0 and 9            f3=f3+c            ex=10*ex+(Asc(c)-48)            j=j+1            Continue While         Endif      Endif      If c=" " Then         j=j+1         Continue While      Endif      If c="-" Then         s=-s         j=j+1         Continue While      Endif      If c="+" Then         j=j+1         Continue While      Endif      If c="." Then         If d=1 Then            j=j+1            Continue While         Endif         d=1      Endif      If (c>"/") And (c<":") Then 'c is digit between 0 and 9         If ((c="0") And (i=0)) Then            If d=0 Then               j=j+1               Continue While            Endif            If (d=1) And (f=0) Then               e=e-1               j=j+1               Continue While            Endif         Endif         If d=0 Then            f1=f1+c            i=i+1         Else            If (c>"0") Then               fp=1            Endif            f2=f2+c            f=f+1         Endif      Endif      If c="E" Or c="D" Then      ep=1      Endif      j=j+1   Wend   If fp=0 Then      f=0      f2=""   Endif   If s=-1 Then s=&h8000 Else s=0   f1=f1+f2   f1=Mid(f1,1,1)+Right(f1,Len(f1)-1)   fln=Len(f1)   If fp=0 Andalso ex=0 Then      ex=es*fln-1   Else      ex=es*ex-1+i+e   End If   If fln>((NUM_BITS+32)*0.3010299956639811) Then      f1=Mid(f1,1,((NUM_BITS+32)*0.3010299956639811)-1)   Endif   Return f1End FunctionFunction str2fp(Byref x As String) As BigFloat_struct   Dim As Long strlen, i, n, ex, sign   Dim As String s2, strin=x   Dim As BigFloat_struct y, z, pw, ten   strin=str2fp_aux(strin, sign, ex)   strlen=Len(strin)   n=1   ten=ui2fp(Valuint("1"+String(n,"0")))   If strlen>n Then      s2=Left(strin, n)      strin=Mid(strin, n+1)   Else      z=ui2fp(Valuint(strin))      If ex<>0 Then         pw=ui2fp(10)         pw=fpipow(pw, strlen-ex-1)         z=fpdiv(z, pw)      End If      Return z   End If   z=ui2fp(Valuint(s2))   While Len(strin)>=n      s2=Left(strin, n)      strin=Mid(strin, n+1)      y=ui2fp(Valuint(s2))      z=fpmul(z, ten)      z=fpadd(z,y)   Wend   pw=ui2fp(10)   pw=fpipow(pw,strlen-ex-1)   z=fpdiv(z, pw)   Return zEnd FunctionFunction fp2str(Byref zz As BigFloat_struct, Byval digits As Long=NUM_DWORDS*9.63) As String   if digits>NUM_DWORDS*9.63 or digits<0 then digits=NUM_DWORDS*9.63   Dim As Long ex, i, n, powten, sign, tn   Dim As Long m   Dim As BigFloat_struct y, ten, z, ep   Dim As String s, s2   If zz.exponent=0 Then      Return "0"   End If   z=zz   z.sign=0   sign=zz.sign   If z.exponent<>0 Then      ex=(z.exponent And &h7FFFFFFFul)-BIAS-1   Else      ex=0   End If   n=NUM_BITS-ex   ep.exponent=-(n+15)+BIAS+1   ep.mantissa[0]=&h10000000ul   z=fpadd(z, ep)   If ex>0 Then      While ex>=3         z=fpdiv_si(z, 10)         powten+=1         ex=(z.exponent And &h7FFFFFFFul)-BIAS-1      Wend   Else      While ex<0         z=fpmul_si(z, 10)         powten-=1         ex=(z.exponent And &h7FFFFFFFul)-BIAS-1      Wend   End If   n=1   m=fp2ui(z)   If m=0 Then      z=fpmul_si(z, 10)      powten-=1      ex=(z.exponent And &h7FFFFFFFul)-BIAS-1      m=fp2ui(z)   End If   s=Trim(Str(m))+"."   If sign<>0 Then      s="-"+s   Else      s=" "+s   End If   ten=ui2fp(Clng("1"+String(n,"0")))   tn=Clng("1"+String(n,"0"))   For i=1 To digits      y=ui2fp(m)      z=fpsub(z,y)      z=fpmul(z,ten)      m=fp2ui(z)      If m<0 Or m>tn Then         m=0      End If      s2=Trim(Str(m))      s2=String(n-len(s2),"0")+s2      s=s+s2   Next   s2=Trim(Str(Abs(powten)))   s2=String(5-len(s2),"0")+s2   If powten<0 Then      s2="-"+s2   Else      s2="+"+s2   End If   s=s+"e"+s2   Return sEnd Function'integer part of numFunction fpfix( Byref num As BigFloat_struct ) As BigFloat_struct   Dim As BigFloat_struct ip   Dim As Long ex, ex2, j, k, c   ex=(num.exponent And &h7FFFFFFFul)-BIAS   If ex<1 Then      Return ip   End If   If ex>=(NUM_BITS) Then      Return num   End If   ex2=ex\32   k=ex2   j=ex Mod 32   While ex2>0      ex2-=1      ip.mantissa[ex2]=num.mantissa[ex2]   Wend   ip.mantissa[k]=num.mantissa[k]   ip.mantissa[k]=shr32(ip.mantissa[k], (32-j)-3, c)   ip.mantissa[k]=shl32(ip.mantissa[k], (32-j)-3, c)   ip.exponent=ex+BIAS   ip.sign=num.sign   Return ipEnd Function'fractional part of numFunction fpfrac( Byref num As BigFloat_struct ) As BigFloat_struct   Dim As BigFloat_struct ip, fp   ip=fpfix(num)   fp=fpsub(num, ip)   Return fpEnd FunctionDim As BigFloat x, x2, x3, y, z, term, fac, pi, epsDim As Long i, n, signDim As String s, spiDim As Double teps.BigNum.exponent=-(NUM_BITS)+BIAS+1eps.BigNum.mantissa[0]=&h10000000ulspi="3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821"pi=spiPrint "Pi+Pi"y=pi + piPrint yPrint "Pi/2"y=pi / 2Print yPrint "integer part of Pi/2"z=fix(y)Print zPrint "fractional part of Pi/2"z=frac(y)Print zPrint "Pi*Pi"z=pi * piPrint zPrint "factorial(10000)"z.BigNum=factorial(10000)Print zPrint "Pi^100"z=pi ^ 100Print zt=Timer'' sin(.5)x=".5"y=xx3=xx2=x * xterm=1fac=termsign=-1i=1While term>eps   x3*= x2   fac/= (i+1)   fac/= (i+2)   term=x3 * fac   If sign<0 Then      y=y - term   Else      y+= term   End If   i+=2   sign=-signWendt=timer-tPrint "sin(.5)"Print yPrint "time to comput sine function ";t;" seconds"PrintPrint "you can specify the number of digits to the toString function"Print pi.toString(16)Printfor j as bigfloat=1 to 20   ? j.toLong, (j/9).toString(16)nextPrint "press return to exit"Sleep`
srvaldez
Posts: 2850
Joined: Sep 25, 2005 21:54

I got the square root working (somewhat) -- not sure that all my routines are kosher because the square root of "1e100" + 1 was not not accurate to 1000 digits
I don't think it's worth pursuing
dodicat
Posts: 7146
Joined: Jan 10, 2006 20:30
Location: Scotland

Thanks srvaldez, nice routines.
I got a strange warning while testing , I got your factorial(10000) to check against my own, I kept your answer inside a hide macro to compare via fbide quick runs, the compiler finds my comment inside the macro and wants to truncate it?

Code: Select all

` #macro hide     '2.846259680917054518906413212119868890148051401702799230794179994274411340 #endmacro 'get warning 'warning 8(2): Literal number too big, truncated sleep  `
srvaldez
Posts: 2850
Joined: Sep 25, 2005 21:54