Anders and Briegel in Python
Vous ne pouvez pas sélectionner plus de 25 sujets Les noms de sujets doivent commencer par une lettre ou un nombre, peuvent contenir des tirets ('-') et peuvent comporter jusqu'à 35 caractères.

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986
  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 act_cz_docs[] = "act_cz(): Do a cz";
  466. static PyObject* act_cz(PyObject* self, PyObject* args)
  467. {
  468. int a, b;
  469. if (!PyArg_ParseTuple(args, "ii", &a, &b)) { return NULL; }
  470. hadamard(q, b);
  471. cnot(q, a, b);
  472. hadamard(q, b);
  473. Py_INCREF(Py_None);
  474. return Py_None;
  475. }
  476. static char get_ket_docs[] = "get_ket(): Get the state vector";
  477. static PyObject* get_ket(PyObject* self, PyObject* noarg)
  478. // Print the state in ket notation (warning: could be huge!)
  479. {
  480. long g; // log_2 of number of nonzero basis states
  481. unsigned long t;
  482. unsigned long t2;
  483. long i;
  484. g = gaussian(q);
  485. PyObject *output = PyDict_New();
  486. if (g > 30) { return output; }
  487. seed(q, g);
  488. dumpbasisstate(q, output, 0);
  489. for (t = 0; t < q->pw[g]-1; t++)
  490. {
  491. t2 = t ^ (t+1);
  492. for (i = 0; i < g; i++)
  493. if (t2 & q->pw[i])
  494. rowmult(q, 2*q->n, q->n + i);
  495. dumpbasisstate(q, output, t+1);
  496. }
  497. return output;
  498. }
  499. void printket(struct QState *q)
  500. // Print the state in ket notation (warning: could be huge!)
  501. {
  502. long g; // log_2 of number of nonzero basis states
  503. unsigned long t;
  504. unsigned long t2;
  505. long i;
  506. g = gaussian(q);
  507. printf("\n2^%ld nonzero basis states", g);
  508. if (g > 31)
  509. {
  510. printf("\nState is WAY too big to print");
  511. return;
  512. }
  513. seed(q, g);
  514. printbasisstate(q);
  515. for (t = 0; t < q->pw[g]-1; t++)
  516. {
  517. t2 = t ^ (t+1);
  518. for (i = 0; i < g; i++)
  519. if (t2 & q->pw[i])
  520. rowmult(q, 2*q->n, q->n + i);
  521. printbasisstate(q);
  522. }
  523. printf("\n");
  524. return;
  525. }
  526. void runprog(struct QProg *h, struct QState *q)
  527. // Simulate the quantum circuit
  528. {
  529. long t;
  530. int m; // measurement result
  531. time_t tp;
  532. double dt;
  533. char mvirgin = 1;
  534. time(&tp);
  535. for (t = 0; t < h->T; t++)
  536. {
  537. if (h->a[t]==CNOT) cnot(q,h->b[t],h->c[t]);
  538. if (h->a[t]==HADAMARD) hadamard(q,h->b[t]);
  539. if (h->a[t]==PHASE) phase(q,h->b[t]);
  540. if (h->a[t]==MEASURE)
  541. {
  542. if (mvirgin && h->DISPTIME)
  543. {
  544. dt = difftime(time(0),tp);
  545. printf("\nGate time: %lf seconds", dt);
  546. printf("\nTime per 10000 gates: %lf seconds", dt*10000.0f/(h->T - h->n));
  547. time(&tp);
  548. }
  549. mvirgin = 0;
  550. m = measure(q,h->b[t],h->SUPPRESSM);
  551. if (!h->SILENT)
  552. {
  553. printf("\nOutcome of measuring qubit %ld: ", h->b[t]);
  554. if (m>1) printf("%d (random)", m-2);
  555. else printf("%d", m);
  556. }
  557. }
  558. if (h->DISPPROG)
  559. {
  560. if (h->a[t]==CNOT) printf("\nCNOT %ld->%ld", h->b[t], h->c[t]);
  561. if (h->a[t]==HADAMARD) printf("\nHadamard %ld", h->b[t]);
  562. if (h->a[t]==PHASE) printf("\nPhase %ld", h->b[t]);
  563. }
  564. }
  565. printf("\n");
  566. if (h->DISPTIME)
  567. {
  568. dt = difftime(time(0),tp);
  569. printf("\nMeasurement time: %lf seconds", dt);
  570. printf("\nTime per 10000 measurements: %lf seconds\n", dt*10000.0f/h->n);
  571. }
  572. if (h->DISPQSTATE)
  573. {
  574. printf("\nFinal state:");
  575. printstate(q);
  576. gaussian(q);
  577. printstate(q);
  578. printket(q);
  579. }
  580. return;
  581. }
  582. void preparestate(struct QState *q, char *s)
  583. // Prepare the initial state's "input"
  584. {
  585. long l;
  586. long b;
  587. l = strlen(s);
  588. for (b = 0; b < l; b++)
  589. {
  590. if (s[b]=='Z')
  591. {
  592. hadamard(q,b);
  593. phase(q,b);
  594. phase(q,b);
  595. hadamard(q,b);
  596. }
  597. if (s[b]=='x') hadamard(q,b);
  598. if (s[b]=='X')
  599. {
  600. hadamard(q,b);
  601. phase(q,b);
  602. phase(q,b);
  603. }
  604. if (s[b]=='y')
  605. {
  606. hadamard(q,b);
  607. phase(q,b);
  608. }
  609. if (s[b]=='Y')
  610. {
  611. hadamard(q,b);
  612. phase(q,b);
  613. phase(q,b);
  614. phase(q,b);
  615. }
  616. }
  617. return;
  618. }
  619. void initstae_(struct QState *q, long n, char *s)
  620. // Initialize state q to have n qubits, and input specified by s
  621. {
  622. long i;
  623. long j;
  624. q->n = n;
  625. q->x = malloc((2*q->n + 1) * sizeof(unsigned long*));
  626. q->z = malloc((2*q->n + 1) * sizeof(unsigned long*));
  627. q->r = malloc((2*q->n + 1) * sizeof(int));
  628. q->over32 = (q->n>>5) + 1;
  629. q->pw[0] = 1;
  630. for (i = 1; i < 32; i++)
  631. q->pw[i] = 2*q->pw[i-1];
  632. for (i = 0; i < 2*q->n + 1; i++)
  633. {
  634. q->x[i] = malloc(q->over32 * sizeof(unsigned long));
  635. q->z[i] = malloc(q->over32 * sizeof(unsigned long));
  636. for (j = 0; j < q->over32; j++)
  637. {
  638. q->x[i][j] = 0;
  639. q->z[i][j] = 0;
  640. }
  641. if (i < q->n)
  642. q->x[i][i>>5] = q->pw[i&31];
  643. else if (i < 2*q->n)
  644. {
  645. j = i-q->n;
  646. q->z[i][j>>5] = q->pw[j&31];
  647. }
  648. q->r[i] = 0;
  649. }
  650. if (s) preparestate(q, s);
  651. return;
  652. }
  653. void readprog(struct QProg *h, char *fn, char *params)
  654. // Read a quantum circuit from filename fn, with optional parameters params
  655. {
  656. long t;
  657. char fn2[255];
  658. FILE *fp;
  659. char c=0;
  660. long val;
  661. long l;
  662. h->DISPQSTATE = 0;
  663. h->DISPTIME = 0;
  664. h->SILENT = 0;
  665. h->DISPPROG = 0;
  666. h->SUPPRESSM = 0;
  667. if (params)
  668. {
  669. l = strlen(params);
  670. for (t = 1; t < l; t++)
  671. {
  672. if ((params[t]=='q')||(params[t]=='Q')) h->DISPQSTATE = 1;
  673. if ((params[t]=='p')||(params[t]=='P')) h->DISPPROG = 1;
  674. if ((params[t]=='t')||(params[t]=='T')) h->DISPTIME = 1;
  675. if ((params[t]=='s')||(params[t]=='S')) h->SILENT = 1;
  676. if ((params[t]=='m')||(params[t]=='M')) h->SUPPRESSM = 1;
  677. }
  678. }
  679. sprintf(fn2, "%s", fn);
  680. fp = fopen(fn2, "r");
  681. if (!fp)
  682. {
  683. sprintf(fn2, "%s.chp", fn);
  684. fp = fopen(fn2, "r");
  685. if (!fp) error(1);
  686. }
  687. while (!feof(fp)&&(c!='#'))
  688. fscanf(fp, "%c", &c);
  689. if (c!='#') error(2);
  690. h->T = 0;
  691. h->n = 0;
  692. while (!feof(fp))
  693. {
  694. fscanf(fp, "%c", &c);
  695. if ((c=='\r')||(c=='\n'))
  696. continue;
  697. fscanf(fp, "%ld", &val);
  698. if (val+1 > h->n) h->n = val+1;
  699. if ((c=='c')||(c=='C'))
  700. {
  701. fscanf(fp, "%ld", &val);
  702. if (val+1 > h->n) h->n = val+1;
  703. }
  704. h->T++;
  705. }
  706. fclose(fp);
  707. h->a = malloc(h->T * sizeof(char));
  708. h->b = malloc(h->T * sizeof(long));
  709. h->c = malloc(h->T * sizeof(long));
  710. fp = fopen(fn2, "r");
  711. while (!feof(fp)&&(c!='#'))
  712. fscanf(fp, "%c", &c);
  713. t=0;
  714. while (!feof(fp))
  715. {
  716. fscanf(fp, "%c", &c);
  717. if ((c=='\r')||(c=='\n'))
  718. continue;
  719. if ((c=='c')||(c=='C')) h->a[t] = CNOT;
  720. if ((c=='h')||(c=='H')) h->a[t] = HADAMARD;
  721. if ((c=='p')||(c=='P')) h->a[t] = PHASE;
  722. if ((c=='m')||(c=='M')) h->a[t] = MEASURE;
  723. fscanf(fp, "%ld", &h->b[t]);
  724. if (h->a[t]==CNOT) fscanf(fp, "%ld", &h->c[t]);
  725. t++;
  726. }
  727. fclose(fp);
  728. return;
  729. }
  730. static char init_docs[] = "init(nqubits): Initialize with n qubits";
  731. static PyObject* init(PyObject* self, PyObject* args)
  732. {
  733. int nqubits;
  734. if (!PyArg_ParseTuple(args, "i", &nqubits)) { return NULL; }
  735. srand(time(0));
  736. h = malloc(sizeof(struct QProg));
  737. q = malloc(sizeof(struct QState));
  738. initstae_(q,nqubits,NULL);
  739. PyObject *response = Py_BuildValue("i", nqubits);
  740. return response;
  741. }
  742. static PyMethodDef chp_funcs[] = {
  743. {"init", (PyCFunction)init, METH_VARARGS, init_docs},
  744. {"get_ket", (PyCFunction)get_ket, METH_NOARGS, get_ket_docs},
  745. {"act_hadamard", (PyCFunction)act_hadamard, METH_VARARGS, act_hadamard_docs},
  746. {"act_cnot", (PyCFunction)act_cnot, METH_VARARGS, act_cnot_docs},
  747. {"act_cz", (PyCFunction)act_cz, METH_VARARGS, act_cz_docs},
  748. {"act_phase", (PyCFunction)act_phase, METH_VARARGS, act_phase_docs},
  749. {NULL}
  750. };
  751. void initchp(int argc, char **argv)
  752. {
  753. Py_InitModule3("chp", chp_funcs, "CHP");
  754. }