' 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