Zwei Zeilen einer Matrix sind orthogonal, wenn ihr inneres Produkt gleich Null ist. Nennen Sie eine Matrix mit allen Zeilen paarweise orthogonal eine orthogonale Matrix . Eine zirkulierende Matrix ist eine Matrix, bei der jeder Zeilenvektor relativ zum vorhergehenden Zeilenvektor um ein Element nach rechts gedreht wird. Wir werden nur an Matrizen interessiert sein, bei denen die Einträge entweder -1oder sind 1.

Aufgabe

Schreiben Sie Code, um in 2 Minuten (für gerade ) so viele verschiedene n/2durch northogonale, zirkulierende Matrizen wie möglich zu zählen n.

Eingang

Der Code hat keine Eingabe. Es kann alle gleichmäßigen Werte ausprobieren, die nes mag. Zum Beispiel könnte der Code alles versuchen n, was Multiplikationen 4von der kleinsten sind, und auch versuchen n = 2.

Ausgabe

Die Anzahl der gefundenen orthogonalen Zirkulationsmatrizen. Ihr Code sollte auch einen einfachen Schalter enthalten, damit er die Matrizen selbst ausgeben kann.

Ergebnis

Die Anzahl der gefundenen zirkulierenden Matrizen.

Hinweise

Orthogonal n/2durch nzirkulierende Matrizen existieren nur, wenn nein Vielfaches von 4oder nkleiner als ist 4.

Eine beispielhafte orthogonale Zirkulationsmatrix ist:

-1  1 -1 -1
-1 -1  1 -1

Tipps für einen naiven Ansatz

Der naivste Ansatz besteht darin, alle möglichen Matrizen zu durchlaufen. Dies kann anhand der folgenden beiden Beobachtungen beschleunigt werden.

• Um die Orthogonalität einer zirkulierenden Matrix zu testen, müssen wir nur jede Zeile mit der ersten vergleichen. Dies ist im Beispielcode implementiert.

• Wir können über Lyndon-Wörter iterieren und dann, wenn wir eine orthogonale Matrix finden, mit der Anzahl möglicher Umdrehungen multiplizieren. Diese Idee ist noch nicht getestet und kann daher fehlerhaft sein.

Beispielcode

Dies ist eine sehr einfache und naive Python-Antwort. Ich habe es mit ausgeführt timeout 120.

import itertools
def check_orthogonal(row):
    for i in xrange(1,int(n/2)):
        if (sum(row[j]*row[(j+i) % n] for j in xrange(n)) != 0):
            return False
    return True

counter = 0
for n in xrange(4,33,4):
    for row in itertools.product([-1,1],repeat = n):
        if check_orthogonal(row):
            counter +=1
            print "Counter is ", counter, ". n = ", n

Richtigkeitsprüfungen

Für n = 4,8,12,16,20,24,28sollte die Anzahl der unterschiedlichen Matrices auf Sie erhalten 12,40,144,128,80,192,560sind.

Ebenen der Großartigkeit

Nach dem Beispielcode zu urteilen, präsentiere ich hiermit zwei Ebenen der Attraktivität, die jede Antwort erreichen kann.

• Die Attraktivität des Silberlevels wird durch eine Punktzahl oder 1156 erreicht .

• Das Goldniveau der Attraktivität ist es, höher zu werden.

Sprachen und Bibliotheken

Sie können jede Sprache oder Bibliothek verwenden, die Sie mögen (die nicht für diese Herausforderung entwickelt wurde). Zum Zwecke der Bewertung werde ich Ihren Code jedoch auf meinem Computer ausführen. Geben Sie daher bitte klare Anweisungen zur Ausführung unter Ubuntu.

Meine Maschine Die Timings werden auf meiner Maschine ausgeführt. Dies ist eine Standard-Ubuntu-Installation auf einem 8 GB AMD FX-8350 Eight-Core-Prozessor. Dies bedeutet auch, dass ich Ihren Code ausführen kann.

Führende Antworten

• 332 von errorr in Octave

• 404 von RT in Python

• 744 durch Probenlösung unter Verwendung von Pypy

• 1156 von Thomas Kwa mit Java . Silberlevel-Attraktivität!

• 1588 von Reimer Behrends in OCaml . Goldlevel-Attraktivität!

OCaml, 1588 (n = 36)

Diese Lösung verwendet den üblichen Bitmusteransatz, um Vektoren von -1s und 1s darzustellen. Das Skalarprodukt wird wie üblich berechnet, indem das xor von zwei Bitvektoren genommen und n / 2 subtrahiert wird. Vektoren sind orthogonal, wenn für ihr xor genau n / 2 Bits gesetzt sind.

Lyndon-Wörter sind an sich nicht als normalisierte Darstellung nützlich, da sie jedes Muster ausschließen, das eine Rotation von sich selbst ist. Sie sind auch relativ teuer zu berechnen. Daher verwendet dieser Code eine etwas einfachere Normalform, die erfordert, dass die längste aufeinanderfolgende Folge von Nullen nach der Drehung (oder eine von ihnen, wenn es mehrere gibt) die höchstwertigen Bits belegen muss. Daraus folgt, dass das niedrigstwertige Bit immer 1 ist.

Beachten Sie auch, dass jeder Kandidatenvektor mindestens n / 4 (und höchstens 3n / 4) haben muss. Wir betrachten daher nur Vektoren mit n / 4 ... n / 2 gesetzten Bits, da wir andere über Komplement und Rotation ableiten können (in der Praxis scheinen alle diese Vektoren zwischen n / 2-2 und n / 2 + 2 zu liegen , aber das scheint auch schwer zu beweisen zu sein).

Wir erstellen diese Normalformen vom niedrigstwertigen Bit an und beachten dabei die Einschränkung, dass alle verbleibenden Nullenläufe (im Code als "Lücken" bezeichnet) unseren normalen Formularanforderungen entsprechen müssen. Insbesondere muss, solange mindestens ein weiteres 1-Bit platziert werden muss, Platz für die aktuelle Lücke vorhanden sein und eine andere, die mindestens so groß ist wie die aktuelle Lücke oder eine andere bisher beobachtete Lücke.

Wir stellen auch fest, dass die Liste der Ergebnisse klein ist. Daher versuchen wir nicht, Duplikate während des Erkennungsprozesses zu vermeiden, sondern zeichnen die Ergebnisse einfach in Worker-Sets auf und berechnen die Vereinigung dieser Sets am Ende.

Es ist erwähnenswert, dass die Laufzeitkosten des Algorithmus immer noch exponentiell und mit einer Geschwindigkeit steigen, die mit der der Brute-Force-Version vergleichbar ist. Dies bringt uns im Wesentlichen eine Reduzierung um einen konstanten Faktor und geht zu Lasten eines Algorithmus, der schwieriger zu parallelisieren ist als die Brute-Force-Version.

Ausgabe für n bis 40:

 4: 12
 8: 40
12: 144
16: 128
20: 80
24: 192
28: 560
32: 0
36: 432
40: 640

Das Programm ist in OCaml geschrieben und kann kompiliert werden mit:

ocamlopt -inline 100 -nodynlink -o orthcirc unix.cmxa bigarray.cmxa orthcirc.ml

Führen Sie aus, um ./orthcirc -helpzu sehen, welche Optionen das Programm unterstützt.

Bei Architekturen, die dies unterstützen, -fno-PICkann dies zu einem kleinen zusätzlichen Leistungsgewinn führen.

Dies ist für OCaml 4.02.3 geschrieben, funktioniert aber möglicherweise auch mit älteren Versionen (sofern diese nicht zu alt sind).

UPDATE: Diese neue Version bietet eine bessere Parallelisierung. Beachten Sie, dass p * (n/4 + 1)pro Instanz des Problems Arbeitsthreads verwendet werden und einige von ihnen immer noch erheblich kürzer als andere ausgeführt werden. Der Wert von pmuss eine Potenz von 2 sein. Die Beschleunigung bei 4-8 Kernen ist minimal (möglicherweise um 10%), lässt sich jedoch bei großen Kernen besser auf eine große Anzahl von Kernen skalieren n.

let max_n = ref 40
let min_n = ref 4
let seq_mode = ref false
let show_res = ref false
let fanout = ref 8

let bitcount16 n =
  let b2 n = match n land 3 with 0 -> 0 | 1 | 2 -> 1 | _ -> 2 in
  let b4 n = (b2 n) + (b2 (n lsr 2)) in
  let b8 n = (b4 n) + (b4 (n lsr 4)) in
  (b8 n) + (b8 (n lsr 8))

let bitcount_data =
  let open Bigarray in
  let tmp = Array1.create int8_signed c_layout 65536 in
    for i = 0 to 65535 do
      Array1.set tmp i (bitcount16 i)
    done;
    tmp

let bitcount n =
  let open Bigarray in
  let bc n = Array1.unsafe_get bitcount_data (n land 65535) in
  (bc n) + (bc (n lsr 16)) + (bc (n lsr 32)) + (bc (n lsr 48))

module IntSet = Set.Make (struct
  type t = int
  let compare = Pervasives.compare
end)

let worker_results = ref IntSet.empty

let test_row vec row mask n =
  bitcount ((vec lxor (vec lsr row) lxor (vec lsl (n-row))) land mask) * 2 = n

let record vec len n =
  let m = (1 lsl n) - 1 in
  let rec test_orth_circ ?(row=2) vec m n =
    if 2 * row >= n then true
    else if not (test_row vec row m n) then false
    else test_orth_circ ~row:(row+1) vec m n
  in if test_row vec 1 m n &&
        test_orth_circ vec m n then
  begin
    for i = 0 to n - 1 do
      let v = ((vec lsr i) lor (vec lsl (n - i))) land m in
      worker_results := IntSet.add v !worker_results;
      worker_results := IntSet.add (v lxor m) !worker_results
    done
  end

let show vec n =
  for i = 0 to n / 2 - 1 do
    let vec' = (vec lsr i) lor (vec lsl (n - i)) in
    for j = 0 to n-1 do
      match (vec' lsr (n-j)) land 1 with
      | 0 -> Printf.printf "  1"
      | _ -> Printf.printf " -1"
    done; Printf.printf "\n"
  done; Printf.printf "\n"; flush stdout

let rec build_normalized ~prefix ~plen ~gap ~maxgap ~maxlen ~bits ~fn =
  if bits = 0 then
    fn prefix plen maxlen
  else begin
    let room = maxlen - gap - plen - bits in
    if room >= gap && room >= maxgap then begin
      build_normalized
        ~prefix:(prefix lor (1 lsl (plen + gap)))
        ~plen:(plen + gap + 1)
        ~gap:0
        ~maxgap:(if gap > maxgap then gap else maxgap)
        ~maxlen
        ~bits:(bits - 1)
        ~fn;
      if room > gap + 1 && room > maxgap then
        build_normalized ~prefix ~plen ~gap:(gap + 1) ~maxgap ~maxlen ~bits ~fn
    end
  end

let rec log2 = function
| 0 -> -1
| n -> 1 + (log2 (n lsr 1))

let rec test_gap n pat =
  if n land pat = 0 then true
  else if pat land 1 = 0 then test_gap n (pat lsr 1)
  else false

let rec test_gaps n maxlen len =
  let fill k = (1 lsl k) -1 in
  if len = 0 then []
  else if test_gap n ((fill maxlen) lxor (fill (maxlen-len))) then
    len :: (test_gaps n maxlen (len-1))
  else test_gaps n maxlen (len-1)

let rec longest_gap n len =
  List.fold_left max 0 (test_gaps n len len)

let start_search low lowbits maxlen bits fn =
  let bits = bits - (bitcount low) in
  let plen = log2 low + 1 in
  let gap = lowbits - plen in
  let maxgap = longest_gap low lowbits in
  worker_results := IntSet.empty;
  if bits >= 0 then
    build_normalized ~prefix:low ~plen ~gap ~maxgap ~maxlen ~bits ~fn;
  !worker_results

let spawn f x =
  let open Unix in
  let safe_fork () = try fork() with _ -> -1 in
  let input, output = pipe () in
  let pid = if !seq_mode then -1 else safe_fork() in
  match pid with
  | -1 -> (* seq_mode selected or fork() failed *)
     close input; close output; (fun () -> f x)
  | 0 -> (* child process *)
    close input;
    let to_parent = out_channel_of_descr output in
    Marshal.to_channel to_parent (f x) [];
    close_out to_parent; exit 0
  | pid -> (* parent process *)
    close output;
    let from_child = in_channel_of_descr input in
    (fun () -> 
      ignore (waitpid [] pid);
      let result = Marshal.from_channel from_child in
      close_in from_child; result)

let worker1 (n, k) =
  start_search 1 1 n k record

let worker2 (n, k, p) =
  start_search (p * 2 + 1) (log2 !fanout + 1) n k record

let spawn_workers n =
  let queue = Queue.create () in
  if n = 4 || n = 8 then begin
    for i = n / 4 to n / 2 do
      Queue.add (spawn worker1 (n, i)) queue
    done
  end else begin
    for i = n / 2 downto n / 4 do
      for p = 0 to !fanout - 1 do
        Queue.add (spawn worker2 (n, i, p)) queue
      done
    done
  end;
  Queue.fold (fun acc w -> IntSet.union acc (w())) IntSet.empty queue

let main () =
  if !max_n > 60 then begin
    print_endline "error: cannot handle n > 60";
    exit 1
  end;
  min_n := max !min_n 4;
  if bitcount !fanout <> 1 then begin
    print_endline "error: number of threads must be a power of 2";
    exit 1;
  end;
  for n = !min_n to !max_n do
    if n mod 4 = 0 then
    let result = spawn_workers n in
    Printf.printf "%2d: %d\n" n (IntSet.cardinal result);
    if !show_res then
      IntSet.iter (fun v -> show v n) result;
    flush stdout
  done

let () =
  let args =[("-m", Arg.Set_int min_n, "min size of the n by n/2 matrix");
             ("-n", Arg.Set_int max_n, "max size of the n by n/2 matrix");
             ("-p", Arg.Set_int fanout, "parallel fanout");
             ("-seq", Arg.Set seq_mode, "run in single-threaded mode");
             ("-show", Arg.Set show_res, "display list of results") ] in
  let usage = ("Usage: " ^
               (Filename.basename Sys.argv.(0)) ^
               " [-n size] [-seq] [-show]") in
  let error _ = Arg.usage args usage; exit 1 in
  Arg.parse args error usage;
  main ()

Java, 1156 Matrizen

Dies verwendet eine ziemlich naive Bitmaske und dauert auf meinem Computer weniger als 15 Sekunden für n = 28.

Zirkulierende Matrizen werden durch ihre ersten Reihen bestimmt. Daher stelle ich die ersten Zeilen der Matrizen als Bitvektoren dar: 1 und 0 stellen -1 und 1 dar. Zwei Zeilen sind orthogonal, wenn die Anzahl der gesetzten Bits, wenn sie zusammen xor'd sind, n / 2 beträgt.

import java.util.Arrays;

class Main {

    static void dispmatrix(long y,int N)
    {
        int[][] arr = new int[N/2][N];
        for(int row=0; row < N/2; row++)
        {
            for(int col=0; col < N; col++)
            {
                arr[row][col] = (int) ((y >>> (N+row-col)) & 1L);
            }
        }
        System.out.println(Arrays.deepToString(arr));
    }

  public static void main(String[] args) {

    int count = 0;
    boolean good;
    boolean display = false;
    long y;
    for(int N=4; N <= 28 ;N += 4)
    {
    long mask = (1L << N) - 1;
    for(long x=0; x < (1L<<N); x++)
    {
        good = true;
        y = x + (x << N);

        for(int row = 1; row < N/2; row++)
        {
            good &= N/2 == Long.bitCount((y ^ (y >>> row)) & mask);
        }

        if(good)
        {
            if(display) dispmatrix(y,N);
            count++;
        }

    }
    System.out.println(count);
    }
  }
}

Ich kann Eclipse momentan nicht zum Laufen bringen, daher wurde dies auf repl.it getestet.

Hier ist die Anzahl der ersten Zeilen, die orthogonal zu den ersten r Zeilen danach für n = 28 sind:

[268435456, 80233200, 23557248, 7060320, 2083424, 640304, 177408, 53088, 14896, 4144, 2128, 1008, 1008, 560]

Optimierungen:

• Da sich die Orthogonalität bei einer zyklischen Drehung beider Zeilen nicht ändert, müssen wir nur überprüfen, ob die erste Zeile orthogonal zum Rest ist.
• Anstatt N / 2-mal manuell eine zyklische N-Bit-Verschiebung durchzuführen, speichere ich die zu verschiebenden Bits oben auf der longund verwende dann eine einzelne andmit den unteren N-Bits, um die benötigten zu extrahieren.

Mögliche weitere Optimierungen:

• Generieren Sie Lyndon-Wörter. Nach meinen Berechnungen ist dies nur dann sinnvoll, wenn Lyndon-Wörter in jeweils weniger als ~ 1000 Zyklen generiert werden können.
• Wenn Lyndon-Wörter zu lange dauern, können wir die Bitvektoren, die beginnen, immer noch nur überprüfen 00und ihre Rotationen (und Rotationen des NOT) verwenden, wenn wir eine orthogonale Matrix finden. Wir können dies tun, weil 0101...01und 1010...10sind keine ersten Zeilen möglich, und alle anderen enthalten entweder a 00oder a 11.
• Verzweigen Sie sich möglicherweise zur Hälfte, wenn die Matrix wahrscheinlich orthogonal ist. Ich weiß jedoch nicht, wie viel Verzweigung kosten wird, also muss ich testen.
• Wenn das oben genannte funktioniert, beginnen Sie mit einer anderen Zeile als Zeile 1?
• Vielleicht gibt es einen mathematischen Weg, um einige Möglichkeiten auszuschließen. Ich weiß nicht, was das sein würde.
• Schreiben Sie natürlich in C ++.

Python (404 Matrizen auf i5-5300U)

Dies wird hauptsächlich als Ausgangspunkt für Verbesserungen anderer verwendet. Dies kann viel aufgeräumt, parallelisiert usw. werden.

import numpy
import itertools
import time

def findCirculantOrthogonalRows(n, t1, timeLimit):
  validRows = []
  testMatrix = numpy.zeros((n//2, n))
  identityMatrixScaled = n*numpy.identity(n//2)
  outOfTime = False
  for startingRowTuple in itertools.product([1, -1], repeat=n):
     for offset in range(n//2):
       for index in range(n):
         testMatrix[offset][index] = startingRowTuple[(index-offset) % n]

     if(numpy.array_equal(identityMatrixScaled, testMatrix.dot(testMatrix.transpose()))):
      validRows.append(startingRowTuple)

    if(time.clock() - t1 >= timeLimit):
      outOfTime = True
      break

  return (validRows, outOfTime)

n = 4
validFirstRows = []
t1 = time.clock()
timeLimit = 120
fullOutput = True

while(True):
  print('calling with', n)
  (moreRows, outOfTime) = findCirculantOrthogonalRows(n, t1, timeLimit)

  if(len(moreRows) > 0):
    validFirstRows.extend(moreRows)

  if(outOfTime == True):
    break

  n += 4

print('Found', len(validFirstRows), 'circulant orthogonal matrices in', timeLimit, 'seconds')

if(fullOutput):
  counter = 1
  for r in validFirstRows:
    n = len(r)
    matrix = numpy.zeros((n//2, n))
    for offset in range(n//2):
      for index in range(n):
        matrix[offset][index] = r[(index-offset) % n]
    print('matrix #', counter, ':\n', matrix)
    counter += 1
    input()

Matlab / Octave, 381/328 Matrizen

Auch nur der naive Ansatz, jede mögliche Kombination auszuprobieren.

counter = 0;
%ok: 2,4,8
%none: 
tic
for n=[2,4,8,12,16,20];
    m=n/2;
    N=2^n-1;
    for k=1:N

        k/N; %remove ; in order to see progress
        v=(dec2bin(k,n)-'0')*2-1;               %first row

        z=toeplitz([v(1) fliplr(v(m+2:n))], v); %create circulante matrix
        w = z * z.';
        if norm(w - diag(diag(w))) < eps
            counter = counter+1;
            %n    %remove % in order to see the matrices
            %z
        end
        if toc>120;
            break;
        end
    end
end
counter