' Walter Trump 2014-04-08
' This is a tutorial program that explains the main features of bit vector programming.
' The number of semi-magic 5x5-squares is calculated in less than two minutes.
' With slightly optimized procedures it would last less than half a minute.

sqWin
sqMain
sqEnd

Proc sqWin
' Open a window - not important
Local Int x = Rand(200), y = Rand(200), w = 320, h = 360
OpenW 1, x, y, w, h
: Me.Caption = "Tutorial - 2014-04-08 - Version 1.07"
: Me.AutoRedraw = 1 : Me.SetFont "Courier New", 11, 1
: Me.ControlBox = True : Me.Visible = True : Me.SetFocus
TopW 1
DoEvents  // Give Windows time to handle certain interactions.

Proc sqConst
' Some values that are constant. The bit vectors are declared as Int32.
' Important: In a program for 8x8-squares you need Uint64 instead of Int32 (replace all Int32)
Public Const n        As Int   = 5                    // order of magic square
Public Const mx       As Int   = n * n                // maximal number (bit)
Public Const mc       As Int   = n * (mx + 1) / 2     // magic constant
Public Const nSeries  As Int   = 1394                 // number of magic series (was calculated in advance)
Public Const sLimit   As Int   = nSeries / (n - 1)    // maximal number of series consisting of a certain entry
Public Const AllBits  As Int32 = 2 ^ mx - 1           // in this bit vector all considered bits (1 to mx) are set

Proc sqDimVar
' Bit vectors. For higher orders use Uint64.
Public Int32 r1, r2, r3, r4, r5      // binary coded rows
Public Int32 c1, c2, c3, c4, c5      // binary coded cols
Public Int32 br0, br1, br2, br3, br4 // brx = all bits of the rows 1 to x (including c1)
Public Int32 bc1, bc2, bc3, bc4      // bcx = all bits of the cols 1 to x (including r1)
Public Int32 OBI(32)                  // integers with only one set bit (One-Bit-Integer)
' General integers
Public Int cnt = 0            // Counter for semi-magic squares
Public Int Trigger = 0        // enables square print in certain intervalls (not important)
Public Int Intervall = 500000 // Trigger-Intervall (not important)
Public Int ir1                // Index of row 1
Public Int ic1                // Index of col 1
Public Int j(n, n)            // Entries of the actual semi-magic square
Public Double StartTime = Timer  // Timer (not important)
' Magic series (bit vectors)
Public Int32   Series(nSeries)           // all series
Public Int32  KSeries(sLimit, mx)        // series with a certain number k
Public Int NumKSeries(mx)                // Number of series with a certain number k
' The following tables contain the remaining possible series
' after reducing the number of series by certain conditions
' Example: tr4s2  means t=table, r4 = series for row 4, s2 = reduction step 2
'         itr4s2 = index of the table or (after reduction is done) number of series in the table
' Reduction step 1 (...s1)
Public Int   itr2s1,        itr3s1,        itr4s1
Public Int32  tr2s1(sLimit), tr3s1(sLimit), tr4s1(sLimit)
Public Int   itc2s1,        itc3s1,        itc4s1
Public Int32  tc2s1(sLimit), tc3s1(sLimit), tc4s1(sLimit)
' Reduction step 2 (...s1)
Public Int   itc2s2,        itc3s2,        itc4s2
Public Int32  tc2s2(sLimit), tc3s2(sLimit), tc4s2(sLimit)
' Reduction step 3 (...s1)
Public Int   itc3s3,        itc4s3
Public Int32  tc3s3(sLimit), tc4s3(sLimit)

Proc sqMain
' Some outputs for the window
Cls
Print " Semi-magic 5x5-Squares ";
' Calling initial procedures
sqConst
sqDimVar
OneBitIntegers
MagicSeries
SplitSeries
Print AT(4, 3); "    row1 =";
Print AT(4, 4); "    col1 =";
Print AT(4, 6); " squares =";
' ------- Start of backtracking
Row1
' -------  End  of backtracking
Print AT(4, 8); " Time:"; Str(Timer - StartTime, 8, 2)

Proc OneBitIntegers
' Fills the array OBI() which consists of all integers with only one set bit.
Local Int    i
Local Int32  b = 1
For i = 1 To 32
OBI(i) = b
b *= 2
End For

Proc MagicSeries
' Fills the array Series()
Local Int i1, i2, i3, i4, i5      // n numbers in the series
Local Int c = 0                   // counter
For i1 = 1 To mx
For i2 = i1 + 1 To mx
For i3 = i2 + 1 To mx
For i4 = i3 + 1 To mx
i5 = mc - i1 - i2 - i3 - i4
If i4 < i5
If i5 <= mx
c++
' Store a bit vector where exactly n bits were set:
Series(c) = OBI(i1) Or OBI(i2) Or OBI(i3) Or OBI(i4) Or OBI(i5)
End If
End If
End For
End For
End For
End For
If c <> nSeries Then Stop    // check if the number of series is correct

Proc SplitSeries
' Fills the arrays KSeries(n,k) which consist of the number k.
' But this certain number k is cleared in the series!
Local Int   i, k, n
Local Int32 b
' Series with number mx = 25
n = 0
For i = 1 To nSeries
b = Series(i)
If (b And OBI(mx))
n++
KSeries(n, mx) = b Xor OBI(mx)   // bit mx is cleared
End If
NumKSeries(mx) = n
End For
' Series without number mx = 25
For k = 1 To mx - 1
n = 0
For i = 1 To nSeries
b = Series(i)
If (b And OBI(mx)) = 0
If (b And OBI(k))
n++
KSeries(n, k) = b Xor OBI(k)   // bit k is cleared
End If
End If
End For
NumKSeries(k) = n
End For

Proc Row1
Local Int x, p                        // p = bit position from 1 to mx
j(1, 1) = mx                          // Largest number in position x=1 and y=1
' Try as r1 all series with number mx = 25:
For ir1 = NumKSeries(mx) DownTo 1
Print AT(14, 3); Str(ir1, 4);       // Show index of row 1
DoEvents
r1 = KSeries(ir1, mx)
br0 = r1 Or OBI(mx)      // bit mx  =25 is set
' Fill j(x,1) for all x, that are the numbers of row 1
p = mx
For x = 2 To n
Repeat : p-- : Until r1 And OBI(p)
j(x, 1) = p
End For
' call procedure for column 1
Col1
End For

Proc Col1
Local Int y, p
' Try all possible bit vectors (series) c1 for the first column.
' To avoid reflected (x-y-interchanged) solutions the index ic1 has to be smaller than ir1.
For ic1 = ir1 - 1 DownTo 1
' mx = 25 has to be an entry of c1.
c1 = KSeries(ic1, mx)
' c1 must not have a common number with r1 (row 1).
If (r1 And c1) = 0
' Get the numbers j(1,y) of the first column from c1
p = mx
For y = 2 To n
Repeat : p-- : Until c1 And OBI(p)
j(1, y) = p
End For
br1 = br0 Or c1
RedRowsStep1
Row2
Print AT(14, 4); Str(ic1, 4);        // Show index of col 1
Print AT(14, 6); cnt;                // number of squares found so far
DoEvents
End If
End For

Proc RedRowsStep1
' Reduction step 1 for rows
' Consider for all rows y all series which contain the number nc1(y) of the first col.
' Check if a series has no common number with r1 and c1 (br1 = r1 or c1).
' Store all correct series (bit vectors) in new arrays: tr2s1(), tr3s1(), ...
Local Int   i
Local Int32 b
' Row 2
itr2s1 = 0
For i = NumKSeries(j(1, 2)) DownTo 1
b = KSeries(i, j(1, 2))
If (b And br1) = 0
itr2s1++
tr2s1(itr2s1) = b
End If
End For
' Row 3
itr3s1 = 0
For i = NumKSeries(j(1, 3)) DownTo 1
b = KSeries(i, j(1, 3))
If (b And br1) = 0
itr3s1++
tr3s1(itr3s1) = b
End If
End For
' Row 4
itr4s1 = 0
For i = NumKSeries(j(1, 4)) DownTo 1
b = KSeries(i, j(1, 4))
If (b And br1) = 0
itr4s1++
tr4s1(itr4s1) = b
End If
End For

Proc Row2
Local Int i
' Try all possible bit vectors r2:
For i = itr2s1 DownTo 1
r2 = tr2s1(i)
br2 = br1 Or r2
RedColsStep1
Row3
End For

Proc RedColsStep1
Local Int   i
Local Int32 b
' ------------ Col 2
itc2s1 = 0
For i = NumKSeries(j(2, 1)) DownTo 1
b = KSeries(i, j(2, 1))
If (b And br1) = 0
If SingleBit(b And r2)
itc2s1++
tc2s1(itc2s1) = b
End If
End If
End For
' ------------ Col 3
itc3s1 = 0
For i = NumKSeries(j(3, 1)) DownTo 1
b = KSeries(i, j(3, 1))
If (b And br1) = 0
If SingleBit(b And r2)
itc3s1++
tc3s1(itc3s1) = b
End If
End If
End For
' ------------ Col 4
itc4s1 = 0
For i = NumKSeries(j(4, 1)) DownTo 1
b = KSeries(i, j(4, 1))
If (b And br1) = 0
If SingleBit(b And r2)
itc4s1++
tc4s1(itc4s1) = b
End If
End If
End For

Proc Row3
Local Int i
For i = itr3s1 DownTo 1
r3 = tr3s1(i)
If (r3 And br2) = 0
br3 = br2 Or r3
RedColsStep2
Row4
End If
End For

Proc RedColsStep2
Local Int   i
' ------------ Col 2 reduce
itc2s2 = 0
For i = itc2s1 DownTo 1
If (tc2s1(i) And r3)
itc2s2++
tc2s2(itc2s2) = tc2s1(i)
End If
End For
' ------------ Col 3 reduce
itc3s2 = 0
For i = itc3s1 DownTo 1
If (tc3s1(i) And r3)
itc3s2++
tc3s2(itc3s2) = tc3s1(i)
End If
End For
' ------------ Col 4 reduce
itc4s2 = 0
For i = itc4s1 DownTo 1
If (tc4s1(i) And r3)
itc4s2++
tc4s2(itc4s2) = tc4s1(i)
End If
End For

Proc Row4
Local Int i
For i = itr4s1 DownTo 1
r4 = tr4s1(i)
If (r4 And br3) = 0
br4 = br3 Or r4
r5 = br4 Xor AllBits
RedColsStep3
Col2
End If
End For

Proc RedColsStep3
Local Int   i
Local Int32 b
' ------------ Col 3 reduce
itc3s3 = 0
For i = itc3s2 DownTo 1
b = tc3s2(i)
If (b And r4)
If (b And r5)
itc3s3++
tc3s3(itc3s3) = b
End If
End If
End For
' ------------ Col 4 reduce
itc4s3 = 0
For i = itc4s2 DownTo 1
b = tc4s2(i)
If (b And r4)
If (b And r5)
itc4s3++
tc4s3(itc4s3) = b
End If
End If
End For

Proc Col2
Local Int i
For i = itc2s2 DownTo 1
c2 = tc2s2(i)
If (c2 And r4)
If (c2 And r5)
bc2 = br1 Or c2
Col3
End If
End If
End For

Proc Col3
Local Int i
For i = itc3s3 DownTo 1
c3 = tc3s3(i)
If (c3 And bc2) = 0
bc3 = bc2 Or c3
Col4
End If
End For

Proc Col4
Local Int i
For i = itc4s3 DownTo 1
c4 = tc4s3(i)
If (c4 And bc3) = 0
cnt++                 // A semi-magic square is found. Increment the counter cnt.
If cnt >= Trigger     // Output of a square only in certain intervalls:
Col5
Trigger += Intervall
End If
End If
End For

Proc Col5
' Calculating the bit vector c5, which automatically is correct:
bc4 = bc3 Or c4
c5 = bc4 Xor AllBits
sqSolution
sqPrint
sqTest
sqSave

Proc sqSolution
' A semi-magic square is found.
' The array j(x,y) is filled with the entries of this square.
' Determine the entries as intersections of rows and columns
j(2, 2) = SingleBit(c2 And r2)
j(2, 3) = SingleBit(c2 And r3)
j(2, 4) = SingleBit(c2 And r4)
j(2, 5) = SingleBit(c2 And r5)
j(3, 2) = SingleBit(c3 And r2)
j(3, 3) = SingleBit(c3 And r3)
j(3, 4) = SingleBit(c3 And r4)
j(3, 5) = SingleBit(c3 And r5)
j(4, 2) = SingleBit(c4 And r2)
j(4, 3) = SingleBit(c4 And r3)
j(4, 4) = SingleBit(c4 And r4)
j(4, 5) = SingleBit(c4 And r5)
j(5, 2) = SingleBit(c5 And r2)
j(5, 3) = SingleBit(c5 And r3)
j(5, 4) = SingleBit(c5 And r4)
j(5, 5) = SingleBit(c5 And r5)

Proc sqTest
Local Int i, x, y, s
Local Int st(64)
' Check if each number from 1 to mx occurs exactly one time:
For i = 1 To mx
st(i) = 0
End For
For y = 1 To n
For x = 1 To n
st(j(x, y))++
End For
End For
For i = 1 To mx
If st(i) <> 1 Then Print "Number Error"; i : Stop
End For
' Check the sum of all rows
For y = 1 To n
s = 0
For x = 1 To n
s += j(x, y)
End For
If s <> mc Then Print "Sum Error in row"; y : Stop
End For
' Check the sum of all cols
For x = 1 To n
s = 0
For y = 1 To n
s += j(x, y)
End For
If s <> mc Then Print "Sum Error in col"; x : Stop
End For

Proc sqPrint
' A semi-magic square is shown in the window:
Local Int x, y
Print AT(1, 10); cnt
Print
For y = 1 To n
For x = 1 To n
Print Str(j(x, y), 4);
End For
Print
End For
DoEvents

Proc sqSave
' A semi-magic square is saved in a file:
Local Int x, y
Open "Semi-magic 5x5-squares.txt" for Append As # 1
DoEvents
Print # 1, cnt
Print # 1
For y = 1 To n
For x = 1 To n
Print # 1, Str(j(x, y), 4);
End For
Print # 1
End For
Print # 1
Close # 1
DoEvents

Function SingleBit(ByVal b As Int32) As Int
' This function returns the number of the set bit (1 to 31)
' of a bit vector where exactly one bit is set.
' Otherwise the function returns 0.
' For 8x8-squares you need uint64 and k = 32 and i = 16
Local Int k = 16, i = 8
If b = 0      Then Return 0
If b = OBI(k) Then Return k
Repeat
If b < OBI(k)
k -= i
Else
k += i
End If
If b = OBI(k) Then Return k
i >>= 1        // shift bits right
Until i = 0
Return 0
End Func

Proc sqEnd
' The calculation is done. Waiting until the window is closed.
Do : Sleep : Until Me Is Nothing : End

Sub Win_1_Close(Cancel?)
' Closing the window will terminate the program.
CloseW 1 : End