--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <time.h>
+#include "complex.c"
+#include "gates.c"
+#include "display.c"
+#include "chacha20.c"
+#define QUBITS_MAX 11
+
+typedef struct
+{
+ char n[128];
+ uint8_t q0, q1, q2;
+ double arg0, arg1, arg2;
+} QInstr;
+
+double qansel_rand()
+{
+ static uint64_t hiddenVariable = 0;
+ static uint32_t blockNumber = 0;
+ if (hiddenVariable == 0)
+ {
+ struct timespec ts;
+ clock_gettime(CLOCK_MONOTONIC, &ts);
+ hiddenVariable = (uint64_t)ts.tv_sec * 1000000000LL + ts.tv_nsec;
+ }
+ uint8_t key[32];
+ uint8_t nonce[12];
+ uint64_t tmpVariable;
+ for (uint8_t i = 0; i < 32; i++)
+ {
+ if (i % 8 == 0) tmpVariable = hiddenVariable;
+ key[i] = tmpVariable & 0xFF;
+ tmpVariable >> 8;
+
+ }
+ for (uint8_t i = 0; i < 12; i++)
+ {
+ nonce[i] = 0;
+ }
+ uint8_t* block = foleo_chacha20(key, nonce, blockNumber++, 4);
+ uint32_t num = 0;
+ for (uint8_t i = 0; i < 4; i++)
+ {
+ num = (num << 8) | block[i];
+ }
+ return ((double)num) / ((double)UINT32_MAX);
+}
+
+void qansel_cnot(cpx_mtx_t* stateVector, uint8_t qubitCount, uint8_t bitA, uint8_t bitB)
+{
+ uint32_t retLen = (uint8_t)pow(2, qubitCount);
+ cpx_mtx_t ret;
+ cpx_mtx_init(&ret, 1, retLen);
+ cpx_t n;
+ for (uint32_t i = 0; i < retLen; i++)
+ {
+ uint8_t bitAVal = (i >> bitA) & 1;
+ uint8_t bitBVal = (i >> bitB) & 1;
+ uint8_t bitBNew = bitAVal ? !bitBVal : bitBVal;
+ uint32_t j = (i & ~(1 << bitB)) | (bitBNew << bitB);
+ cpx_mtx_get(stateVector, 0, i, &n);
+ cpx_mtx_set(&ret, 0, j, &n);
+ }
+ cpx_mtx_free(stateVector);
+ stateVector->ptr = ret.ptr;
+ stateVector->rows = ret.rows;
+ stateVector->cols = ret.cols;
+}
+
+void qansel_swap(cpx_mtx_t* stateVector, uint8_t qubitCount, uint8_t bitA, uint8_t bitB)
+{
+ uint32_t retLen = (uint8_t)pow(2, qubitCount);
+ cpx_mtx_t ret;
+ cpx_mtx_init(&ret, 1, retLen);
+ cpx_t n;
+ for (uint32_t i = 0; i < retLen; i++)
+ {
+ uint8_t bitAVal = (i >> bitA) & 1;
+ uint8_t bitBVal = (i >> bitB) & 1;
+ uint8_t bitANew = bitBVal;
+ uint8_t bitBNew = bitAVal;
+ uint32_t j = (i & ~((1 << bitA) | (1 << bitB))) | ((bitANew << bitA) | (bitBNew << bitB));
+ cpx_mtx_get(stateVector, 0, i, &n);
+ cpx_mtx_set(&ret, 0, j, &n);
+ }
+ cpx_mtx_free(stateVector);
+ stateVector->ptr = ret.ptr;
+ stateVector->rows = ret.rows;
+ stateVector->cols = ret.cols;
+}
+
+void qansel_fredkin(cpx_mtx_t* stateVector, uint8_t qubitCount, uint8_t bitA, uint8_t bitB, uint8_t bitC)
+{
+ uint32_t retLen = (uint8_t)pow(2, qubitCount);
+ cpx_mtx_t ret;
+ cpx_mtx_init(&ret, 1, retLen);
+ cpx_t n;
+ for (uint32_t i = 0; i < retLen; i++)
+ {
+ uint8_t bitAVal = (i >> bitA) & 1;
+ uint8_t bitBVal = (i >> bitB) & 1;
+ uint8_t bitCVal = (i >> bitC) & 1;
+ uint8_t bitBNew = bitAVal ? bitCVal : bitBVal;
+ uint8_t bitCNew = bitAVal ? bitBVal : bitCVal;
+ uint32_t j = (i & ~((1 << bitB) | (1 << bitC))) | ((bitBNew << bitB) | (bitCNew << bitC));
+ cpx_mtx_get(stateVector, 0, i, &n);
+ cpx_mtx_set(&ret, 0, j, &n);
+ }
+ cpx_mtx_free(stateVector);
+ stateVector->ptr = ret.ptr;
+ stateVector->rows = ret.rows;
+ stateVector->cols = ret.cols;
+}
+
+
+void qansel_toffoli(cpx_mtx_t* stateVector, uint8_t qubitCount, uint8_t bitA, uint8_t bitB, uint8_t bitC)
+{
+ uint32_t retLen = (uint8_t)pow(2, qubitCount);
+ cpx_mtx_t ret;
+ cpx_mtx_init(&ret, 1, retLen);
+ cpx_t n;
+ for (uint32_t i = 0; i < retLen; i++)
+ {
+ uint8_t bitAVal = (i >> bitA) & 1;
+ uint8_t bitBVal = (i >> bitB) & 1;
+ uint8_t bitCVal = (i >> bitC) & 1;
+ uint8_t bitCNew = (bitAVal && bitBVal) ? !bitCVal : bitCVal;
+ uint32_t j = (i & ~(1 << bitC)) | (bitCNew << bitC);
+ cpx_mtx_get(stateVector, 0, i, &n);
+ cpx_mtx_set(&ret, 0, j, &n);
+ }
+ cpx_mtx_free(stateVector);
+ stateVector->ptr = ret.ptr;
+ stateVector->rows = ret.rows;
+ stateVector->cols = ret.cols;
+}
+
+double* qansel_unitary(double theta, double phi, double lambda)
+{
+ cpx_mtx_t m;
+ cpx_t a, b, c, d;
+ a.real = cos(theta/2.0);
+ a.imaginary = 0;
+ b.real = -cos(lambda) * sin(theta/2.0);
+ b.imaginary = sin(lambda) * sin(theta/2.0);
+ c.real = cos(phi) * sin(theta/2.0);
+ c.imaginary = sin(phi) * sin(theta/2.0);
+ d.real = cos(phi + lambda) * cos(theta/2.0);
+ d.imaginary = sin(phi + lambda) * cos(theta/2.0);
+ cpx_mtx_init(&m, 2, 2);
+ cpx_mtx_set(&m, 0, 0, &a);
+ cpx_mtx_set(&m, 0, 1, &b);
+ cpx_mtx_set(&m, 1, 0, &c);
+ cpx_mtx_set(&m, 1, 1, &d);
+ return m.ptr;
+}
+
+void qansel_instruction(cpx_mtx_t* stateVector, uint8_t qubitCount, QInstr* instr)
+{
+ cpx_mtx_t tmp;
+ cpx_mtx_t gate;
+ gate.rows = 2;
+ gate.cols = 2;
+ double* gate_ptr;
+ switch (instr->n[0])
+ {
+ case 'h': gate_ptr = Hadamard; break;
+ case 'x': gate_ptr = PauliX; break;
+ case 'y': gate_ptr = PauliY; break;
+ case 'z': gate_ptr = PauliZ; break;
+ case 's': gate_ptr = PhaseS; break;
+ case 't': gate_ptr = PhaseT; break;
+ case 'u':
+ gate_ptr = qansel_unitary(instr->arg0, instr->arg1, instr->arg2);
+ break;
+ default: gate_ptr = Identity; break;
+ }
+
+ cpx_t n;
+ cpx_mtx_t filter;
+ cpx_mtx_init(&filter, 2, 2);
+ uint8_t qubit = qubitCount - (instr->q0) - 1;
+ if (qubit == 0)
+ {
+ memcpy(filter.ptr, gate_ptr, 16 * sizeof(double));
+ }
+ else
+ {
+ memcpy(filter.ptr, Identity, 16 * sizeof(double));
+ }
+
+ for (uint8_t i = 1; i < qubitCount; i++)
+ {
+ cpx_mtx_init(&tmp, filter.rows * 2, filter.cols * 2);
+ if (qubit == i)
+ {
+ gate.ptr = gate_ptr;
+ }
+ else
+ {
+ gate.ptr = Identity;
+ }
+ cpx_mtx_knk(&tmp, &filter, &gate);
+ free(filter.ptr);
+ filter.ptr = tmp.ptr;
+ filter.rows = tmp.rows;
+ filter.cols = tmp.cols;
+ }
+
+ cpx_mtx_init(&tmp, stateVector->rows, stateVector->cols);
+ cpx_mtx_mul(&tmp, stateVector, &filter);
+ free(stateVector->ptr);
+ stateVector->ptr = tmp.ptr;
+ free(filter.ptr);
+ if (instr->n[0] == 'u') free(gate_ptr);
+}
+
+
+uint8_t qansel_measure(cpx_mtx_t* stateVector, uint8_t qubitCount, uint8_t qubit)
+{
+ uint32_t qubitCountPow2 = (uint32_t)pow(2, qubitCount);
+ cpx_t n;
+ double prob0 = 0;
+ for (uint32_t i = 0; i < qubitCountPow2; i++)
+ {
+ uint8_t bit = (i >> qubit) & 1;
+ cpx_mtx_get(stateVector, 0, i, &n);
+ if (bit == 0) prob0 += cpx_magsqr(&n);
+ }
+
+ double r = qansel_rand();
+ uint8_t newBit = r < prob0 ? 0 : 1;
+ double probTot = 0;
+ for (uint32_t i = 0; i < qubitCountPow2; i++)
+ {
+ uint8_t bit = (i >> qubit) & 1;
+ cpx_mtx_get(stateVector, 0, i, &n);
+ if (bit != newBit)
+ {
+ n.real = 0;
+ n.imaginary = 0;
+ }
+ else
+ {
+ probTot += cpx_magsqr(&n);
+ }
+ cpx_mtx_set(stateVector, 0, i, &n);
+ }
+
+ double multiplier = sqrt(1 / probTot);
+ for (uint32_t i = 0; i < qubitCountPow2; i++)
+ {
+ uint8_t bit = (i >> qubit) & 1;
+ cpx_mtx_get(stateVector, 0, i, &n);
+ if (bit == newBit)
+ {
+ n.real *= multiplier;
+ n.imaginary *= multiplier;
+ }
+ cpx_mtx_set(stateVector, 0, i, &n);
+ }
+
+ return newBit;
+}
+
+void qansel_run(uint8_t qubitCount, uint8_t bitCount, QInstr* instr, uint32_t instrLen, uint8_t gfx)
+{
+ uint32_t qubitCountPow2 = (uint32_t)pow(2, qubitCount);
+
+ uint8_t bitVector[bitCount];
+ cpx_mtx_t stateVector;
+ cpx_mtx_init(&stateVector, 1, qubitCountPow2);
+ cpx_mtx_set2(&stateVector, 0, 0, 1, 0);
+ if (gfx) display(&stateVector, qubitCount);
+ uint8_t flags;
+
+ for (uint8_t i = 0; i < bitCount; i++) bitVector[i] = 0;
+ for (uint32_t i = 0; i < instrLen; i++)
+ {
+ if (strcmp(instr[i].n, "measure") == 0)
+ {
+ bitVector[instr[i].q1] = qansel_measure(&stateVector, qubitCount, instr[i].q0);
+ }
+ else if (strcmp(instr[i].n, "cswap") == 0)
+ {
+ qansel_fredkin(&stateVector, qubitCount, instr[i].q0, instr[i].q1, instr[i].q2);
+ }
+ else if (strcmp(instr[i].n, "ccx") == 0)
+ {
+ qansel_toffoli(&stateVector, qubitCount, instr[i].q0, instr[i].q1, instr[i].q2);
+ }
+ else if (strcmp(instr[i].n, "cx") == 0)
+ {
+ qansel_cnot(&stateVector, qubitCount, instr[i].q0, instr[i].q1);
+ }
+ else if (strcmp(instr[i].n, "swap") == 0)
+ {
+ qansel_swap(&stateVector, qubitCount, instr[i].q0, instr[i].q1);
+ }
+ else if (strcmp(instr[i].n, "if_all") == 0)
+ {
+ uint8_t val = 0;
+ for (int32_t j = bitCount - 1; j >= 0; j--)
+ {
+ val = (val << 1) | bitVector[j];
+ }
+ if (val != instr[i].q0) i++;
+ }
+ else if (strcmp(instr[i].n, "if") == 0)
+ {
+ if (bitVector[instr[i].q0] != instr[i].q1) i++;
+ }
+ else if (strcmp(instr[i].n, "printq_all") == 0)
+ {
+ printf("[ "); cpx_mtx_print(&stateVector); printf(" ]\n");
+ }
+ else if (strcmp(instr[i].n, "printc_all") == 0)
+ {
+ for (int32_t j = bitCount - 1; j >= 0; j--)
+ {
+ putchar('0' + bitVector[j]);
+ }
+ putchar('\n');
+ }
+ else if (strcmp(instr[i].n, "printq") == 0)
+ {
+ cpx_mtx_t tmp;
+ cpx_mtx_init(&tmp, 1, 2);
+ for (uint32_t j = 0; j < qubitCountPow2; j++)
+ {
+ if ((j >> instr[i].q0) & 1)
+ {
+ cpx_t a, b;
+ cpx_mtx_get(&tmp, 0, 1, &a);
+ cpx_mtx_get(&stateVector, 0, j, &b);
+ a.real += b.real;
+ a.imaginary += b.imaginary;
+ cpx_mtx_set(&tmp, 0, 1, &a);
+ }
+ else
+ {
+ cpx_t a, b;
+ cpx_mtx_get(&tmp, 0, 0, &a);
+ cpx_mtx_get(&stateVector, 0, j, &b);
+ a.real += b.real;
+ a.imaginary += b.imaginary;
+ cpx_mtx_set(&tmp, 0, 0, &a);
+ }
+ }
+ double multiplier = 0;
+ cpx_t n;
+ cpx_mtx_get(&tmp, 0, 0, &n);
+ multiplier += cpx_magsqr(&n);
+ cpx_mtx_get(&tmp, 0, 1, &n);
+ multiplier += cpx_magsqr(&n);
+ multiplier = sqrt(1 / multiplier);
+ n.real *= multiplier;
+ n.imaginary *= multiplier;
+ cpx_mtx_set(&tmp, 0, 1, &n);
+ cpx_mtx_get(&tmp, 0, 0, &n);
+ n.real *= multiplier;
+ n.imaginary *= multiplier;
+ cpx_mtx_set(&tmp, 0, 0, &n);
+
+ printf("[ "); cpx_mtx_print(&tmp); printf(" ]\n");
+ cpx_mtx_free(&tmp);
+ }
+ else if (strcmp(instr[i].n, "sample_all") == 0)
+ {
+ for (uint32_t j = 0; j < qubitCountPow2; j++)
+ {
+ uint32_t tmp = j;
+ for (uint8_t k = 0; k < qubitCount; k++)
+ {
+ putchar('0' + (tmp >> (qubitCount - 1) & 1));
+ tmp <<= 1;
+ }
+ cpx_t n;
+ cpx_mtx_get(&stateVector, 0, j, &n);
+ printf(": %.00f%%\n", cpx_magsqr(&n) * 100);
+ }
+ }
+ else if (strcmp(instr[i].n, "sample") == 0)
+ {
+ double prob = 0;
+ for (uint32_t j = 0; j < qubitCountPow2; j++)
+ {
+ cpx_t n;
+ cpx_mtx_get(&stateVector, 0, j, &n);
+ if ((j >> instr[i].q0) & 1)
+ {
+ prob += cpx_magsqr(&n);
+ }
+ }
+ printf("%.00f%%\n", prob * 100.0);
+ }
+ else if (strcmp(instr[i].n, "reset") == 0)
+ {
+ cpx_mtx_set2(&stateVector, 0, 0, 1, 0);
+ for (uint32_t j = 1; j < qubitCountPow2; j++)
+ {
+ cpx_mtx_set2(&stateVector, 0, j, 0, 0);
+ }
+ for (uint8_t j = 0; j < bitCount; j++)
+ {
+ bitVector[j] = 0;
+ }
+ }
+ else
+ {
+ qansel_instruction(&stateVector, qubitCount, instr + i);
+ }
+
+ if (gfx) display(&stateVector, qubitCount);
+ }
+ display(NULL, -1);
+
+ cpx_mtx_free(&stateVector);
+}
+
+void main()
+{
+ char** lines = malloc(0);
+ uint32_t* lineIDs = malloc(0);
+ char* text = malloc(0);
+ uint32_t textLen = 0;
+ uint32_t linesLen = 0;
+ int c;
+ int pc = -1;
+ uint8_t comment = 0;
+ uint8_t commentM = 0;
+ uint32_t lineID = 1;
+ uint8_t skipSpaces = 1;
+ uint8_t inGate = 0;
+ while ( (c = getchar()) != EOF )
+ {
+ if (c == '/' && commentM == 0 && comment == 0)
+ {
+ commentM = 1;
+ }
+ else if (c == '/' && commentM == 1 && comment == 0)
+ {
+ comment = 1;
+ commentM = 0;
+ }
+ else if (c == '\n')
+ {
+ comment = 0;
+ commentM = 0;
+ lineID += 1;
+ }
+ else if (comment || (c == ' ' && skipSpaces)) {}
+ else if (c != '\n' && c != '\t' && c != ';' && (c != ')' || inGate))
+ {
+ skipSpaces = 0;
+ if (c >= 'A' && c <= 'Z') c += 'a' - 'A';
+ if (c == 'u') inGate = 1;
+ text = realloc(text, textLen + 1);
+ text[textLen++] = c;
+ pc = c;
+ }
+ else if (c == ';' || (c == ')' && !inGate))
+ {
+ inGate = 0;
+ skipSpaces = 1;
+ if (c == ')')
+ {
+ text = realloc(text, textLen + 1);
+ text[textLen++] = ')';
+ }
+ text = realloc(text, textLen + 1);
+ text[textLen++] = 0;
+ lineIDs = realloc(lineIDs, (linesLen + 1) * sizeof(uint32_t));
+ lineIDs[linesLen] = lineID;
+ lines = realloc(lines, (linesLen + 1) * sizeof(char*));
+ lines[linesLen] = malloc(strlen(text) + 1);
+ strcpy(lines[linesLen++], text);
+ text = realloc(text, 0);
+ textLen = 0;
+ pc = ';';
+ }
+ }
+ text = realloc(text, textLen + 1);
+ text[textLen++] = 0;
+ if (strlen(text) > 0)
+ {
+ free(text);
+ fprintf(stderr, "QAnsel: Invalid trailing text.\n");
+ exit(1);
+ }
+ free(text);
+
+ uint8_t qubitCount = 0xFF;
+ uint8_t bitCount = 0xFF;
+ QInstr* instr = malloc(0);
+ uint32_t instrLen = 0;
+
+ uint8_t errFound = 0;
+ for (uint32_t i = 0; i < linesLen; i++)
+ {
+ lineID = i;
+
+ char g;
+ int q0, q1, q2;
+ float a0, a1, a2;
+
+ if (sscanf(lines[i], "qreg q[%i]", &q0) == 1)
+ {
+ if (qubitCount == 0xFF)
+ {
+ qubitCount = q0;
+ if (qubitCount < 1 || qubitCount > QUBITS_MAX)
+ {
+ fprintf(stderr, "QAnsel: Invalid count");
+ errFound = 1;
+ break;
+ }
+ }
+ else
+ {
+ fprintf(stderr, "QAnsel: Repeated initialization");
+ errFound = 1;
+ break;
+ }
+ }
+ else if (sscanf(lines[i], "creg c[%i]", &q0) == 1)
+ {
+ if (bitCount == 0xFF)
+ {
+ bitCount = q0;
+ if (bitCount < 1 || bitCount > QUBITS_MAX)
+ {
+ fprintf(stderr, "QAnsel: Invalid count");
+ errFound = 1;
+ break;
+ }
+ }
+ else
+ {
+ fprintf(stderr, "QAnsel: Repeated initialization");
+ errFound = 1;
+ break;
+ }
+ }
+ else if (sscanf(lines[i], "u(%f,%f,%f) q[%i]", &a0, &a1, &a2, &q0) == 4)
+ {
+ if (qubitCount == 0xFF)
+ {
+ fprintf(stderr, "QAnsel: Instruction before initialization.\n");
+ errFound = 1;
+ break;
+ }
+ if (q0 < 0 || q0 >= qubitCount)
+ {
+ fprintf(stderr, "QAnsel: Invalid index ");
+ errFound = 1;
+ break;
+ }
+ instr = realloc(instr, (instrLen + 1) * sizeof(QInstr));
+ instr[instrLen].n[0] = 'u';
+ instr[instrLen].n[1] = 0;
+ instr[instrLen].q0 = q0;
+ instr[instrLen].arg0 = a0;
+ instr[instrLen].arg1 = a1;
+ instr[instrLen++].arg2 = a2;
+ }
+ else if
+ (
+ sscanf(lines[i], "h q[%i]", &q0) == 1
+ || sscanf(lines[i], "x q[%i]", &q0) == 1
+ || sscanf(lines[i], "y q[%i]", &q0) == 1
+ || sscanf(lines[i], "z q[%i]", &q0) == 1
+ || sscanf(lines[i], "t q[%i]", &q0) == 1
+ || sscanf(lines[i], "z q[%i]", &q0) == 1
+ )
+ {
+ g = lines[i][0];
+ if (qubitCount == 0xFF)
+ {
+ fprintf(stderr, "QAnsel: Instruction before initialization.\n");
+ errFound = 1;
+ break;
+ }
+ if (q0 < 0 || q0 >= qubitCount)
+ {
+ fprintf(stderr, "QAnsel: Invalid index ");
+ errFound = 1;
+ break;
+ }
+ instr = realloc(instr, (instrLen + 1) * sizeof(QInstr));
+ instr[instrLen].n[0] = g;
+ instr[instrLen].n[1] = 0;
+ instr[instrLen++].q0 = (uint8_t)q0;
+ }
+ else if (sscanf(lines[i], "cx q[%i], q[%i]", &q0, &q1) == 2)
+ {
+ if (qubitCount == 0xFF)
+ {
+ fprintf(stderr, "QAnsel: Instruction before initialization.\n");
+ errFound = 1;
+ break;
+ }
+ if (q0 > qubitCount || q1 > qubitCount | q0 < 0 || q1 < 0 || q0 == q1)
+ {
+ fprintf(stderr, "QAnsel: Invalid index");
+ errFound = 1;
+ break;
+ }
+ instr = realloc(instr, (instrLen + 1) * sizeof(QInstr));
+ strcpy(instr[instrLen].n, "cx");
+ instr[instrLen].q0 = q0;
+ instr[instrLen++].q1 = q1;
+ }
+ else if (sscanf(lines[i], "swap q[%i], q[%i]", &q0, &q1) == 2)
+ {
+ if (qubitCount == 0xFF)
+ {
+ fprintf(stderr, "QAnsel: Instruction before initialization.\n");
+ errFound = 1;
+ break;
+ }
+ if (q0 > qubitCount || q1 > qubitCount | q0 < 0 || q1 < 0 || q0 == q1)
+ {
+ fprintf(stderr, "QAnsel: Invalid index");
+ errFound = 1;
+ break;
+ }
+ instr = realloc(instr, (instrLen + 1) * sizeof(QInstr));
+ strcpy(instr[instrLen].n, "swap");
+ instr[instrLen].q0 = q0;
+ instr[instrLen++].q1 = q1;
+ }
+ else if
+ (
+ sscanf(lines[i], "cswap q[%i], q[%i], q[%i]", &q0, &q1, &q2) == 3
+ || sscanf(lines[i], "fredkin q[%i], q[%i], q[%i]", &q0, &q1, &q2) == 3
+ )
+ {
+ if (qubitCount == 0xFF)
+ {
+ fprintf(stderr, "QAnsel: Instruction before initialization.\n");
+ errFound = 1;
+ break;
+ }
+ if (qubitCount < 3)
+ {
+ fprintf(stderr, "QAnsel: Three qubit gate used with insufficient qubits initialized.\n");
+ errFound = 1;
+ break;
+ }
+ if (q0 >= qubitCount || q1 >= qubitCount || q2 >= qubitCount || q0 < 0 || q1< 0 || q2 < 0 || q0 == q1 || q1 == q2 || q0 == q2)
+ {
+ fprintf(stderr, "QAnsel: Invalid index");
+ errFound = 1;
+ break;
+ }
+ instr = realloc(instr, (instrLen + 1) * sizeof(QInstr));
+ strcpy(instr[instrLen].n, "cswap");
+ instr[instrLen].q0 = q0;
+ instr[instrLen].q1 = q1;
+ instr[instrLen++].q2 = q2;
+ }
+ else if
+ (
+ sscanf(lines[i], "ccx q[%i], q[%i], q[%i]", &q0, &q1, &q2) == 3
+ || sscanf(lines[i], "toffoli q[%i], q[%i], q[%i]", &q0, &q1, &q2) == 3
+ )
+ {
+ if (qubitCount == 0xFF)
+ {
+ fprintf(stderr, "QAnsel: Instruction before initialization.\n");
+ errFound = 1;
+ break;
+ }
+ if (qubitCount < 3)
+ {
+ fprintf(stderr, "QAnsel: Three qubit gate used with insufficient qubits initialized.\n");
+ errFound = 1;
+ break;
+ }
+ if (q0 >= qubitCount || q1 >= qubitCount || q2 >= qubitCount || q0 < 0 || q1< 0 || q2 < 0 || q0 == q1 || q1 == q2 || q0 == q2)
+ {
+ fprintf(stderr, "QAnsel: Invalid index");
+ errFound = 1;
+ break;
+ }
+ instr = realloc(instr, (instrLen + 1) * sizeof(QInstr));
+ strcpy(instr[instrLen].n, "ccx");
+ instr[instrLen].q0 = q0;
+ instr[instrLen].q1 = q1;
+ instr[instrLen++].q2 = q2;
+ }
+ else if (sscanf(lines[i], "measure q[%i] -> c[%i]", &q0, &q1) == 2)
+ {
+ if (bitCount == 0xFF)
+ {
+ fprintf(stderr, "QAnsel: Measure instruction used before bit initialization.\n");
+ errFound = 1;
+ break;
+ }
+ if (q0 >= qubitCount || q1 >= bitCount || q0 < 0 || q1 < 0)
+ {
+ fprintf(stderr, "QAnsel: Invalid index");
+ errFound = 1;
+ break;
+ }
+ instr = realloc(instr, (instrLen + 1) * sizeof(QInstr));
+ strcpy(instr[instrLen].n, "measure");
+ instr[instrLen].q0 = q0;
+ instr[instrLen++].q1 = q1;
+ }
+ else if (sscanf(lines[i], "if(c==%i)", &q0) == 1)
+ {
+ if (bitCount == 0xFF)
+ {
+ fprintf(stderr, "QAnsel: If instruction used before bit initialization.\n");
+ errFound = 1;
+ break;
+ }
+ if (q0 < 0)
+ {
+ fprintf(stderr, "QAnsel: Invalid comparison");
+ errFound = 1;
+ break;
+ }
+ instr = realloc(instr, (instrLen + 1) * sizeof(QInstr));
+ strcpy(instr[instrLen].n, "if_all");
+ instr[instrLen++].q0 = q0;
+ }
+ else if (sscanf(lines[i], "if(c[%i]==%i)", &q0, &q1) == 2)
+ {
+ if (bitCount == 0xFF)
+ {
+ fprintf(stderr, "QAnsel: If instruction used before bit initialization.\n");
+ errFound = 1;
+ break;
+ }
+ if (q0 < 0 || q1 < 0 || q0 > bitCount)
+ {
+ fprintf(stderr, "QAnsel: Invalid comparison");
+ errFound = 1;
+ break;
+ }
+ instr = realloc(instr, (instrLen + 1) * sizeof(QInstr));
+ strcpy(instr[instrLen].n, "if");
+ instr[instrLen].q0 = q0;
+ instr[instrLen++].q1 = q1;
+ }
+ else if (strcmp(lines[i], "print q") == 0)
+ {
+ if (qubitCount == 0xFF)
+ {
+ fprintf(stderr, "QAnsel: Qubit instruction used before initialization.\n");
+ errFound = 1;
+ break;
+ }
+ instr = realloc(instr, (instrLen + 1) * sizeof(QInstr));
+ strcpy(instr[instrLen].n, "printq_all");
+ instrLen++;
+ }
+ else if (strcmp(lines[i], "print c") == 0)
+ {
+ if (bitCount == 0xFF)
+ {
+ fprintf(stderr, "QAnsel: Bit instruction used before initialization.\n");
+ errFound = 1;
+ break;
+ }
+ instr = realloc(instr, (instrLen + 1) * sizeof(QInstr));
+ strcpy(instr[instrLen].n, "printc_all");
+ instrLen++;
+ }
+ else if (sscanf(lines[i], "print q[%i]", &q0) == 1)
+ {
+ if (qubitCount == 0xFF)
+ {
+ fprintf(stderr, "QAnsel: Qubit instruction used before initialization.\n");
+ errFound = 1;
+ break;
+ }
+ if (q0 >= qubitCount || q0 < 0)
+ {
+ fprintf(stderr, "QAnsel: Invalid index");
+ errFound = 1;
+ break;
+ }
+ instr = realloc(instr, (instrLen + 1) * sizeof(QInstr));
+ strcpy(instr[instrLen].n, "printq");
+ instr[instrLen++].q0 = q0;
+ }
+ else if (sscanf(lines[i], "print c[%i]", &q0) == 1)
+ {
+ if (bitCount == 0xFF)
+ {
+ fprintf(stderr, "QAnsel: Bit instruction used before initialization.\n");
+ errFound = 1;
+ break;
+ }
+ if (q0 >= bitCount || q0 < 0)
+ {
+ fprintf(stderr, "QAnsel: Invalid index");
+ errFound = 1;
+ break;
+ }
+ instr = realloc(instr, (instrLen + 1) * sizeof(QInstr));
+ strcpy(instr[instrLen].n, "printc");
+ instr[instrLen++].q0 = q0;
+ }
+ else if (strcmp(lines[i], "sample") == 0)
+ {
+ if (qubitCount == 0xFF)
+ {
+ fprintf(stderr, "QAnsel: Qubit instruction used before initialization.\n");
+ errFound = 1;
+ break;
+ }
+ instr = realloc(instr, (instrLen + 1) * sizeof(QInstr));
+ strcpy(instr[instrLen].n, "sample_all");
+ instrLen++;
+ }
+ else if (sscanf(lines[i], "sample q[%i]", &q0) == 1)
+ {
+ if (qubitCount == 0xFF)
+ {
+ fprintf(stderr, "QAnsel: Qubit instruction used before initialization.\n");
+ errFound = 1;
+ break;
+ }
+ if (q0 >= qubitCount || q0 < 0)
+ {
+ fprintf(stderr, "QAnsel: Invalid index");
+ errFound = 1;
+ break;
+ }
+ instr = realloc(instr, (instrLen + 1) * sizeof(QInstr));
+ strcpy(instr[instrLen].n, "sample");
+ instr[instrLen++].q0 = q0;
+ }
+ else if (strcmp(lines[i], "reset") == 0)
+ {
+ if (qubitCount == 0xFF || bitCount == 0xFF)
+ {
+ fprintf(stderr, "QAnsel: Instruction used before initialization.\n");
+ errFound = 1;
+ break;
+ }
+ instr = realloc(instr, (instrLen + 1) * sizeof(QInstr));
+ strcpy(instr[instrLen].n, "reset");
+ instrLen++;
+ }
+ else
+ {
+ fprintf(stderr, "QAnsel: Syntax error");
+ errFound = 1;
+ break;
+ }
+ }
+
+ for (uint32_t i = 0; i < linesLen; i++) free(lines[i]);
+ free(lines);
+
+ if (errFound)
+ {
+ printf(" on line %i.\n", lineIDs[lineID]);
+ free(lineIDs);
+ free(instr);
+ exit(1);
+ }
+
+ qansel_run(qubitCount, bitCount, instr, instrLen, 0);
+ free(instr);
+ free(lineIDs);
+}
--- /dev/null
+#ifndef __cpx__
+#define __cpx__
+#include <stdint.h>
+#include <stddef.h>
+#include <math.h>
+
+typedef struct
+{
+ double real, imaginary;
+} cpx_t;
+
+typedef struct
+{
+ double *ptr;
+ size_t rows, cols;
+} cpx_mtx_t;
+
+uint8_t* cpx_str(cpx_t* n)
+{
+ uint8_t* r;
+ size_t z;
+
+ double rl = n->real;
+ double ig = n->imaginary >= 0 ? n->imaginary : -(n->imaginary);
+ if (ig == 0)
+ {
+ z = snprintf(NULL, 0, "%f", rl);
+ r = malloc(z + 1);
+ sprintf(r, "%f", rl);
+
+ }
+ else
+ {
+ uint8_t op = n->imaginary >= 0 ? '+' : '-';
+ z = snprintf(NULL, 0, "%f %c %fi", rl, op, ig);
+ r = malloc(z + 1);
+ sprintf(r, "%f %c %fi", rl, op, ig);
+ }
+ return r;
+}
+
+cpx_t cpx_new(double r, double i)
+{
+ cpx_t n;
+ n.real = r;
+ n.imaginary = i;
+ return n;
+}
+
+void cpx_add(cpx_t* r, cpx_t* a, cpx_t* b)
+{
+ r->real = a->real + b->real;
+ r->imaginary = a->imaginary + b->imaginary;
+}
+
+void cpx_sub(cpx_t* r, cpx_t* a, cpx_t* b)
+{
+ r->real = a->real - b->real;
+ r->imaginary = a->imaginary - b->imaginary;
+}
+
+void cpx_mul(cpx_t* r, cpx_t* a, cpx_t* b)
+{
+ //FOIL
+ double first = a->real * b->real; //real
+ double outer = a->real * b->imaginary; //imaginary
+ double inner = a->imaginary * b->real; //imaginary
+ double last = -(a->imaginary * b->imaginary); //real
+ r->real = first + last;
+ r->imaginary = outer + inner;
+}
+
+//non-complex matrix multiply
+// shared = colsA = rowsB
+void cpx_ncpx_mmul(double* ptrR, double* ptrA, double* ptrB, size_t rowsA, size_t colsB, size_t shared)
+{
+ size_t colsA = shared;
+ size_t rowsB = shared;
+ size_t rowsR = rowsA;
+ size_t colsR = colsB;
+ for (size_t rowR = 0; rowR < rowsR; rowR++)
+ {
+ for (size_t colR = 0; colR < colsR; colR++)
+ {
+ size_t posR = colR + rowR * colsR;
+ size_t rowA = rowR;
+ size_t colB = colR;
+ ptrR[posR] = 0;
+ for (size_t i = 0; i < shared; i++)
+ {
+ size_t posA = i + rowA * colsA;
+ size_t posB = colB + i * colsB;
+ ptrR[posR] += ptrA[posA] * ptrB[posB];
+ }
+ }
+ }
+}
+
+//non-complex kronecker product
+void cpx_ncpx_mknk(double* ptrR, double* ptrA, double* ptrB, size_t rowsA, size_t colsA, size_t rowsB, size_t colsB)
+{
+ size_t rowsR = rowsA * rowsB;
+ size_t colsR = colsA * colsB;
+ for (size_t rowR = 0; rowR < rowsR; rowR++)
+ {
+ for (size_t colR = 0; colR < colsR; colR++)
+ {
+ size_t rowA = rowR / rowsB;
+ size_t colA = colR / colsB;
+ size_t rowB = rowR % rowsB;
+ size_t colB = colR % colsB;
+ printf(">%i,%i|%i<\n", colR, rowR, colR + rowR * colsR);
+ ptrR[colR + rowR * colsR] =
+ ptrA[colA + rowA * colsA]
+ * ptrB[colB + rowB * colsB]
+ ;
+ }
+ }
+}
+
+double cpx_magsqr(cpx_t* n)
+{
+ return (n->real * n->real) + (n->imaginary * n->imaginary);
+}
+
+double cpx_mag(cpx_t* n)
+{
+ return sqrt((n->real * n->real) + (n->imaginary * n->imaginary));
+}
+
+void cpx_mtx_mul(cpx_mtx_t* r, cpx_mtx_t* a, cpx_mtx_t* b)
+{
+ r->rows = a->rows;
+ r->cols = b->cols;
+ cpx_ncpx_mmul(r->ptr, a->ptr, b->ptr, a->rows * 2, b->cols * 2, a->cols * 2);
+}
+
+void cpx_mtx_set(cpx_mtx_t* m, size_t row, size_t col, cpx_t* n)
+{
+ row *= 2;
+ col *= 2;
+ size_t cols = m->cols * 2;
+ m->ptr[col + row * cols] = n->real;
+ m->ptr[(col + 1) + row * cols] = -(n->imaginary);
+ m->ptr[col + (row + 1) * cols] = n->imaginary;
+ m->ptr[(col + 1) + (row + 1) * cols] = n->real;
+}
+
+void cpx_mtx_set2(cpx_mtx_t* m, size_t row, size_t col, double real, double imaginary)
+{
+ row *= 2;
+ col *= 2;
+ size_t cols = m->cols * 2;
+ m->ptr[col + row * cols] = real;
+ m->ptr[(col + 1) + row * cols] = -imaginary;
+ m->ptr[col + (row + 1) * cols] = imaginary;
+ m->ptr[(col + 1) + (row + 1) * cols] = real;
+}
+
+void cpx_mtx_get(cpx_mtx_t* m, size_t row, size_t col, cpx_t* n)
+{
+ row *= 2;
+ col *= 2;
+ size_t cols = m->cols * 2;
+ n->real = m->ptr[(col + 1) + (row + 1) * cols];
+ n->imaginary = m->ptr[col + (row + 1) * cols];
+}
+
+double cpx_mtx_get_real(cpx_mtx_t* m, size_t row, size_t col)
+{
+ row *= 2;
+ col *= 2;
+ size_t cols = m->cols * 2;
+ return m->ptr[(col + 1) + (row + 1) * cols];
+}
+
+double cpx_mtx_get_imaginary(cpx_mtx_t* m, size_t row, size_t col)
+{
+ row *= 2;
+ col *= 2;
+ size_t cols = m->cols * 2;
+ return m->ptr[col + (row + 1) * cols];
+}
+
+void cpx_mtx_init(cpx_mtx_t* m, size_t rows, size_t cols)
+{
+ m->rows = rows;
+ m->cols = cols;
+ size_t z = (rows * 2) * (cols * 2);
+ m->ptr = malloc(z * sizeof(double));
+ for (size_t i = 0; i < z; i++) m->ptr[i] = 0;
+}
+
+void cpx_mtx_expand_row(cpx_mtx_t* m)
+{
+ m->rows += 1;
+ size_t z = (m->rows * 2) * (m->cols * 2);
+ m->ptr = realloc(m->ptr, z * sizeof(double));
+}
+
+
+void cpx_mtx_free(cpx_mtx_t* m)
+{
+ if (m->ptr != NULL) free(m->ptr);
+ m->rows = 0;
+ m->cols = 0;
+}
+
+void cpx_mtx_knk2
+(
+ double* ptrR,
+ size_t rowsR,
+ size_t colsR,
+ double* ptrA,
+ size_t rowsA,
+ size_t colsA,
+ double* ptrB,
+ size_t rowsB,
+ size_t colsB
+)
+{
+ for (size_t rowR = 0; rowR < rowsR; rowR++)
+ {
+ for (size_t colR = 0; colR < colsR; colR++)
+ {
+ size_t rowA = rowR / rowsB;
+ size_t colA = colR / colsB;
+ size_t rowB = rowR % rowsB;
+ size_t colB = colR % colsB;
+
+ double r1 = ptrA[((colA * 2) + 1) + ((rowA * 2) + 1) * (colsA * 2)];
+ double i1 = ptrA[(colA * 2) + ((rowA * 2) + 1) * (colsA * 2)];
+ double r2 = ptrB[((colB * 2) + 1) + ((rowB * 2) + 1) * (colsB * 2)];
+ double i2 = ptrB[(colB * 2) + ((rowB * 2) + 1) * (colsB * 2)];
+
+ double first = r1 * r2; //real
+ double outer = r1 * i2; //imaginary
+ double inner = i1 * r2; //imaginary
+ double last = -(i1 * i2); //real
+ r1 = first + last;
+ i1 = outer + inner;
+
+ ptrR[(colR * 2) + (rowR * 2) * (colsR * 2)] = r1;
+ ptrR[((colR * 2) + 1) + (rowR * 2) * (colsR * 2)] = -i1;
+ ptrR[(colR * 2) + ((rowR * 2) + 1) * (colsR * 2)] = i1;
+ ptrR[((colR * 2) + 1) + ((rowR * 2) + 1) * (colsR * 2)] = r1;
+
+ }
+ }
+}
+
+void cpx_mtx_knk(cpx_mtx_t* r, cpx_mtx_t* a, cpx_mtx_t* b)
+{
+ size_t rowsA = a->rows;
+ size_t colsA = a->cols;
+ size_t rowsB = b->rows;
+ size_t colsB = b->cols;
+ size_t rowsR = rowsA * rowsB;
+ size_t colsR = colsA * colsB;
+ cpx_mtx_init(r, rowsR, colsR);
+ for (size_t rowR = 0; rowR < rowsR; rowR++)
+ {
+ for (size_t colR = 0; colR < colsR; colR++)
+ {
+ size_t rowA = rowR / rowsB;
+ size_t colA = colR / colsB;
+ size_t rowB = rowR % rowsB;
+ size_t colB = colR % colsB;
+ cpx_t n1, n2;
+
+ cpx_mtx_get(a, rowA, colA, &n1);
+ cpx_mtx_get(b, rowB, colB, &n2);
+ cpx_mul(&n1, &n1, &n2);
+
+ cpx_mtx_set(r, rowR, colR, &n1);
+ }
+ }
+}
+
+
+void cpx_mtx_(cpx_mtx_t* m)
+{
+ for (size_t r = 0; r < m->rows * 2; r++)
+ {
+ for (size_t c = 0; c < m->cols * 2; c++)
+ {
+ if (c > 0) printf(", ");
+ printf("%f", m->ptr[c + r * m->cols * 2]);
+ }
+ printf("\n");
+ }
+}
+
+void cpx_mtx_print(cpx_mtx_t* m)
+{
+ for (size_t r = 0; r < m->rows; r++)
+ {
+ if (r > 0) printf("\n");
+ for (size_t c = 0; c < m->cols; c++)
+ {
+ cpx_t n;
+ cpx_mtx_get(m, r, c, &n);
+ uint8_t* s = cpx_str(&n);
+ if (c > 0) printf(", ");
+ printf("%s", s);
+ free(s);
+ }
+ }
+}
+
+#endif
Foleosoft / QAnsel.git / commitdiff