about bignum

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

about bignum

Postby srvaldez » Nov 22, 2021 0:23

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

Re: about bignum

Postby dodicat » Nov 24, 2021 19:20

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
             next


Print
Print shiftl("2021",2021)
print "done"
sleep


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

Re: about bignum

Postby srvaldez » Nov 24, 2021 19:51

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

Re: about bignum

Postby srvaldez » Nov 27, 2021 21:52

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

Re: about bignum

Postby srvaldez » Dec 07, 2021 4:13

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 Boolean
End Extern

Const NUMBER_OF_DIGITS = 96
Const 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
#Endif

Const NUM_DWORDS    =   NUM_BITS\32
Const 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 error

Type 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 mantissa
End Type

Constructor 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=0
End Constructor

Constructor 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]
   Next
End Constructor

Operator 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]
   Next
End Operator

Destructor bigfloat_struct ( )
   If this.mantissa=0 Then
      Print "unable to de-allocate memory"
      End(-4)
   Else
      Deallocate(this.mantissa)
   End If
   this.mantissa=0
End Destructor

Dim Shared As bigfloat_struct Log102
Log102.sign=0
Log102.exponent=(-2)+BIAS+1
If NUM_DWORDS>=0 Then Log102.mantissa[0]=&h13441350
If NUM_DWORDS>=1 Then Log102.mantissa[1]=&h9F79FEF3
If NUM_DWORDS>=2 Then Log102.mantissa[2]=&h11F12B35
If NUM_DWORDS>=3 Then Log102.mantissa[3]=&h816F922F

Private 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 Function

Private 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 Function

Function 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 z
End Function

Function 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 n
End Function

Function 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 n
End Function

Sub 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 If
End Sub

Sub 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 If
End Sub

Function 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 If
End Function

Function 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 If
End Function

Function 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 fac1
End Function

Function 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 fac1
End Function

Sub 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 Sub

Sub 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 Sub

Function 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 fac1
End Function

Function 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 fac1
End Function

Function 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 fac1
End Function

Function 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 fac1
End Function

Function 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 fac1
End Function


Function 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 fac1
End Function

Function 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 z
End Function

Function 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 bf
End Function

Function 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*sign
End Function

Function 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 dbl
End Function

Function 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 f
End Function

Function 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 f
End Function

Function 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 Function

Function 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 Function

Function 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 f1
End Function

Function 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 z
End Function

Function 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 s
End Function

'integer part of num
Function 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 ip
End Function

'fractional part of num
Function fpfrac( Byref num As bigfloat_struct ) As bigfloat_struct
   Dim As bigfloat_struct ip, fp
   ip=fpfix(num)
   fp=fpsub(num, ip)
   Return fp
End Function

Dim As bigfloat_struct x, x2, x3, y, z, term, fac, pi, eps
Dim As Long i, n, sign
Dim As String s, spi
Dim As Double t

eps.exponent=-(NUM_BITS)+BIAS+1
eps.mantissa[0]=&h10000000ul
spi="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=x
x3=x
x2=fpmul(x, x)
term=ui2fp(1)
fac=term
sign=-1
i=1
While 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=-sign
Wend
t=timer-t
Print "sin(.5)"
Print fp2str(y)
Print "time to comput sine function ";t;" seconds"
Print
Print "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

Re: about bignum

Postby srvaldez » Dec 07, 2021 19:52

with overloading

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 Boolean
End Extern

Const NUMBER_OF_DIGITS = 96
Const 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
#Endif

Const NUM_DWORDS    =   NUM_BITS\32
Const 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 error

Type 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 mantissa
End Type

Constructor 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=0
End Constructor

Constructor 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]
   Next
End Constructor

Operator 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]
   Next
End Operator

Destructor BigFloat_struct ( )
   If this.mantissa=0 Then
      Print "unable to de-allocate memory"
      End(-4)
   Else
      Deallocate(this.mantissa)
   End If
   this.mantissa=0
End Destructor

Type 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_struct

End Type

Declare Function fpcmp(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As Long
Declare Function fpadd(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fpsub(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fpmul(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fpmul_si(Byref x As BigFloat_struct, Byval y As Long, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fpdiv_si(Byref num As BigFloat_struct, Byval den As Long, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fpdiv(Byref x As BigFloat_struct, Byref y As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fpipow(Byref x As BigFloat_struct, Byval e As Longint, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function dbl2fp(Byval x As Double) As BigFloat_struct
Declare Function fp2dbl(Byref x As BigFloat_struct) As Double
Declare Function ui2fp(Byval m As Ulong, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fp2ui(Byref x As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As Ulong
Declare Function si2fp(Byval m As Long, Byval dwords As Long=NUM_DWORDS) As BigFloat_struct
Declare Function fp2si(Byref x As BigFloat_struct, Byval dwords As Long=NUM_DWORDS) As Long
Declare Function str2fp(Byref x As String) As BigFloat_struct
Declare Function fp2str(Byref zz As BigFloat_struct, Byval digits As Long=NUM_DWORDS*9.63) As String
Declare Function fpfix( Byref num As BigFloat_struct ) As BigFloat_struct
Declare Function fpfrac( Byref num As BigFloat_struct ) As BigFloat_struct

Function Bigfloat.toString( byval places as long=NUM_DWORDS*9.63 ) As String
   Function = fp2str( this.BigNum, places )
End Function

Function Bigfloat.toLong ( ) As Long
   dim as double x
   x=fp2dbl(this.BigNum)
   return clng(x)
End Function


Function Bigfloat.toDouble ( ) As Double
   Function = fp2dbl(this.BigNum)
End Function

Constructor Bigfloat ( )
   this.BigNum=si2fp(0, NUM_DWORDS)
End Constructor

Constructor Bigfloat ( Byval rhs As Long )
   this.BigNum=si2fp( rhs, NUM_DWORDS )
End Constructor

Constructor Bigfloat ( Byref rhs As String )
   this.BigNum=str2fp( rhs )
End Constructor

Constructor 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]
   next
End Constructor

Constructor Bigfloat ( Byval rhs As Double )
   this.BigNum=str2fp( str( rhs ) )
End Constructor

Operator 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]
   next
End Operator

Operator Bigfloat.let ( Byval rhs As Long )
   this.BigNum=si2fp( rhs, NUM_DWORDS )
End Operator

Operator Bigfloat.let ( Byref rhs As String )
   this.BigNum=str2fp( rhs )
End Operator

Operator Bigfloat.Let ( Byval rhs As Double )
   this.BigNum=str2fp( str(rhs) )
End Operator

Operator Bigfloat.cast ( ) As String
   Operator = fp2str(this.BigNum)
End Operator

Operator Bigfloat.cast ( ) As Long
   dim as double x
   x=fp2dbl(this.BigNum)
   Operator = clng(x)
End Operator

Operator Bigfloat.cast ( ) As Double
   Operator = fp2dbl(this.BigNum)
End Operator

Destructor 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)<=0
End 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 If
End 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 = result
End Operator

Operator + ( 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 = result
End Operator

Operator + ( 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 = result
End Operator

Operator + (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 = result
End Operator

Operator - ( Byref rhs As Bigfloat ) As Bigfloat
   Dim As Bigfloat result = rhs
   result.BigNum.sign = result.BigNum.sign xor &h8000
   Operator = result
End Operator

Operator - ( Byref lhs As Bigfloat, Byval rhs As Bigfloat ) As Bigfloat
   Dim As Bigfloat result
   result.BigNum= fpsub(lhs.BigNum, rhs.BigNum, NUM_DWORDS)
   Operator = result
End Operator

Operator - ( 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 = result
End Operator

Operator - ( 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 = result
End Operator


Operator - ( 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 = result
End Operator

Operator - ( 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 = result
End Operator

Operator - (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 = result
End Operator

Operator * ( Byref lhs As Bigfloat, Byval rhs As Bigfloat ) As Bigfloat
   Dim As Bigfloat result
   result.BigNum= fpmul(lhs.BigNum, rhs.BigNum, NUM_DWORDS)
   Operator = result
End Operator

Operator * ( 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 = result
End Operator

Operator * ( 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 = result
End Operator

Operator * ( 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 = result
End Operator

Operator * ( 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 = result
End Operator

Operator * (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 = result
End Operator

Operator / ( Byref lhs As Bigfloat, Byval rhs As Bigfloat ) As Bigfloat
   Dim As Bigfloat result
   result.BigNum= fpdiv(lhs.BigNum, rhs.BigNum, NUM_DWORDS)
   Operator = result
End Operator

Operator / ( 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 = result
End Operator

Operator / ( 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 = result
End Operator

Operator / ( 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 = result
End Operator

Operator / ( 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 = result
End Operator

Operator / (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 = result
End Operator

Operator ^ ( Byref lhs As Bigfloat, Byval rhs As Longint ) As Bigfloat
   Dim As Bigfloat lhs2
   lhs2.BigNum=fpipow(lhs.BigNum, rhs, NUM_DWORDS)
   Operator = lhs2
end operator

Operator BigFloat.+= (Byref rhs As BigFloat)
   this.BigNum = fpadd(this.BigNum, rhs.BigNum, NUM_DWORDS)
End Operator

Operator 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 Operator

Operator BigFloat.+= (Byval rhs As Double)
   Dim As BigFloat result
   result.BigNum=str2fp( str(rhs) )
   this.BigNum = fpadd(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator

Operator BigFloat.+= (Byref rhs As String)
   Dim As BigFloat result
   result.BigNum=str2fp(rhs)
   this.BigNum = fpadd(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator

Operator BigFloat.-= (Byref rhs As BigFloat)
   this.BigNum = fpsub(this.BigNum, rhs.BigNum, NUM_DWORDS)
End Operator

Operator 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 Operator

Operator BigFloat.-= (Byval rhs As Double)
   Dim As BigFloat result
   result.BigNum=str2fp( str(rhs) )
   this.BigNum = fpsub(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator

Operator BigFloat.-= (Byref rhs As String)
   Dim As BigFloat result
   result.BigNum=str2fp(rhs)
   this.BigNum = fpsub(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator

Operator BigFloat.*= (Byref rhs As BigFloat)
   this.BigNum = fpmul(this.BigNum, rhs.BigNum, NUM_DWORDS)
End Operator

Operator 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 Operator

Operator BigFloat.*= (Byval rhs As Double)
   Dim As BigFloat result
   result.BigNum=str2fp( str(rhs) )
   this.BigNum = fpmul(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator

Operator BigFloat.*= (Byref rhs As String)
   Dim As BigFloat result
   result.BigNum=str2fp(rhs)
   this.BigNum = fpmul(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator

Operator BigFloat./= (Byref rhs As BigFloat)
   this.BigNum = fpdiv(this.BigNum, rhs.BigNum, NUM_DWORDS)
End Operator

Operator 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 Operator

Operator BigFloat./= (Byval rhs As Double)
   Dim As BigFloat result
   result.BigNum=str2fp( str(rhs) )
   this.BigNum = fpdiv(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator

Operator BigFloat./= (Byref rhs As String)
   Dim As BigFloat result
   result.BigNum=str2fp(rhs)
   this.BigNum = fpdiv(this.BigNum, result.BigNum, NUM_DWORDS)
End Operator

Operator = ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long
   Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)=0
End Operator

Operator < ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long
   Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)<0
End Operator

Operator > ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long
   Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)>0
End Operator

Operator <= ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long
   Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)<=0
End Operator

Operator >= ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long
   Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)>=0
End Operator

Operator <> ( Byref lhs As BigFloat, Byref rhs As BigFloat ) As Long
   Operator = fpcmp(lhs.BigNum, rhs.BigNum, NUM_DWORDS)<>0
End Operator

Operator Abs(Byref rhs As BigFloat) As BigFloat
   Dim As BigFloat result
   result.BigNum=rhs.BigNum
   result.BigNum.sign = 0
   Operator = result
End Operator

Operator  Fix (Byref x As BigFloat) As BigFloat
   Dim As BigFloat result
   result.BigNum =fpfix(x.BigNum)
   Operator = result
End Operator

Operator Frac(Byref x As BigFloat) As BigFloat
   Dim As BigFloat result
   result.BigNum=fpfrac(x.BigNum)
   Operator = result
End Operator

Dim Shared As BigFloat_struct Log102
Log102.sign=0
Log102.exponent=(-2)+BIAS+1
If NUM_DWORDS>=0 Then Log102.mantissa[0]=&h13441350
If NUM_DWORDS>=1 Then Log102.mantissa[1]=&h9F79FEF3
If NUM_DWORDS>=2 Then Log102.mantissa[2]=&h11F12B35
If NUM_DWORDS>=3 Then Log102.mantissa[3]=&h816F922F

Private 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 Function

Private 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 Function

Function 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 z
End Function

Function 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 n
End Function

Function 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 n
End Function

Sub 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 If
End Sub

Sub 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 If
End Sub

Function 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 If
End Function

Function 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 If
End Function

Function 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 fac1
End Function

Function 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 fac1
End Function

Sub 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 Sub

Sub 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 Sub

Function 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 fac1
End Function

Function 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 fac1
End Function

Function 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 fac1
End Function

Function 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 fac1
End Function

Function 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 fac1
End Function


Function 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 fac1
End Function

Function 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 z
End Function

Function 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 bf
End Function

Function 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*sign
End Function

Function 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 f
End Function

Function 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 f
End Function

Function 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 Function

Function 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 Function

Function 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 f1
End Function

Function 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 z
End Function

Function 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 s
End Function

'integer part of num
Function 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 ip
End Function

'fractional part of num
Function fpfrac( Byref num As BigFloat_struct ) As BigFloat_struct
   Dim As BigFloat_struct ip, fp
   ip=fpfix(num)
   fp=fpsub(num, ip)
   Return fp
End Function

Dim As BigFloat x, x2, x3, y, z, term, fac, pi, eps
Dim As Long i, n, sign
Dim As String s, spi
Dim As Double t

eps.BigNum.exponent=-(NUM_BITS)+BIAS+1
eps.BigNum.mantissa[0]=&h10000000ul
spi="3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821"

pi=spi
Print "Pi+Pi"
y=pi + pi
Print y
Print "Pi/2"
y=pi / 2
Print y
Print "integer part of Pi/2"
z=fix(y)
Print z
Print "fractional part of Pi/2"
z=frac(y)
Print z
Print "Pi*Pi"
z=pi * pi
Print z
Print "factorial(10000)"
z.BigNum=factorial(10000)
Print z
Print "Pi^100"
z=pi ^ 100
Print z

t=Timer
'' sin(.5)
x=".5"
y=x
x3=x
x2=x * x
term=1
fac=term
sign=-1
i=1
While 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=-sign
Wend
t=timer-t
Print "sin(.5)"
Print y
Print "time to comput sine function ";t;" seconds"
Print
Print "you can specify the number of digits to the toString function"
Print pi.toString(16)
Print
for j as bigfloat=1 to 20
   ? j.toLong, (j/9).toString(16)
next
Print "press return to exit"
Sleep
srvaldez
Posts: 2850
Joined: Sep 25, 2005 21:54

Re: about bignum

Postby srvaldez » Dec 08, 2021 20:52

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

Re: about bignum

Postby dodicat » Dec 08, 2021 22:10

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

Re: about bignum

Postby srvaldez » Dec 08, 2021 22:49

hi dodicat
thanks for testing, about your macro warning it's strange to say the lest
about the math routines being kosher, I just printed out the value of sqr(1e100+1) in hex and it's ok, it's my bigfloat to string that's at fault
I did a torture test of sorts by calculating Pi to 715000 places using the Brent_Salamin algorithm and they agree with what Maple produces but the same calculation is at least 60% faster using my decfloat routines

Return to “General”

Who is online

Users browsing this forum: No registered users and 5 guests