Anders and Briegel in Python
Вы не можете выбрать более 25 тем Темы должны начинаться с буквы или цифры, могут содержать дефисы(-) и должны содержать не более 35 символов.

8 лет назад
8 лет назад
8 лет назад
8 лет назад
8 лет назад
8 лет назад
8 лет назад
8 лет назад
8 лет назад
8 лет назад
8 лет назад
8 лет назад
8 лет назад
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. }