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

Postby jdebord » Jul 15, 2014 7:41

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

Postby dodicat » Jul 15, 2014 15:50

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.

I downloaded this build from MOD'S post in this page (fifth post down):

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 = MatOk
END 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 matrix
End Type

'matrix multiply(NOT USED HERE)
Operator *(m1 As matrix,m2 As matrix) As matrix
Dim rows As Integer=Ubound(m1.element,1)
Dim columns As Integer=Ubound(m2.element,2)
If Ubound(m1.element,2)<> Ubound(m2.element,1)Then
Print "Can't do"
Exit Operator
End If
Dim As matrix ans
Redim ans.element(1 To rows,1 To columns)
Dim rxc As Double
For 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 c
Next r
Operator= ans
End Operator

'rounding function
Function 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 r
End Function

'matrix printer
Operator matrix.cast() As String
Dim As String ans,comma
For 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 a
Operator= ans
End Operator

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

function 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 ans
End function
'=============  TEST ======================================================
screen 20


do
    cls
dim as integer d=1+int(rnd*15)'set a dimension

dim as double a(1 to d,1 to d+1),det
Dim As matrix m,rhs

Redim m.element(1 To d,1 To d)
Redim rhs.element(1 To d,1 To 1)


'FILL UP EVERYTHING
for 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 c
next r


GaussJordan(a(),det)'  FBMATH

print "SYSTEM DIMENSION ";d
print"FBMATH determinant "; det
print"   UDT determinant "; m.determinant
print
print "FBMATH solution "
for n as integer=1 to d
    print str(a(n,d+1))
next n
print
print "UDT solution "
print m.GaussJordan(rhs)
print
print "Press a key or esc"
sleep
loop until inkey=chr(27)
sleep


 
jdebord
Posts: 529
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

Re: FBMath Update

Postby jdebord » Jul 17, 2014 7:56

Thank you dodicat. I have downloaded the new compiler.

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

Postby fxm » Jul 17, 2014 13:57

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

Postby dodicat » Jul 17, 2014 14:12

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

Operator Vector.[](I as integer) ByRef As double
    operator= element(I)
End Operator

Operator vector.cast() As String
Dim As String ans
For a As Integer=lbound(element) To Ubound(element)
        ans+=str(this[a])+chr(10)
Next a
Operator= ans
End Operator

dim as vector Vct
redim Vct.element(1 to 10)

for n as integer=1 to 10
    Vct[n]=n
next n

print Vct

'=========================================
print
redim preserve Vct.element(1 to 15)
for n as integer=11 to 15
    Vct[n]=n*n+.1
next n
print Vct

sleep


 


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 __ .element
end type

dim as matrix m
redim 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

Postby jdebord » Jul 18, 2014 7:33

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

Postby fxm » Jul 18, 2014 9:01

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 Type
Constructor Pt (Byval x0 As Integer = 0, Byval y0 As Integer = 0)
  This.x = x0 : This.y = y0
End Constructor

Type Vector
  Dim As Double element(Any, Any)
  Declare Operator [](Byref N As Pt) Byref As Double
End Type
Operator Vector.[](Byref N As Pt) Byref As Double
  Operator = This.element(N.x, N.y)
End Operator


Dim as Vector v
Redim v.element(0 To 9, 0 To 9)
v[Pt(1, 2)] = 12
Print v[Pt(1, 2)]
fxm
Posts: 9994
Joined: Apr 22, 2009 12:46
Location: Paris suburbs, FRANCE

Re: FBMath Update

Postby fxm » Jul 18, 2014 9:45

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 Double
End Type
Operator Vector1.[](Byval index As Integer) Byref As Double
  Operator = This.element1(index)
End Operator

Type Vector2
  Dim As Vector1 element2(Any)
  Declare Operator [](Byval index As Integer) Byref As Vector1
End Type
Operator Vector2.[](Byval index As Integer) Byref As Vector1
  Operator = This.element2(index)
End Operator


Dim As Vector2 v
Redim v.element2(0 To 9)
For I As Integer = 0 To 9
  Redim (v.element2(I).element1)(0 To 9)
Next I
v[3][4] = 34   '' equivalent to: v.element2(3).element1(4) = 34
Print 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

Postby jdebord » Jul 19, 2014 8:10

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

Postby dodicat » Jul 19, 2014 23:46

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 String
End Type
Operator vector.[](Byval index As Integer) Byref As Double
  Operator = Columns(index)
End Operator

Type 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 Double
End Type
Operator matrix.[](Byval index As Integer) Byref As vector
  Operator = Rows(index)
End Operator


operator vector.cast() as string
dim as string ans,comma
For a As Integer=lbound(Columns) To Ubound(Columns)
    If a=Ubound(Columns) Then comma="" Else comma=" , "
    ans=ans+str(this[a])+comma
next a
operator=ans
end operator

'matrix printer
Operator matrix.cast() As String
Dim As String ans,comma
For a As Integer=lbound(Rows) To Ubound(Rows)
    print Rows(a)
Next a
Operator= ans
End Operator

Function 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-1
Redim M.Rows(r)
For I as integer  = r
  Redim (M.Rows(I).Columns)(c)
Next I
end scope
#endmacro

'=================================================
Dim As matrix MAT

RedimMatrix(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 b
next a

print "MATRIX"
print MAT
print "Determinant"
print MAT.determinant
sleep
Roland Chastain
Posts: 948
Joined: Nov 24, 2011 19:49
Location: France
Contact:

Re: FBMath Update

Postby Roland Chastain » Aug 13, 2014 3:49

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

Postby dodicat » Aug 15, 2014 22:36

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

CONST Asize = 97279 '                  number array size
CONST ErrTrap = 0 '                    optionally stop on error

CONST LB = 15, MB = &H8000& '          binary storage base
CONST LD = 4,  MD = 10000 '            decimal string base < MB
CONST t0 = -1, t1 = -2 '               pointers to swap-registers
CONST t2 = -3, t3 = -4

CONST MH = MB \ 2, M1 = &H7FFF% '      bitmask MB - 1
CONST M2 = M1 * MB, SB = NOT M1 '      borrow, sign bit
DIM SHARED AS SHORT n(Asize), uj
'n() is the largeint number array,
'signs and lengths are at positions -1,
'followed by (length + 1) largeint-words
DIM SHARED i() AS INTEGER
'the i()'s are number indices
DIM SHARED AS INTEGER Bezsw, Errsw '   switches
DIM SHARED AS INTEGER Lognr, Prmnr '   filehandles

SUB 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 IF
END 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, length
END SUB

SUB 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) = 0
END SUB

FUNCTION LibErr AS INTEGER
   LibErr = Errsw: Errsw = 0
END FUNCTION

SUB Work (ByRef f AS STRING)
   f = TRIM$(ENVIRON$("LARGEINT"))
   IF f = "" THEN f = CURDIR$
   IF RIGHT$(f, 1) <> "\" THEN f += "\"
END SUB

FUNCTION LargeInit (ByVal ix AS SHORT, ByRef f AS STRING) AS INTEGER
DIM s AS STRING, t AS INTEGER
LargeInit = 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
   '
   InitMTbyTime
LargeInit = uj: uj = uj - 2
END FUNCTION

SUB Term
   CLOSE 'all files, then terminate
END SUB

' ***************************************************************************************
'                                    Direct access
' ***************************************************************************************
FUNCTION Getl (ByVal p AS SHORT) AS INTEGER
   Getl = n(i(p) - 1) AND M1
END FUNCTION

FUNCTION Gets (ByVal p AS SHORT) AS INTEGER
   Gets = n(i(p) - 1) AND SB
END FUNCTION

FUNCTION Getw (ByVal p AS SHORT, ByVal j AS SHORT) AS INTEGER
   Getw = n(i(p) + j)
END FUNCTION

SUB Setl (ByVal p AS SHORT, ByVal a AS SHORT)
   n(i(p) - 1) = (n(i(p) - 1) AND SB) OR ABS(a)
END SUB

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

SUB Setw (ByVal p AS SHORT, ByVal j AS SHORT, ByVal a AS SHORT)
   n(i(p) + j) = ABS(a)
END SUB

SUB 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 IF
END SUB

FUNCTION ExQ (ByVal p AS SHORT, ByRef k AS SHORT, ByVal q AS SHORT, _
 ByRef j AS SHORT) AS INTEGER
DIM 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 IF
END 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 j
END SUB

SUB Swp (ByVal p AS SHORT, ByVal q AS SHORT)
   SWAP i(p), i(q)
END SUB

FUNCTION Cmp (ByVal p AS SHORT, ByVal q AS SHORT) AS INTEGER
DIM 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 IF
Cmp = x
END FUNCTION

FUNCTION Isf (ByVal p AS SHORT, ByVal a AS SHORT) AS INTEGER
DIM 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 1
Isf = r = s AND n(i(p)) = ABS(a)
END FUNCTION

SUB 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 t
END 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) - 1
DIM AS INTEGER im = i(q), lm = (n(im - 1) AND M1) - 1
DIM 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 + 1
END SUB

SUB Add (ByVal p AS SHORT, ByVal q AS SHORT)
   n(i(q) - 1) XOR= SB
   Subt p, q
   n(i(q) - 1) XOR= SB
END SUB

SUB Inc (ByVal p AS SHORT, ByVal a AS SHORT)
DIM AS INTEGER j, b = SGN(Gets(p) OR 1) * a
DIM 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 IF
END SUB

FUNCTION Ismin1 (ByVal p AS SHORT, ByVal m AS SHORT) AS INTEGER
   Inc p, 1
   Ismin1 = Cmp(p, m) = 0
   Inc p, -1
END FUNCTION

FUNCTION Divint (ByVal p AS SHORT, ByVal a AS SHORT) AS INTEGER
DIM 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, lp
Divint = c
END FUNCTION

FUNCTION EstQ (ByVal t AS INTEGER, ByVal c AS INTEGER) AS INTEGER
DIM 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 = M1
EstQ = c
END FUNCTION

SUB Divd (ByVal p AS SHORT, ByVal d AS SHORT, ByVal q AS SHORT)
DIM AS INTEGER iq = i(q) - 1, ip = i(p) - 1
DIM AS INTEGER it = i(d) - 1, lt = n(it) AND M1
DIM 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 IF
END SUB

SUB Ldiv (ByVal p AS SHORT, ByVal d AS SHORT)
   Divd p, d, t1: SWAP i(p), i(t1)
END SUB

SUB 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) - 1
DIM AS INTEGER iq = i(q), lq = (n(iq - 1) AND M1) - 1
DIM 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 + 1
END SUB

SUB Lmul (ByVal p AS SHORT, ByVal q AS SHORT)
   Mult p, q, t2: SWAP i(p), i(t2)
END SUB

SUB 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 + 1
END SUB

SUB Chs (ByVal p AS SHORT)
   n(i(p) - 1) XOR= SB '               signed magnitude
END SUB

FUNCTION Absf (ByVal p AS SHORT) AS INTEGER
DIM AS INTEGER s = n(i(p) - 1)
   n(i(p) - 1) = s AND M1
Absf = s AND SB
END 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 3
DIM 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, lx
END SUB

SUB Lsft (ByVal p AS SHORT, ByVal r AS SHORT)
DIM AS INTEGER ip = i(p), lp = n(ip - 1) AND M1
DIM 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, lp
END SUB

SUB Rsft (ByVal p AS SHORT, ByVal r AS SHORT)
DIM AS INTEGER ip = i(p), lp = (n(ip - 1) AND M1) - 1
DIM 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 + 1
END SUB

FUNCTION Odd (ByVal p AS SHORT) AS INTEGER
DIM AS INTEGER a, lp, ip = i(p), t = 0
Odd = 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 IF
Odd = t
END FUNCTION

FUNCTION Bitl (ByVal p AS SHORT) AS INTEGER
DIM AS INTEGER a, t = Getl(p) - 1
   a = n(i(p) + t): t *= LB
   WHILE a
      a SHR= 1: t += 1
   WEND
Bitl = t
END 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 IF
END SUB

SUB 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 IF
END SUB

FUNCTION Mp2 (ByVal p AS SHORT, ByVal a AS SHORT) AS INTEGER
DIM AS INTEGER m = n(i(p)) AND a - 1
   IF Gets(p) THEN m = a - m AND a - 1
Mp2 = m
END FUNCTION

SUB 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, t1
END SUB

SUB Modsqu (ByVal p AS SHORT, ByVal m AS SHORT)
   Squ p, t2: SWAP i(p), i(t2)
   Divd p, m, t1
END SUB

SUB 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, s
END SUB

FUNCTION Hp2 (ByVal p AS SHORT) AS INTEGER
DIM AS INTEGER a = n(i(p) + Getl(p) - 1), a0 = MH
   DO UNTIL a AND a0
      a0 SHR= 1: LOOP '                bitscan L=>R
Hp2 = a0
END FUNCTION

SUB 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 IF
END SUB

SUB Powr (ByVal p AS SHORT, ByVal c AS INTEGER)
   Letf t2, ABS(c): Letf t0, 0
   Modpwr p, t2, t0
END 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 > a
END SUB

SUB 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, 0
END SUB

SUB 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 IF
END SUB

SUB 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 =                / q
END SUB

SUB 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 k
END SUB

FUNCTION IsSqr (ByVal p AS SHORT, ByVal r AS SHORT) AS INTEGER STATIC
DIM AS INTEGER q64(63), q63(62), q55(54), sw, t, k
IF 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 = -1
END 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 IF
END FUNCTION

FUNCTION Kronec (ByVal p AS SHORT, ByVal q AS SHORT) AS INTEGER
DIM 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
   LOOP
Kronec = x
END FUNCTION

FUNCTION Nxtprm (ByRef sw AS SHORT) AS INTEGER STATIC
DIM AS INTEGER c(2), cp(2), cv(2), d(2), b(2)
DIM AS INTEGER fl, r, clp, t, ini '    3-row stack
DIM b4 AS STRING * 4
IF NOT ini THEN
   clp = LOF(Prmnr): t = 0: ini = -1
END IF
IF 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 IF
END IF
END FUNCTION

FUNCTION PrimCeil AS INTEGER
DIM b4 AS STRING * 4
DIM 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 IF
PrimCeil = c
END FUNCTION

FUNCTION IsPPrm (ByVal p AS SHORT, ByVal d AS SHORT, ByVal w AS SHORT) AS INTEGER
DIM AS INTEGER j, k, a, x = 0, lp = Getl(p)
DIM AS SHORT r = 2
IsPPrm = 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 r
IsPPrm = x
END FUNCTION

SUB 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, s
END 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 - 1
END SUB

#ifndef USE_PREFIX
FUNCTION Logf (ByVal p AS SHORT) AS DOUBLE
#else
FUNCTION fbLogf (ByVal p AS SHORT) AS DOUBLE
#endif
STATIC lMB AS DOUBLE, sw AS INTEGER
DIM AS INTEGER j, t = Getl(p) - 1, k = t - 3
DIM AS DOUBLE p3 = n(i(p) + t)
IF NOT sw THEN
   lMB = LOG(MB): sw = -1
END 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 IF
END FUNCTION

FUNCTION Bufl (ByVal p AS SHORT) AS INTEGER
Bufl = 2
   IF Isf(p, 0) THEN EXIT FUNCTION
   #ifndef USE_PREFIX
   Bufl = 2 + INT(Logf(p) * .4343)
   #else
   Bufl = 2 + INT(fbLogf(p) * .4343)
   #endif
END FUNCTION

SUB CnvSt (ByRef gp AS ZSTRING PTR, ByVal p AS SHORT)
DIM AS STRING f, g = *gp
DIM 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 = g
END SUB

SUB RatCnv (ByRef gp AS ZSTRING PTR, ByVal p AS SHORT, ByVal q AS SHORT, _
 ByVal b AS SHORT)
DIM AS STRING g = *gp, h
DIM 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 + h
END 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 SELECT
END SUB

SUB 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, sw
END SUB

SUB 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 + 1
DIM g AS STRING
   IF t < 10 THEN t = 10
   g = STRING$(t, "0")
   RatCnv g, p, q, 10
   Printf f, RTRIM$(g), "", sw AND 1
END SUB

SUB 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 SELECT
END SUB
'=================================================
'end of largeinteger
'=================================================
'Example:
CONST p0 = 0, p1 = 1        ' largeint pointers

DIM AS INTEGER m = LargeInit(p1, "fibonacc")
DIM AS DOUBLE  tim,tim2
DIM AS INTEGER nx = 15000
DIM AS STRING  g

CLS
LOCATE 1, 2           
tim=timer
Letf p0, 0: Letf p1, 1
FOR m = 1 TO nx
   Add p0, p1: Swp p0, p1            ' F_n = F_n-2 + F_n-1
NEXT m
tim2=timer
print "Fib ";15000
Printn p0, g, "", 2
print tim2-tim
sleep

 


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 Function

Function  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 =term
End Function

dim as double t1=timer,t2
var ans= fibonacci(15000)
t2=timer
print "fib ";15000
print ans
print t2-t1
sleep
 
jdebord
Posts: 529
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

Re: FBMath Update

Postby jdebord » Aug 19, 2014 7:53

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

Postby Roland Chastain » Aug 19, 2014 9:37

jdebord wrote:I apologize for the inconvenience.


Thank you for your answer.

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

Postby frisian » Aug 19, 2014 9:53

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.

Return to “Projects”

Who is online

Users browsing this forum: No registered users and 15 guests