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
FBMath Update
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.
I downloaded this build from MOD'S post in this page (fifth post down):
http://www.freebasic.net/forum/viewtopi ... &start=165
I've (as an exercise) compared Your GaussJordan and a UDT GaussJordan of my own.
Using random Matrix elements:
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/viewtopi ... &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
Re: FBMath Update
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) ?
Is it possible to overload the () operator, so as to write A(i,j) instead of A.element(i,j) ?
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).
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).
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)
As a silly workaround, just define some characters to substitute .element, or to save typing just make .element .e
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
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
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!
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 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)]
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 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)
Re: FBMath Update
Thanks again. A notation like A[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^^[j] which is not very handy!
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
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: 1006
- Joined: Nov 24, 2011 19:49
- Location: France
- Contact:
Re: FBMath Update
Hello!
I have compilation errors with the largint demos:
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'
...
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:
here is a second method for fibonacci as a check:
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
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
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.
In the meantime, dodicat's solution seems to be the best choice.
I apologize for the inconvenience.
-
Roland Chastain
- Posts: 1006
- Joined: Nov 24, 2011 19:49
- Location: France
- Contact:
Re: FBMath Update
Thank you for your answer.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.
Re: FBMath Update
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.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'
...
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.
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.