diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml new file mode 100755 index 0000000..136322d --- /dev/null +++ b/.github/workflows/publish.yml @@ -0,0 +1,28 @@ +name: Publish + +on: + push: + branches: + - main + +jobs: + publish: + runs-on: ubuntu-latest + steps: + - name: Check out + uses: actions/checkout@v4 + + - name: Install Emacs + run: | + sudo apt-add-repository ppa:ubuntu-elisp/ppa + sudo apt update + sudo apt-get install emacs-snapshot + + - name: Build the site + run: ./publi.sh + + - name: Publish generated content to GitHub Pages + uses: JamesIves/github-pages-deploy-action@v4 + with: + branch: gh-pages + folder: public diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..fb23c93 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +.packages \ No newline at end of file diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..a1a4447 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2024 Dr. Ramón L. Panadés-Barrueta + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/README.org b/README.org new file mode 100644 index 0000000..5740135 --- /dev/null +++ b/README.org @@ -0,0 +1,6 @@ +#+TITLE: A natural sciences and array programming blog + +Welcome to the source files from my blog [[https://panadestein.github.io/tac/][The Alien's Cuneiform]]. +Here, I delve into solving computational problems across various +disciplines using my preferred [[https://mlochbaum.github.io/BQN/][tool of thought]]. Every post is actually +an =org-mode= literate programming computational notebook. diff --git a/assets/BQN386.ttf b/assets/BQN386.ttf new file mode 100644 index 0000000..062ab33 Binary files /dev/null and b/assets/BQN386.ttf differ diff --git a/assets/style.css b/assets/style.css new file mode 100644 index 0000000..1c93771 --- /dev/null +++ b/assets/style.css @@ -0,0 +1,60 @@ +/* BQN386 Font Face */ +@font-face { + font-family: 'BQN'; + src: url('BQN386.ttf') format('truetype'); +} + +/* Blog header settings */ + +.container { + max-width: 500px; /* Set your desired maximum width */ + margin: 0 auto; /* Center the container */ + padding: 10px; +} + +.header { + display: flex; + justify-content: space-between; + align-items: center; + padding: 10px 0; +} + +.header-left { + font-size: 1.5em; +} + +/* Applying the BQN386 Font */ +body, input, select, textarea, pre, code, kbd { + font-family: 'BQN', monospace; +} + +body { + text-align: justify; +} + +/* Basic styling for code elements */ +code, pre { + background-color: #e3e7e7; + color: #292929; + font-size: 0.94em; +} + +pre, pre.src { + overflow-x: auto; + border: 1px solid #ccc; + border-radius: 3px; + box-shadow: 3px 3px 3px #eee; + background-color: #e3e7e7; +} + +pre.src { + position: relative; + margin: 0 2em 0 2em; + padding-top: 0.7em; +} + +/* Hyperlinks */ +a:link, a:visited { + color: #267CB9; + text-decoration: none; +} diff --git a/build-site.el b/build-site.el new file mode 100755 index 0000000..534cebd --- /dev/null +++ b/build-site.el @@ -0,0 +1,79 @@ +;;; package --- Build Org website + +;;; Commentary: +;; Build website from Org-mode source files + +;;; Code: + +;; Set a package installation directory to avoid conflicts + +(require 'package) +(setq package-user-dir (expand-file-name "./.packages")) +(setq package-archives '(("melpa" . "https://melpa.org/packages/") + ("elpa" . "https://elpa.gnu.org/packages/"))) + +;; Initialize the package system + +(package-initialize) +(unless package-archive-contents + (package-refresh-contents)) + +;; Install packages needed for HTML export + +(package-install 'htmlize) +(package-install 'reformatter) +(package-install 'color-theme-modern) + +(require 'htmlize) +(require 'ox-publish) +(require 'font-lock) + +;; Using this library is a work-around to get color in HTML exports. +;; Otherwise Emacs in batch mode cannot get the correct faces + +(load-theme 'greiner t t) +(enable-theme 'greiner) + +;; Set some variables for the export + +(global-font-lock-mode t) +(setq org-html-validation-link nil + org-html-head-include-scripts nil + org-html-include-default-style nil + org-confirm-babel-evaluate nil + org-src-fontify-natively t) + +;; Define the project to be published + +(setq org-publish-project-alist + (list + (list "assets" + :base-directory "./assets" + :base-extension "css\\|ttf" + :publishing-directory "./public/assets" + :publishing-function 'org-publish-attachment + :recursive t) + (list "supp" + :base-directory "./supp" + :base-extension 'any + :publishing-directory "./public/supp" + :publishing-function 'org-publish-attachment + :recursive t) + (list "tac" + :recursive t + :base-directory "./src" + :publishing-directory "./public" + :publishing-function 'org-html-publish-to-html + :with-author nil + :with-creator nil + :with-toc nil + :section-numbers nil + :time-stamp-file nil))) + +;; Generate site + +(org-publish-all t) +(message "Build completed") + +(provide 'build-site) +;;; build-site.el ends here diff --git a/publi.sh b/publi.sh new file mode 100755 index 0000000..ec4f9e1 --- /dev/null +++ b/publi.sh @@ -0,0 +1,5 @@ +#!/usr/bin/env bash + +# Run Emacs Lisp publishing script + +emacs -Q --script build-site.el diff --git a/src/hf.org b/src/hf.org new file mode 100644 index 0000000..cd508ac --- /dev/null +++ b/src/hf.org @@ -0,0 +1,3 @@ +# -*- eval: (face-remap-add-relative 'default '(:family "BQN386 Unicode" :height 180)); -*- +#+TITLE: A classical Hartree-Fock program +#+HTML_HEAD: diff --git a/src/index.org b/src/index.org new file mode 100644 index 0000000..40720cf --- /dev/null +++ b/src/index.org @@ -0,0 +1,13 @@ +#+HTML_HEAD: + +#+HTML:
+#+HTML:
+#+HTML:
The Alien's Cuneiform
+#+HTML: +#+HTML:
+ +1. [[./qbqn.org][BQN's Quantum Noise]] +2. [[./spodat.org][Songs to pave the seasons]] +3. [[./hf.org][A classical Hartree-Fock program]] + +#+HTML:
diff --git a/src/qbqn.org b/src/qbqn.org new file mode 100644 index 0000000..d3cdce0 --- /dev/null +++ b/src/qbqn.org @@ -0,0 +1,200 @@ +# -*- eval: (face-remap-add-relative 'default '(:family "BQN386 Unicode" :height 180)); -*- +#+TITLE: BQN's Quantum Noise +#+HTML_HEAD: + +** Preamble + +We will implement and test a compact quantum interpreter in the BQN[fn:1] programming language. +Initially, we import the necessary system functions and define a 1-modifier for handling +complex matrix products. Next, we define a namespace containing various quantum gates: + +#+name: preamble +#+begin_src bqn :exports code :results none :tangle ../supp/perf_qi/q.bqn + Sin‿Cos‿GCD ← •math + U ← •rand.Range + _cp ← {(-´𝔽¨)⋈(+´𝔽¨)⟜⌽} + + g ← { + IM ← (⊢⋈≢⥊⟜0⊢)¨ + x‿id‿h ⇐ (⌽‿⊢{𝕎𝕩}¨<=⌜˜↕2) ∾○IM <-⌾(1‿1⊸⊑)2‿2⥊÷√2 + swap‿cnot ⇐ IM ⟨1‿2, 2‿3⟩ {⌽⌾(𝕨⊸⊏)𝕩}¨ <=⌜˜↕4 + P ⇐ {⟨=⌜˜↕4, 4‿4⥊0⟩ {𝕨⌾(3‿3⊸⊑) 𝕩}¨˜ Sin⊸⋈⟜Cos 𝕩} + } +#+end_src + +** Interpreter + +The (call_count-chars() chars[fn:2]) quantum interpreter is based on references [[https://arxiv.org/abs/1711.02086][arXiv:1711.02086]] +and [[https://arxiv.org/abs/1608.03355][arXiv:1608.03355]]. For simplicity, we always measure at the end of the execution: + +#+name: interpreter +#+begin_src bqn :exports code :results none :tangle ../supp/perf_qi/q.bqn + Q ← {𝕊qb‿sc‿r: + wf ← (1⌾⊑⋈⊢)⥊⟜0 2⋆qb + M‿K ← ⟨+˝∘×⎉1‿∞ _cp, {1𝕊𝕩:𝕩; 𝕨𝕊1:𝕨; 𝕨∾∘×⟜<_cp𝕩}⟩ + E ← {0𝕊𝕩:1; K⍟(𝕨-1)˜𝕩} + L ← {K´ ⟨(qb-𝕨+⌊2⋆⁼≠⊑𝕩) E g.id, 𝕩, 𝕨 E g.id⟩} + T ← ∾↕∘≠{a←𝕩, i←𝕨{𝕩⊑a}•_while_{𝕩<𝕨}𝕨⊑a, 𝕨<◶⟨⟩‿{(⊢∾1⊸↓∘⌽)𝕨+↕𝕩-𝕨}i}¨< + A ← {qs‿gs𝕊v: + 1⊸=◶{𝕊𝕩: ui ← 0 L gs + v M˜ {𝕊⟨⟩:ui; (M˜´M⟜ui⟜M˜M´) L⟜g.swap¨ 𝕩} T qs (⌽∘⊢∾¬∘∊/⊣)˜ ↕qb + }‿(v M˜ gs L˜ ⊑qs) ≠qs} + »⊸<∨` 0>r-`>+○(ט)˝ wf A´ ⌽sc + } +#+end_src + +** Shor's algorithm + +As a test case, we employ the quantum circuit of Shor's algorithm +for the number fifteen and base eleven, following references +[[https://arxiv.org/abs/1804.03719][arXiv:1804.03719]] and [[https://arxiv.org/abs/2306.09122][arXiv:2306.09122]]. The resulting compiled circuit +uses five qubits, three of which serve as control. To enhance +statistical accuracy, the experiment is repeated multiple times. +Additionally, we define a classical post-processing function: + +#+name: test +#+begin_src bqn :exports code :results none :tangle ../supp/perf_qi/q.bqn + n‿a‿qb‿r ← ⟨15, 11, 5, 0 U˜ 2⋆3⟩ + + sc ← ⟨ + ⟨0⟩‿g.h ⋄ ⟨1⟩‿g.h ⋄ ⟨2⟩‿g.h + ⟨2, 3⟩‿g.cnot ⋄ ⟨2, 4⟩‿g.cnot ⋄ ⟨1⟩‿g.h + ⟨⟨1, 0⟩, g.P π÷2⟩ ⋄ ⟨0⟩‿g.h + ⟨⟨1, 2⟩, g.P π÷4⟩ ⋄ ⟨⟨0, 2⟩, g.P π÷2⟩ ⋄ ⟨2⟩‿g.h + ⟩ + + C ← {n (⊣≡×´∘GCD) +‿-{𝕩𝕎1}¨ +˝{Q qb‿sc‿𝕩}¨ r +#+end_src + +#+RESULTS: run +: 1 + +Compare the result with that from a real [[./supp/ibm_eagle/shor_factorize_fifteen.html][quantum computer]]. + +** Epilogue + +Why the hieroglyphs, you may ask? The tacit and functional style, coupled with numerous combinators, +makes programming feel like solving a fun algebraic puzzle rather than drafting a manifesto. +BQN is the epitome of minimalism's power: + +#+begin_export html +
+Primitive's stats +#+end_export + +The src_bqn[:exports code]{prog} string contains the full source code. We used: + +#+begin_src bqn :noweb yes :noweb-prefix no :exports none :tangle no :results none + prog ← "<><><><>" +#+end_src + +#+begin_src bqn :noweb yes :noweb-prefix no :exports both :tangle no :wrap example + prog (+´⊸≍⟜≠∊)˜ ⊑¨•primitives +#+end_src + +#+RESULTS: +#+begin_example +⟨ 44 64 ⟩ +#+end_example + +With this distribution: + +#+begin_src bqn :noweb yes :noweb-prefix no :exports both :tangle no :wrap example + ⍉>(⍷∾≠)¨∘(⊐⊸⊔∊/⊣)⟜(⊑¨•primitives)˜ prog +#+end_src + +#+RESULTS: +#+begin_example +┌─ +╵ '-' '´' '¨' '⋈' '+' '⟜' '⌽' '⊢' '≢' '⥊' '<' '=' '⌜' '˜' '↕' '∾' '○' '⌾' '⊸' '⊑' '÷' '√' '⊏' '⋆' '˝' '∘' '×' '⎉' '≡' '⊣' '⌊' '⁼' '≠' '⍟' '◶' '↓' '¬' '∊' '/' '»' '∨' '`' '>' '⍒' + 8 8 10 5 8 3 6 7 1 5 9 6 3 12 6 5 2 5 7 9 5 1 1 5 4 8 5 1 3 3 1 1 5 1 2 1 1 1 1 1 1 2 3 1 + ┘ +#+end_example + +#+begin_export html +
+#+end_export + +BQN is also fast: + + +#+begin_export html +
+Benchmarks +#+end_export + +While the interpreter's performance is not particularly optimized, here is a comparison with the equivalent Common Lisp code: + +#+begin_src bash :exports results :tangle no :results raw :wrap example + hyperfine --runs 5 'cbqn -f ./perf/q.bqn' 'sbcl --script ./perf/q.lisp' +#+end_src + +#+RESULTS: +#+begin_example +Benchmark 1: cbqn -f ./perf/q.bqn + Time (mean ± σ): 5.468 s ± 0.077 s [User: 5.427 s, System: 0.005 s] + Range (min … max): 5.358 s … 5.535 s 5 runs + +Benchmark 2: sbcl --script ./perf/q.lisp + Time (mean ± σ): 37.114 s ± 0.893 s [User: 37.544 s, System: 0.207 s] + Range (min … max): 36.457 s … 38.634 s 5 runs + +Summary + cbqn -f ./perf/q.bqn ran + 6.79 ± 0.19 times faster than sbcl --script ./perf/q.lisp +#+end_example + +And here is a full program's profile. All time is spent in the Kronecker and matrix products: + +#+begin_src bqn :exports both :tangle no :results raw :wrap example + )profile C >+˝{Q qb‿sc‿𝕩}¨ r +#+end_src + +#+RESULTS: +#+begin_example +Got 25361 samples +(REPL): 25361 samples: + 2│ Q ← {𝕊qb‿sc‿r: + 1│ wf ← (1⌾⊑⋈⊢)⥊⟜0 2⋆qb + 2471│ M‿K ← ⟨+˝∘×⎉1‿∞ _cp, {1𝕊𝕩:𝕩; 𝕨𝕊1:𝕨; 𝕨∾∘×⟜<_cp𝕩}⟩ + 26│ E ← {0𝕊𝕩:1; K⍟(𝕨-1)˜𝕩} + 39│ L ← {K´ ⟨(qb-𝕨+⌊2⋆⁼≠⊑𝕩) E g.id, 𝕩, 𝕨 E g.id⟩} + 16│ T ← ∾↕∘≠{a←𝕩, i←𝕨{𝕩⊑a}•_while_{𝕩<𝕨}𝕨⊑a, 𝕨<◶⟨⟩‿{(⊢∾1⊸↓∘⌽)𝕨+↕𝕩-𝕨}i}¨< + 1│ A ← {qs‿gs𝕊v: + 4│ 1⊸=◶{𝕊𝕩: ui ← 0 L gs + 22430│ v M˜ {𝕊⟨⟩:ui; (M˜´M⟜ui⟜M˜M´) L⟜g.swap¨ 𝕩} T qs (⌽∘⊢∾¬∘∊/⊣)˜ ↕qb + 366│ }‿(v M˜ gs L˜ ⊑qs) ≠qs} + 5│ »⊸<∨` 0>r-`>+○(ט)˝ wf A´ ⌽sc + │ } +#+end_example + +#+begin_export html +
+#+end_export + +Try running the simulation in the call_generate-bqn-link() and explore it! The complete source code is hosted in a GitHub [[https://github.com/Panadestein/qbqn][repository]]. +Additionally, I recommend reading Stylewarning's [[https://www.stylewarning.com/posts/quantum-interpreter/][blog post]], +which served as an inspiration for this one and kindly lists it as a re-implementation. + +#+name: generate-bqn-link +#+begin_src emacs-lisp :noweb yes :noweb-prefix no :exports none :results raw :tangle no + (let* ((bqn-code (concat "<>\n\n" "<>\n\n" "<>\n\n" "<>")) + (encoded (base64-encode-string (encode-coding-string bqn-code 'utf-8) t))) + (concat "[[https://mlochbaum.github.io/BQN/try.html#code=" encoded "][BQN repl]]")) +#+end_src + +#+name: count-chars +#+begin_src emacs-lisp :noweb yes :noweb-prefix no :exports none :results raw :tangle no + (- (length "<>") 4) +#+end_src + +[fn:1] This post's title is a playful recursive acronym that employs quantum computing terminology, without any specific significance beyond that. +[fn:2] I optimized it up to this number, but I wasn't targeting the Kolmogorov complexity. +[fn:3] Hilbert's [[https://maa.org/press/periodicals/convergence/david-hilberts-radio-address-english-translation][radio address]] in 1930. diff --git a/src/spodat.org b/src/spodat.org new file mode 100644 index 0000000..e75d092 --- /dev/null +++ b/src/spodat.org @@ -0,0 +1,92 @@ +# -*- eval: (face-remap-add-relative 'default '(:family "BQN386 Unicode" :height 180)); -*- +#+TITLE: Songs to pave the seasons +#+HTML_HEAD: + +Analysis of my [[https://support.spotify.com/us/article/understanding-my-data/][Spotify data]] for the period 2016-2024. + +** Technical details + +This is a suitable task for an array language, so I rely on [[https://mlochbaum.github.io/BQN/index.html][BQN]] which is my +favorite one: + +#+begin_src bqn :results none + ⟨P⇐Parse⟩ ← •Import "../bqn-libs/json.bqn" + @ ⊣ spd ← P •file.Chars "data.json" +#+end_src + +Additional structural information is needed, namely index of the songs, +artists, elapsed time, and the length of the desired result. + +#+begin_src bqn :results none + s‿a‿m‿l ← ⟨7, 8, 3, 1+↕20⟩ +#+end_src + +Define a (50 chars) dyadic block function to execute the queries: + +#+begin_src bqn :results none + Q ← {l≍˘∾(⍷𝕨⊸⊏˘)¨ l⊏ ((⍒(+´m⊸⊏˘)¨)⊸⊏ 𝕨⊸⊏˘⊐⊸⊔⊢) >1⊏¨ 𝕩} +#+end_src + +** Top songs + +#+begin_src bqn :exports both + s Q spd +#+end_src + +#+RESULTS: +#+begin_example +┌─ +╵ 1 "Countless Skies" + 2 "Divertimento I, K.136: Allegro" + 3 "The Numbers" + 4 "Autre temps" + 5 "Ghost of Perdition" + 6 "Crossing the Road Material" + 7 "Hoppípolla" + 8 "Ether" + 9 "Colossus" + 10 "River" + 11 "El Tete" + 12 "Will o the Wisp" + 13 "Pakumba" + 14 "Damned Rope" + 15 "Eternal Rains Will Come" + 16 "La femme d'argent" + 17 "Bajanda" + 18 "Nimrodel - Medley" + 19 "Breathe (In The Air) - 2011 Remastered Version" + 20 "In The Shadow Of Our Pale Companion" + ┘ +#+end_example + +** Top artists + +#+begin_src bqn :exports both + a Q spd +#+end_src + +#+RESULTS: +#+begin_example +┌─ +╵ 1 "Opeth" + 2 "Wolfgang Amadeus Mozart" + 3 "Pink Floyd" + 4 "Be'lakor" + 5 "Sigur Rós" + 6 "Coldplay" + 7 "Mogwai" + 8 "Chocolate Mc" + 9 "Radiohead" + 10 "Joaquín Sabina" + 11 "Iron & Wine" + 12 "Rammstein" + 13 "Alcest" + 14 "Buena Fe" + 15 "Silvio Rodríguez" + 16 "Amon Amarth" + 17 "In Mourning" + 18 "Lamb of God" + 19 "Camel" + 20 "Omnium Gatherum" + ┘ +#+end_example diff --git a/supp/hf_zo/README.org b/supp/hf_zo/README.org new file mode 100644 index 0000000..13976c9 --- /dev/null +++ b/supp/hf_zo/README.org @@ -0,0 +1,10 @@ +#+TITLE: The Hartree-Fock implementation from Szabo-Ostlund + +This is pretty much the original code listed in Appendix B of the legendary [[https://store.doverpublications.com/products/9780486691862][book]]. It has some +minimal modifications introduced [[http://www.ccl.net/cca/software/SOURCES/FORTRAN/szabo/index.html][here]]. Compiling it is not so easy anymore, but it is still +possible with legacy options: + +#+begin_src bash + gfortran -o hfzo hf_so.f -std=legacy +#+end_src + diff --git a/supp/hf_zo/hf_so.f b/supp/hf_zo/hf_so.f new file mode 100644 index 0000000..d963068 --- /dev/null +++ b/supp/hf_zo/hf_so.f @@ -0,0 +1,570 @@ +C******************************************************************** +C +C MINIMAL BASIS STO-3G CALCULATION ON HEH+ +C +C THIS IS A LITTLE DUMMY MAIN PROGRAM WHICH CALLS HFCALC +C +C APPENDIX B: TWO-ELECTRON SELF-CONSISTENT-FIELD PROGRAM +C OF MODERN QUANTUM CHEMISTRY by +C Attila Szabo and Neil S. Ostlund +C Ed. 2nd (1989) Dover Publications INC. +C +C Labourly Typed by Michael Zitolo (Feb., 2005) +C Edited and Compiled by Michael Zitolo and Xihua Chen +C +C Cleaned up and debugged again by Andrew Long (2012) +C and Daniele (kalium) Dondi (2013) +C********************************************************************* + + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + IOP=2 + N=3 + R=1.4632D0 + ZETA1=2.0925D0 + ZETA2=1.24D0 + ZA=2.0D0 + ZB=1.0D0 + CALL HFCALC(IOP,N,R,ZETA1,ZETA2,ZA,ZB) + END + +C********************************************************************* + SUBROUTINE HFCALC(IOP,N,R,ZETA1,ZETA2,ZA,ZB) +C +C DOES A HARTREE-FOCK CALCULATION FOR A TWO-ELECTRON DIATOMIC +C USING THE 1S MINIMAL STO-NG BASIS SET +C MINIMAL BASIS SET HAS BASIS FUNCTIONS 1 AND 2 ON NUCLEI A AND B +C +C IOP=0 NO PRINTING WHATSOEVER (TO OPTIMIZE EXPONENTS, SAY) +C IOP=1 PRINT ONLY CONVERGED RESULTS +C IOP=2 PRINT EVERY ITERATION +C N STO-NG CALCULATION (N=1,2 OR 3) +C R BONDLENGTH (AU) +C ZETA1 SLATER ORBITAL EXPONENT (FUNCTION 1) +C ZETA2 SLATER ORBITAL EXPONENT (FUNCTION 2) +C ZA ATOMIC NUMBER (ATOM A) +C ZB ATOMIC NUMBER (ATOM B) +C +C********************************************************************* + + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + IF (IOP.EQ.0) GO TO 20 + PRINT 10,N,ZA,ZB + 10 FORMAT(' ',2X,'STO-',I1,'G FOR ATOMIC NUMBERS ',F5.2,' AND ', + $ F5.2//) + 20 CONTINUE +C CALCULATE ALL THE ONE AND TWO ELECTRON INTEGRALS + CALL INTGRL(IOP,N,R,ZETA1,ZETA2,ZA,ZB) +C BE INEFFICIENT AND PUT ALL INTEGRALS IN PRETTY ARRAYS + CALL COLECT(IOP,N,R,ZETA1,ZETA2,ZA,ZB) +C PERFORM THE SCF CALCULATION + CALL SCF(IOP,N,R,ZETA1,ZETA2,ZA,ZB) + RETURN + END + +C********************************************************************* + SUBROUTINE INTGRL(IOP,N,R,ZETA1,ZETA2,ZA,ZB) +C +C CALCULATES ALL THE BASIC INTEGRALS NEEDED FOR SCF CALCULATION +C +C********************************************************************* + + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + COMMON/INT/S12,T11,T12,T22,V11A,V12A,V22A,V11B,V12B,V22B, + $ V1111,V2111,V2121,V2211,V2221,V2222 + DIMENSION COEF(3,3),EXPON(3,3),D1(3),A1(3),D2(3),A2(3) + DATA PI/3.1415926535898D0/ +C THESE ARE THE CONTRACTION COEFFICIENTS AND EXPONENTS FOR +C A NORMALIZED SLATER ORBITAL WITH EXPONENT 1.0 IN TERMS OF +C NORMALIZED 1S PRIMITIVE GAUSSIANS + DATA COEF,EXPON/1.0D0,2*0.0D0,0.678914D0,0.430129D0,0.0D0, + $ 0.444635D0,0.535328D0,0.154329D0,0.270950D0,2*0.0D0,0.151623D0, + $ 0.851819D0,0.0D0,0.109818D0,0.405771D0,2.22766D0/ + R2=R*R +C SCALE THE EXPONENTS (A) OF PRIMITIVE GAUSSIANS +C INCLUDE NORMALIZATION IN CONTRACTION COEFFICIENTS (D) + DO 10 I=1,N + A1(I)=EXPON(I,N)*(ZETA1**2) + D1(I)=COEF(I,N)*((2.0D0*A1(I)/PI)**0.75D0) + A2(I)=EXPON(I,N)*(ZETA2**2) + D2(I)=COEF(I,N)*((2.0D0*A2(I)/PI)**0.75D0) + 10 CONTINUE +C D AND A ARE NOW THE CONTRACTION COEFFICIENTS AND EXPONENTS +C IN TERMS OF UNNORMALIZED PRIMITIVE GAUSSIANS + S12=0.0D0 + T11=0.0D0 + T12=0.0D0 + T22=0.0D0 + V11A=0.0D0 + V12A=0.0D0 + V22A=0.0D0 + V11B=0.0D0 + V12B=0.0D0 + V22B=0.0D0 + V1111=0.0D0 + V2111=0.0D0 + V2121=0.0D0 + V2211=0.0D0 + V2221=0.0D0 + V2222=0.0D0 +C CALCULATE ONE-ELECTRON INTEGRALS +C CENTER A IS FIRST ATOM, CETER B IS SECOND ATOM +C ORIGIN IS ON CENTER A +C V12A = OFF-DIAGONAL NUCLEAR ATTRACTION TO CENTER A, ETC. + DO 20 I=1,N + DO 20 J=1,N +C RAP2 = SQUARED DISTANCE BETWEEN CENTER A AND CENTER P, ETC. + RAP=A2(J)*R/(A1(I)+A2(J)) + RAP2=RAP**2 + RBP2=(R-RAP)**2 + S12=S12+S(A1(I),A2(J),R2)*D1(I)*D2(J) + T11=T11+T(A1(I),A1(J),0.0D0)*D1(I)*D1(J) + T12=T12+T(A1(I),A2(J),R2)*D1(I)*D2(J) + T22=T22+T(A2(I),A2(J),0.0D0)*D2(I)*D2(J) + V11A=V11A+V(A1(I),A1(J),0.0D0,0.0D0,ZA)*D1(I)*D1(J) + V12A=V12A+V(A1(I),A2(J),R2,RAP2,ZA)*D1(I)*D2(J) + V22A=V22A+V(A2(I),A2(J),0.0D0,R2,ZA)*D2(I)*D2(J) + V11B=V11B+V(A1(I),A1(J),0.0D0,R2,ZB)*D1(I)*D1(J) + V12B=V12B+V(A1(I),A2(J),R2,RBP2,ZB)*D1(I)*D2(J) + V22B=V22B+V(A2(I),A2(J),0.0D0,0.0D0,ZB)*D2(I)*D2(J) + 20 CONTINUE +C CALCULATE TWO-ELECTRON INTEGRALS + DO 30 I=1,N + DO 30 J=1,N + DO 30 K=1,N + DO 30 L=1,N + RAP=A2(I)*R/(A2(I)+A1(J)) + RBP=R-RAP + RAQ=A2(K)*R/(A2(K)+A1(L)) + RBQ=R-RAQ + RPQ=RAP-RAQ + RAP2=RAP*RAP + RBP2=RBP*RBP + RAQ2=RAQ*RAQ + RBQ2=RBQ*RBQ + RPQ2=RPQ*RPQ + V1111=V1111+TWOE(A1(I),A1(J),A1(K),A1(L),0.0D0,0.0D0,0.0D0) + $ *D1(I)*D1(J)*D1(K)*D1(L) + V2111=V2111+TWOE(A2(I),A1(J),A1(K),A1(L),R2,0.0D0,RAP2) + $ *D2(I)*D1(J)*D1(K)*D1(L) + V2121=V2121+TWOE(A2(I),A1(J),A2(K),A1(L),R2,R2,RPQ2) + $ *D2(I)*D1(J)*D2(K)*D1(L) + V2211=V2211+TWOE(A2(I),A2(J),A1(K),A1(L),0.0D0,0.0D0,R2) + $ *D2(I)*D2(J)*D1(K)*D1(L) + V2221=V2221+TWOE(A2(I),A2(J),A2(K),A1(L),0.0D0,R2,RBQ2) + $ *D2(I)*D2(J)*D2(K)*D1(L) + V2222=V2222+TWOE(A2(I),A2(J),A2(K),A2(L),0.0D0,0.0D0,0.0D0) + $ *D2(I)*D2(J)*D2(K)*D2(L) + 30 CONTINUE + IF (IOP.EQ.0) GO TO 90 + PRINT 40 + 40 FORMAT(3X,'R',10X,'ZETA1',6X,'ZETA2',6X,'S12',8X,'T11'/) + PRINT 50, R,ZETA1,ZETA2,S12,T11 + 50 FORMAT(5F11.6//) + PRINT 60 + 60 FORMAT(3X,'T12',8X,'T22',8X,'V11A',7X,'V12A',7X,'V22A'/) + PRINT 50, T12,T22,V11A,V12A,V22A + PRINT 70 + 70 FORMAT(3X,4HV11B,7X,4HV12B,7X,4HV22B,7X,'V1111',6X,'V2111'/) + PRINT 50, V11B,V12B,V22B,V1111,V2111 + PRINT 80 + 80 FORMAT(3X,5HV2121,6X,5HV2211,6X,5HV2221,6X,5HV2222/) + PRINT 50, V2121,V2211,V2221,V2222 + 90 RETURN + END + +C********************************************************************* + FUNCTION F0(ARG) +C +C CALCULATES THE F FUNCTION +C FO ONLY (S-TYPE ORBITALS) +C +C********************************************************************* + + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + DATA PI/3.1415926535898D0/ + IF (ARG.LT.1.0D-6) GO TO 10 +C F0 IN TERMS OF THE ERROR FUNCTION + F0=DSQRT(PI/ARG)*DERFOTHER(DSQRT(ARG))/2.0D0 + GO TO 20 +C ASYMPTOTIC VALUE FOR SMALL ARGUMENTS + 10 F0=1.0D0-ARG/3.0D0 + 20 CONTINUE + RETURN + END + +C********************************************************************* + FUNCTION DERFOTHER(ARG) +C +C CALCULATES THE ERROR FUNCTION ACCORDING TO A RATIONAL +C APPROXIMATION FROM M. ARBRAMOWITZ AND I.A. STEGUN, +C HANDBOOK OF MATHEMATICAL FUNCTIONS, DOVER. +C ABSOLUTE ERROR IS LESS THAN 1.5*10**(-7) +C CAN BE REPLACED BY A BUILT-IN FUNCTION ON SOME MACHINES +C +C********************************************************************* + + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + DIMENSION A(5) + DATA P/0.3275911D0/ + DATA A/0.254829592D0,-0.284496736D0,1.421413741D0, + $ -1.453152027D0,1.061405429D0/ + T=1.0D0/(1.0D0+P*ARG) + TN=T + POLY=A(1)*TN + DO 10 I=2,5 + TN=TN*T + POLY=POLY+A(I)*TN + 10 CONTINUE + DERFOTHER=1.0D0-POLY*DEXP(-ARG*ARG) + RETURN + END + +C********************************************************************* + FUNCTION S(A,B,RAB2) +C +C CALCULATES OVERLAPS FOR UN-NORMALIZED PRIMITIVES +C +C********************************************************************* + + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + DATA PI/3.1415926535898D0/ + S=(PI/(A+B))**1.5D0*DEXP(-A*B*RAB2/(A+B)) + RETURN + END + +C********************************************************************* + FUNCTION T(A,B,RAB2) +C +C CALCULATES KINETIC ENERGY INTEGRALS FOR UN-NORMALIZED PRIMITIVES +C +C********************************************************************* + + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + DATA PI/3.1415926535898D0/ + T=A*B/(A+B)*(3.0D0-2.0D0*A*B*RAB2/(A+B))*(PI/(A+B))**1.5D0 + $ *DEXP(-A*B*RAB2/(A+B)) + RETURN + END + +C********************************************************************* + FUNCTION V(A,B,RAB2,RCP2,ZC) +C +C CALCULATES UN-NORMALIZED NUCLEAR ATTRACTION INTEGRALS +C +C********************************************************************* + + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + DATA PI/3.1415926535898D0/ + V=2.0D0*PI/(A+B)*F0((A+B)*RCP2)*DEXP(-A*B*RAB2/(A+B)) + V=-V*ZC + RETURN + END + +C********************************************************************* + FUNCTION TWOE(A,B,C,D,RAB2,RCD2,RPQ2) +C +C CALCULATES TWO-ELECTRON INTEGRALS FOR UN-NORMALIZED PRIMITIVES +C A,B,C,D ARE THE EXPONENTS ALPHA, BETA, ETC. +C RAB2 EQUALS SQUARED DISTANCE BETWEEN CENTER A AND CENTER B, ETC. +C********************************************************************* + + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + DATA PI/3.1415926535898D0/ + TWOE=2.0D0*(PI**2.5D0)/((A+B)*(C+D)*DSQRT(A+B+C+D)) + $ *F0((A+B)*(C+D)*RPQ2/(A+B+C+D)) + $ *DEXP(-A*B*RAB2/(A+B)-C*D*RCD2/(C+D)) + RETURN + END + +C********************************************************************* + SUBROUTINE COLECT(IOP,N,R,ZETA1,ZETA2,ZA,ZB) +C +C THIS TAKES THE BASIC INTEGRALS FROM COMMON AND ASSEMBLES THE +C RELEVENT MATRICES, THAT IS S,H,X,XT, AND TWO-ELECTRON INTEGRALS +C +C********************************************************************* + + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + COMMON/MATRIX/S(2,2),X(2,2),XT(2,2),H(2,2),F(2,2),G(2,2),C(2,2), + $ FPRIME(2,2),CPRIME(2,2),P(2,2),OLDP(2,2),TT(2,2,2,2),E(2,2) + COMMON/INT/S12,T11,T12,T22,V11A,V12A,V22A,V11B,V12B,V22B, + $ V1111,V2111,V2121,V2211,V2221,V2222 +C FORM CORE HAMILTONIAN + H(1,1)=T11+V11A+V11B + H(1,2)=T12+V12A+V12B + H(2,1)=H(1,2) + H(2,2)=T22+V22A+V22B +C FORM OVERLAP MATRIX + S(1,1)=1.0D0 + S(1,2)=S12 + S(2,1)=S(1,2) + S(2,2)=1.0D0 +C USE CANONICAL ORTHOGONALIZATION + X(1,1)=1.0D0/DSQRT(2.0D0*(1.0D0+S12)) + X(2,1)=X(1,1) + X(1,2)=1.0D0/DSQRT(2.0D0*(1.0D0-S12)) + X(2,2)=-X(1,2) +C TRANSPOSE OF TRANSFORMATION MATRIX + XT(1,1)=X(1,1) + XT(1,2)=X(2,1) + XT(2,1)=X(1,2) + XT(2,2)=X(2,2) +C MATRIX OF TWO-ELE�CTRON INTEGRALS + TT(1,1,1,1)=V1111 + TT(2,1,1,1)=V2111 + TT(1,2,1,1)=V2111 + TT(1,1,2,1)=V2111 + TT(1,1,1,2)=V2111 + TT(2,1,2,1)=V2121 + TT(1,2,2,1)=V2121 + TT(2,1,1,2)=V2121 + TT(1,2,1,2)=V2121 + TT(2,2,1,1)=V2211 + TT(1,1,2,2)=V2211 + TT(2,2,2,1)=V2221 + TT(2,2,1,2)=V2221 + TT(2,1,2,2)=V2221 + TT(1,2,2,2)=V2221 + TT(2,2,2,2)=V2222 + IF (IOP.EQ.0) GO TO 40 + CALL MATOUT(S,2,2,2,2,4HS ) + CALL MATOUT(X,2,2,2,2,4HX ) + CALL MATOUT(H,2,2,2,2,4HH ) + PRINT 10 + 10 FORMAT(//) + DO 30 I=1,2 + DO 30 J=1,2 + DO 30 K=1,2 + DO 30 L=1,2 + PRINT 20, I,J,K,L,TT(I,J,K,L) + 20 FORMAT(3X,1H(,4I2,2H ),F10.6) + 30 CONTINUE + 40 RETURN + END + +C********************************************************************* + SUBROUTINE SCF(IOP,N,R,ZETA1,ZETA2,ZA,ZB) +C +C PERFORMS THE SCF ITERATIONS +C +C********************************************************************* + + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + COMMON/MATRIX/S(2,2),X(2,2),XT(2,2),H(2,2),F(2,2),G(2,2),C(2,2), + $ FPRIME(2,2),CPRIME(2,2),P(2,2),OLDP(2,2),TT(2,2,2,2),E(2,2) + DATA PI/3.1415926535898D0/ +C CONVERGENCE CRITERION FOR DENSITY MATRIX + DATA CRIT/1.0D-4/ +C MAXIMUM NUMBER OF ITERATIONS + DATA MAXIT/25/ +C ITERATION NUMBER + ITER=0 +C USE CORE-HAMILTONIAN FOR INITIAL GUESS AT F, I.E. (P=0) + DO 10 I=1,2 + DO 10 J=1,2 + 10 P(I,J)=0.0D0 + IF (IOP.LT.2) GO TO 20 + CALL MATOUT(P,2,2,2,2,4HP ) +C START OF ITERATION LOOP + 20 ITER=ITER+1 + IF (IOP.LT.2) GO TO 40 + PRINT 30, ITER + 30 FORMAT(/,4X,28HSTART OF ITERATION NUMBER = ,I2) + 40 CONTINUE +C FORM TWO-ELECTRON PART OF FOCK MATRIX FROM P + CALL FORMG + IF (IOP.LT.2) GO TO 50 + CALL MATOUT(G,2,2,2,2,4HG ) + 50 CONTINUE +C ADD CORE HAMILTONIAN TO GET FOCK MATRIX + DO 60 I=1,2 + DO 60 J=1,2 + F(I,J) = H(I,J)+G(I,J) + 60 CONTINUE +C CALCULATE ELECTRONIC ENERGY + EN=0.0D0 + DO 70 I=1,2 + DO 70 J=1,2 + EN=EN+0.5D0*P(I,J)*(H(I,J)+F(I,J)) + 70 CONTINUE + IF (IOP.LT.2) GO TO 90 + CALL MATOUT(F,2,2,2,2,4HF ) + PRINT 80, EN + 80 FORMAT(///,4X,20HELECTRONIC ENERGY = ,D20.12) + 90 CONTINUE +C TRANSFORM FOCK MATRIX USING G FOR TEMPORARY STORAGE + CALL MULT(F,X,G,2,2) + CALL MULT(XT,G,FPRIME,2,2) +C DIAGONALIZE TRANSFORMED FOCK MATRIX + CALL DIAG(FPRIME,CPRIME,E) +C TRANSFORM EIGENVECTORS TO GET MATRIX C + CALL MULT(X,CPRIME,C,2,2) +C FORM NEW DENSITY MATRIX + DO 100 I=1,2 + DO 100 J=1,2 +C SAVE PRESENT DENSITY MATRIX +C BEFORE CREATING NEW ONE + OLDP(I,J)=P(I,J) + P(I,J)=0.0D0 + DO 100 K=1,1 + P(I,J)=P(I,J)+2.0D0*C(I,K)*C(J,K) + 100 CONTINUE + IF (IOP.LT.2) GO TO 110 + CALL MATOUT(FPRIME,2,2,2,2,"F' ") + CALL MATOUT(CPRIME,2,2,2,2,"C' ") + CALL MATOUT(E,2,2,2,2,'E ') + CALL MATOUT(C,2,2,2,2,'C ') + CALL MATOUT(P,2,2,2,2,'P ') + 110 CONTINUE +C CALCULATE DELTA + DELTA=0.0D0 + DO 120 I=1,2 + DO 120 J=1,2 + DELTA=DELTA+(P(I,J)-OLDP(I,J))**2 + 120 CONTINUE + DELTA=DSQRT(DELTA/4.0D0) + IF (IOP.EQ.0) GO TO 140 + PRINT 130, DELTA + 130 FORMAT(/,4X,39HDELTA(CONVERGENCE OF DENSITY MATRIX) = + $F10.6,/) + 140 CONTINUE +C CHECK FOR CONVERGENCE + IF (DELTA.LT.CRIT) GO TO 160 +C NOT YET CONVERGED +C TEST FOR MAXIMUM NUMBER OF ITERATIONS +C IF MAXIMUM NUMBER NOT YET REACHED +C GO BACK FOR ANOTHER ITERATION + IF(ITER.LT.MAXIT) GO TO 20 +C SOMETHING WRONG HERE + PRINT 150 + 150 FORMAT(4X,21HNO CONVERGENCE IN SCF) + STOP + 160 CONTINUE +C CALCULATION CONVERGED IF IT GOT HERE +C ADD NUCLEAR REPULSION TO GET TOTAL ENERGY + ENT=EN+ZA*ZB/R + IF (IOP.EQ.0) GO TO 180 + PRINT 170, EN, ENT + 170 FORMAT(//,4X,21HCALCULATION CONVERGED,//, + $4X,20HELECTRONIC ENERGY = ,D20.12,//, + $4X,20HTOTAL ENERGY = ,D20.12 ) + 180 CONTINUE + IF (IOP.NE.1) GO TO 190 +C PRINT OUT THE FINAL RESULTS IF +C HAVE NOT DONE SO ALREADY + CALL MATOUT(G,2,2,2,2,4HG ) + CALL MATOUT(F,2,2,2,2,4HF ) + CALL MATOUT(E,2,2,2,2,4HE ) + CALL MATOUT(C,2,2,2,2,4HC ) + CALL MATOUT(P,2,2,2,2,4HP ) + 190 CONTINUE +C PS MATRIX HAS MULLIKEN POPULATIONS + CALL MULT(P,S,OLDP,2,2) + IF(IOP.EQ.0) GO TO 200 + CALL MATOUT(OLDP,2,2,2,2,4HPS ) + 200 CONTINUE + RETURN + END + +C********************************************************************* + SUBROUTINE FORMG +C +C CALCULATES THE G MATRIX FROM THE DENSITY MATRIX +C AND TWO-ELECTRON INTEGRALS +C +C********************************************************************* + + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + COMMON/MATRIX/S(2,2),X(2,2),XT(2,2),H(2,2),F(2,2),G(2,2),C(2,2), + $FPRIME(2,2),CPRIME(2,2),P(2,2),OLDP(2,2),TT(2,2,2,2),E(2,2) + DO 10 I=1,2 + DO 10 J=1,2 + G(I,J)=0.0D0 + DO 10 K=1,2 + DO 10 L=1,2 + G(I,J)=G(I,J)+P(K,L)*(TT(I,J,K,L)-0.5D0*TT(I,L,K,J)) + 10 CONTINUE + RETURN + END + +C********************************************************************* + SUBROUTINE DIAG(F,C,E) +C +C DIAGONALIZES F TO GIVE EIGENVECTORS IN C AND EIGENVALUES IN E +C THETA IS THE ANGLE DESCRIBING SOLUTION +C +C********************************************************************* + + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + DIMENSION F(2,2),C(2,2),E(2,2) + DATA PI/3.1415926535898D0/ + IF (DABS(F(1,1)-F(2,2)).GT.1.0D-20) GO TO 10 +C HERE IS SYMMETRY DETERMINED SOLUTION (HOMONUCLEAR DIATOMIC) + THETA=PI/4.0D0 + GO TO 20 + 10 CONTINUE +C SOLUTION FOR HETERONUCLEAR DIATOMIC + THETA=0.5D0*DATAN(2.0D0*F(1,2)/(F(1,1)-F(2,2))) + 20 CONTINUE + C(1,1)=DCOS(THETA) + C(2,1)=DSIN(THETA) + C(1,2)=DSIN(THETA) + C(2,2)=-DCOS(THETA) + E(1,1)=F(1,1)*DCOS(THETA)**2+F(2,2)*DSIN(THETA)**2 + $ +F(1,2)*DSIN(2.0D0*THETA) + E(2,2)=F(2,2)*DCOS(THETA)**2+F(1,1)*DSIN(THETA)**2 + $ -F(1,2)*DSIN(2.0D0*THETA) + E(2,1)=0.0D0 + E(1,2)=0.0D0 +C ORDER EIGENVALUES AND EIGENVECTORS + IF (E(2,2).GT.E(1,1)) GO TO 30 + TEMP=E(2,2) + E(2,2)=E(1,1) + E(1,1)=TEMP + TEMP=C(1,2) + C(1,2)=C(1,1) + C(1,1)=TEMP + TEMP=C(2,2) + C(2,2)=C(2,1) + C(2,1)=TEMP + 30 RETURN + END + +C********************************************************************* + SUBROUTINE MULT(A,B,C,IM,M) +C +C MULTIPLIES TWO SQUARE MATRICES A AND B TO GET C +C +C********************************************************************* + + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + DIMENSION A(IM,IM),B(IM,IM),C(IM,IM) + DO 10 I=1,M + DO 10 J=1,M + C(I,J)=0.0D0 + DO 10 K=1,M + 10 C(I,J)=C(I,J)+A(I,K)*B(K,J) + RETURN + END + +C********************************************************************* + SUBROUTINE MATOUT(A,IM,IN,M,N,LABEL) +C +C PRINT MATRICES OF SIZE M BY N +C +C********************************************************************* + + IMPLICIT DOUBLE PRECISION(A-H,O-Z) + DIMENSION A(IM,IN) + IHIGH=0 + 10 LOW=IHIGH+1 + IHIGH=IHIGH+5 + IHIGH=MIN(IHIGH,N) + PRINT 20, LABEL,(I,I=LOW,IHIGH) + 20 FORMAT(///,3X,5H THE ,A4,6H ARRAY,/,15X,5(10X,I3,6X)//) + DO 30 I=1,M + 30 PRINT 40, I,(A(I,J),J=LOW,IHIGH) + 40 FORMAT(I10,5X,5(1X,D18.10)) + IF (N-IHIGH) 50,50,10 + 50 RETURN + END diff --git a/supp/ibm_eagle/shor_factorize_fifteen.html b/supp/ibm_eagle/shor_factorize_fifteen.html new file mode 100644 index 0000000..fef5007 --- /dev/null +++ b/supp/ibm_eagle/shor_factorize_fifteen.html @@ -0,0 +1,7808 @@ + + + + + +shor_factorize_fifteen + + + + + + + + + + + + +
+ + + + + + + + + +
+ + diff --git a/supp/perf_qi/q.bqn b/supp/perf_qi/q.bqn new file mode 100644 index 0000000..3684deb --- /dev/null +++ b/supp/perf_qi/q.bqn @@ -0,0 +1,36 @@ +Sin‿Cos‿GCD ← •math +U ← •rand.Range +_cp ← {(-´𝔽¨)⋈(+´𝔽¨)⟜⌽} + +g ← { + IM ← (⊢⋈≢⥊⟜0⊢)¨ + x‿id‿h ⇐ (⌽‿⊢{𝕎𝕩}¨<=⌜˜↕2) ∾○IM <-⌾(1‿1⊸⊑)2‿2⥊÷√2 + swap‿cnot ⇐ IM ⟨1‿2, 2‿3⟩ {⌽⌾(𝕨⊸⊏)𝕩}¨ <=⌜˜↕4 + P ⇐ {⟨=⌜˜↕4, 4‿4⥊0⟩ {𝕨⌾(3‿3⊸⊑) 𝕩}¨˜ Sin⊸⋈⟜Cos 𝕩} +} + +Q ← {𝕊qb‿sc‿r: + wf ← (1⌾⊑⋈⊢)⥊⟜0 2⋆qb + M‿K ← ⟨+˝∘×⎉1‿∞ _cp, {1𝕊𝕩:𝕩; 𝕨𝕊1:𝕨; 𝕨∾∘×⟜<_cp𝕩}⟩ + E ← {0𝕊𝕩:1; K⍟(𝕨-1)˜𝕩} + L ← {K´ ⟨(qb-𝕨+⌊2⋆⁼≠⊑𝕩) E g.id, 𝕩, 𝕨 E g.id⟩} + T ← ∾↕∘≠{a←𝕩, i←𝕨{𝕩⊑a}•_while_{𝕩<𝕨}𝕨⊑a, 𝕨<◶⟨⟩‿{(⊢∾1⊸↓∘⌽)𝕨+↕𝕩-𝕨}i}¨< + A ← {qs‿gs𝕊v: + 1⊸=◶{𝕊𝕩: ui ← 0 L gs + v M˜ {𝕊⟨⟩:ui; (M˜´M⟜ui⟜M˜M´) L⟜g.swap¨ 𝕩} T qs (⌽∘⊢∾¬∘∊/⊣)˜ ↕qb + }‿(v M˜ gs L˜ ⊑qs) ≠qs} + »⊸<∨` 0>r-`>+○(ט)˝ wf A´ ⌽sc +} + +n‿a‿qb‿r ← ⟨15, 11, 5, 0 U˜ 2⋆3⟩ + +sc ← ⟨ + ⟨0⟩‿g.h ⋄ ⟨1⟩‿g.h ⋄ ⟨2⟩‿g.h + ⟨2, 3⟩‿g.cnot ⋄ ⟨2, 4⟩‿g.cnot ⋄ ⟨1⟩‿g.h + ⟨⟨1, 0⟩, g.P π÷2⟩ ⋄ ⟨0⟩‿g.h + ⟨⟨1, 2⟩, g.P π÷4⟩ ⋄ ⟨⟨0, 2⟩, g.P π÷2⟩ ⋄ ⟨2⟩‿g.h +⟩ + +C ← {n (⊣≡×´∘GCD) +‿-{𝕩𝕎1}¨ +˝{Q qb‿sc‿𝕩}¨ r diff --git a/supp/perf_qi/q.lisp b/supp/perf_qi/q.lisp new file mode 100644 index 0000000..b9d4f3b --- /dev/null +++ b/supp/perf_qi/q.lisp @@ -0,0 +1,229 @@ +(defstruct machine + quantum-state + measurement-register) + +(defun make-quantum-state (n) + (let ((s (make-array (expt 2 n) :initial-element 0.0d0))) + (setf (aref s 0) 1.0d0) + s)) + +(defun dimension-qubits (d) + (- (integer-length d) 1)) + +(defun apply-operator (unitary-operator ket) + (let* ((oper-size (array-dimension unitary-operator 0)) + (new-ket (make-array oper-size :initial-element 0.0d0))) + (dotimes (i oper-size) + (let ((element 0.0d0)) + (dotimes (j oper-size) + (incf element (* (aref unitary-operator i j) (aref ket j)))) + (setf (aref new-ket i) element))) + (replace ket new-ket))) + +(defun compose-operator (unitary-left unitary-right) + (let* ((m (array-dimension unitary-left 0)) + (n (array-dimension unitary-left 1)) + (p (array-dimension unitary-right 1)) + (unitary-new (make-array (list m p) :initial-element 0.0d0))) + (dotimes (i m unitary-new) + (dotimes (j p) + (let ((dot-prod 0.0d0)) + (dotimes (k n) + (setf dot-prod (+ dot-prod + (* (aref unitary-left i k) + (aref unitary-right k j))))) + (setf (aref unitary-new i j) dot-prod)))))) + +(defun observe (machine) + (let ((b (sample (machine-quantum-state machine)))) + (collapse (machine-quantum-state machine) b) + (setf (machine-measurement-register machine) b) + machine)) + +(defun sample (state) + (let ((u (random 1.0d0))) + (dotimes (i (length state)) + (decf u (expt (abs (aref state i)) 2)) + (when (minusp u) (return i))))) + +(defun collapse (state qubit) + (fill state 0.0d0) + (setf (aref state qubit) 1.0d0)) + +(defparameter +I+ #2A((1 0) + (0 1)) + "The identity gate.") + +(defun apply-gate (state unitary-operator qubits) + (assert (= (length qubits) + (dimension-qubits (array-dimension unitary-operator 0)))) + (if (= (length qubits) 1) + (apply-singleq-gate state unitary-operator (first qubits)) + (apply-multiq-gate state unitary-operator qubits))) + +(defun kronecker-product (unitary-A unitary-B) + (destructuring-bind (m n) (array-dimensions unitary-A) + (destructuring-bind (p q) (array-dimensions unitary-B) + (let ((kron-dot-unitary (make-array (list (* m p) (* n q))))) + (dotimes (i m kron-dot-unitary) + (dotimes (j n) + (let ((unitary-A-ij (aref unitary-A i j)) + (pointer-row (* i p)) + (pointer-col (* j q))) + (dotimes (k p) + (dotimes (l q) + (setf (aref kron-dot-unitary + (+ pointer-row k) (+ pointer-col l)) + (* unitary-A-ij (aref unitary-B k l)))))))))))) + +(defun kronecker-expt (unitary n) + (cond + ((< n 1) #2A((1))) + ((= n 1) unitary) + (t (kronecker-product (kronecker-expt unitary (- n 1)) unitary)))) + +(defun lift (unitary i n) + (let ((left-factors (kronecker-expt +I+ (- n i (dimension-qubits + (array-dimension unitary 0))))) + (right-factors (kronecker-expt +I+ i))) + (kronecker-product left-factors + (kronecker-product unitary right-factors)))) + +(defparameter +SWAP+ #2A((1 0 0 0) + (0 0 1 0) + (0 1 0 0) + (0 0 0 1)) + "The SWAP gate") + +(defun perm-as-trans (permutation) + (let ((transpositions nil)) + (dotimes (natural (length permutation) (nreverse transpositions)) + (let ((permuted (elt permutation natural))) + (loop :while (< permuted natural) :do + (setf permuted (elt permutation permuted))) + (when (> permuted natural) + (push (cons natural permuted) transpositions)))))) + +(defun trans-as-adjacent (transpositions) + (mapcan (lambda (trans) + (let ((start (car trans)) + (end (cdr trans))) + (if (= end (1+ start)) + (list start) + (nconc (loop :for i :from start :below (1- end) :collect i) + (loop :for i :from (1- end) :downto start :collect i))))) + transpositions)) + +(defun apply-singleq-gate (state unitary q) + (apply-operator (lift unitary q (dimension-qubits (length state))) + state)) + +(defun apply-multiq-gate (state unitary qubits) + (let ((n (dimension-qubits (length state)))) + (let* ((unitary-init (lift unitary 0 n)) + (perm (append (reverse qubits) + (remove-if (lambda (i) (member i qubits)) + (loop for i :below n :collect i)))) + (trans (trans-as-adjacent (perm-as-trans perm)))) + (if trans + (let ((to->from (reduce #'compose-operator trans + :key (lambda (i) (lift +swap+ i n)))) + (from->to (reduce #'compose-operator (reverse trans) + :key (lambda (i) (lift +swap+ i n))))) + (apply-operator + (compose-operator to->from + (compose-operator unitary-init from->to)) + state)) + (apply-operator unitary-init state))))) + +(defun clq (qprog machine) + (loop :for (instruction . parameters) :in qprog + :do (ecase instruction + ((GATE) + (destructuring-bind (gate &rest qubits) parameters + (apply-gate (machine-quantum-state machine) gate qubits))) + ((MEASURE) + (observe machine))) + :finally (return machine))) + +(defparameter +X+ #2A((0 1) + (1 0)) + "Pauli X (NOT) gate") + +(defparameter +CNOT+ #2A((1 0 0 0) + (0 1 0 0) + (0 0 0 1) + (0 0 1 0)) + "Controlled NOT (CNOT) gate") + +(defparameter +H+ + (make-array '(2 2) + :initial-contents + (let ((s (/ (sqrt 2)))) + (list (list s s) + (list s (- s))))) + "Hadamard gate") + +(defun cphase (angle) + "Controlled phase (CPHASE) gate" + (make-array '(4 4) :initial-contents `((1 0 0 0) + (0 1 0 0) + (0 0 1 0) + (0 0 0 ,(cis angle))))) + +(defparameter +TOFFOLI+ + (let ((matrix (make-array '(8 8) + :element-type 'double-float + :initial-element 0.0d0))) + (dotimes (i 8) + (setf (aref matrix i i) 1d0)) + (setf (aref matrix 6 6) 0.0d0) + (setf (aref matrix 6 7) 1.0d0) + (setf (aref matrix 7 7) 0.0d0) + (setf (aref matrix 7 6) 1.0d0) + matrix) + "Controlled-controlled NOT (CCNOT) or Toffoli gate") + +(defun shor-fifteen () + `((GATE ,+H+ 0) + (GATE ,+H+ 1) + (GATE ,+H+ 2) + (GATE ,+CNOT+ 2 3) + (GATE ,+CNOT+ 2 4) + (GATE ,+H+ 1) + (GATE ,(cphase (/ pi 2)) 1 0) + (GATE ,+H+ 0) + (GATE ,(cphase (/ pi 4)) 1 2) + (GATE ,(cphase (/ pi 2)) 0 2) + (GATE ,+H+ 2) + (MEASURE) + )) + +(defun counts (realizations qubits q-circuit &rest args) + (let* ((hilbert (expt 2 qubits)) + (results (make-array hilbert :initial-element 0))) + (dotimes (i realizations results) + (let ((state (machine-quantum-state + (clq (apply q-circuit args) + (make-machine :quantum-state (make-quantum-state qubits) + :measurement-register 0))))) + (dotimes (j hilbert) + (when (> (aref state j) 0) (incf (aref results j)))))))) + +(defun counts-control (cnt n-control) + (let* ((bits-control (expt 2 n-control)) + (result (make-array bits-control :initial-element 0))) + (loop :for i :from 0 :below (length cnt) :by bits-control + :do + (loop :for j :from 0 :below bits-control + :do (incf (aref result j) (aref cnt (+ i j))))) + result)) + +(defun factorization (control-state N a n-control) + (let* ((counts (cdr control-state)) + (max-val (1+ (position (reduce 'max counts) counts))) + (r (/ (expt 2 n-control) max-val))) + (list (gcd (1+ (expt a (/ r 2))) N) + (gcd (1- (expt a (/ r 2))) N)))) + +(factorization (coerce (counts-control (counts 1024 5 #'shor-fifteen) 3) 'list) 15 11 3)