## FBMath Update

User projects written in or related to FreeBASIC.
jdebord
Posts: 529
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

### FBMath Update

The FreeBASIC Math library (FBMath) has been updated.

* Added: Two hypergeometric functions, 1F1 and 2F1
* Added: Numerical inversion of Laplace transforms

The new version is here:

http://sourceforge.net/projects/fbmath/ ... mat065.zip
dodicat
Posts: 6727
Joined: Jan 10, 2006 20:30
Location: Scotland

### Re: FBMath Update

Thanks for the update.

Bye the way jdebord, in case you are not aware, the latest freebasic builds (unofficial yet) can redim UDT arrays.
This means you can make up general Matrix (or vector) types very easily.

http://www.freebasic.net/forum/viewtopic.php?f=17&t=19095&start=165

I've (as an exercise) compared Your GaussJordan and a UDT GaussJordan of my own.
Using random Matrix elements:

Code: Select all

`' ******************************************************************' Solution of a system of linear equations by Gauss-Jordan method' ******************************************************************' ------------------------------------------------------------------' Error codes' ------------------------------------------------------------------#DEFINE MatOk 0  ' No error#DEFINE MatSing -2  ' Quasi-singular matrix#DEFINE MatErrDim -3 ' Non-compatible dimensions' ------------------------------------------------------------------' Machine-dependent constant' ------------------------------------------------------------------#DEFINE MachEp 2.220446049250313D-16  ' Floating point precision: 2^(-52)' ------------------------------------------------------------------' Global variable' ------------------------------------------------------------------COMMON SHARED ErrCode AS INTEGER' Error code from the latest function evaluation' ******************************************************************SUB GaussJordan (A() AS DOUBLE, BYREF Det AS DOUBLE)' ------------------------------------------------------------------' Gauss-Jordan algorithm for a matrix A(L..N, L..M) with M >= N' ------------------------------------------------------------------' On input:'   * The submatrix A(L..N, L..N) contains the system matrix'   * The submatrix A(L..N, (N+1)..M) contains the constant vector(s)'' On output:'   * The submatrix A(L..N, L..N) contains the inverse matrix'   * The submatrix A(L..N, (N+1)..M) contains the solution vector(s)'   * The determinant of the system matrix is returned in Det'   * The error code is returned in the global variable ErrCode:'       ErrCode = MatOk     ==> no error'       ErrCode = MatErrDim ==> non-compatible dimensions (N < M)'       ErrCode = MatSing   ==> quasi-singular matrix' ------------------------------------------------------------------  DIM AS INTEGER L, N, M  ' Bounds of A  DIM AS INTEGER I, J, K  ' Loop variables  DIM AS INTEGER Ik, Jk   ' Pivot coordinates  DIM AS DOUBLE  Pvt      ' Pivot  DIM AS DOUBLE  T        ' Auxiliary variable  L = LBOUND(A, 1)  N = UBOUND(A, 1)  M = UBOUND(A, 2)    IF N > M THEN    ErrCode = MatErrDim    EXIT SUB  END IF  DIM AS INTEGER PRow(L TO N)  ' Stores line of pivot  DIM AS INTEGER PCol(L TO N)  ' Stores column of pivot  DIM AS DOUBLE  MCol(L TO N)  ' Stores a column of the matrix  Det = 1  K = L  DO WHILE K <= N    ' Search for largest pivot in submatrix A[K..N, K..N]    Pvt = A(K, K)    Ik = K    Jk = K    FOR I = K TO N      FOR J = K TO N        IF ABS(A(I, J)) > ABS(Pvt) THEN          Pvt = A(I, J)          Ik = I          Jk = J        END IF      NEXT J    NEXT I    ' Pivot too small ==> quasi-singular matrix    IF ABS(Pvt) < MachEp THEN      Det = 0      ErrCode = MatSing      EXIT SUB    END IF    ' Save pivot position    PRow(K) = Ik    PCol(K) = Jk    ' Update determinant    Det = Det * Pvt    IF Ik <> K THEN Det = -Det    IF Jk <> K THEN Det = -Det    ' Exchange current row (K) with pivot row (Ik)    IF Ik <> K THEN      FOR J = L TO M        SWAP A(K, J), A(Ik, J)      NEXT J    END IF    ' Exchange current column (K) with pivot column (Jk)    IF Jk <> K THEN      FOR I = L TO N        SWAP A(I, K), A(I, Jk)      NEXT I    END IF    ' Store col. K of A into MCol and set this col. to 0    FOR I = L TO N      IF I <> K THEN        MCol(I) = A(I, K)        A(I, K) = 0      ELSE        MCol(I) = 0        A(I, K) = 1      END IF    NEXT I    ' Transform pivot row    FOR J = L TO M      A(K, J) = A(K, J) / Pvt    NEXT J    ' Transform other rows    FOR I = L TO N      IF I <> K THEN        T = MCol(I)        FOR J = L TO M          A(I, J) = A(I, J) - T * A(K, J)        NEXT J      END IF    NEXT I    K = K + 1  LOOP  ' Exchange lines of whole matrix  FOR I = N TO L STEP -1    Ik = PCol(I)    IF Ik <> I THEN      FOR J = L TO M        SWAP A(I, J), A(Ik, J)      NEXT J    END IF  NEXT I  ' Exchange columns of inverse matrix  FOR J = N TO L STEP -1    Jk = PRow(J)    IF Jk <> J THEN      FOR I = L TO N        SWAP A(I, J), A(I, Jk)      NEXT I    END IF  NEXT J  ErrCode = MatOkEND SUB'=======================  JDEBORD END ========================'======================= DODICAT START =======================Type matrix     Dim As Double element(Any,Any)    Declare Operator Cast() As String    Declare Function determinant() As Double    Declare Function GaussJordan(rhs As matrix ) As matrixEnd Type'matrix multiply(NOT USED HERE)Operator *(m1 As matrix,m2 As matrix) As matrixDim rows As Integer=Ubound(m1.element,1)Dim columns As Integer=Ubound(m2.element,2)If Ubound(m1.element,2)<> Ubound(m2.element,1)ThenPrint "Can't do"Exit OperatorEnd IfDim As matrix ansRedim ans.element(1 To rows,1 To columns)Dim rxc As DoubleFor r As Integer=1 To rows    For c As Integer=1 To columns        rxc=0        For k As Integer = 1 To Ubound(m1.element,2)            rxc=rxc+m1.element(r,k)*m2.element(k,c)        Next k        ans.element(r,c)=rxc    Next cNext rOperator= ansEnd Operator'rounding functionFunction round (a As Double,b As Integer) As Double     Var y = (Abs(a)-Int(Abs(a))) * (10 ^ b),i=Int(y):y-=i    If y >= .5 Then i+= 1     i /= (10 ^ b)    Var r = Int(Abs(a))+i    If a < 0 Then r = -r    Return rEnd Function'matrix printerOperator matrix.cast() As StringDim As String ans,commaFor a As Integer=1 To Ubound(this.element,1)    For b As Integer=1 To Ubound(this.element,2)        If b=Ubound(element,2) Then comma="" Else comma=" , "        'ans=ans+Str(round(element(a,b),12))+comma        ans=ans+str(element(a,b))    Next b    ans=ans+Chr(10)Next aOperator= ansEnd OperatorFunction matrix.determinant() As Double    Dim As Double det=1    Dim As Integer n=Ubound(this.element),sign=1    Dim As matrix b=(This)    #macro pivot(num)    For p1 As Integer  = num To n - 1        For p2 As Integer  = p1 + 1 To n              If Abs(b.element(p1,num))<Abs(b.element(p2,num)) Then                sign=-sign                For g As Integer=1 To n                    Swap b.element(p1,g),b.element(p2,g)                Next g            End If        Next p2    Next p1    #endmacro    For k As Integer=1 To n-1        pivot(k)        For row As Integer =k To n-1            If b.element(row+1,k)=0 Then Exit For            Var f=b.element(k,k)/b.element(row+1,k)            For g As Integer=1 To n                b.element((row+1),g)=(b.element((row+1),g)*f-b.element(k,g))/f            Next g        Next row        det=det*b.element(k,k)    Next k    Return sign*det*b.element(n,n)End Functionfunction matrix.GaussJordan(rhs As matrix) as matrix    Dim As Integer n=ubound(element)    Dim As matrix  ans=rhs    dim as matrix b=this,r=rhs    dim as matrix detb=this    #macro pivot(num)    For p1 As Integer  = num To n - 1        For p2 As Integer  = p1 + 1 To n              If Abs(b.element(p1,num))<Abs(b.element(p2,num)) Then                Swap r.element(p1,1),r.element(p2,1)                For g As Integer=1 To n                    Swap b.element(p1,g),b.element(p2,g)                Next g            End If        Next p2    Next p1    #endmacro    dim as double f    For k As Integer=1 To n-1        pivot(k)                      For row As Integer =k To n-1            If b.element(row+1,k)=0 Then Exit For            f=b.element(k,k)/b.element(row+1,k)            r.element(row+1,1)=r.element(row+1,1)*f-r.element(k,1)            For g As Integer=1 To n                b.element((row+1),g)=b.element((row+1),g)*f-b.element(k,g)            Next g        Next row    Next k    'back substitute     For z As Integer=n To 1 Step -1        ans.element(z,1)=r.element(z,1)/b.element(z,z)        For j As Integer = n To z+1 Step -1            ans.element(z,1)=ans.element(z,1)-(b.element(z,j)*ans.element(j,1)/b.element(z,z))        Next j    Next z    return ansEnd function'=============  TEST ======================================================screen 20do    clsdim as integer d=1+int(rnd*15)'set a dimensiondim as double a(1 to d,1 to d+1),detDim As matrix m,rhsRedim m.element(1 To d,1 To d)Redim rhs.element(1 To d,1 To 1)'FILL UP EVERYTHINGfor r as integer=1 to d    a(r,d+1)=rnd*5-rnd*5 'jdebord's     rhs.element(r,1)=a(r,d+1)'mine    for c as integer=1 to d        a(r,c)=rnd*50-rnd*50'jdebord's        m.element(r,c)=a(r,c)'mine    next cnext rGaussJordan(a(),det)'  FBMATHprint "SYSTEM DIMENSION ";dprint"FBMATH determinant "; detprint"   UDT determinant "; m.determinantprintprint "FBMATH solution "for n as integer=1 to d    print str(a(n,d+1))next nprintprint "UDT solution "print m.GaussJordan(rhs)printprint "Press a key or esc"sleeploop until inkey=chr(27)sleep `
jdebord
Posts: 529
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

### Re: FBMath Update

Is it possible to overload the () operator, so as to write A(i,j) instead of A.element(i,j) ?
fxm
Posts: 9994
Joined: Apr 22, 2009 12:46
Location: Paris suburbs, FRANCE

### Re: FBMath Update

At the opposite of C++, the operator '()' (parenthesis) cannot be overloaded in FB.

It is a pity because in C++, operator '()' is commonly overloaded with two parameters to index multidimensional arrays, or to retrieve a subset of a one dimensional array (returning all the elements from parameter 1 to parameter 2).
dodicat
Posts: 6727
Joined: Jan 10, 2006 20:30
Location: Scotland

### Re: FBMath Update

In the latest build the operator {} can be overloaded (see the changelog).
With this and some code by Gothon a while back, you could do a C type array but only in one dimension (say a vector)

Code: Select all

`Type Vector    element(any) As double    Declare Operator [](As integer) ByRef As double    Declare Operator Cast() As StringEnd TypeOperator Vector.[](I as integer) ByRef As double    operator= element(I)End OperatorOperator vector.cast() As StringDim As String ansFor a As Integer=lbound(element) To Ubound(element)        ans+=str(this[a])+chr(10)Next aOperator= ansEnd Operatordim as vector Vctredim Vct.element(1 to 10)for n as integer=1 to 10    Vct[n]=nnext nprint Vct'=========================================printredim preserve Vct.element(1 to 15)for n as integer=11 to 15    Vct[n]=n*n+.1next nprint Vctsleep `

As a silly workaround, just define some characters to substitute .element, or to save typing just make .element .e

Code: Select all

` type matrix    as double element(any,any)    #define __ .elementend typedim as matrix mredim m __(1 to 3,1 to 3) m __(1,2)=9 print m __(1,2) 'test out some __'s to see if there is a clash! print __fb_version__  var __=90 print __ 'seems OK sleep `
jdebord
Posts: 529
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

### Re: FBMath Update

Thanks. I think I will wait for the next evolution of FB. I don't mind using [] instead of () but I will miss the 2-dimensional arrays!
fxm
Posts: 9994
Joined: Apr 22, 2009 12:46
Location: Paris suburbs, FRANCE

### Re: FBMath Update

The parameter of operator '[]' can be not only numeric type but also any user-defined type:

Code: Select all

`Type Pt  Dim As Integer x, y  Declare Constructor (Byval x0 As Integer = 0, Byval y0 As Integer = 0)End TypeConstructor Pt (Byval x0 As Integer = 0, Byval y0 As Integer = 0)  This.x = x0 : This.y = y0End ConstructorType Vector  Dim As Double element(Any, Any)  Declare Operator [](Byref N As Pt) Byref As DoubleEnd TypeOperator Vector.[](Byref N As Pt) Byref As Double  Operator = This.element(N.x, N.y)End OperatorDim as Vector vRedim v.element(0 To 9, 0 To 9)v[Pt(1, 2)] = 12Print v[Pt(1, 2)]`
fxm
Posts: 9994
Joined: Apr 22, 2009 12:46
Location: Paris suburbs, FRANCE

### Re: FBMath Update

Another solution by calling operators '[]' in cascade on one-dimension arrays:

Code: Select all

`Type Vector1  Dim As Double element1(Any)  Declare Operator [](Byval index As Integer) Byref As DoubleEnd TypeOperator Vector1.[](Byval index As Integer) Byref As Double  Operator = This.element1(index)End OperatorType Vector2  Dim As Vector1 element2(Any)  Declare Operator [](Byval index As Integer) Byref As Vector1End TypeOperator Vector2.[](Byval index As Integer) Byref As Vector1  Operator = This.element2(index)End OperatorDim As Vector2 vRedim v.element2(0 To 9)For I As Integer = 0 To 9  Redim (v.element2(I).element1)(0 To 9)Next Iv[3][4] = 34   '' equivalent to: v.element2(3).element1(4) = 34Print v[3][4]  '' equivalent to: Print v.element2(3).element1(4)`
jdebord
Posts: 529
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

### Re: FBMath Update

Thanks again. A notation like A[i][j] for a matrix element could be convenient. It reminds me of what I did with Turbo Pascal, using pointers to arrays of pointers, ending with notations such as A^[i]^[j] which is not very handy!
dodicat
Posts: 6727
Joined: Jan 10, 2006 20:30
Location: Scotland

### Re: FBMath Update

I've messed around with the [][] notation for UDT matrices.
I have only processed a determinant.
Doesn't look much like the basic language anymore!!

YOU NEED THE GIT FB BUILD

Code: Select all

` Type vector  Dim As Double Columns(Any)  Declare Operator [](Byval index As Integer) Byref As Double  Declare Operator Cast() As StringEnd TypeOperator vector.[](Byval index As Integer) Byref As Double  Operator = Columns(index)End OperatorType matrix   Dim As vector Rows(Any)  as integer rs,cs  Declare Operator [](Byval index As Integer) Byref As vector  Declare Operator Cast() As String   Declare Function determinant() As DoubleEnd TypeOperator matrix.[](Byval index As Integer) Byref As vector  Operator = Rows(index)End Operatoroperator vector.cast() as stringdim as string ans,commaFor a As Integer=lbound(Columns) To Ubound(Columns)    If a=Ubound(Columns) Then comma="" Else comma=" , "    ans=ans+str(this[a])+commanext aoperator=ansend operator'matrix printerOperator matrix.cast() As StringDim As String ans,commaFor a As Integer=lbound(Rows) To Ubound(Rows)    print Rows(a)Next aOperator= ansEnd OperatorFunction matrix.determinant() As Double    Dim As Double det=1    if rs<>cs then print "Matrix must be square":exit function    Dim As Integer n=rs,sign=1    Dim As matrix b=This    #macro pivot(num)    For p1 As Integer  = num To n - 1        For p2 As Integer  = p1 + 1 To n              If Abs(b[p1][num])<Abs(b[p2][num]) Then                sign=-sign                For g As Integer=1 To n                    Swap b[p1][g],b[p2][g]                Next g            End If        Next p2    Next p1    #endmacro    For k As Integer=1 To n-1        pivot(k)        For row As Integer =k To n-1            If b[row+1][k]=0 Then Exit For            Var f=(b[k][k])/(b[row+1][k])            For g As Integer=1 To n                b[row+1][g]=((b[row+1][g])*f-(b[k][g]))/f            Next g        Next row        det=det*b[k][k]    Next k    Return sign*det*b[n][n]End Function#macro RedimMatrix(M,r,c)scope    dim as integer z    for z  = c:next z: M.cs=z-1    for z  = r:next z: M.rs=z-1Redim M.Rows(r)For I as integer  = r  Redim (M.Rows(I).Columns)(c)Next Iend scope#endmacro'=================================================Dim As matrix MATRedimMatrix(MAT,1 to 4,1 to 4)for a as integer=1 to MAT.rs    for b as integer=1 to MAT.cs        MAT[a][b]=rnd*5    next bnext aprint "MATRIX"print MATprint "Determinant"print MAT.determinantsleep`
Roland Chastain
Posts: 948
Joined: Nov 24, 2011 19:49
Location: France
Contact:

### Re: FBMath Update

Hello!

I have compilation errors with the largint demos:

C:\FreeBASIC\fbmath\demo\largeint\bench.o:fake:(.text+0x56): undefined reference to `LargeInit@8'
C:\FreeBASIC\fbmath\demo\largeint\bench.o:fake:(.text+0xef): undefined reference to `Prints@8'
C:\FreeBASIC\fbmath\demo\largeint\bench.o:fake:(.text+0x101): undefined reference to `Letf@8'
...
dodicat
Posts: 6727
Joined: Jan 10, 2006 20:30
Location: Scotland

### Re: FBMath Update

Hi Roland.
jdebord will probably tell you better the reason.
However I have noticed that the first sub and function of largeinteger.bas are incomplete.
I have simply filled them in.
Now largeinteger.bas can stand alone, well, it does here anyway.
Fibonacci example:

Code: Select all

`'Subject: Basic library for large-integer arithmetic.'Author : Sjoerd.J.Schaper'Update : 13-02-2013'Code   : FreeBasic 0.24.0' ref: ..\random\rng.bas SUB InitMTbyTime ()     randomize timer end sub ' Initializes the generator from the value of TIMER FUNCTION RanGen2 () AS DOUBLE     return rnd     end function' Returns a random number on [0,1)-real-intervalCONST Asize = 97279 '                  number array sizeCONST ErrTrap = 0 '                    optionally stop on errorCONST LB = 15, MB = &H8000& '          binary storage baseCONST LD = 4,  MD = 10000 '            decimal string base < MBCONST t0 = -1, t1 = -2 '               pointers to swap-registersCONST t2 = -3, t3 = -4CONST MH = MB \ 2, M1 = &H7FFF% '      bitmask MB - 1CONST M2 = M1 * MB, SB = NOT M1 '      borrow, sign bitDIM SHARED AS SHORT n(Asize), uj'n() is the largeint number array,'signs and lengths are at positions -1,'followed by (length + 1) largeint-wordsDIM SHARED i() AS INTEGER'the i()'s are number indicesDIM SHARED AS INTEGER Bezsw, Errsw '   switchesDIM SHARED AS INTEGER Lognr, Prmnr '   filehandlesSUB Errorh (ByRef f AS STRING, ByVal sw AS SHORT)   IF Errsw = 0 THEN      Errsw = sw      PRINT " Error: "; f      IF Lognr THEN PRINT #Lognr, " Error: "; f      IF ErrTrap THEN ERROR sw   END IFEND SUB' ***************************************************************************************'                            Initialization and assignment' ***************************************************************************************SUB Letf (ByVal p AS SHORT, ByVal c AS INTEGER)DIM AS INTEGER j = 0, c0 = ABS(c)   DO      n(i(p) + j) = c0 AND M1      c0 SHR= LB: j += 1 '             split DWord c   LOOP WHILE c0: n(i(p) + j) = 0   n(i(p) - 1) = (SGN(c) AND SB) OR j ' sign, lengthEND SUBSUB Rndf (ByVal p AS SHORT, ByVal k AS SHORT)DIM AS INTEGER j, r = ABS(k)DIM AS INTEGER t = r \ LB, m = r - t * LB   IF m THEN      t += 1   ELSE      m = LB   END IF   IF t > uj THEN      Errorh "overflow Sub Rndf", 6      EXIT SUB   END IF   FOR j = i(p) TO i(p) + t - 2      n(j) = INT(RanGen2 * MB)   NEXT j: r = 1   FOR j = 2 TO m: r SHL= 1: NEXT   n(i(p) + t - 1) = r + INT(RanGen2 * r)   IF t = 0 THEN t = 1: n(i(p)) = 0   n(i(p) - 1) = t: n(i(p) + t) = 0END SUBFUNCTION LibErr AS INTEGER   LibErr = Errsw: Errsw = 0END FUNCTIONSUB Work (ByRef f AS STRING)   f = TRIM\$(ENVIRON\$("LARGEINT"))   IF f = "" THEN f = CURDIR\$   IF RIGHT\$(f, 1) <> "\" THEN f += "\"END SUBFUNCTION LargeInit (ByVal ix AS SHORT, ByRef f AS STRING) AS INTEGERDIM s AS STRING, t AS INTEGERLargeInit = 0   Bezsw = 0: Errsw = 0   '   t = Asize \ (5 + ix) '       "5" for four swap-rows + row zero   IF t < 5 OR t > M1 THEN      Errorh "out of memory LargeInit fail", 4      EXIT FUNCTION   END IF   uj = t   REDIM i(-4 TO ix) AS INTEGER   FOR t = -4 TO ix      i(t) = 1 + (t + 4) * uj '        initialize indices      Letf t, 0 '                      set numbers to zero   NEXT t   '   s = SPACE\$(64): Work s   Lognr = 0: Prmnr = FREEFILE   OPEN s + "PrimFlgs.bin" FOR BINARY SHARED AS Prmnr   IF VARPTR(f) THEN      IF LEN(f) THEN         Lognr = FREEFILE         OPEN s + f + ".log" FOR OUTPUT AS Lognr      END IF   END IF   '   InitMTbyTimeLargeInit = uj: uj = uj - 2END FUNCTIONSUB Term   CLOSE 'all files, then terminateEND SUB' ***************************************************************************************'                                    Direct access' ***************************************************************************************FUNCTION Getl (ByVal p AS SHORT) AS INTEGER   Getl = n(i(p) - 1) AND M1END FUNCTIONFUNCTION Gets (ByVal p AS SHORT) AS INTEGER   Gets = n(i(p) - 1) AND SBEND FUNCTIONFUNCTION Getw (ByVal p AS SHORT, ByVal j AS SHORT) AS INTEGER   Getw = n(i(p) + j)END FUNCTIONSUB Setl (ByVal p AS SHORT, ByVal a AS SHORT)   n(i(p) - 1) = (n(i(p) - 1) AND SB) OR ABS(a)END SUBSUB Sets (ByVal p AS SHORT, ByVal a AS SHORT)   n(i(p) - 1) = (n(i(p) - 1) AND M1) OR (SGN(a) AND SB)END SUBSUB Setw (ByVal p AS SHORT, ByVal j AS SHORT, ByVal a AS SHORT)   n(i(p) + j) = ABS(a)END SUBSUB EnQ (ByVal q AS SHORT, ByVal p AS SHORT, ByVal k AS SHORT)DIM AS INTEGER j, o, lq = Getl(q), lx = lq + Getl(p)   IF lx < uj THEN      o = i(p) - lq      FOR j = lq TO lx - 1         n(i(q) + j) = n(o + j)      NEXT j      IF k = 0 THEN k = 1      n(i(q) + lx) = -ABS(k) '         end of record      n(i(q) - 1) = lx + 1   ELSE      Errorh "list overflow Sub EnQ", 6   END IFEND SUBFUNCTION ExQ (ByVal p AS SHORT, ByRef k AS SHORT, ByVal q AS SHORT, _ ByRef j AS SHORT) AS INTEGERDIM AS INTEGER s, o, lq = Getl(q)   IF j < lq THEN      ExQ = -1: o = i(p) - j - 1      DO WHILE j < lq         s = n(i(q) + j): j += 1         IF s < 0 THEN EXIT DO '       end of record         n(o + j) = s      LOOP      k = ABS(s): n(o + j) = 0      n(i(p) - 1) = o + j - i(p)   ELSE      ExQ = 0   END IFEND FUNCTION' ***************************************************************************************'                                 Copying and comparison' ***************************************************************************************SUB Dup (ByVal p AS SHORT, ByVal q AS SHORT)DIM AS INTEGER j, o = i(p) - i(q)   FOR j = i(q) - 1 TO i(q) + Getl(p)      n(j) = n(o + j)   NEXT jEND SUBSUB Swp (ByVal p AS SHORT, ByVal q AS SHORT)   SWAP i(p), i(q)END SUBFUNCTION Cmp (ByVal p AS SHORT, ByVal q AS SHORT) AS INTEGERDIM AS INTEGER lp = n(i(p) - 1), lq = n(i(q) - 1)DIM AS INTEGER j, o, s, x = SGN((lp AND SB) OR 1)   '   IF ((lp XOR lq) AND SB) = 0 THEN '  equal signs      lp AND= M1      s = lp - (lq AND M1)      IF s = 0 THEN '                  equal lengths         o = i(q) - i(p)         FOR j = i(p) + lp - 1 TO i(p) STEP -1            s = n(j) - n(o + j) '      L=>R check            IF s THEN EXIT FOR         NEXT j         IF s = 0 THEN x = 0      END IF      IF s < 0 THEN x = -x   END IFCmp = xEND FUNCTIONFUNCTION Isf (ByVal p AS SHORT, ByVal a AS SHORT) AS INTEGERDIM AS INTEGER s, r = n(i(p) - 1)   IF a = 0 THEN r = r AND M1 '        mask off sign(0)   s = (SGN(a) AND SB) OR 1Isf = r = s AND n(i(p)) = ABS(a)END FUNCTIONSUB Lftj (ByVal p AS SHORT, ByVal k AS SHORT)DIM AS INTEGER j, t, ip = i(p)   FOR j = k TO 1 STEP -1      IF n(ip + j) THEN EXIT FOR   NEXT j   t = j + 1: n(ip + t) = 0   n(ip - 1) = (n(ip - 1) AND SB) OR tEND SUB' ***************************************************************************************'                              Basic arithmetic functions' ***************************************************************************************SUB Subt (ByVal p AS SHORT, ByVal q AS SHORT)DIM AS INTEGER ix = i(p), lx = (n(ix - 1) AND M1) - 1DIM AS INTEGER im = i(q), lm = (n(im - 1) AND M1) - 1DIM AS INTEGER j, o, c, s = n(ix - 1) XOR n(im - 1)   '   IF s AND SB THEN '                  distinct signs: add      IF lx < lm THEN SWAP lm, lx: im = ix      IF lx >= uj THEN         Errorh "overflow Sub Add", 6         EXIT SUB      END IF      FOR j = im + lm + 2 TO im + lx        n(j) = 0: NEXT j      '      o = i(q) - i(p)      c = 0: FOR j = i(p) TO i(p) + lx         c = c + n(o + j) + n(j)         n(j) = c AND M1: c SHR= LB      NEXT j: n(j) = c      '   ELSE '                              subtract      s = lx - lm      IF s = 0 THEN '                  equal lengths         o = im - ix         FOR j = ix + lx TO ix STEP -1            s = n(j) - n(o + j)            IF s THEN lx = j - ix: EXIT FOR         NEXT j         IF s = 0 THEN            n(ix - 1) = 1: n(ix) = 0 ' unsigned zero p            n(ix + 1) = 0: EXIT SUB         END IF         lm = lx      END IF      '      IF s < 0 THEN '                  p:= -(q - p)         n(ix - 1) XOR= SB         SWAP lm, lx: SWAP im, ix      END IF      FOR j = im + lm + 2 TO im + lx        n(j) = 0: NEXT j      '      im -= i(p): ix -= i(p)      c = MB: FOR j = i(p) TO i(p) + lx         c += M2 - n(im + j) + n(ix + j)         n(j) = c AND M1: c SHR= LB      NEXT j      n(j) = c AND M1   END IF   Lftj p, lx + 1END SUBSUB Add (ByVal p AS SHORT, ByVal q AS SHORT)   n(i(q) - 1) XOR= SB   Subt p, q   n(i(q) - 1) XOR= SBEND SUBSUB Inc (ByVal p AS SHORT, ByVal a AS SHORT)DIM AS INTEGER j, b = SGN(Gets(p) OR 1) * aDIM AS INTEGER lp, c = CINT(n(i(p))) + b   IF c >= 0 AND c < MB THEN      n(i(p)) = c   ELSE      c = MB + b: lp = Getl(p)      FOR j = i(p) TO i(p) + lp '      add single word         c += M2 + n(j)         n(j) = c AND M1: c SHR= LB      NEXT j      IF c = MB THEN         Lftj p, lp      ELSE         n(i(p) - 1) XOR= SB '         negate         n(i(p)) = MB - n(i(p))         n(i(p) + 1) = 0      END IF   END IFEND SUBFUNCTION Ismin1 (ByVal p AS SHORT, ByVal m AS SHORT) AS INTEGER   Inc p, 1   Ismin1 = Cmp(p, m) = 0   Inc p, -1END FUNCTIONFUNCTION Divint (ByVal p AS SHORT, ByVal a AS SHORT) AS INTEGERDIM AS INTEGER j, c, c0, c2, lp = Getl(p) - 1   c = 0: c2 = ABS(a)   FOR j = i(p) + lp TO i(p) STEP -1      c0 = c SHL LB OR n(j) '          paste carry w/dividend      n(j) = c0 \ c2 '                 store quotient      c = c0 - n(j) * c2 '             remainder   NEXT j   Lftj p, lpDivint = cEND FUNCTIONFUNCTION EstQ (ByVal t AS INTEGER, ByVal c AS INTEGER) AS INTEGERDIM AS LONGINT p3 = CLNGINT(n(t)) SHL LB   p3 = (p3 + n(t - 1)) SHL LB + n(t - 2)   c = p3 \ c: IF c > M1 THEN c = M1EstQ = cEND FUNCTIONSUB Divd (ByVal p AS SHORT, ByVal d AS SHORT, ByVal q AS SHORT)DIM AS INTEGER iq = i(q) - 1, ip = i(p) - 1DIM AS INTEGER it = i(d) - 1, lt = n(it) AND M1DIM AS INTEGER l2, j, k, t, c, c0, c2   '   IF lt = 1 THEN '                    single-word divisor      IF Isf(d, 0) THEN         Errorh "div. by zero Sub Divd", 1         EXIT SUB      END IF      '      n(iq) = (n(ip) AND SB) OR 1      n(iq + 1) = Divint(p, n(it + 1))      n(iq + 2) = 0: SWAP i(q), i(p)      n(ip) XOR= n(it) AND SB   ELSE      '      n(iq) = n(ip) XOR n(it)      l2 = (n(ip) AND M1) - lt      IF l2 < 0 THEN         n(iq) = (n(iq) AND SB) OR 1 ' signed zero q         n(iq + 1) = 0: n(iq + 2) = 0         EXIT SUB      END IF      '      it += 1: lt -= 1 '               divisor prefix{2}      c2 = CINT(n(it + lt)) SHL LB + n(it + lt - 1)      iq += 1      '      ip += 1 - i(t2)      FOR j = i(t2) - 1 TO i(t2) + l2        n(j) = n(ip + j): NEXT j '     initial remainder      '      FOR k = l2 TO 0 STEP -1 '        long division         t = i(p) + k + lt + 1         c0 = EstQ(t, c2) '            estimate quotient digit k         '         ip = i(p) - it + k         t = i(t2) - it + k         DO WHILE c0            c = MB            FOR j = it TO it + lt '    subtract scaled divisor,               c += M2 - c0 * n(j) + n(ip + j)               n(t + j) = c AND M1: c SHR= LB            NEXT j            IF c + n(ip + j) = MB THEN               SWAP i(p), i(t2): EXIT DO            END IF '                   swap dividend and remainder            c0 -= 1 '                  or decrease quotient         LOOP         n(iq + k) = c0      NEXT k      Lftj p, lt      Lftj q, l2   END IFEND SUBSUB Ldiv (ByVal p AS SHORT, ByVal d AS SHORT)   Divd p, d, t1: SWAP i(p), i(t1)END SUBSUB Mult (ByVal p AS SHORT, ByVal q AS SHORT, ByVal r AS SHORT)   DIM AS INTEGER ip = i(p), lp = (n(ip - 1) AND M1) - 1DIM AS INTEGER iq = i(q), lq = (n(iq - 1) AND M1) - 1DIM AS INTEGER j, k, c, c0, ir = i(r), lx = lp + lq   IF lx >= uj THEN      Errorh "overflow Sub Mult", 6      EXIT SUB   END IF   IF lp < lq THEN SWAP ip, iq: SWAP lp, lq   n(ir - 1) = n(ip - 1) XOR n(iq - 1)   '   ip -= ir: c0 = n(iq)   c = 0: FOR j = ir TO ir + lp '      initialize destination      c += c0 * n(ip + j)      n(j) = c AND M1: c SHR= LB   NEXT j: n(j) = c   '   iq -= ir   FOR k = ir + 1 TO ir + lq '         multiply and add      ip -= 1: c0 = n(iq + k)      c = 0: FOR j = k TO k + lp         c += c0 * n(ip + j) + n(j)         n(j) = c AND M1: c SHR= LB      NEXT j: n(j) = c   NEXT k   Lftj r, lx + 1END SUBSUB Lmul (ByVal p AS SHORT, ByVal q AS SHORT)   Mult p, q, t2: SWAP i(p), i(t2)END SUBSUB Squ (ByVal p AS SHORT, ByVal q AS SHORT)DIM AS INTEGER j, o, k, c, c2, ip = i(p), iq = i(q)DIM AS INTEGER lp = (n(ip - 1) AND M1) - 1, lx = lp + lp   IF lx >= uj THEN      Errorh "overflow Sub Squ", 6      EXIT SUB   END IF   '   j = iq: n(iq - 1) = 0   FOR k = ip TO ip + lp '             initialize destination       c = CINT(n(k)) * n(k)       n(j) = c AND M1: c SHR= LB       n(j + 1) = c: j += 2   NEXT k   '   o = ip - iq   FOR k = 1 TO lp '                   add cross terms      o -= 1: iq += 1      c2 = CINT(n(ip + k)) SHL 1      c = 0: FOR j = iq TO iq + k - 1         c += c2 * n(o + j) + n(j)         n(j) = c AND M1: c SHR= LB      NEXT j      c += n(j)      n(j) = c AND M1      n(j + 1) += c SHR LB   NEXT k   Lftj q, lx + 1END SUBSUB Chs (ByVal p AS SHORT)   n(i(p) - 1) XOR= SB '               signed magnitudeEND SUBFUNCTION Absf (ByVal p AS SHORT) AS INTEGERDIM AS INTEGER s = n(i(p) - 1)   n(i(p) - 1) = s AND M1Absf = s AND SBEND FUNCTION' ***************************************************************************************'                                   Bit manipulation' ***************************************************************************************SUB Boolf (ByVal p AS SHORT, ByVal q AS SHORT, ByVal k AS SHORT)DIM AS INTEGER lx = Getl(p) - 1, t = k AND 3DIM AS INTEGER lm, im, j, o   IF t THEN      lm = Getl(q) - 1: im = i(q)      IF lx < lm THEN SWAP lx, lm: im = i(p)      IF t > 1 THEN         FOR j = im + lm + 2 TO im + lx           n(j) = 0: NEXT j      ELSE         lx = lm      END IF   END IF   '   o = i(q) - i(p)   SELECT CASE t   CASE 1      FOR j = i(p) TO i(p) + lx         n(j) AND= n(o + j)      NEXT j   CASE 2      FOR j = i(p) TO i(p) + lx         n(j) XOR= n(o + j)      NEXT j   CASE 3      FOR j = i(p) TO i(p) + lx         n(j) OR= n(o + j)      NEXT j   CASE ELSE      FOR j = i(p) TO i(p) + lx '      one's complement         n(j) = M1 AND NOT n(j) '      excepting sign bits      NEXT j   END SELECT   Lftj p, lxEND SUBSUB Lsft (ByVal p AS SHORT, ByVal r AS SHORT)DIM AS INTEGER ip = i(p), lp = n(ip - 1) AND M1DIM AS INTEGER k, m, j, c   IF n(ip) = 0 AND lp = 1 THEN EXIT SUB   '   IF r > 0 THEN      k = r \ LB: m = r - k * LB      IF m THEN         c = 0         FOR j = ip TO ip + lp            c SHR= LB            c OR= CINT(n(j)) SHL m            n(j) = c AND M1 '          store product         NEXT j         IF n(ip + lp) THEN lp += 1         n(ip + lp) = 0      END IF   ELSE      k = ABS(r)   END IF   '   IF k THEN      lp += k      IF lp > uj THEN         Errorh "overflow Sub Lsft", 6         EXIT SUB      END IF      FOR j = ip + lp TO ip + k STEP -1        n(j) = n(j - k): NEXT j '      ShL words      FOR k = j TO ip STEP -1        n(k) = 0: NEXT k   END IF   Setl p, lpEND SUBSUB Rsft (ByVal p AS SHORT, ByVal r AS SHORT)DIM AS INTEGER ip = i(p), lp = (n(ip - 1) AND M1) - 1DIM AS INTEGER k, m, j, c, c0   '   IF r > 0 THEN      k = r \ LB: m = r - k * LB      IF m THEN         c = 0         FOR j = ip + lp TO ip STEP -1            c0 = c SHL LB OR n(j)            c = c0 AND M1            n(j) = c0 SHR m AND M1 '   store quotient         NEXT j         IF n(ip + lp) = 0 AND lp THEN lp -= 1      END IF   ELSE      k = ABS(r)   END IF   '   IF k THEN      lp -= k      IF lp < 0 THEN Letf p, 0: EXIT SUB      FOR j = ip TO ip + lp        n(j) = n(j + k): NEXT j '      ShR words   END IF   n(ip + lp + 1) = 0   Setl p, lp + 1END SUBFUNCTION Odd (ByVal p AS SHORT) AS INTEGERDIM AS INTEGER a, lp, ip = i(p), t = 0Odd = 0   lp = n(ip - 1) AND M1   IF lp = 1 THEN      a = n(ip)      IF a = 0 THEN a = 1   ELSE      DO UNTIL n(ip + t)         t += 1      LOOP      a = n(ip + t): t *= LB   END IF   DO UNTIL a AND 1      a SHR= 1: t += 1   LOOP   IF t = 0 THEN EXIT FUNCTION   IF lp = 1 THEN      n(ip) = a   ELSE      Rsft p, t   END IFOdd = tEND FUNCTIONFUNCTION Bitl (ByVal p AS SHORT) AS INTEGERDIM AS INTEGER a, t = Getl(p) - 1   a = n(i(p) + t): t *= LB   WHILE a      a SHR= 1: t += 1   WENDBitl = tEND FUNCTION' ***************************************************************************************'                             Modular arithmetic functions' ***************************************************************************************SUB Moddiv (ByVal p AS SHORT, ByVal m AS SHORT)DIM s AS INTEGER   Divd p, m, t1   IF Gets(p) THEN '                   make positive residue      IF NOT Isf(p, 0) THEN         s = n(i(m) - 1)         n(i(m) - 1) = s OR SB '      -Abs(m)         Subt p, m: Sets m, s      ELSE         Sets p, 0      END IF   END IFEND SUBSUB Modbal (ByVal p AS SHORT, ByVal m AS SHORT)DIM AS INTEGER s, z, r   Divd p, m, t1   IF NOT Isf(p, 0) THEN      s = Absf(p): z = Absf(m)      Dup m, t1: Subt t1, p      r = Cmp(p, t1)      IF r = 1 THEN '                  Abs(2p) > m    SWAP i(p), i(t1) '            balance p mod m    s XOR= SB      ELSEIF r = 0 THEN    s = 0      END IF      Sets p, s: Sets m, z   ELSE      Sets p, 0   END IFEND SUBFUNCTION Mp2 (ByVal p AS SHORT, ByVal a AS SHORT) AS INTEGERDIM AS INTEGER m = n(i(p)) AND a - 1   IF Gets(p) THEN m = a - m AND a - 1Mp2 = mEND FUNCTIONSUB Modmlt (ByVal p AS SHORT, ByVal q AS SHORT, ByVal m AS SHORT)   Mult p, q, t2: SWAP i(p), i(t2)   Divd p, m, t1END SUBSUB Modsqu (ByVal p AS SHORT, ByVal m AS SHORT)   Squ p, t2: SWAP i(p), i(t2)   Divd p, m, t1END SUBSUB Euclid (ByVal p AS SHORT, ByVal q AS SHORT, ByVal d AS SHORT)DIM AS INTEGER s = Absf(q)   IF Isf(q, 0) THEN Dup p, d: EXIT SUB   IF Gets(p) THEN Moddiv p, q   '   Letf d, 1: Letf t0, 0   DO UNTIL Isf(p, 0)      SWAP i(p), i(q)      SWAP i(d), i(t0)      Divd p, q, t1 '                  Euclidean division & inversion      Mult t0, t1, t2: Subt d, t2 '    r_t = r_t-2 - q_t * r_t-1   LOOP   SWAP i(p), i(t0)   SWAP i(d), i(q)   '   IF Gets(p) AND NOT Bezsw THEN Add p, q   Sets q, sEND SUBFUNCTION Hp2 (ByVal p AS SHORT) AS INTEGERDIM AS INTEGER a = n(i(p) + Getl(p) - 1), a0 = MH   DO UNTIL a AND a0      a0 SHR= 1: LOOP '                bitscan L=>RHp2 = a0END FUNCTIONSUB Modpwr (ByVal p AS SHORT, ByVal q AS SHORT, ByVal m AS SHORT)DIM AS INTEGER j, k, s, a, a0, sw = NOT Isf(m, 0)   IF Isf(q, 0) THEN Letf p, 1: EXIT SUB   '   IF sw THEN '                        enable reduction mod m      s = Absf(m)      IF Gets(q) THEN         Euclid p, m, t3 '             handle negative exponent         IF NOT Isf(t3, 1) THEN            Errorh "impossible inverse Sub Modpwr", 1            EXIT SUB         END IF      ELSE         Moddiv p, m '                 initial reduction      END IF   ELSE      s = Absf(p)      IF (n(i(q)) AND 1) = 0 THEN s = 0   END IF   '   k = i(q) + Getl(q) - 1   FOR j = k TO i(q) STEP -1 '         L=>R binary exponentiation      a = n(j): a0 = MH '              unsigned bitvector q      IF j = k THEN         Dup p, t0         a0 = Hp2(q) SHR 1 '           handle highest set bit(q)      END IF      WHILE a0         Squ t0, t1         SWAP i(t0), i(t1) '           square t0         IF sw THEN Divd t0, m, t1 '   reduce Mod m         IF a AND a0 THEN '            read bit            Mult t0, p, t1 '           t0 times base p            SWAP i(t0), i(t1)            IF sw THEN Divd t0, m, t1         END IF         a0 SHR= 1      WEND   NEXT j   SWAP i(p), i(t0)   IF sw THEN      Sets m, s   ELSE      Sets p, s   END IFEND SUBSUB Powr (ByVal p AS SHORT, ByVal c AS INTEGER)   Letf t2, ABS(c): Letf t0, 0   Modpwr p, t2, t0END SUB' ***************************************************************************************'                             Number theoretic functions' ***************************************************************************************SUB Fctl (ByVal p AS SHORT, ByVal a AS SHORT)DIM AS INTEGER r, c, c0   Letf p, 1: r = 1   IF a < 3 THEN      IF a = 2 THEN Inc p, 1      EXIT SUB   END IF   DO      c = r      DO: r += 1         c0 = c: c *= r      LOOP UNTIL c > MB OR r > a      Letf t1, c0: Lmul p, t1 '        partial product   LOOP UNTIL r > aEND SUBSUB Gcd (ByVal p AS SHORT, ByVal q AS SHORT, ByVal d AS SHORT)   Dup p, d: Dup q, t0   DO UNTIL Isf(t0, 0)      SWAP i(t0), i(d)      Divd t0, d, t1   LOOP   Sets d, 0END SUBSUB Lcm (ByVal p AS SHORT, ByVal q AS SHORT, ByVal d AS SHORT)   Gcd p, q, d   IF Isf(d, 1) THEN      Mult p, q, d   ELSE      IF Isf(d, 0) THEN EXIT SUB      Dup p, t0: Divd t0, d, t1      Mult t1, q, d   END IFEND SUBSUB Bezout (ByVal p AS SHORT, ByVal q AS SHORT, ByVal d AS SHORT)   Dup p, t3: Bezsw = -1   Euclid p, q, d: Bezsw = 0   IF NOT Isf(d, 1) THEN Ldiv t3, d   Mult t3, p, t0 '                                p * p^-1)   Inc t0, -1: Chs t0 '                       (1 -   Divd t0, q, t3: SWAP i(q), i(t3) '  q^-1 =                / qEND SUBSUB Isqrt (ByVal p AS SHORT, ByVal q AS SHORT, ByVal r AS SHORT)DIM AS INTEGER t, iq, ir = i(r), ip = i(p)DIM AS INTEGER j, k, c, cx, c0, c2, lq, lr, lp = n(ip - 1)   IF lp AND SB THEN      Errorh "illegal argument Sub Isqrt", 5      EXIT SUB   END IF     IF (lp AND 1) = 0 THEN lp -= 1 '    even pwordlength   c = CINT(n(ip + lp)) SHL LB + n(ip + lp - 1)   c2 = INT(SQR(c)): Letf r, c2 '      initial root and   Letf q, c - c2 * c2 '              -quadratic residue   lr = 1: c2 SHL= LB + 1 '            double root prefix{2}   FOR k = lp - 2 TO 1 STEP -2 '       expand (r + x)^2 <= p      IF lr = 2 THEN c2 OR= cx SHL 1      Lsft r, -1: lr += 1 '            root * MB      n(ir + lr + 1) = 0      Lsft q, -2: iq = i(q) '          q * MB^2      n(iq + 1) = n(ip + k) '          shift in pwords k, k-1      n(iq) = n(ip + k - 1)      lq = n(iq - 1): Lftj q, lq      FOR j = lq TO lr   n(iq + j + 1) = 0: NEXT '      adjust current residue      iq += lr: c = c2      IF n(iq + 1) THEN    iq += 1: c SHR= LB      END IF      cx = EstQ(iq, c) '               appr.x = q{3} / 2r{2}      DO    n(ir) = cx '                  shift x in root    IF cx = 0 THEN EXIT DO    iq = i(t0): t = ir - iq + 1    c = cx * cx: c0 = cx SHL 1    FOR j = iq TO iq + lr '       t0:= (2r + x) * x       n(j) = c AND M1: c SHR= LB       c += c0 * n(t + j)    NEXT j: n(j) = c    n(iq - 1) = SB    Lftj t0, lr + 1: Add t0, q    IF Gets(t0) = 0 THEN '        verify t0 <= p - r^2       SWAP i(q), i(t0): EXIT DO    END IF    Divd t0, r, t1 '              correction step    cx -= 1 + n(i(t1)) SHR 1 '    ceil[t0 / 2(r + x)]      LOOP   NEXT kEND SUBFUNCTION IsSqr (ByVal p AS SHORT, ByVal r AS SHORT) AS INTEGER STATICDIM AS INTEGER q64(63), q63(62), q55(54), sw, t, kIF NOT sw THEN   FOR t = 0 TO 31      k = t * t      q64(k AND 63) = -1 '             tabulate quadratic residues      q63(k MOD 63) = -1      q55(k MOD 55) = -1   NEXT t: sw = -1END IF   IsSqr = 0   IF q64(n(i(p)) AND 63) THEN      Dup p, t0: k = Divint(t0, 3465)      IF q63(k MOD 63) THEN    IF q55(k MOD 55) THEN '       >98% non-squares filtered       IF Gets(p) THEN EXIT FUNCTION       Isqrt p, t3, r       IsSqr = Isf(t3, 0) '       check possible square    END IF      END IF   END IFEND FUNCTIONFUNCTION Kronec (ByVal p AS SHORT, ByVal q AS SHORT) AS INTEGERDIM AS INTEGER r, x = 1   '   Dup p, t0: Dup q, t1   IF (n(i(t1)) AND 1) = 0 THEN '      even(N)      IF Isf(t1, 0) THEN '             (a/0)         Sets t0, 0         IF NOT Isf(t0, 1) THEN x = 0         Kronec = x: EXIT FUNCTION      END IF      IF (n(i(t0)) AND 1) = 0 THEN '   both even         Kronec = 0: EXIT FUNCTION      END IF      '      IF Odd(t1) AND 1 THEN '          make odd(N)         r = n(i(t0)) AND 7         IF r = 3 OR r = 5 THEN x = -x      END IF   END IF   '   IF Gets(t0) THEN '                  negative(a)      IF Mp2(t1, 4) = 3 THEN x = -x '  N mod 4 = 3      Sets t0, 0 '                     Abs(a)   END IF   Sets t1, 0 '                        Abs(N)   '   DO      IF (n(i(t0)) AND 1) = 0 THEN '   even(a)         IF Isf(t0, 0) THEN            IF NOT Isf(t1, 1) THEN x = 0            EXIT DO         END IF         '         IF Odd(t0) AND 1 THEN '       make odd(a)            r = n(i(t1)) AND 7            IF r = 3 OR r = 5 THEN x = -x         END IF      END IF      r = n(i(t0)) AND n(i(t1))      IF r AND 2 THEN x = -x '         reciprocity      SWAP i(t0), i(t1)      Divd t0, t1, t3   LOOPKronec = xEND FUNCTIONFUNCTION Nxtprm (ByRef sw AS SHORT) AS INTEGER STATICDIM AS INTEGER c(2), cp(2), cv(2), d(2), b(2)DIM AS INTEGER fl, r, clp, t, ini '    3-row stackDIM b4 AS STRING * 4IF NOT ini THEN   clp = LOF(Prmnr): t = 0: ini = -1END IFIF sw AND 1 THEN   IF fl THEN      DO         d(t) XOR= 6: c(t) += d(t) '   skip multiples of 2 and 3         IF b(t) = 0 THEN            b(t) = 30            IF cp(t) < clp THEN               GET #Prmnr, cp(t), b4 ' next bitvector in PrimFlgs.bin               cv(t) = CVL(b4)               cp(t) += 4            ELSE               cv(t) = &H1BF6FDBF '    5-folds excluded            END IF         END IF         r = cv(t) AND 1 '             read bit         b(t) -= 1: cv(t) SHR= 1      LOOP UNTIL r   ELSE      c(t) += d(t)      d(t) SHL= 1: fl = c(t) = 5   END IF   Nxtprm = c(t)ELSE   '   IF sw < 0 THEN '                    pop      IF t THEN t -= 1: fl = c(t) > 4   ELSE      IF sw AND t < 2 THEN t += 1 '    push      c(t) = 2: cp(t) = 1: fl = 0      d(t) = 1: b(t) = 0: sw = -1 '    initialize   Nxtprm = 2   END IFEND IFEND FUNCTIONFUNCTION PrimCeil AS INTEGERDIM b4 AS STRING * 4DIM AS INTEGER c = 0   IF LOF(Prmnr) > 3 THEN      GET #Prmnr, 1, b4      IF CVL(b4) = &HB76BDBF THEN '    valid primelist         c = 5 + (LOF(Prmnr) \ 4) * 90      END IF   END IFPrimCeil = cEND FUNCTIONFUNCTION IsPPrm (ByVal p AS SHORT, ByVal d AS SHORT, ByVal w AS SHORT) AS INTEGERDIM AS INTEGER j, k, a, x = 0, lp = Getl(p)DIM AS SHORT r = 2IsPPrm = 0   IF Gets(p) OR Isf(p, 1) THEN EXIT FUNCTION   '   IF lp = 1 OR w >= 0 THEN '          trial divisions      k = lp * 54      IF k > 669 THEN k = 669 '        pr(669) = 4999      n(i(d) - 1) = 1: n(i(d) + 1) = 0      FOR j = 1 TO k         Dup p, t1: n(i(d)) = Nxtprm(r)         IF Divint(t1, n(i(d))) = 0 THEN            x = Isf(t1, 1): EXIT FOR ' d | N         END IF         x = Cmp(d, t1) = 1         IF x THEN EXIT FOR '          d > sqrt(N)      NEXT j      r = Nxtprm(-2)      IF j <= k THEN         IF x THEN x = 1         IsPPrm = x: EXIT FUNCTION      END IF   END IF   '   a = 2 + INT(RanGen2 * 255) '        1-byte random witness   Dup p, d   Inc d, -1: k = Odd(d) '             determine N - 1 = (2^ k) * odd(m)   '   FOR r = 1 TO ABS(w) '               Miller-Rabin test      Letf t3, a      Modpwr t3, d, p '                if N is prime then      x = Isf(t3, 1) OR Ismin1(t3, p) ' a^ m Mod N = 1 or      IF x = 0 THEN '         (a^ m)^ (2^ j) Mod N = -1         FOR j = 1 TO k - 1 '    for some j in [0, k - 1]            Modsqu t3, p            x = Ismin1(t3, p)            IF x THEN EXIT FOR         NEXT j         IF x = 0 THEN EXIT FOR      END IF      a += 1   NEXT rIsPPrm = xEND FUNCTIONSUB Triald (ByVal q AS SHORT, ByVal p AS SHORT, ByVal c AS INTEGER)DIM AS INTEGER cp, k, r, s = n(i(p) - 1)DIM AS SHORT sw = 2   n(i(q) - 1) = 0 '                   clear list {q}   n(i(p) - 1) = s AND M1 '            Abs(p)   r = n(i(p)) < 2 AND n(i(p) - 1) = 1   IF r THEN EXIT SUB   DO      cp = Nxtprm(sw) '                read primes      IF cp > c THEN EXIT DO      Letf t3, cp      FOR k = 0 TO 32766         Dup p, t0: Divd t0, t3, t1 '  trial divide p         IF NOT Isf(t0, 0) THEN EXIT FOR ' nonzero remainder         SWAP i(p), i(t1) '            running quotient p:= p/dr      NEXT k      IF k THEN         IF k < 32767 THEN            EnQ q, t3, k '             store pfactor         ELSE            Errorh "exp. overflow Sub Triald", 6            EXIT SUB         END IF      END IF      r = Cmp(t3, t1) = 1 '            divisor > Isqrt(p)   LOOP UNTIL r   IF r AND NOT Isf(p, 1) THEN      EnQ q, p, 1: Letf p, 1   END IF   cp = Nxtprm(-2)   Sets p, sEND SUB' ***************************************************************************************'                                 Conversion functions' ***************************************************************************************SUB Readst (ByVal p AS SHORT, ByRef gp AS ZSTRING PTR)DIM AS STRING g = RTRIM\$(*gp), h = "&H"DIM AS INTEGER sg = INSTR(g, "-")DIM AS INTEGER d, j, k, m, lp   j = INSTR(g, "0x") OR INSTR(g, h)   IF sg OR j THEN      m = j + 1: IF m < sg THEN m = sg      g = MID\$(g, m + 1)   ELSE      g = LTRIM\$(g)   END IF   Letf p, 0: k = LEN(g)   IF g = "0" OR k = 0 THEN EXIT SUB   '   IF j = 0 THEN      Letf t1, MD: d = LD: h = "" '    decimal   ELSE      Letf t1, &H1000: d = 3 '         or hex base   END IF   lp = k \ d: m = k - lp * d   IF m = 0 THEN m = d: lp -= 1   '   k = 1: FOR j = 0 TO lp      IF j THEN Lmul p, t1 '           p:= p * base      Inc p, VAL(h + MID\$(g, k, m)) '  add base-digit      k += m: m = d   NEXT j   k = n(i(p) - 1)   IF sg THEN n(i(p) - 1) = SB   Lftj p, k - 1END SUB#ifndef USE_PREFIXFUNCTION Logf (ByVal p AS SHORT) AS DOUBLE#elseFUNCTION fbLogf (ByVal p AS SHORT) AS DOUBLE#endifSTATIC lMB AS DOUBLE, sw AS INTEGERDIM AS INTEGER j, t = Getl(p) - 1, k = t - 3DIM AS DOUBLE p3 = n(i(p) + t)IF NOT sw THEN   lMB = LOG(MB): sw = -1END IF   IF k < 0 THEN k = 0   FOR j = t - 1 TO k STEP -1      p3 = p3 * MB + n(i(p) + j)   NEXT j   IF p3 THEN      RETURN LOG(p3) + k * lMB   ELSE      Errorh "zero argument Logf", 11      RETURN SB   END IFEND FUNCTIONFUNCTION Bufl (ByVal p AS SHORT) AS INTEGERBufl = 2   IF Isf(p, 0) THEN EXIT FUNCTION   #ifndef USE_PREFIX   Bufl = 2 + INT(Logf(p) * .4343)   #else   Bufl = 2 + INT(fbLogf(p) * .4343)   #endifEND FUNCTIONSUB CnvSt (ByRef gp AS ZSTRING PTR, ByVal p AS SHORT)DIM AS STRING f, g = *gpDIM AS INTEGER k = LEN(g), sg, t   IF k < 2 THEN      Errorh "no buffer Sub CnvSt", 14      EXIT SUB   END IF   IF Isf(p, 0) THEN *gp = " 0": EXIT SUB   '   Dup p, t2: sg = Absf(t2)   FOR t = k + 1 TO 2 STEP -LD      f = LTRIM\$(STR\$(Divint(t2, MD)))      MID\$(g, t - LEN(f)) = f '        stuff   NEXT t   t = 0: WHILE g[t] = 48 '            "0"      g[t] = 32: t += 1 '              " "   WEND   IF sg THEN g[t - 1] = 45 '          "-"   *gp = gEND SUBSUB RatCnv (ByRef gp AS ZSTRING PTR, ByVal p AS SHORT, ByVal q AS SHORT, _ ByVal b AS SHORT)DIM AS STRING g = *gp, hDIM AS INTEGER k = LEN(g), r, s, sg, e, t, m, j, d   *gp = STRING\$(k, CHR\$(0))   IF Isf(q, 0) THEN      g = "oo": k = 2      IF Gets(p) THEN g = "-oo"   END IF   IF Isf(p, 0) THEN g = "0": k = 1   IF k < 10 THEN *gp = g: EXIT SUB   '   IF b < 2 OR b > 17 THEN b = 10   r = Absf(p): s = Absf(q): sg = r XOR s   #ifndef USE_PREFIX   e = 1 + INT((Logf(p) - Logf(q)) / LOG(b))   #else   e = 1 + INT((fbLogf(p) - fbLogf(q)) / LOG(b))   #endif   '   t = k - e - 1 '                     scaling factor   Letf t3, b   IF ABS(t) > 1 THEN Powr t3, t   IF t > 0 THEN      Mult p, t3, t0   ELSE      Dup p, t0      IF t < 0 THEN Ldiv t0, t3   END IF   Divd t0, q, t1 '                    integral quotient   Sets p, r: Sets q, s   '   IF b = 10 THEN      CnvSt g, t1   ELSE      s = -1: t = 1      WHILE t < MB         s += 1: m = t: t *= b '       max. s with b^s < MB      WEND      h = "0123456789ABCDEFG"      FOR j = k TO 1 STEP -s '         base conversion         t = j: r = Divint(t1, m)         WHILE r            d = r \ b: t -= 1            g[t] = h[r - d * b]: r = d         WEND      NEXT j      t = 0: WHILE g[t] = 48 '         "0"         g[t] = 32: t += 1 '           " "      WEND   END IF   g = LTRIM\$(g)   s = LEN(g): e -= k - s '            exponent   '   IF Isf(t0, 0) THEN '                q divides scaled p      r = 2: IF e > 0 THEN r = 2 + e      FOR t = s TO r STEP -1         IF g[t - 1] - 48 THEN EXIT FOR      NEXT t      g = LEFT\$(g, t) '                strip trailing zeroes   END IF   '   IF sg THEN k -= 1   t = 1: s = b <> 10   DO      IF e OR s THEN '                 make suffix h    IF e THEN h = LTRIM\$(STR\$(ABS(e)))    SELECT CASE e    CASE IS < 0: h = "E-" + h    CASE IS > 0: h = "E+" + h    CASE ELSE: h = "E0"    END SELECT    IF s THEN h = " " + LTRIM\$(STR\$(b)) + h    d = k - LEN(h)    '    SELECT CASE e    CASE -4 TO -1 '               fixed point < 1       IF s AND (e < 7 - k) THEN EXIT DO       g = STRING\$(-e, "0") + g: e = 0    CASE 1 TO k - 2 '             fp > 1       IF s AND (e > k - 6) THEN EXIT DO       t = 1 + e: e = 0    CASE ELSE: EXIT DO    END SELECT      ELSE    h = "": d = k: EXIT DO      END IF   LOOP   '   IF LEN(g) > t THEN      g = LEFT\$(LEFT\$(g, t) + "." + MID\$(g, t + 1), d)   END IF   IF sg THEN g = "-" + g   *gp = g + hEND SUB' ***************************************************************************************'                                   Obtaining output' ***************************************************************************************SUB Printf (ByRef f AS STRING, ByRef g AS STRING, ByRef h AS STRING, ByVal sw AS SHORT)DIM s AS STRING, k AS INTEGER   SELECT CASE sw   CASE 0      PRINT f; g; h;      IF Lognr THEN PRINT #Lognr, f; g; h;   CASE 1      PRINT f; g; h      IF Lognr THEN PRINT #Lognr, f; g; h   CASE ELSE      k = LEN(g)      IF LEFT\$(g, 1) = "-" THEN k -= 1      s = "  [" + LTRIM\$(STR\$(k)) + "]"      PRINT f; g; h; s      IF Lognr THEN PRINT #Lognr, f; g; h; s   END SELECTEND SUBSUB Printn (ByVal p AS SHORT, ByRef f AS STRING, ByRef h AS STRING, _ ByVal sw AS SHORT)DIM AS STRING g = STRING\$(Bufl(p), "0")   CnvSt g, p   Printf f, LTRIM\$(g), h, swEND SUBSUB Printr (ByVal p AS SHORT, ByVal q AS SHORT, ByRef f AS STRING, ByVal k AS SHORT, _ ByVal sw AS SHORT)DIM AS INTEGER t = k + 1DIM g AS STRING   IF t < 10 THEN t = 10   g = STRING\$(t, "0")   RatCnv g, p, q, 10   Printf f, RTRIM\$(g), "", sw AND 1END SUBSUB Prints (ByRef f AS STRING, ByVal sw AS SHORT)   SELECT CASE sw   CASE 0      PRINT f;      IF Lognr THEN PRINT #Lognr, f;   CASE 1      PRINT f      IF Lognr THEN PRINT #Lognr, f   CASE ELSE      PRINT f: PRINT      IF Lognr THEN PRINT #Lognr, f: PRINT #Lognr,   END SELECTEND SUB'================================================='end of largeinteger'================================================='Example:CONST p0 = 0, p1 = 1        ' largeint pointersDIM AS INTEGER m = LargeInit(p1, "fibonacc")DIM AS DOUBLE  tim,tim2DIM AS INTEGER nx = 15000DIM AS STRING  gCLSLOCATE 1, 2           tim=timerLetf p0, 0: Letf p1, 1FOR m = 1 TO nx   Add p0, p1: Swp p0, p1            ' F_n = F_n-2 + F_n-1NEXT mtim2=timerprint "Fib ";15000Printn p0, g, "", 2print tim2-timsleep `

here is a second method for fibonacci as a check:

Code: Select all

`Dim Shared As String * 20 _addQmod="01234567890123456789"Dim Shared As String * 20 _addbool_addbool=Chr(0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1)Function plusINT(num1 As String,num2 As String) As String    Var n_=0,flag=0    Dim As Ubyte addup=Any,addcarry=Any    #macro finish()    answer=Ltrim(answer,"0")    If flag Then Swap num1,num2    Return answer    #endmacro    If Len(num2)>Len(num1) Then         Swap num2,num1        flag=1    End If    Var diff=Len(num1)-Len(num2)    Var k=Sgn(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()    End If    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 FunctionFunction  fibonacci(n As Integer) As String    Dim As String sl,l,term    sl="1": l="1"    If n=1 Then Return "1"    If n=2 Then Return "1"    n=n-2    For x As Integer= 1 To n        term=plusINT(l,sl)        sl=l        l=term    Next x    Function =termEnd Functiondim as double t1=timer,t2var ans= fibonacci(15000)t2=timerprint "fib ";15000print ansprint t2-t1sleep `
jdebord
Posts: 529
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

### Re: FBMath Update

I have checked that the compiled module largeint.o is present in the library file libfbmath.a and that the demo programs are compiled correctly. So, the problem is with the linker. Unfortunately, I have not been able to fix it at this time. I will check this with the developer of the largeint library.

In the meantime, dodicat's solution seems to be the best choice.

I apologize for the inconvenience.
Roland Chastain
Posts: 948
Joined: Nov 24, 2011 19:49
Location: France
Contact:

### Re: FBMath Update

jdebord wrote:I apologize for the inconvenience.

I like FBMath. As I am not a mathematician, it's not always easy for me to understand the examples, but I like the coding style and the cleanliness of the package and of the documentation.
frisian
Posts: 249
Joined: Oct 08, 2009 17:25

### Re: FBMath Update

Roland Chastain wrote:Hello!

I have compilation errors with the largint demos:

C:\FreeBASIC\fbmath\demo\largeint\bench.o:fake:(.text+0x56): undefined reference to `LargeInit@8'
C:\FreeBASIC\fbmath\demo\largeint\bench.o:fake:(.text+0xef): undefined reference to `Prints@8'
C:\FreeBASIC\fbmath\demo\largeint\bench.o:fake:(.text+0x101): undefined reference to `Letf@8'
...

I had that same problem, I tried some things and found that if you remove the alias "..." in the large integer section of FBMath.bi all works as it should be.

Example:

DECLARE FUNCTION LargeInit alias "LargeInit" (ByVal k AS SHORT, ByRef f AS STRING) AS INTEGER

into

DECLARE FUNCTION LargeInit (ByVal k AS SHORT, ByRef f AS STRING) AS INTEGER

Best to rebuild the libFBMath.a . with altered FBMath.bi files.

I include a link to my modified FBMath.bi.

Updated FBMath.bi one alias had escaped the removal.

EDIT: removed link, no need for that, FBmath has been updated.
Last edited by frisian on Oct 07, 2014 14:22, edited 1 time in total.

### Who is online

Users browsing this forum: No registered users and 15 guests