Anders and Briegel in Python
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

971 lines
21KB

  1. // CHP: CNOT-Hadamard-Phase
  2. // Stabilizer Quantum Computer Simulator
  3. // by Scott Aaronson
  4. // Last modified June 30, 2004
  5. // Thanks to Simon Anders and Andrew Cross for bugfixes
  6. #include <stdio.h>
  7. #include <stdlib.h>
  8. #include <string.h>
  9. #include <Python.h>
  10. #include <time.h>
  11. #define CNOT 0
  12. #define HADAMARD 1
  13. #define PHASE 2
  14. #define MEASURE 3
  15. struct QProg
  16. // Quantum circuit
  17. {
  18. long n; // # of qubits
  19. long T; // # of gates
  20. char *a; // Instruction opcode
  21. long *b; // Qubit 1
  22. long *c; // Qubit 2 (target for CNOT)
  23. int DISPQSTATE; // whether to print the state (q for final state only, Q for every iteration)
  24. int DISPTIME; // whether to print the execution time
  25. int SILENT; // whether NOT to print measurement results
  26. int DISPPROG; // whether to print instructions being executed as they're executed
  27. int SUPPRESSM; // whether to suppress actual computation of determinate measurement results
  28. };
  29. struct QState
  30. // Quantum state
  31. {
  32. // To save memory and increase speed, the bits are packed 32 to an unsigned long
  33. long n; // # of qubits
  34. unsigned long **x; // (2n+1)*n matrix for stabilizer/destabilizer x bits (there's one "scratch row" at
  35. unsigned long **z; // (2n+1)*n matrix for z bits the bottom)
  36. int *r; // Phase bits: 0 for +1, 1 for i, 2 for -1, 3 for -i. Normally either 0 or 2.
  37. unsigned long pw[32]; // pw[i] = 2^i
  38. long over32; // floor(n/8)+1
  39. };
  40. // Globals
  41. struct QProg *h;
  42. struct QState *q;
  43. void error(int k)
  44. {
  45. if (k==0) printf("\nSyntax: chp [-options] <filename> [input]\n");
  46. if (k==1) printf("\nFile not found\n");
  47. exit(0);
  48. }
  49. void cnot(struct QState *q, long b, long c)
  50. // Apply a CNOT gate with control b and target c
  51. {
  52. long i;
  53. long b5;
  54. long c5;
  55. unsigned long pwb;
  56. unsigned long pwc;
  57. b5 = b>>5;
  58. c5 = c>>5;
  59. pwb = q->pw[b&31];
  60. pwc = q->pw[c&31];
  61. for (i = 0; i < 2*q->n; i++)
  62. {
  63. if (q->x[i][b5]&pwb) q->x[i][c5] ^= pwc;
  64. if (q->z[i][c5]&pwc) q->z[i][b5] ^= pwb;
  65. if ((q->x[i][b5]&pwb) && (q->z[i][c5]&pwc) &&
  66. (q->x[i][c5]&pwc) && (q->z[i][b5]&pwb))
  67. q->r[i] = (q->r[i]+2)%4;
  68. if ((q->x[i][b5]&pwb) && (q->z[i][c5]&pwc) &&
  69. !(q->x[i][c5]&pwc) && !(q->z[i][b5]&pwb))
  70. q->r[i] = (q->r[i]+2)%4;
  71. }
  72. return;
  73. }
  74. void hadamard(struct QState *q, long b)
  75. // Apply a Hadamard gate to qubit b
  76. {
  77. long i;
  78. unsigned long tmp;
  79. long b5;
  80. unsigned long pw;
  81. b5 = b>>5;
  82. pw = q->pw[b&31];
  83. for (i = 0; i < 2*q->n; i++)
  84. {
  85. tmp = q->x[i][b5];
  86. q->x[i][b5] ^= (q->x[i][b5] ^ q->z[i][b5]) & pw;
  87. q->z[i][b5] ^= (q->z[i][b5] ^ tmp) & pw;
  88. if ((q->x[i][b5]&pw) && (q->z[i][b5]&pw)) q->r[i] = (q->r[i]+2)%4;
  89. }
  90. return;
  91. }
  92. void phase(struct QState *q, long b)
  93. // Apply a phase gate (|0>->|0>, |1>->i|1>) to qubit b
  94. {
  95. long i;
  96. long b5;
  97. unsigned long pw;
  98. b5 = b>>5;
  99. pw = q->pw[b&31];
  100. for (i = 0; i < 2*q->n; i++)
  101. {
  102. if ((q->x[i][b5]&pw) && (q->z[i][b5]&pw)) q->r[i] = (q->r[i]+2)%4;
  103. q->z[i][b5] ^= q->x[i][b5]&pw;
  104. }
  105. return;
  106. }
  107. void rowcopy(struct QState *q, long i, long k)
  108. // Sets row i equal to row k
  109. {
  110. long j;
  111. for (j = 0; j < q->over32; j++)
  112. {
  113. q->x[i][j] = q->x[k][j];
  114. q->z[i][j] = q->z[k][j];
  115. }
  116. q->r[i] = q->r[k];
  117. return;
  118. }
  119. void rowswap(struct QState *q, long i, long k)
  120. // Swaps row i and row k
  121. {
  122. rowcopy(q, 2*q->n, k);
  123. rowcopy(q, k, i);
  124. rowcopy(q, i, 2*q->n);
  125. return;
  126. }
  127. void rowset(struct QState *q, long i, long b)
  128. // Sets row i equal to the bth observable (X_1,...X_n,Z_1,...,Z_n)
  129. {
  130. long j;
  131. long b5;
  132. unsigned long b31;
  133. for (j = 0; j < q->over32; j++)
  134. {
  135. q->x[i][j] = 0;
  136. q->z[i][j] = 0;
  137. }
  138. q->r[i] = 0;
  139. if (b < q->n)
  140. {
  141. b5 = b>>5;
  142. b31 = b&31;
  143. q->x[i][b5] = q->pw[b31];
  144. }
  145. else
  146. {
  147. b5 = (b - q->n)>>5;
  148. b31 = (b - q->n)&31;
  149. q->z[i][b5] = q->pw[b31];
  150. }
  151. return;
  152. }
  153. int clifford(struct QState *q, long i, long k)
  154. // Return the phase (0,1,2,3) when row i is LEFT-multiplied by row k
  155. {
  156. long j;
  157. long l;
  158. unsigned long pw;
  159. long e=0; // Power to which i is raised
  160. for (j = 0; j < q->over32; j++)
  161. for (l = 0; l < 32; l++)
  162. {
  163. pw = q->pw[l];
  164. if ((q->x[k][j]&pw) && (!(q->z[k][j]&pw))) // X
  165. {
  166. if ((q->x[i][j]&pw) && (q->z[i][j]&pw)) e++; // XY=iZ
  167. if ((!(q->x[i][j]&pw)) && (q->z[i][j]&pw)) e--; // XZ=-iY
  168. }
  169. if ((q->x[k][j]&pw) && (q->z[k][j]&pw)) // Y
  170. {
  171. if ((!(q->x[i][j]&pw)) && (q->z[i][j]&pw)) e++; // YZ=iX
  172. if ((q->x[i][j]&pw) && (!(q->z[i][j]&pw))) e--; // YX=-iZ
  173. }
  174. if ((!(q->x[k][j]&pw)) && (q->z[k][j]&pw)) // Z
  175. {
  176. if ((q->x[i][j]&pw) && (!(q->z[i][j]&pw))) e++; // ZX=iY
  177. if ((q->x[i][j]&pw) && (q->z[i][j]&pw)) e--; // ZY=-iX
  178. }
  179. }
  180. e = (e+q->r[i]+q->r[k])%4;
  181. if (e>=0) return e;
  182. else return e+4;
  183. }
  184. void rowmult(struct QState *q, long i, long k)
  185. // Left-multiply row i by row k
  186. {
  187. long j;
  188. q->r[i] = clifford(q,i,k);
  189. for (j = 0; j < q->over32; j++)
  190. {
  191. q->x[i][j] ^= q->x[k][j];
  192. q->z[i][j] ^= q->z[k][j];
  193. }
  194. return;
  195. }
  196. void printstate(struct QState *q)
  197. // Print the destabilizer and stabilizer for state q
  198. {
  199. long i;
  200. long j;
  201. long j5;
  202. unsigned long pw;
  203. for (i = 0; i < 2*q->n; i++)
  204. {
  205. if (i == q->n)
  206. {
  207. printf("\n");
  208. for (j = 0; j < q->n+1; j++)
  209. printf("-");
  210. }
  211. if (q->r[i]==2) printf("\n-");
  212. else printf("\n+");
  213. for (j = 0; j < q->n; j++)
  214. {
  215. j5 = j>>5;
  216. pw = q->pw[j&31];
  217. if ((!(q->x[i][j5]&pw)) && (!(q->z[i][j5]&pw))) printf("I");
  218. if ((q->x[i][j5]&pw) && (!(q->z[i][j5]&pw))) printf("X");
  219. if ((q->x[i][j5]&pw) && (q->z[i][j5]&pw)) printf("Y");
  220. if ((!(q->x[i][j5]&pw)) && (q->z[i][j5]&pw)) printf("Z");
  221. }
  222. }
  223. printf("\n");
  224. return;
  225. }
  226. int measure(struct QState *q, long b, int sup)
  227. // Measure qubit b
  228. // Return 0 if outcome would always be 0
  229. // 1 if outcome would always be 1
  230. // 2 if outcome was random and 0 was chosen
  231. // 3 if outcome was random and 1 was chosen
  232. // sup: 1 if determinate measurement results should be suppressed, 0 otherwise
  233. {
  234. int ran = 0;
  235. long i;
  236. long p; // pivot row in stabilizer
  237. long m; // pivot row in destabilizer
  238. long b5;
  239. unsigned long pw;
  240. b5 = b>>5;
  241. pw = q->pw[b&31];
  242. for (p = 0; p < q->n; p++) // loop over stabilizer generators
  243. {
  244. if (q->x[p+q->n][b5]&pw) ran = 1; // if a Zbar does NOT commute with Z_b (the
  245. if (ran) break; // operator being measured), then outcome is random
  246. }
  247. // If outcome is indeterminate
  248. if (ran)
  249. {
  250. rowcopy(q, p, p + q->n); // Set Xbar_p := Zbar_p
  251. rowset(q, p + q->n, b + q->n); // Set Zbar_p := Z_b
  252. q->r[p + q->n] = 2*(rand()%2); // moment of quantum randomness
  253. for (i = 0; i < 2*q->n; i++) // Now update the Xbar's and Zbar's that don't commute with
  254. if ((i!=p) && (q->x[i][b5]&pw)) // Z_b
  255. rowmult(q, i, p);
  256. if (q->r[p + q->n]) return 3;
  257. else return 2;
  258. }
  259. // If outcome is determinate
  260. if ((!ran) && (!sup))
  261. {
  262. for (m = 0; m < q->n; m++) // Before we were checking if stabilizer generators commute
  263. if (q->x[m][b5]&pw) break; // with Z_b; now we're checking destabilizer generators
  264. rowcopy(q, 2*q->n, m + q->n);
  265. for (i = m+1; i < q->n; i++)
  266. if (q->x[i][b5]&pw)
  267. rowmult(q, 2*q->n, i + q->n);
  268. if (q->r[2*q->n]) return 1;
  269. else return 0;
  270. /*for (i = m+1; i < q->n; i++)
  271. if (q->x[i][b5]&pw)
  272. {
  273. rowmult(q, m + q->n, i + q->n);
  274. rowmult(q, i, m);
  275. }
  276. return (int)q->r[m + q->n];*/
  277. }
  278. return 0;
  279. }
  280. long gaussian(struct QState *q)
  281. // Do Gaussian elimination to put the stabilizer generators in the following form:
  282. // At the top, a minimal set of generators containing X's and Y's, in "quasi-upper-triangular" form.
  283. // (Return value = number of such generators = log_2 of number of nonzero basis states)
  284. // At the bottom, generators containing Z's only in quasi-upper-triangular form.
  285. {
  286. long i = q->n;
  287. long k;
  288. long k2;
  289. long j;
  290. long j5;
  291. long g; // Return value
  292. unsigned long pw;
  293. for (j = 0; j < q->n; j++)
  294. {
  295. j5 = j>>5;
  296. pw = q->pw[j&31];
  297. for (k = i; k < 2*q->n; k++) // Find a generator containing X in jth column
  298. if (q->x[k][j5]&pw) break;
  299. if (k < 2*q->n)
  300. {
  301. rowswap(q, i, k);
  302. rowswap(q, i-q->n, k-q->n);
  303. for (k2 = i+1; k2 < 2*q->n; k2++)
  304. if (q->x[k2][j5]&pw)
  305. {
  306. rowmult(q, k2, i); // Gaussian elimination step
  307. rowmult(q, i-q->n, k2-q->n);
  308. }
  309. i++;
  310. }
  311. }
  312. g = i - q->n;
  313. for (j = 0; j < q->n; j++)
  314. {
  315. j5 = j>>5;
  316. pw = q->pw[j&31];
  317. for (k = i; k < 2*q->n; k++) // Find a generator containing Z in jth column
  318. if (q->z[k][j5]&pw) break;
  319. if (k < 2*q->n)
  320. {
  321. rowswap(q, i, k);
  322. rowswap(q, i-q->n, k-q->n);
  323. for (k2 = i+1; k2 < 2*q->n; k2++)
  324. if (q->z[k2][j5]&pw)
  325. {
  326. rowmult(q, k2, i);
  327. rowmult(q, i-q->n, k2-q->n);
  328. }
  329. i++;
  330. }
  331. }
  332. return g;
  333. }
  334. long innerprod(struct QState *q1, struct QState *q2)
  335. // Returns -1 if q1 and q2 are orthogonal
  336. // Otherwise, returns a nonnegative integer s such that the inner product is (1/sqrt(2))^s
  337. {
  338. return 0;
  339. }
  340. void printbasisstate(struct QState *q)
  341. // Prints the result of applying the Pauli operator in the "scratch space" of q to |0...0>
  342. {
  343. long j;
  344. long j5;
  345. unsigned long pw;
  346. int e = q->r[2*q->n];
  347. for (j = 0; j < q->n; j++)
  348. {
  349. j5 = j>>5;
  350. pw = q->pw[j&31];
  351. if ((q->x[2*q->n][j5]&pw) && (q->z[2*q->n][j5]&pw)) // Pauli operator is "Y"
  352. e = (e+1)%4;
  353. }
  354. if (e==0) printf("\n +|");
  355. if (e==1) printf("\n+i|");
  356. if (e==2) printf("\n -|");
  357. if (e==3) printf("\n-i|");
  358. for (j = 0; j < q->n; j++)
  359. {
  360. j5 = j>>5;
  361. pw = q->pw[j&31];
  362. if (q->x[2*q->n][j5]&pw) printf("1");
  363. else printf("0");
  364. }
  365. printf(">");
  366. return;
  367. }
  368. void dumpbasisstate(struct QState *q, PyObject *output, int index)
  369. // Dumps the result of applying the Pauli operator in the "scratch space" of q to |0...0>
  370. {
  371. long j;
  372. long j5;
  373. unsigned long pw;
  374. int e = q->r[2*q->n];
  375. char buffer[1024];
  376. strcpy(buffer, "");
  377. for (j = 0; j < q->n; j++)
  378. {
  379. j5 = j>>5;
  380. pw = q->pw[j&31];
  381. if ((q->x[2*q->n][j5]&pw) && (q->z[2*q->n][j5]&pw)) // Pauli operator is "Y"
  382. e = (e+1)%4;
  383. }
  384. if (e==0) strcat(buffer, "+|");
  385. if (e==1) strcat(buffer, "+i|");
  386. if (e==2) strcat(buffer, "-|");
  387. if (e==3) strcat(buffer, "-i|");
  388. int key;
  389. key = 0;
  390. for (j = 0; j < q->n; j++)
  391. {
  392. j5 = j>>5;
  393. pw = q->pw[j&31];
  394. if (q->x[2*q->n][j5]&pw) key = key ^ (1 << j);
  395. }
  396. PyDict_SetItem(output, Py_BuildValue("i", key), Py_BuildValue("i", e));
  397. return;
  398. }
  399. void seed(struct QState *q, long g)
  400. // Finds a Pauli operator P such that the basis state P|0...0> occurs with nonzero amplitude in q, and
  401. // writes P to the scratch space of q. For this to work, Gaussian elimination must already have been
  402. // performed on q. g is the return value from gaussian(q).
  403. {
  404. long i;
  405. long j;
  406. long j5;
  407. unsigned long pw;
  408. int f;
  409. long min=0;
  410. q->r[2*q->n] = 0;
  411. for (j = 0; j < q->over32; j++)
  412. {
  413. q->x[2*q->n][j] = 0; // Wipe the scratch space clean
  414. q->z[2*q->n][j] = 0;
  415. }
  416. for (i = 2*q->n - 1; i >= q->n + g; i--)
  417. {
  418. f = q->r[i];
  419. for (j = q->n - 1; j >= 0; j--)
  420. {
  421. j5 = j>>5;
  422. pw = q->pw[j&31];
  423. if (q->z[i][j5]&pw)
  424. {
  425. min = j;
  426. if (q->x[2*q->n][j5]&pw) f = (f+2)%4;
  427. }
  428. }
  429. if (f==2)
  430. {
  431. j5 = min>>5;
  432. pw = q->pw[min&31];
  433. q->x[2*q->n][j5] ^= pw; // Make the seed consistent with the ith equation
  434. }
  435. }
  436. return;
  437. }
  438. static char act_hadamard_docs[] = "act_hadamard(): Do a Hadamard";
  439. static PyObject* act_hadamard(PyObject* self, PyObject* args)
  440. {
  441. int which;
  442. if (!PyArg_ParseTuple(args, "i", &which)) { return NULL; }
  443. hadamard(q, which);
  444. Py_INCREF(Py_None);
  445. return Py_None;
  446. }
  447. static char act_phase_docs[] = "act_phase(): Do a phase";
  448. static PyObject* act_phase(PyObject* self, PyObject* args)
  449. {
  450. int which;
  451. if (!PyArg_ParseTuple(args, "i", &which)) { return NULL; }
  452. phase(q, which);
  453. Py_INCREF(Py_None);
  454. return Py_None;
  455. }
  456. static char act_cnot_docs[] = "act_cnot(): Do a CNOT";
  457. static PyObject* act_cnot(PyObject* self, PyObject* args)
  458. {
  459. int a, b;
  460. if (!PyArg_ParseTuple(args, "ii", &a, &b)) { return NULL; }
  461. cnot(q, a, b);
  462. Py_INCREF(Py_None);
  463. return Py_None;
  464. }
  465. static char get_ket_docs[] = "get_ket(): Get the state vector";
  466. static PyObject* get_ket(PyObject* self, PyObject* noarg)
  467. // Print the state in ket notation (warning: could be huge!)
  468. {
  469. long g; // log_2 of number of nonzero basis states
  470. unsigned long t;
  471. unsigned long t2;
  472. long i;
  473. g = gaussian(q);
  474. PyObject *output = PyDict_New();
  475. if (g > 30) { return output; }
  476. seed(q, g);
  477. dumpbasisstate(q, output, 0);
  478. for (t = 0; t < q->pw[g]-1; t++)
  479. {
  480. t2 = t ^ (t+1);
  481. for (i = 0; i < g; i++)
  482. if (t2 & q->pw[i])
  483. rowmult(q, 2*q->n, q->n + i);
  484. dumpbasisstate(q, output, t+1);
  485. }
  486. return output;
  487. }
  488. void printket(struct QState *q)
  489. // Print the state in ket notation (warning: could be huge!)
  490. {
  491. long g; // log_2 of number of nonzero basis states
  492. unsigned long t;
  493. unsigned long t2;
  494. long i;
  495. g = gaussian(q);
  496. printf("\n2^%ld nonzero basis states", g);
  497. if (g > 31)
  498. {
  499. printf("\nState is WAY too big to print");
  500. return;
  501. }
  502. seed(q, g);
  503. printbasisstate(q);
  504. for (t = 0; t < q->pw[g]-1; t++)
  505. {
  506. t2 = t ^ (t+1);
  507. for (i = 0; i < g; i++)
  508. if (t2 & q->pw[i])
  509. rowmult(q, 2*q->n, q->n + i);
  510. printbasisstate(q);
  511. }
  512. printf("\n");
  513. return;
  514. }
  515. void runprog(struct QProg *h, struct QState *q)
  516. // Simulate the quantum circuit
  517. {
  518. long t;
  519. int m; // measurement result
  520. time_t tp;
  521. double dt;
  522. char mvirgin = 1;
  523. time(&tp);
  524. for (t = 0; t < h->T; t++)
  525. {
  526. if (h->a[t]==CNOT) cnot(q,h->b[t],h->c[t]);
  527. if (h->a[t]==HADAMARD) hadamard(q,h->b[t]);
  528. if (h->a[t]==PHASE) phase(q,h->b[t]);
  529. if (h->a[t]==MEASURE)
  530. {
  531. if (mvirgin && h->DISPTIME)
  532. {
  533. dt = difftime(time(0),tp);
  534. printf("\nGate time: %lf seconds", dt);
  535. printf("\nTime per 10000 gates: %lf seconds", dt*10000.0f/(h->T - h->n));
  536. time(&tp);
  537. }
  538. mvirgin = 0;
  539. m = measure(q,h->b[t],h->SUPPRESSM);
  540. if (!h->SILENT)
  541. {
  542. printf("\nOutcome of measuring qubit %ld: ", h->b[t]);
  543. if (m>1) printf("%d (random)", m-2);
  544. else printf("%d", m);
  545. }
  546. }
  547. if (h->DISPPROG)
  548. {
  549. if (h->a[t]==CNOT) printf("\nCNOT %ld->%ld", h->b[t], h->c[t]);
  550. if (h->a[t]==HADAMARD) printf("\nHadamard %ld", h->b[t]);
  551. if (h->a[t]==PHASE) printf("\nPhase %ld", h->b[t]);
  552. }
  553. }
  554. printf("\n");
  555. if (h->DISPTIME)
  556. {
  557. dt = difftime(time(0),tp);
  558. printf("\nMeasurement time: %lf seconds", dt);
  559. printf("\nTime per 10000 measurements: %lf seconds\n", dt*10000.0f/h->n);
  560. }
  561. if (h->DISPQSTATE)
  562. {
  563. printf("\nFinal state:");
  564. printstate(q);
  565. gaussian(q);
  566. printstate(q);
  567. printket(q);
  568. }
  569. return;
  570. }
  571. void preparestate(struct QState *q, char *s)
  572. // Prepare the initial state's "input"
  573. {
  574. long l;
  575. long b;
  576. l = strlen(s);
  577. for (b = 0; b < l; b++)
  578. {
  579. if (s[b]=='Z')
  580. {
  581. hadamard(q,b);
  582. phase(q,b);
  583. phase(q,b);
  584. hadamard(q,b);
  585. }
  586. if (s[b]=='x') hadamard(q,b);
  587. if (s[b]=='X')
  588. {
  589. hadamard(q,b);
  590. phase(q,b);
  591. phase(q,b);
  592. }
  593. if (s[b]=='y')
  594. {
  595. hadamard(q,b);
  596. phase(q,b);
  597. }
  598. if (s[b]=='Y')
  599. {
  600. hadamard(q,b);
  601. phase(q,b);
  602. phase(q,b);
  603. phase(q,b);
  604. }
  605. }
  606. return;
  607. }
  608. void initstae_(struct QState *q, long n, char *s)
  609. // Initialize state q to have n qubits, and input specified by s
  610. {
  611. long i;
  612. long j;
  613. q->n = n;
  614. q->x = malloc((2*q->n + 1) * sizeof(unsigned long*));
  615. q->z = malloc((2*q->n + 1) * sizeof(unsigned long*));
  616. q->r = malloc((2*q->n + 1) * sizeof(int));
  617. q->over32 = (q->n>>5) + 1;
  618. q->pw[0] = 1;
  619. for (i = 1; i < 32; i++)
  620. q->pw[i] = 2*q->pw[i-1];
  621. for (i = 0; i < 2*q->n + 1; i++)
  622. {
  623. q->x[i] = malloc(q->over32 * sizeof(unsigned long));
  624. q->z[i] = malloc(q->over32 * sizeof(unsigned long));
  625. for (j = 0; j < q->over32; j++)
  626. {
  627. q->x[i][j] = 0;
  628. q->z[i][j] = 0;
  629. }
  630. if (i < q->n)
  631. q->x[i][i>>5] = q->pw[i&31];
  632. else if (i < 2*q->n)
  633. {
  634. j = i-q->n;
  635. q->z[i][j>>5] = q->pw[j&31];
  636. }
  637. q->r[i] = 0;
  638. }
  639. if (s) preparestate(q, s);
  640. return;
  641. }
  642. void readprog(struct QProg *h, char *fn, char *params)
  643. // Read a quantum circuit from filename fn, with optional parameters params
  644. {
  645. long t;
  646. char fn2[255];
  647. FILE *fp;
  648. char c=0;
  649. long val;
  650. long l;
  651. h->DISPQSTATE = 0;
  652. h->DISPTIME = 0;
  653. h->SILENT = 0;
  654. h->DISPPROG = 0;
  655. h->SUPPRESSM = 0;
  656. if (params)
  657. {
  658. l = strlen(params);
  659. for (t = 1; t < l; t++)
  660. {
  661. if ((params[t]=='q')||(params[t]=='Q')) h->DISPQSTATE = 1;
  662. if ((params[t]=='p')||(params[t]=='P')) h->DISPPROG = 1;
  663. if ((params[t]=='t')||(params[t]=='T')) h->DISPTIME = 1;
  664. if ((params[t]=='s')||(params[t]=='S')) h->SILENT = 1;
  665. if ((params[t]=='m')||(params[t]=='M')) h->SUPPRESSM = 1;
  666. }
  667. }
  668. sprintf(fn2, "%s", fn);
  669. fp = fopen(fn2, "r");
  670. if (!fp)
  671. {
  672. sprintf(fn2, "%s.chp", fn);
  673. fp = fopen(fn2, "r");
  674. if (!fp) error(1);
  675. }
  676. while (!feof(fp)&&(c!='#'))
  677. fscanf(fp, "%c", &c);
  678. if (c!='#') error(2);
  679. h->T = 0;
  680. h->n = 0;
  681. while (!feof(fp))
  682. {
  683. fscanf(fp, "%c", &c);
  684. if ((c=='\r')||(c=='\n'))
  685. continue;
  686. fscanf(fp, "%ld", &val);
  687. if (val+1 > h->n) h->n = val+1;
  688. if ((c=='c')||(c=='C'))
  689. {
  690. fscanf(fp, "%ld", &val);
  691. if (val+1 > h->n) h->n = val+1;
  692. }
  693. h->T++;
  694. }
  695. fclose(fp);
  696. h->a = malloc(h->T * sizeof(char));
  697. h->b = malloc(h->T * sizeof(long));
  698. h->c = malloc(h->T * sizeof(long));
  699. fp = fopen(fn2, "r");
  700. while (!feof(fp)&&(c!='#'))
  701. fscanf(fp, "%c", &c);
  702. t=0;
  703. while (!feof(fp))
  704. {
  705. fscanf(fp, "%c", &c);
  706. if ((c=='\r')||(c=='\n'))
  707. continue;
  708. if ((c=='c')||(c=='C')) h->a[t] = CNOT;
  709. if ((c=='h')||(c=='H')) h->a[t] = HADAMARD;
  710. if ((c=='p')||(c=='P')) h->a[t] = PHASE;
  711. if ((c=='m')||(c=='M')) h->a[t] = MEASURE;
  712. fscanf(fp, "%ld", &h->b[t]);
  713. if (h->a[t]==CNOT) fscanf(fp, "%ld", &h->c[t]);
  714. t++;
  715. }
  716. fclose(fp);
  717. return;
  718. }
  719. static char init_docs[] = "init(nqubits): Initialize with n qubits";
  720. static PyObject* init(PyObject* self, PyObject* args)
  721. {
  722. int nqubits;
  723. if (!PyArg_ParseTuple(args, "i", &nqubits)) { return NULL; }
  724. srand(time(0));
  725. h = malloc(sizeof(struct QProg));
  726. q = malloc(sizeof(struct QState));
  727. initstae_(q,nqubits,NULL);
  728. PyObject *response = Py_BuildValue("i", nqubits);
  729. return response;
  730. }
  731. static PyMethodDef chp_funcs[] = {
  732. {"init", (PyCFunction)init, METH_VARARGS, init_docs},
  733. {"get_ket", (PyCFunction)get_ket, METH_NOARGS, get_ket_docs},
  734. {"act_hadamard", (PyCFunction)act_hadamard, METH_VARARGS, act_hadamard_docs},
  735. {"act_cnot", (PyCFunction)act_cnot, METH_VARARGS, act_cnot_docs},
  736. {"act_phase", (PyCFunction)act_phase, METH_VARARGS, act_phase_docs},
  737. {NULL}
  738. };
  739. void initchp(int argc, char **argv)
  740. {
  741. Py_InitModule3("chp", chp_funcs, "CHP");
  742. }