packages feed

newsynth (empty) → 0.1.0.0

raw patch · 23 files changed

+5759/−0 lines, 23 filesdep +basedep +fixedprecdep +newsynthbuild-type:Customsetup-changedbinary-added

Dependencies added: base, fixedprec, newsynth, random, superdoc, time

Files

+ ChangeLog view
@@ -0,0 +1,31 @@+ChangeLog++v0.1.0.0 2013/12/14+	(2013/12/13) PS1 - added new rings DComplex and QComplex.+	Improvements to the ring QOmega.+	(2013/12/13) PS1 - more uniform naming of rings. Old names+	DInteger, DReal, EReal, DComplex, EComplex have become ZRootTwo,+	DRootTwo, QRootTwo, DRComplex, and QRComplex, respectively.+	(2013/12/13) PS1 - adjusted output syntax to remove Unicode and+	make output parseable.+	(2013/12/13) PS1 - refactored as a Cabal package.+	(2013/12/11) PS1 - removed erroneous Adjoint2 instances for Double+	and Float.+	(2013/12/11) PS1 - removed dependency on numbers package, removed+	Random instance for FixedPrec (this is now in fixedprec package).+	(2013/12/11) PS1 - added EmptyDataDecls pragma, to keep GHC happy.+	(2013/12/10) PS1 - removed some unnecessary type class+	dependencies.+	(2013/09/25) PS1 - added "alternate" version of multi-qubit+	synthesis algorithm, using only generators of determinant 1 if+	possible.+	(2013/09/25) PS1 - renamed some constructors and deconstructors+	for matrices.++Release 2013/09/02+	Released as part of Quipper 0.5.+	(2013/07/08) PS1 - moved definitions of U2 and SO3 to Matrix.hs.+	(2013/07/05) PS1 - added RotationDecomposition module.++Release 2013/06/19+	Initial public release, as part of Quipper 0.4.
+ LICENSE view
@@ -0,0 +1,675 @@+              GNU GENERAL PUBLIC LICENSE+                Version 3, 29 June 2007++ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>+ Everyone is permitted to copy and distribute verbatim copies+ of this license document, but changing it is not allowed.++                     Preamble++  The GNU General Public License is a free, copyleft license for+software and other kinds of works.++  The licenses for most software and other practical works are designed+to take away your freedom to share and change the works.  By contrast,+the GNU General Public License is intended to guarantee your freedom to+share and change all versions of a program--to make sure it remains free+software for all its users.  We, the Free Software Foundation, use the+GNU General Public License for most of our software; it applies also to+any other work released this way by its authors.  You can apply it to+your programs, too.++  When we speak of free software, we are referring to freedom, not+price.  Our General Public Licenses are designed to make sure that you+have the freedom to distribute copies of free software (and charge for+them if you wish), that you receive source code or can get it if you+want it, that you can change the software or use pieces of it in new+free programs, and that you know you can do these things.++  To protect your rights, we need to prevent others from denying you+these rights or asking you to surrender the rights.  Therefore, you have+certain responsibilities if you distribute copies of the software, or if+you modify it: responsibilities to respect the freedom of others.++  For example, if you distribute copies of such a program, whether+gratis or for a fee, you must pass on to the recipients the same+freedoms that you received.  You must make sure that they, too, receive+or can get the source code.  And you must show them these terms so they+know their rights.++  Developers that use the GNU GPL protect your rights with two steps:+(1) assert copyright on the software, and (2) offer you this License+giving you legal permission to copy, distribute and/or modify it.++  For the developers' and authors' protection, the GPL clearly explains+that there is no warranty for this free software.  For both users' and+authors' sake, the GPL requires that modified versions be marked as+changed, so that their problems will not be attributed erroneously to+authors of previous versions.++  Some devices are designed to deny users access to install or run+modified versions of the software inside them, although the manufacturer+can do so.  This is fundamentally incompatible with the aim of+protecting users' freedom to change the software.  The systematic+pattern of such abuse occurs in the area of products for individuals to+use, which is precisely where it is most unacceptable.  Therefore, we+have designed this version of the GPL to prohibit the practice for those+products.  If such problems arise substantially in other domains, we+stand ready to extend this provision to those domains in future versions+of the GPL, as needed to protect the freedom of users.++  Finally, every program is threatened constantly by software patents.+States should not allow patents to restrict development and use of+software on general-purpose computers, but in those that do, we wish to+avoid the special danger that patents applied to a free program could+make it effectively proprietary.  To prevent this, the GPL assures that+patents cannot be used to render the program non-free.++  The precise terms and conditions for copying, distribution and+modification follow.++                TERMS AND CONDITIONS++  0. Definitions.++  "This License" refers to version 3 of the GNU General Public License.++  "Copyright" also means copyright-like laws that apply to other kinds of+works, such as semiconductor masks.+ +  "The Program" refers to any copyrightable work licensed under this+License.  Each licensee is addressed as "you".  "Licensees" and+"recipients" may be individuals or organizations.++  To "modify" a work means to copy from or adapt all or part of the work+in a fashion requiring copyright permission, other than the making of an+exact copy.  The resulting work is called a "modified version" of the+earlier work or a work "based on" the earlier work.++  A "covered work" means either the unmodified Program or a work based+on the Program.++  To "propagate" a work means to do anything with it that, without+permission, would make you directly or secondarily liable for+infringement under applicable copyright law, except executing it on a+computer or modifying a private copy.  Propagation includes copying,+distribution (with or without modification), making available to the+public, and in some countries other activities as well.++  To "convey" a work means any kind of propagation that enables other+parties to make or receive copies.  Mere interaction with a user through+a computer network, with no transfer of a copy, is not conveying.++  An interactive user interface displays "Appropriate Legal Notices"+to the extent that it includes a convenient and prominently visible+feature that (1) displays an appropriate copyright notice, and (2)+tells the user that there is no warranty for the work (except to the+extent that warranties are provided), that licensees may convey the+work under this License, and how to view a copy of this License.  If+the interface presents a list of user commands or options, such as a+menu, a prominent item in the list meets this criterion.++  1. Source Code.++  The "source code" for a work means the preferred form of the work+for making modifications to it.  "Object code" means any non-source+form of a work.++  A "Standard Interface" means an interface that either is an official+standard defined by a recognized standards body, or, in the case of+interfaces specified for a particular programming language, one that+is widely used among developers working in that language.++  The "System Libraries" of an executable work include anything, other+than the work as a whole, that (a) is included in the normal form of+packaging a Major Component, but which is not part of that Major+Component, and (b) serves only to enable use of the work with that+Major Component, or to implement a Standard Interface for which an+implementation is available to the public in source code form.  A+"Major Component", in this context, means a major essential component+(kernel, window system, and so on) of the specific operating system+(if any) on which the executable work runs, or a compiler used to+produce the work, or an object code interpreter used to run it.++  The "Corresponding Source" for a work in object code form means all+the source code needed to generate, install, and (for an executable+work) run the object code and to modify the work, including scripts to+control those activities.  However, it does not include the work's+System Libraries, or general-purpose tools or generally available free+programs which are used unmodified in performing those activities but+which are not part of the work.  For example, Corresponding Source+includes interface definition files associated with source files for+the work, and the source code for shared libraries and dynamically+linked subprograms that the work is specifically designed to require,+such as by intimate data communication or control flow between those+subprograms and other parts of the work.++  The Corresponding Source need not include anything that users+can regenerate automatically from other parts of the Corresponding+Source.++  The Corresponding Source for a work in source code form is that+same work.++  2. Basic Permissions.++  All rights granted under this License are granted for the term of+copyright on the Program, and are irrevocable provided the stated+conditions are met.  This License explicitly affirms your unlimited+permission to run the unmodified Program.  The output from running a+covered work is covered by this License only if the output, given its+content, constitutes a covered work.  This License acknowledges your+rights of fair use or other equivalent, as provided by copyright law.++  You may make, run and propagate covered works that you do not+convey, without conditions so long as your license otherwise remains+in force.  You may convey covered works to others for the sole purpose+of having them make modifications exclusively for you, or provide you+with facilities for running those works, provided that you comply with+the terms of this License in conveying all material for which you do+not control copyright.  Those thus making or running the covered works+for you must do so exclusively on your behalf, under your direction+and control, on terms that prohibit them from making any copies of+your copyrighted material outside their relationship with you.++  Conveying under any other circumstances is permitted solely under+the conditions stated below.  Sublicensing is not allowed; section 10+makes it unnecessary.++  3. Protecting Users' Legal Rights From Anti-Circumvention Law.++  No covered work shall be deemed part of an effective technological+measure under any applicable law fulfilling obligations under article+11 of the WIPO copyright treaty adopted on 20 December 1996, or+similar laws prohibiting or restricting circumvention of such+measures.++  When you convey a covered work, you waive any legal power to forbid+circumvention of technological measures to the extent such circumvention+is effected by exercising rights under this License with respect to+the covered work, and you disclaim any intention to limit operation or+modification of the work as a means of enforcing, against the work's+users, your or third parties' legal rights to forbid circumvention of+technological measures.++  4. Conveying Verbatim Copies.++  You may convey verbatim copies of the Program's source code as you+receive it, in any medium, provided that you conspicuously and+appropriately publish on each copy an appropriate copyright notice;+keep intact all notices stating that this License and any+non-permissive terms added in accord with section 7 apply to the code;+keep intact all notices of the absence of any warranty; and give all+recipients a copy of this License along with the Program.++  You may charge any price or no price for each copy that you convey,+and you may offer support or warranty protection for a fee.++  5. Conveying Modified Source Versions.++  You may convey a work based on the Program, or the modifications to+produce it from the Program, in the form of source code under the+terms of section 4, provided that you also meet all of these conditions:++    a) The work must carry prominent notices stating that you modified+    it, and giving a relevant date.++    b) The work must carry prominent notices stating that it is+    released under this License and any conditions added under section+    7.  This requirement modifies the requirement in section 4 to+    "keep intact all notices".++    c) You must license the entire work, as a whole, under this+    License to anyone who comes into possession of a copy.  This+    License will therefore apply, along with any applicable section 7+    additional terms, to the whole of the work, and all its parts,+    regardless of how they are packaged.  This License gives no+    permission to license the work in any other way, but it does not+    invalidate such permission if you have separately received it.++    d) If the work has interactive user interfaces, each must display+    Appropriate Legal Notices; however, if the Program has interactive+    interfaces that do not display Appropriate Legal Notices, your+    work need not make them do so.++  A compilation of a covered work with other separate and independent+works, which are not by their nature extensions of the covered work,+and which are not combined with it such as to form a larger program,+in or on a volume of a storage or distribution medium, is called an+"aggregate" if the compilation and its resulting copyright are not+used to limit the access or legal rights of the compilation's users+beyond what the individual works permit.  Inclusion of a covered work+in an aggregate does not cause this License to apply to the other+parts of the aggregate.++  6. Conveying Non-Source Forms.++  You may convey a covered work in object code form under the terms+of sections 4 and 5, provided that you also convey the+machine-readable Corresponding Source under the terms of this License,+in one of these ways:++    a) Convey the object code in, or embodied in, a physical product+    (including a physical distribution medium), accompanied by the+    Corresponding Source fixed on a durable physical medium+    customarily used for software interchange.++    b) Convey the object code in, or embodied in, a physical product+    (including a physical distribution medium), accompanied by a+    written offer, valid for at least three years and valid for as+    long as you offer spare parts or customer support for that product+    model, to give anyone who possesses the object code either (1) a+    copy of the Corresponding Source for all the software in the+    product that is covered by this License, on a durable physical+    medium customarily used for software interchange, for a price no+    more than your reasonable cost of physically performing this+    conveying of source, or (2) access to copy the+    Corresponding Source from a network server at no charge.++    c) Convey individual copies of the object code with a copy of the+    written offer to provide the Corresponding Source.  This+    alternative is allowed only occasionally and noncommercially, and+    only if you received the object code with such an offer, in accord+    with subsection 6b.++    d) Convey the object code by offering access from a designated+    place (gratis or for a charge), and offer equivalent access to the+    Corresponding Source in the same way through the same place at no+    further charge.  You need not require recipients to copy the+    Corresponding Source along with the object code.  If the place to+    copy the object code is a network server, the Corresponding Source+    may be on a different server (operated by you or a third party)+    that supports equivalent copying facilities, provided you maintain+    clear directions next to the object code saying where to find the+    Corresponding Source.  Regardless of what server hosts the+    Corresponding Source, you remain obligated to ensure that it is+    available for as long as needed to satisfy these requirements.++    e) Convey the object code using peer-to-peer transmission, provided+    you inform other peers where the object code and Corresponding+    Source of the work are being offered to the general public at no+    charge under subsection 6d.++  A separable portion of the object code, whose source code is excluded+from the Corresponding Source as a System Library, need not be+included in conveying the object code work.++  A "User Product" is either (1) a "consumer product", which means any+tangible personal property which is normally used for personal, family,+or household purposes, or (2) anything designed or sold for incorporation+into a dwelling.  In determining whether a product is a consumer product,+doubtful cases shall be resolved in favor of coverage.  For a particular+product received by a particular user, "normally used" refers to a+typical or common use of that class of product, regardless of the status+of the particular user or of the way in which the particular user+actually uses, or expects or is expected to use, the product.  A product+is a consumer product regardless of whether the product has substantial+commercial, industrial or non-consumer uses, unless such uses represent+the only significant mode of use of the product.++  "Installation Information" for a User Product means any methods,+procedures, authorization keys, or other information required to install+and execute modified versions of a covered work in that User Product from+a modified version of its Corresponding Source.  The information must+suffice to ensure that the continued functioning of the modified object+code is in no case prevented or interfered with solely because+modification has been made.++  If you convey an object code work under this section in, or with, or+specifically for use in, a User Product, and the conveying occurs as+part of a transaction in which the right of possession and use of the+User Product is transferred to the recipient in perpetuity or for a+fixed term (regardless of how the transaction is characterized), the+Corresponding Source conveyed under this section must be accompanied+by the Installation Information.  But this requirement does not apply+if neither you nor any third party retains the ability to install+modified object code on the User Product (for example, the work has+been installed in ROM).++  The requirement to provide Installation Information does not include a+requirement to continue to provide support service, warranty, or updates+for a work that has been modified or installed by the recipient, or for+the User Product in which it has been modified or installed.  Access to a+network may be denied when the modification itself materially and+adversely affects the operation of the network or violates the rules and+protocols for communication across the network.++  Corresponding Source conveyed, and Installation Information provided,+in accord with this section must be in a format that is publicly+documented (and with an implementation available to the public in+source code form), and must require no special password or key for+unpacking, reading or copying.++  7. Additional Terms.++  "Additional permissions" are terms that supplement the terms of this+License by making exceptions from one or more of its conditions.+Additional permissions that are applicable to the entire Program shall+be treated as though they were included in this License, to the extent+that they are valid under applicable law.  If additional permissions+apply only to part of the Program, that part may be used separately+under those permissions, but the entire Program remains governed by+this License without regard to the additional permissions.++  When you convey a copy of a covered work, you may at your option+remove any additional permissions from that copy, or from any part of+it.  (Additional permissions may be written to require their own+removal in certain cases when you modify the work.)  You may place+additional permissions on material, added by you to a covered work,+for which you have or can give appropriate copyright permission.++  Notwithstanding any other provision of this License, for material you+add to a covered work, you may (if authorized by the copyright holders of+that material) supplement the terms of this License with terms:++    a) Disclaiming warranty or limiting liability differently from the+    terms of sections 15 and 16 of this License; or++    b) Requiring preservation of specified reasonable legal notices or+    author attributions in that material or in the Appropriate Legal+    Notices displayed by works containing it; or++    c) Prohibiting misrepresentation of the origin of that material, or+    requiring that modified versions of such material be marked in+    reasonable ways as different from the original version; or++    d) Limiting the use for publicity purposes of names of licensors or+    authors of the material; or++    e) Declining to grant rights under trademark law for use of some+    trade names, trademarks, or service marks; or++    f) Requiring indemnification of licensors and authors of that+    material by anyone who conveys the material (or modified versions of+    it) with contractual assumptions of liability to the recipient, for+    any liability that these contractual assumptions directly impose on+    those licensors and authors.++  All other non-permissive additional terms are considered "further+restrictions" within the meaning of section 10.  If the Program as you+received it, or any part of it, contains a notice stating that it is+governed by this License along with a term that is a further+restriction, you may remove that term.  If a license document contains+a further restriction but permits relicensing or conveying under this+License, you may add to a covered work material governed by the terms+of that license document, provided that the further restriction does+not survive such relicensing or conveying.++  If you add terms to a covered work in accord with this section, you+must place, in the relevant source files, a statement of the+additional terms that apply to those files, or a notice indicating+where to find the applicable terms.++  Additional terms, permissive or non-permissive, may be stated in the+form of a separately written license, or stated as exceptions;+the above requirements apply either way.++  8. Termination.++  You may not propagate or modify a covered work except as expressly+provided under this License.  Any attempt otherwise to propagate or+modify it is void, and will automatically terminate your rights under+this License (including any patent licenses granted under the third+paragraph of section 11).++  However, if you cease all violation of this License, then your+license from a particular copyright holder is reinstated (a)+provisionally, unless and until the copyright holder explicitly and+finally terminates your license, and (b) permanently, if the copyright+holder fails to notify you of the violation by some reasonable means+prior to 60 days after the cessation.++  Moreover, your license from a particular copyright holder is+reinstated permanently if the copyright holder notifies you of the+violation by some reasonable means, this is the first time you have+received notice of violation of this License (for any work) from that+copyright holder, and you cure the violation prior to 30 days after+your receipt of the notice.++  Termination of your rights under this section does not terminate the+licenses of parties who have received copies or rights from you under+this License.  If your rights have been terminated and not permanently+reinstated, you do not qualify to receive new licenses for the same+material under section 10.++  9. Acceptance Not Required for Having Copies.++  You are not required to accept this License in order to receive or+run a copy of the Program.  Ancillary propagation of a covered work+occurring solely as a consequence of using peer-to-peer transmission+to receive a copy likewise does not require acceptance.  However,+nothing other than this License grants you permission to propagate or+modify any covered work.  These actions infringe copyright if you do+not accept this License.  Therefore, by modifying or propagating a+covered work, you indicate your acceptance of this License to do so.++  10. Automatic Licensing of Downstream Recipients.++  Each time you convey a covered work, the recipient automatically+receives a license from the original licensors, to run, modify and+propagate that work, subject to this License.  You are not responsible+for enforcing compliance by third parties with this License.++  An "entity transaction" is a transaction transferring control of an+organization, or substantially all assets of one, or subdividing an+organization, or merging organizations.  If propagation of a covered+work results from an entity transaction, each party to that+transaction who receives a copy of the work also receives whatever+licenses to the work the party's predecessor in interest had or could+give under the previous paragraph, plus a right to possession of the+Corresponding Source of the work from the predecessor in interest, if+the predecessor has it or can get it with reasonable efforts.++  You may not impose any further restrictions on the exercise of the+rights granted or affirmed under this License.  For example, you may+not impose a license fee, royalty, or other charge for exercise of+rights granted under this License, and you may not initiate litigation+(including a cross-claim or counterclaim in a lawsuit) alleging that+any patent claim is infringed by making, using, selling, offering for+sale, or importing the Program or any portion of it.++  11. Patents.++  A "contributor" is a copyright holder who authorizes use under this+License of the Program or a work on which the Program is based.  The+work thus licensed is called the contributor's "contributor version".++  A contributor's "essential patent claims" are all patent claims+owned or controlled by the contributor, whether already acquired or+hereafter acquired, that would be infringed by some manner, permitted+by this License, of making, using, or selling its contributor version,+but do not include claims that would be infringed only as a+consequence of further modification of the contributor version.  For+purposes of this definition, "control" includes the right to grant+patent sublicenses in a manner consistent with the requirements of+this License.++  Each contributor grants you a non-exclusive, worldwide, royalty-free+patent license under the contributor's essential patent claims, to+make, use, sell, offer for sale, import and otherwise run, modify and+propagate the contents of its contributor version.++  In the following three paragraphs, a "patent license" is any express+agreement or commitment, however denominated, not to enforce a patent+(such as an express permission to practice a patent or covenant not to+sue for patent infringement).  To "grant" such a patent license to a+party means to make such an agreement or commitment not to enforce a+patent against the party.++  If you convey a covered work, knowingly relying on a patent license,+and the Corresponding Source of the work is not available for anyone+to copy, free of charge and under the terms of this License, through a+publicly available network server or other readily accessible means,+then you must either (1) cause the Corresponding Source to be so+available, or (2) arrange to deprive yourself of the benefit of the+patent license for this particular work, or (3) arrange, in a manner+consistent with the requirements of this License, to extend the patent+license to downstream recipients.  "Knowingly relying" means you have+actual knowledge that, but for the patent license, your conveying the+covered work in a country, or your recipient's use of the covered work+in a country, would infringe one or more identifiable patents in that+country that you have reason to believe are valid.+  +  If, pursuant to or in connection with a single transaction or+arrangement, you convey, or propagate by procuring conveyance of, a+covered work, and grant a patent license to some of the parties+receiving the covered work authorizing them to use, propagate, modify+or convey a specific copy of the covered work, then the patent license+you grant is automatically extended to all recipients of the covered+work and works based on it.++  A patent license is "discriminatory" if it does not include within+the scope of its coverage, prohibits the exercise of, or is+conditioned on the non-exercise of one or more of the rights that are+specifically granted under this License.  You may not convey a covered+work if you are a party to an arrangement with a third party that is+in the business of distributing software, under which you make payment+to the third party based on the extent of your activity of conveying+the work, and under which the third party grants, to any of the+parties who would receive the covered work from you, a discriminatory+patent license (a) in connection with copies of the covered work+conveyed by you (or copies made from those copies), or (b) primarily+for and in connection with specific products or compilations that+contain the covered work, unless you entered into that arrangement,+or that patent license was granted, prior to 28 March 2007.++  Nothing in this License shall be construed as excluding or limiting+any implied license or other defenses to infringement that may+otherwise be available to you under applicable patent law.++  12. No Surrender of Others' Freedom.++  If conditions are imposed on you (whether by court order, agreement or+otherwise) that contradict the conditions of this License, they do not+excuse you from the conditions of this License.  If you cannot convey a+covered work so as to satisfy simultaneously your obligations under this+License and any other pertinent obligations, then as a consequence you may+not convey it at all.  For example, if you agree to terms that obligate you+to collect a royalty for further conveying from those to whom you convey+the Program, the only way you could satisfy both those terms and this+License would be to refrain entirely from conveying the Program.++  13. Use with the GNU Affero General Public License.++  Notwithstanding any other provision of this License, you have+permission to link or combine any covered work with a work licensed+under version 3 of the GNU Affero General Public License into a single+combined work, and to convey the resulting work.  The terms of this+License will continue to apply to the part which is the covered work,+but the special requirements of the GNU Affero General Public License,+section 13, concerning interaction through a network will apply to the+combination as such.++  14. Revised Versions of this License.++  The Free Software Foundation may publish revised and/or new versions of+the GNU General Public License from time to time.  Such new versions will+be similar in spirit to the present version, but may differ in detail to+address new problems or concerns.++  Each version is given a distinguishing version number.  If the+Program specifies that a certain numbered version of the GNU General+Public License "or any later version" applies to it, you have the+option of following the terms and conditions either of that numbered+version or of any later version published by the Free Software+Foundation.  If the Program does not specify a version number of the+GNU General Public License, you may choose any version ever published+by the Free Software Foundation.++  If the Program specifies that a proxy can decide which future+versions of the GNU General Public License can be used, that proxy's+public statement of acceptance of a version permanently authorizes you+to choose that version for the Program.++  Later license versions may give you additional or different+permissions.  However, no additional obligations are imposed on any+author or copyright holder as a result of your choosing to follow a+later version.++  15. Disclaimer of Warranty.++  THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY+APPLICABLE LAW.  EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR+PURPOSE.  THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM+IS WITH YOU.  SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.++  16. Limitation of Liability.++  IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF+SUCH DAMAGES.++  17. Interpretation of Sections 15 and 16.++  If the disclaimer of warranty and limitation of liability provided+above cannot be given local legal effect according to their terms,+reviewing courts shall apply local law that most closely approximates+an absolute waiver of all civil liability in connection with the+Program, unless a warranty or assumption of liability accompanies a+copy of the Program in return for a fee.++              END OF TERMS AND CONDITIONS++     How to Apply These Terms to Your New Programs++  If you develop a new program, and you want it to be of the greatest+possible use to the public, the best way to achieve this is to make it+free software which everyone can redistribute and change under these terms.++  To do so, attach the following notices to the program.  It is safest+to attach them to the start of each source file to most effectively+state the exclusion of warranty; and each file should have at least+the "copyright" line and a pointer to where the full notice is found.++    <one line to give the program's name and a brief idea of what it does.>+    Copyright (C) <year>  <name of author>++    This program is free software: you can redistribute it and/or modify+    it under the terms of the GNU General Public License as published by+    the Free Software Foundation, either version 3 of the License, or+    (at your option) any later version.++    This program is distributed in the hope that it will be useful,+    but WITHOUT ANY WARRANTY; without even the implied warranty of+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the+    GNU General Public License for more details.++    You should have received a copy of the GNU General Public License+    along with this program.  If not, see <http://www.gnu.org/licenses/>.++Also add information on how to contact you by electronic and paper mail.++  If the program does terminal interaction, make it output a short+notice like this when it starts in an interactive mode:++    <program>  Copyright (C) <year>  <name of author>+    This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.+    This is free software, and you are welcome to redistribute it+    under certain conditions; type `show c' for details.++The hypothetical commands `show w' and `show c' should show the appropriate+parts of the General Public License.  Of course, your program's commands+might be different; for a GUI interface, you would use an "about box".++  You should also get your employer (if you work as a programmer) or school,+if any, to sign a "copyright disclaimer" for the program, if necessary.+For more information on this, and how to apply and follow the GNU GPL, see+<http://www.gnu.org/licenses/>.++  The GNU General Public License does not permit incorporating your program+into proprietary programs.  If your program is a subroutine library, you+may consider it more useful to permit linking proprietary applications with+the library.  If this is what you want to do, use the GNU Lesser General+Public License instead of this License.  But first, please read+<http://www.gnu.org/philosophy/why-not-lgpl.html>.+
+ Quantum/Synthesis/ArcTan2.hs view
@@ -0,0 +1,36 @@+-- | This module provides a replacement for Haskell's 'atan2'. The+-- problem is that Haskell's standard implementation of 'atan2'+-- depends on the 'RealFloat' class, which limits its applicability.+-- So we provide a new 'ArcTan2' class with an 'arctan2' function.+-- +-- Unlike Haskell's 'atan2', the 'arctan2' function may not take+-- signed zeros and signed infinities into account. But it works at+-- fixed-precision types such as 'FixedPrec'.++module Quantum.Synthesis.ArcTan2 where++import Data.Number.FixedPrec++-- ----------------------------------------------------------------------+-- * The arctan2 function++-- | We provide a replacement for Haskell's 'atan2', because the+-- latter depends on the 'RealFloat' class, which limits its+-- applicability.+class ArcTan2 a where+  arctan2 :: a -> a -> a+  +instance ArcTan2 Double where+  arctan2 = atan2++instance ArcTan2 Float where+  arctan2 = atan2++instance (Precision e) => ArcTan2 (FixedPrec e) where+  arctan2 y x+    | x == 0 && y == 0 = 0+    | abs y <= x       = atan (y/x)+    | abs x <= y       = pi/2 - atan (x/y)+    | abs x <= -y      = -pi/2 - atan (x/y)+    | y >= 0           = pi + atan (y/x)+    | otherwise        = -pi + atan (y/x)
+ Quantum/Synthesis/Clifford.hs view
@@ -0,0 +1,304 @@+{-# LANGUAGE OverlappingInstances #-}++{-# OPTIONS_GHC -fno-warn-incomplete-patterns #-}++-- | This module provides an efficient symbolic representation of the+-- Clifford group on one qubit. This group is generated by /S/, /H/,+-- and the scalar ω = [exp /i/π\/4]. It has 192 elements. ++module Quantum.Synthesis.Clifford (+  -- * The Clifford group+  Clifford,+  +  -- ** Constructors+  clifford_X,+  clifford_Y,+  clifford_Z,+  clifford_H,+  clifford_S,+  clifford_SH,+  clifford_E,+  clifford_W,+  ToClifford(to_clifford),+  +  -- ** Deconstructors+  clifford_decompose,+  Axis(..),+  clifford_decompose_coset,+  +  -- ** Group operations+  clifford_id,+  clifford_mult,+  clifford_inv,+  +  -- ** Conjugation by /T/+  clifford_tconj  +  ) where++-- ----------------------------------------------------------------------+-- * The Clifford group++-- $ We could, in principle, implement the Clifford group as an+-- enumerated type with 192 elements, and a large 192×192 lookup+-- table for the group multiplication. Instead, we take advantage of+-- some of the internal structure of the group to reduce the size of+-- the lookup tables. The resulting implementation is still very+-- efficient.++-- | A type representing single-qubit Clifford operators.+data Clifford = Clifford Int Int Int Int+                deriving (Eq, Ord)++instance Show Clifford where+  show (Clifford a b c d) = "C" ++ show a ++ show b ++ show c ++ show d++-- ----------------------------------------------------------------------+-- ** Constructors++-- | The Pauli /X/-gate as a Clifford operator.+clifford_X :: Clifford+clifford_X = Clifford 0 1 0 0++-- | The Pauli /Y/-gate as a Clifford operator.+clifford_Y :: Clifford+clifford_Y = Clifford 0 1 2 2++-- | The Pauli /Z/-gate as a Clifford operator.+clifford_Z :: Clifford+clifford_Z = Clifford 0 0 2 0++-- | The Hadamard gate as a Clifford operator.+clifford_H :: Clifford+clifford_H = Clifford 1 0 1 5++-- | The Clifford operator /S/.+clifford_S :: Clifford+clifford_S = Clifford 0 0 1 0++-- | The Clifford operator /SH/.+clifford_SH :: Clifford+clifford_SH = clifford_S `clifford_mult` clifford_H++-- | The Clifford operator /E/ = /H//S/[sup 3]ω[sup 3]. This operator is+-- uniquely determined by the properties /E/³ = /I/, +-- /EXE/⁻¹ = /Y/, /EYE/⁻¹ = /Z/, and /EZE/⁻¹ = /X/.+-- +-- \[image E.png]+clifford_E :: Clifford+clifford_E = Clifford 1 0 0 0++-- | The Clifford operator ω = [exp /i/π\/4].+clifford_W :: Clifford+clifford_W = Clifford 0 0 0 1++-- | A type class for things that can be exactly converted to a+-- Clifford operator. One particular instance of this is 'String', so+-- that Clifford operators can be denoted, e.g.,+-- +-- > to_clifford "-iX"+-- +-- The valid characters for such string conversions are @\"XYZHSEIWi-\"@.+class ToClifford a where+  -- | Convert any suitable thing to a Clifford operator.+  to_clifford :: a -> Clifford+  +instance ToClifford Clifford where+  to_clifford = id+  +instance ToClifford Char where+  to_clifford 'E' = clifford_E+  to_clifford 'X' = clifford_X+  to_clifford 'S' = clifford_S+  to_clifford 'W' = clifford_W+  to_clifford 'I' = clifford_id+  to_clifford 'i' = Clifford 0 0 0 2+  to_clifford '-' = Clifford 0 0 0 4+  to_clifford 'H' = clifford_H+  to_clifford 'Y' = clifford_Y+  to_clifford 'Z' = clifford_Z+  to_clifford x = error $ "ToClifford Char: unknown gate " ++ show x++instance ToClifford a => ToClifford [a] where+  to_clifford [] = clifford_id+  to_clifford (h:t) = to_clifford h `clifford_mult` to_clifford t++-- ----------------------------------------------------------------------+-- ** Deconstructors++-- | Given a Clifford operator /U/, return (/a/, /b/, /c/, /d/) such that+-- +-- * /U/ = /E/[sup /a/]/X/[sup /b/]/S/[sup /c/]ω[sup /d/],+-- +-- * /a/ ∈ {0, 1, 2}, /b/ ∈ {0, 1}, /c/ ∈ {0, …, 3}, and /d/ ∈ {0, …,+-- 7}.+-- +-- Here, /E/ = /H//S/[sup 3]ω[sup 3]. Note that /E/, /X/, /S/, and ω have order+-- 3, 2, 4, and 8, respectively. Moreover, each Clifford operator can+-- be uniquely represented as above.+clifford_decompose :: (ToClifford a) => a -> (Int, Int, Int, Int)+clifford_decompose m = (a,b,c,d) where+  Clifford a b c d = to_clifford m++-- | A axis is either /I/, /H/, or /SH/.+data Axis = Axis_I | Axis_H | Axis_SH+           deriving (Eq, Show)++instance ToClifford Axis where+  to_clifford Axis_I = to_clifford "I"+  to_clifford Axis_H = to_clifford "H"+  to_clifford Axis_SH = to_clifford "SH"++-- | Given a Clifford operator /U/, return (/K/, /b/, /c/, /d/) such that+-- +-- * /U/ = /K//X/[sup /b/]/S/[sup /c/]ω[sup /d/],+-- +-- * /K/ ∈ {/I/, /H/, /SH/}, /b/ ∈ {0, 1}, /c/ ∈ {0, …, 3}, and /d/ ∈ {0, …,+-- 7}.+clifford_decompose_coset :: (ToClifford a) => a -> (Axis, Int, Int, Int)+clifford_decompose_coset u = case op of+  Clifford 0 b c d -> (Axis_I, b, c, d)+  Clifford 1 b c d -> (Axis_H, b', c', d') where+    Clifford 0 b' c' d' = clifford_inv "H" `clifford_mult` op+  Clifford 2 b c d -> (Axis_SH, b', c', d') where+    Clifford 0 b' c' d' = clifford_inv "SH" `clifford_mult` op+  where+    op = to_clifford u+  +-- ----------------------------------------------------------------------+-- ** Group operations++-- | The identity Clifford operator.+clifford_id :: Clifford+clifford_id = Clifford 0 0 0 0++-- | Clifford multiplication.+clifford_mult :: Clifford -> Clifford -> Clifford+clifford_mult u1 u2 = u where+  -- U = U1 U2   +  --   = A1 B1 C1 D1 A2 B2 C2 D2+  --   = A1 (B1 C1 A2) B2 C2 D1 D2+  --   = A1 (A3 B3 C3 D3) B2 C2 D1 D2+  --   = A1 A3 B3 (C3 B2) C2 D3 D1 D2+  --   = A1 A3 B3 (B2 C4 D4) C2 D3 D1 D2+  --   = (A1 A3) (B3 B2) (C4 C2) (D4 D3 D1 D2)+  --   = A B C D+  Clifford a1 b1 c1 d1 = u1+  Clifford a2 b2 c2 d2 = u2+  (a3, b3, c3, d3) = conj3 b1 c1 a2+  (c4, d4) = conj2 c3 b2+  a = (a1 + a3) `mod` 3+  b = (b3 + b2) `mod` 2+  c = (c4 + c2) `mod` 4+  d = (d4 + d3 + d1 + d2) `mod` 8+  u = Clifford a b c d++-- | Clifford inverse.+clifford_inv :: (ToClifford a) => a -> Clifford+clifford_inv op = Clifford a2 b2 c2 d3 where+  -- U⁻¹ = (A B C)⁻¹ D⁻¹ = (A2 B2 C2 D2) D⁻¹+  Clifford a b c d = to_clifford op+  (a2, b2, c2, d2) = cinv a b c+  d3 = (d2 - d) `mod` 8++-- ----------------------------------------------------------------------+-- ** Conjugation by /T/++-- | Given a Clifford gate /C/, return an axis /K/ ∈ {/I/, /H/, /SH/}+-- and a Clifford gate /C'/ such that+-- +-- * /C//T/ = /K//T//C/'.+clifford_tconj ::  Clifford -> (Axis, Clifford)+clifford_tconj u = (k, v) where+  -- U T = A1 B1 C1 D1 T+  --     = (A1 B1 T) C1 D1+  --     = (K T B1 C2 D2) C1 D1+  --     = K T B1 (C2 C1) (D2 D1)+  Clifford a1 b1 c1 d1 = u+  (k, c2, d2) = tconj a1 b1+  c = (c2 + c1) `mod` 4+  d = (d2 + d1) `mod` 8+  v = Clifford 0 b1 c d++-- ----------------------------------------------------------------------+-- ** Lookup tables++-- | 'conj2' /c/ /b/ returns (/c/', /d/') such that+-- +-- * /S/[sup /c/]/X/[sup /b/] = /X/[sup /b/]/S/[sup /c/']ω[sup /d/'].+conj2 :: Int -> Int -> (Int, Int)+conj2 0 0 = (0,0)+conj2 0 1 = (0,0)+conj2 1 0 = (1,0)+conj2 1 1 = (3,2)+conj2 2 0 = (2,0)+conj2 2 1 = (2,4)+conj2 3 0 = (3,0)+conj2 3 1 = (1,6)++-- | 'conj3' /b/ /c/ /a/ returns (/a/', /b/', /c/', /d/') such that+-- +-- * /X/[sup /b/]/S/[sup /c/]/E/[sup /a/] = /E/[sup /a/']/X/[sup /b/']/S/[sup /c/']ω[sup /d/'].+conj3 :: Int -> Int -> Int -> (Int, Int, Int, Int)+conj3 0 0 0 = (0,0,0,0)+conj3 0 0 1 = (1,0,0,0)+conj3 0 0 2 = (2,0,0,0)+conj3 0 1 0 = (0,0,1,0)+conj3 0 1 1 = (2,0,3,6)+conj3 0 1 2 = (1,1,3,4)+conj3 0 2 0 = (0,0,2,0)+conj3 0 2 1 = (1,1,2,2)+conj3 0 2 2 = (2,1,0,0)+conj3 0 3 0 = (0,0,3,0)+conj3 0 3 1 = (2,1,3,6)+conj3 0 3 2 = (1,0,1,2)+conj3 1 0 0 = (0,1,0,0)+conj3 1 0 1 = (1,0,2,0)+conj3 1 0 2 = (2,1,2,2)+conj3 1 1 0 = (0,1,1,0)+conj3 1 1 1 = (2,1,1,0)+conj3 1 1 2 = (1,1,1,0)+conj3 1 2 0 = (0,1,2,0)+conj3 1 2 1 = (1,1,0,6)+conj3 1 2 2 = (2,0,2,6)+conj3 1 3 0 = (0,1,3,0)+conj3 1 3 1 = (2,0,1,4)+conj3 1 3 2 = (1,0,3,2)++-- | 'cinv' /a/ /b/ /c/ returns (/a/', /b/', /c/', /d/') such that+-- +-- * (/E/[sup /a/]/X/[sup /b/]/S/[sup /c/])⁻¹ = /E/[sup /a/']/X/[sup /b/']/S/[sup /c/']ω[sup /d/'].+cinv :: Int -> Int -> Int -> (Int, Int, Int, Int)+cinv 0 0 0 = (0,0,0,0)+cinv 0 0 1 = (0,0,3,0)+cinv 0 0 2 = (0,0,2,0)+cinv 0 0 3 = (0,0,1,0)+cinv 0 1 0 = (0,1,0,0)+cinv 0 1 1 = (0,1,1,6)+cinv 0 1 2 = (0,1,2,4)+cinv 0 1 3 = (0,1,3,2)+cinv 1 0 0 = (2,0,0,0)+cinv 1 0 1 = (1,0,1,2)+cinv 1 0 2 = (2,1,0,0)+cinv 1 0 3 = (1,1,3,4)+cinv 1 1 0 = (2,1,2,2)+cinv 1 1 1 = (1,1,1,6)+cinv 1 1 2 = (2,0,2,2)+cinv 1 1 3 = (1,0,3,4)+cinv 2 0 0 = (1,0,0,0)+cinv 2 0 1 = (2,1,3,6)+cinv 2 0 2 = (1,1,2,2)+cinv 2 0 3 = (2,0,3,6)+cinv 2 1 0 = (1,0,2,0)+cinv 2 1 1 = (2,1,1,6)+cinv 2 1 2 = (1,1,0,2)+cinv 2 1 3 = (2,0,1,6)++-- | 'tconj2' /a/ /b/ returns (/K/, /c/, /d/) such that+-- +-- * /E/[sup /a/]/X/[sup /b/]/T/ = /K//T//X/[sup /b/]/S/[sup /c/]ω[sup /d/].+tconj 0 0 = (Axis_I,  0, 0)+tconj 0 1 = (Axis_I,  1, 7)+tconj 1 0 = (Axis_H,  3, 3)+tconj 1 1 = (Axis_H,  2, 0)+tconj 2 0 = (Axis_SH, 0, 5)+tconj 2 1 = (Axis_SH, 1, 4)
+ Quantum/Synthesis/CliffordT.hs view
@@ -0,0 +1,606 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE FlexibleContexts #-}++-- | This module provides a representation of the single-qubit+-- Clifford+/T/ operators, Matsumoto-Amano normal forms, and functions+-- for the exact synthesis of single-qubit Clifford+/T/ operators.+--+-- Matsumoto-Amano normal forms and the Matsumoto-Amano exact+-- synthesis algorithm are described in the paper:+--+-- * Ken Matsumoto, Kazuyuki Amano. Representation of Quantum Circuits+-- with Clifford and π\/8 Gates. <http://arxiv.org/abs/0806.3834>.++module Quantum.Synthesis.CliffordT where++import Quantum.Synthesis.Ring+import Quantum.Synthesis.Matrix+import Quantum.Synthesis.Clifford+import Quantum.Synthesis.MultiQubitSynthesis++import Data.List+import Data.Bits++-- ----------------------------------------------------------------------+-- * Clifford+/T/ interchange format++-- $ It is convenient to have a simple but exact \"interchange+-- format\" for operators in the single-qubit Clifford+/T/+-- group. Different operator representations can be converted to and+-- from this format.+--+-- Our format is simply a list of gates from /X/, /Y/, /Z/, /H/, /S/,+-- /T/, and /E/ = /H//S/[sup 3]ω[sup 3], with the obvious+-- interpretation as a matrix product. We also include the global+-- phase gate /W/ = ω = [exp /i/π\/4]. The /W/ gate is ignored when+-- converting to or from representations that cannot represent global+-- phase (such as the Bloch sphere representation).++-- | An enumeration type to represent symbolic basic gates (/X/, /Y/,+-- /Z/, /H/, /S/, /T/, /W/, /E/).+-- +-- Note: when we use a list of 'Gate's to express a sequence of+-- operators, the operators are meant to be applied right-to-left,+-- i.e., as in the mathematical notation for matrix multiplication.+-- This is the opposite of the quantum circuit notation.+data Gate = X | Y | Z | H | S | T | E | W+          deriving (Show, Eq)++-- | A type class for all things that can be exactly converted to a+-- list of gates. These are the exact representations of the+-- single-qubit Clifford+/T/ group.+class ToGates a where+  -- | Convert any suitable thing to a list of gates.+  to_gates :: a -> [Gate]++instance ToGates Gate where+  to_gates x = [x]++instance (ToGates a) => ToGates [a] where+  to_gates x = concat [ to_gates y | y <- x ]++instance ToGates Char where+  to_gates 'X' = [X]+  to_gates 'Y' = [Y]+  to_gates 'Z' = [Z]+  to_gates 'H' = [H]+  to_gates 'S' = [S]+  to_gates 'T' = [T]+  to_gates 'E' = [E]+  to_gates 'W' = [W]+  to_gates 'I' = []+  to_gates '-' = [W,W,W,W]+  to_gates 'i' = [W,W]+  to_gates x = error $ "to_gates[Char]: undefined -- " ++ (show x)++instance ToGates Axis where+  to_gates Axis_I = []+  to_gates Axis_H = [H]+  to_gates Axis_SH = [S,H]++instance ToGates Clifford where+  to_gates op = as ++ xs ++ ss ++ ws where+    (k, b, c, d) = clifford_decompose_coset op+    as = to_gates k+    xs  = replicate b X+    ss  = replicate c S+    ws  = replicate d W++-- | A type class for all things that a list of gates can be converted+-- to. For example, a list of gates can be converted to an element of+-- /U/(2) or an element of /SO/(3), using various (exact or+-- approximate) representations of the matrix entries.+class FromGates a where+  -- | Convert a list of gates to any suitable type.+  from_gates :: [Gate] -> a++instance FromGates String where+  from_gates = concat . map show++instance FromGates [Gate] where+  from_gates = id++-- | Invert a gate list.+invert_gates :: [Gate] -> [Gate]+invert_gates gs = aux [] gs where+  aux acc [] = acc+  aux acc (X:t) = aux (X:acc) t+  aux acc (Y:t) = aux (Y:acc) t+  aux acc (Z:t) = aux (Z:acc) t+  aux acc (H:t) = aux (H:acc) t+  aux acc (S:t) = aux (Z:S:acc) t+  aux acc (T:t) = aux (Z:S:T:acc) t+  aux acc (E:t) = aux (E:E:acc) t+  aux acc (W:t) = aux (W:W:W:W:W:W:W:acc) t++-- | Convert any precise format to any format.+convert :: (ToGates a, FromGates b) => a -> b+convert = from_gates . to_gates++-- ----------------------------------------------------------------------+-- * Matrices in /U/(2) and /SO/(3)++-- ----------------------------------------------------------------------+-- ** Matrices in /U/(2)++-- | The Pauli /X/ operator.+u2_X :: (Ring a) => U2 a+u2_X = matrix2x2 (0, 1)+                 (1, 0)++-- | The Pauli /Y/ operator.+u2_Y :: (ComplexRing a) => U2 a+u2_Y = matrix2x2 (0, -i)+                 (i,  0)++-- | The Pauli /Z/ operator.+u2_Z :: (Ring a) => U2 a+u2_Z = matrix2x2 (1,  0)+                 (0, -1)++-- | The Hadamard operator.+u2_H :: (RootHalfRing a) => U2 a+u2_H = roothalf * matrix2x2 (1,  1)+                            (1, -1)++-- | The /S/ operator.+u2_S :: (ComplexRing a) => U2 a+u2_S = matrix2x2 (1, 0)+                 (0, i)++-- | The /T/ operator.+u2_T :: (OmegaRing a) => U2 a+u2_T = matrix2x2 (1,     0)+                 (0, omega)++-- | The /E/ operator.+u2_E :: (OmegaRing a, RootHalfRing a) => U2 a+u2_E = roothalf * matrix2x2 (omega^3,  omega)+                            (omega^3, -omega)++-- | The /W/ = [exp /i/π\/4] global phase operator.+u2_W :: (OmegaRing a) => U2 a+u2_W = matrix2x2 (omega,     0)+                 (0,     omega)++-- | Convert a symbolic gate to the corresponding operator.+u2_of_gate :: (RootHalfRing a, ComplexRing a) => Gate -> U2 a+u2_of_gate X = u2_X+u2_of_gate Y = u2_Y+u2_of_gate Z = u2_Z+u2_of_gate H = u2_H+u2_of_gate S = u2_S+u2_of_gate T = u2_T+u2_of_gate E = u2_E+u2_of_gate W = u2_W++instance (RootHalfRing a, ComplexRing a) => FromGates (U2 a) where+  from_gates = product' . map u2_of_gate where+    product' = foldl' (*) 1++-- ----------------------------------------------------------------------+-- ** Matrices in /SO/(3)++-- $ This is the Bloch sphere representation of single qubit+-- operators.++-- | The Pauli /X/ operator.+so3_X :: (Ring a) => SO3 a+so3_X = matrix3x3 (1,  0,  0)+                  (0, -1,  0)+                  (0,  0, -1)++-- | The Pauli /Y/ operator.+so3_Y :: (Ring a) => SO3 a+so3_Y = matrix3x3 (-1, 0,  0)+                  ( 0, 1,  0)+                  ( 0, 0, -1)++-- | The Pauli /Z/ operator.+so3_Z :: (Ring a) => SO3 a+so3_Z = matrix3x3 (-1,  0, 0)+                  ( 0, -1, 0)+                  ( 0,  0, 1)++-- | The Hadamard operator.+so3_H :: (Ring a) => SO3 a+so3_H = matrix3x3 (0,  0, 1)+                  (0, -1, 0)+                  (1,  0, 0)++-- | The operator /S/.+so3_S :: (Ring a) => SO3 a+so3_S = matrix3x3 (0, -1, 0)+                  (1,  0, 0)+                  (0,  0, 1)++-- | The operator /E/.+so3_E :: (Ring a) => SO3 a+so3_E = matrix3x3 (0, 0, 1)+                  (1, 0, 0)+                  (0, 1, 0)++-- | The /T/ operator.+so3_T :: (RootHalfRing a) => SO3 a+so3_T = matrix3x3 (r, -r,  0)+                  (r,  r,  0)+                  (0,  0,  1)+  where r = roothalf++-- | Convert a symbolic gate to the corresponding Bloch sphere+-- operator.+so3_of_gate :: (RootHalfRing a) => Gate -> SO3 a+so3_of_gate X = so3_X+so3_of_gate Y = so3_Y+so3_of_gate Z = so3_Z+so3_of_gate H = so3_H+so3_of_gate S = so3_S+so3_of_gate T = so3_T+so3_of_gate E = so3_E+so3_of_gate W = 1++instance (RootHalfRing a) => FromGates (SO3 a) where+  from_gates = product . map so3_of_gate++-- ----------------------------------------------------------------------+-- ** Conversions++-- | Conversion from /U/(2) to /SO/(3).+so3_of_u2 :: (Adjoint a, ComplexRing a, RealPart a b, HalfRing b) => U2 a -> SO3 b+so3_of_u2 u = matrix_of_function f where+  f i j = half * (real $ tr (sigma i * u * sigma j * adj u))+  sigma 0 = u2_X+  sigma 1 = u2_Y+  sigma 2 = u2_Z+  sigma _ = error "so3_of_u2" -- not reached++-- | Convert a Clifford operator to a matrix in /SO/(3).+so3_of_clifford :: (ToClifford a, Ring b) => a -> SO3 b+so3_of_clifford m = so3_E^a * so3_X^b * so3_S^c where+  (a,b,c,d) = clifford_decompose m++-- | Convert a matrix in /SO/(3) to a Clifford gate. Throw an error if+-- the matrix isn't Clifford.+clifford_of_so3 :: (Ring a, Eq a, Adjoint a) => SO3 a -> Clifford+clifford_of_so3 m = case columns_of_matrix m of+  [_, _, [ 1, 0, 0]] -> with "H"+  [_, _, [-1, 0, 0]] -> with "HX"+  [_, _, [ 0, 1, 0]] -> with "SH"+  [_, _, [ 0,-1, 0]] -> with "SHX"+  [_, _, [ 0, 0,-1]] -> with "X"+  [_, [-1, 0, 0], _] -> with "S"+  [_, [ 0,-1, 0], _] -> with "SS"+  [_, [ 1, 0, 0], _] -> with "SSS"+  [[1, 0, 0], [0, 1, 0], [0, 0, 1]] -> clifford_id+  _ -> error "clifford_of_so3: not a Clifford operator"+  where+    with s = op `clifford_mult` op1 where+      op = to_clifford s+      m1 = adj (so3_of_clifford op) * m+      op1 = clifford_of_so3 m1++instance (Ring a, Eq a, Adjoint a) => ToClifford (SO3 a) where+  to_clifford = clifford_of_so3++-- ----------------------------------------------------------------------+-- * Matsumoto-Amano normal forms++-- $ A Matsumoto-Amano normal form is a sequence of Clifford+/T/+-- operators that is of the form+--+-- * (ε | /T/) (/HT/ | /SHT/)[sup *] /C/.+--+-- Here, ε is the empty sequence, /C/ is any Clifford operator, and+-- the meanings of @\"|\"@ and @\"*\"@ are as for regular+-- expressions. Every single-qubit Clifford+/T/ operator has a unique+-- Matsumoto-Amano normal form.++-- ----------------------------------------------------------------------+-- ** Representation of normal forms++-- | A representation of normal forms, optimized for right+-- multiplication.+data NormalForm = NormalForm Syllables Clifford+                  deriving (Eq)++-- | Syllables is a circuit of the form (ε|/T/) (/HT/|/SHT/)[sup *].+data Syllables =+  S_I                  -- ^ The empty sequence ε.+  | S_T                -- ^ The sequence /T/.+  | SApp_HT Syllables  -- ^ A sequence of the form …/HT/.+  | SApp_SHT Syllables -- ^ A sequence of the form …/SHT/.+            deriving (Eq, Show)++instance ToGates NormalForm where+  to_gates (NormalForm ts c) = to_gates ts ++ to_gates c++instance ToGates Syllables where+  to_gates S_I = []+  to_gates S_T = [T]+  to_gates (SApp_HT ts) = to_gates ts ++ [H, T]+  to_gates (SApp_SHT ts) = to_gates ts ++ [S, H, T]++instance Show NormalForm where+  show x = case to_gates x of+    [] -> "I"+    gs -> concat $ map show gs++-- | Right-multiply the given normal form by a gate.+normalform_append :: NormalForm -> Gate -> NormalForm+normalform_append (NormalForm ts c) X =+  NormalForm ts (c `clifford_mult` clifford_X)+normalform_append (NormalForm ts c) Y =+  NormalForm ts (c `clifford_mult` clifford_Y)+normalform_append (NormalForm ts c) Z =+  NormalForm ts (c `clifford_mult` clifford_Z)+normalform_append (NormalForm ts c) H =+  NormalForm ts (c `clifford_mult` clifford_H)+normalform_append (NormalForm ts c) S =+  NormalForm ts (c `clifford_mult` clifford_S)+normalform_append (NormalForm ts c) E =+  NormalForm ts (c `clifford_mult` clifford_E)+normalform_append (NormalForm ts c) W =+  NormalForm ts (c `clifford_mult` clifford_W)+normalform_append (NormalForm ts c) T+  | k == Axis_H = NormalForm (SApp_HT ts) c'+  | k == Axis_SH = NormalForm (SApp_SHT ts) c'+  | otherwise = case ts of+      S_I -> NormalForm S_T c'+      S_T -> NormalForm S_I (clifford_S `clifford_mult` c')+      SApp_HT ts' -> NormalForm ts' (clifford_HS `clifford_mult` c')+      SApp_SHT ts' -> NormalForm ts' (clifford_SHS `clifford_mult` c')+  where+    (k, c') = clifford_tconj c+    clifford_HS = to_clifford "HS"+    clifford_SHS = to_clifford "SHS"++-- ----------------------------------------------------------------------+-- ** Group operations on normal forms++-- | The identity as a normal form.+nf_id :: NormalForm+nf_id = NormalForm S_I clifford_id++-- | Multiply two normal forms. The right factor can be any+-- 'ToGates'.+nf_mult :: (ToGates b) => NormalForm -> b -> NormalForm+nf_mult a b = foldl' normalform_append a (to_gates b)++-- | Invert a normal form. The input can be any 'ToGates'.+nf_inv :: (ToGates a) => a -> NormalForm+nf_inv = from_gates . invert_gates . to_gates++-- ----------------------------------------------------------------------+-- ** Conversion to normal form++-- | Convert any 'ToGates' list to a 'NormalForm', thereby normalizing it.+normalize :: (ToGates a) => a -> NormalForm+normalize = nf_mult nf_id++instance FromGates NormalForm where+  from_gates = normalize++-- ----------------------------------------------------------------------+-- * Exact synthesis++-- ----------------------------------------------------------------------+-- ** Synthesis from /SO/(3)++-- | Input an exact matrix in /SO/(3), and output the corresponding+-- Clifford+/T/ normal form. It is an error if the given matrix is not+-- an element of /SO/(3), i.e., orthogonal with determinant 1.+--+-- This implementation uses the Matsumoto-Amano algorithm.+-- +-- Note: the list of gates will be returned in right-to-left order,+-- i.e., as in the mathematical notation for matrix multiplication.+-- This is the opposite of the quantum circuit notation.+synthesis_bloch :: SO3 DRootTwo -> [Gate]+synthesis_bloch m = aux m1 k+  where+    (m1, k) = denomexp_decompose m++    aux :: SO3 ZRootTwo -> Integer -> [Gate]    +    aux m 0 = to_gates (clifford_of_so3 m)+    aux m k = to_gates axis ++ [T] ++ aux m4 (k-1)+      where+        Matrix p = matrix_map parity m+        v1 = vector_head p+        v2 = vector_head (vector_tail p)+        v = list_of_vector $ vector_zipwith (\x y -> x + y - x*y) v1 v2+        axis = case v of+          [1, 1, 0] -> Axis_I+          [0, 1, 1] -> Axis_H+          [1, 0, 1] -> Axis_SH+          _ -> error "synthesis_bloch: not unitary"+        m2 = adj (so3_of_clifford axis) * m+        m3 = adj sqrt2T * m2+        m4 = matrix_map half_ZRootTwo m3+    sqrt2T = matrix3x3 (1, -1, 0) (1, 1, 0) (0, 0, roottwo)++    -- Divide a 'ZRootTwo' of the form 2/a/ + 2/b/√2 by 2, or throw an+    -- error if it is not of the required form.+    half_ZRootTwo :: ZRootTwo -> ZRootTwo+    half_ZRootTwo (RootTwo a b)+      | even a && even b = RootTwo a' b'+      | otherwise = error "synthesis_bloch: not unitary"+      where+        a' = a `div` 2+        b' = b `div` 2++instance (ToQOmega a) => ToGates (SO3 a) where+  to_gates = synthesis_bloch . matrix_map (to_dyadic . to_real . toQOmega)+    where+      to_real :: QOmega -> QRootTwo+      to_real x = case fromQOmega x :: QRComplex of +        (Cplx a 0) -> a+        _ -> error "to_gates: not a real number"++-- ----------------------------------------------------------------------+-- ** Synthesis from /U/(2)++instance ToGates TwoLevel where+  to_gates (TL_X 0 1) = [X]+  to_gates (TL_X 1 0) = [X]+  to_gates (TL_H 0 1) = [H]+  to_gates (TL_H 1 0) = [X,H,X]+  to_gates (TL_T k 0 1)+    | k `mod` 2 == 1 = [T] ++ to_gates (TL_T (k-1) 0 1)+    | k `mod` 4 == 2 = [S] ++ to_gates (TL_T (k-2) 0 1)+    | k `mod` 8 == 4 = [Z]+    | otherwise = []+  to_gates (TL_T k 1 0) = [X] ++ to_gates (TL_T k 0 1) ++ [X]+  to_gates (TL_omega k 1) = to_gates (TL_T k 0 1)+  to_gates (TL_omega k 0) = to_gates (TL_T k 1 0)+  to_gates _ = error $ "ToGates TwoLevel: invalid gate"++-- | Input an exact matrix in /U/(2), and output the corresponding+-- Clifford+/T/ normal form. The behavior is undefined if the given+-- matrix is not an element of /U/(2), i.e., unitary with determinant+-- 1.+--+-- We use a variant of the Kliuchnikov-Maslov-Mosca algorithm, as+-- implemented in "Quantum.Synthesis.MultiQubitSynthesis".+-- +-- Note: the list of gates will be returned in right-to-left order,+-- i.e., as in the mathematical notation for matrix multiplication.+-- This is the opposite of the quantum circuit notation.+synthesis_u2 :: U2 DOmega -> [Gate]+synthesis_u2 = to_gates . normalize . synthesis_nqubit++instance (ToQOmega a) => ToGates (U2 a) where+  to_gates = synthesis_u2 . matrix_map (fromDOmega . to_dyadic . toQOmega)++-- ----------------------------------------------------------------------+-- * Compact representation of normal forms++-- $ It is sometimes useful to store Clifford+/T/ operators in a file;+-- for this purpose, we provide a very succinct encoding of+-- Clifford+/T/ operators as bit strings, which are in turns+-- represented as integers.+--+-- Our bitwise encoding is as follows. The first regular expression+-- represents the set of Matsumoto-Amano normal forms (with a+-- particular presentation of the rightmost Clifford operator). The+-- second regular expression, which has the same form, defines the+-- corresponding bit string encoding.+--+-- * (ε|/T/) (/HT/|/SHT/)[sup *] (ε|/H/|/SH/) (ε|/X/) (ε|/S²/) (ε|/S/) (ε|ω⁴) (ε|ω²) (ε|ω)+--+-- * (10|11) (0|1)[sup *] (00|01|10) (0|1) (0|1) (0|1) (0|1) (0|1) (0|1)+--+-- As a special case, the leading bits 10 are omitted in case the+-- encoded operator is a Clifford operator. This ensures that the+-- encoding of a Clifford operator is an integer from 0 to 191.+-- +-- This format has the property that the encoded Clifford+/T/+-- operator can, in principle, be read off directly from the hexadecimal+-- representation of the bit string, with the following decoding:+--+-- Leftmost one or two hexadecimal digits:+--+-- >  0 = n/a             4 = HT              8 = HTHT            c = THTHT+-- >  1 = see below       5 = SHT             9 = HTSHT           d = THTSHT+-- >  2 = ε               6 = THT             a = SHTHT           e = TSHTHT+-- >  3 = T               7 = TSHT            b = SHTSHT          f = TSHTSHT+-- >+-- >  10 = HTHTHT         14 = SHTHTHT        18 = THTHTHT        1c = TSHTHTHT+-- >  11 = HTHTSHT        15 = SHTHTSHT       19 = THTHTSHT       1d = TSHTHTSHT+-- >  12 = HTSHTHT        16 = SHTSHTHT       1a = THTSHTHT       1e = TSHTSHTHT+-- >  13 = HTSHTSHT       17 = SHTSHTSHT      1b = THTSHTSHT      1f = TSHTSHTSHT+--+-- Central hexadecimal digit:+--+-- >  0 = HTHTHTHT        4 = HTSHTHTHT       8 = SHTHTHTHT       c = SHTSHTHTHT+-- >  1 = HTHTHTSHT       5 = HTSHTHTSHT      9 = SHTHTHTSHT      d = SHTSHTHTSHT+-- >  2 = HTHTSHTHT       6 = HTSHTSHTHT      a = SHTHTSHTHT      e = SHTSHTSHTHT+-- >  3 = HTHTSHTSHT      7 = HTSHTSHTSHT     b = SHTHTSHTSHT     f = SHTSHTSHTSHT+--+-- Second-to-rightmost hexadecimal digit:+--+-- >  0 = ε               4 = H               8 = SH              c = n/a+-- >  1 = SS              5 = HSS             9 = SHSS            d = n/a+-- >  2 = X               6 = HX              a = SHX             e = n/a+-- >  3 = XSS             7 = HXSS            b = SHXSS           f = n/a+--+-- Rightmost hexadecimal digit:+--+-- >  0 = ε               4 = ω⁴              8 = S               c = Sω⁴+-- >  1 = ω               5 = ω⁵              9 = Sω              d = Sω⁵+-- >  2 = ω²              6 = ω⁶              a = Sω²             e = Sω⁶+-- >  3 = ω³              7 = ω⁷              b = Sω³             f = Sω⁷+--+-- For example, the hexadecimal integer+--+-- > 6bf723e31+--+-- encodes the Clifford+/T/ operator+--+-- > THT SHTHTSHTSHT SHTSHTSHTSHT HTSHTSHTSHT HTHTSHTHT HTHTSHTSHT SHTSHTSHTHT XSS ω.++-- | Compactly encode a 'NormalForm' as an 'Integer'.+normalform_pack :: NormalForm -> Integer+normalform_pack (NormalForm S_I op) = clifford_pack op+normalform_pack (NormalForm s op) = 256 * syllables_pack s + clifford_pack op+  where+    syllables_pack :: Syllables -> Integer+    syllables_pack S_I = 2+    syllables_pack S_T = 3+    syllables_pack (SApp_HT s) = (syllables_pack s `shiftL` 1) + 0+    syllables_pack (SApp_SHT s) = (syllables_pack s `shiftL` 1) + 1++-- | Decode a 'NormalForm' from its 'Integer' encoding. This is the+-- inverse of 'normalform_pack'.+normalform_unpack :: Integer -> NormalForm+normalform_unpack n +  | n < 0 = error "normalform_unpack: invalid encoding"+  | n < 192 = NormalForm S_I op+  | n < 768 = error "normalform_unpack: invalid encoding"+  | otherwise = NormalForm s op +  where+    s = syllables_unpack (n `shiftR` 8)+    op = clifford_unpack (n .&. 0xff)++    syllables_unpack :: Integer -> Syllables+    syllables_unpack 0 = error "normalform_unpack: invalid encoding"+    syllables_unpack 1 = error "normalform_unpack: invalid encoding"+    syllables_unpack 2 = S_I+    syllables_unpack 3 = S_T+    syllables_unpack n+      | even n     = SApp_HT s+      | otherwise  = SApp_SHT s+      where+        s = syllables_unpack (n `shiftR` 1)++-- | Encode a Clifford operator as an integer in the range 0−191.+clifford_pack :: Clifford -> Integer+clifford_pack op = toInteger (64 * encode k + 32*b + 8*c + d)+  where+    (k, b, c, d) = clifford_decompose_coset op+    encode Axis_I = 0+    encode Axis_H = 1+    encode Axis_SH = 2++-- | Decode a Clifford operator from its integer encoding. This is the+-- inverse of 'clifford_pack'+clifford_unpack :: Integer -> Clifford+clifford_unpack n +  | n < 0 || n > 191 = error "clifford_unpack: invalid encoding"+  | otherwise = decode k * (clifford_X^b) * (clifford_S^c) * (clifford_W^d)+  where+    d = n .&. 0x7+    c = (n `shiftR` 3) .&. 0x3+    b = (n `shiftR` 5) .&. 0x1+    k = (n `shiftR` 6) .&. 0x3+    decode 0 = clifford_id+    decode 1 = clifford_H+    decode _ = clifford_SH+    (*) = clifford_mult+    (^) x n = foldl (*) clifford_id (genericReplicate n x)++instance ToGates Integer where+  to_gates = to_gates . normalform_unpack++instance FromGates Integer where+  from_gates = normalform_pack . from_gates
+ Quantum/Synthesis/EuclideanDomain.hs view
@@ -0,0 +1,153 @@+{-# LANGUAGE FlexibleInstances #-}++-- | This module provides a type class for Euclidean domains. A+-- Euclidean domain is a ring with a notion of division with+-- remainder, and therefore greatest common divisors.++module Quantum.Synthesis.EuclideanDomain where++import Quantum.Synthesis.Ring+import Data.Maybe++-- ----------------------------------------------------------------------+-- * Euclidean domains++-- ----------------------------------------------------------------------+-- ** Definition++-- | A type class for Euclidean domains. A Euclidean domain is a ring+-- with a Euclidean function and a division with remainder.+class (Eq a, Ring a) => EuclideanDomain a where+  -- | The Euclidean function for the Euclidean domain. This is a+  -- function /rank/ : /R/\\{0} → ℕ such that:+  -- +  -- * for all nonzero /a/, /b/ ∈ /R/, /rank/(/a/) ≤ /rank/(/ab/);+  -- +  -- * if /b/ ≠ 0 and (/q/,/r/) = /a/ `divmod` /b/, then either /r/ =+  -- 0 or /rank/(/r/) < /rank/(/b/).+  rank :: a -> Integer+  -- | Given /a/ and /b/≠0, return a quotient and remainder for+  -- division of /a/ by /b/. Specifically, return (/q/,/r/) such that+  -- /a/ = /qb/ + /r/, and such that /r/ = 0 or /rank/(/r/) < /rank/(/b/).+  divmod :: a -> a -> (a,a)++-- ----------------------------------------------------------------------+-- Particular Euclidean domains++instance EuclideanDomain Integer where+  rank x = x+  divmod x y = divMod x y++instance EuclideanDomain ZComplex where+  rank x = abs (norm x)+  divmod x y = (q, r) where+    (Cplx l m) = x * adj y+    k = norm y+    q1 = l `rounddiv` k+    q2 = m `rounddiv` k+    q = Cplx q1 q2+    r = x - y * q++instance EuclideanDomain ZRootTwo where+  rank x = abs (norm x)+  divmod x y@(RootTwo c d) = (q, r) where+    (RootTwo l m) = x * adj2 y+    k = norm y+    q1 = l `rounddiv` k+    q2 = m `rounddiv` k+    q = RootTwo q1 q2+    r = x - y * q+  +instance EuclideanDomain ZOmega where+  rank x = abs (norm x)+  divmod x y = (q, r) where+    (Omega a' b' c' d') = x * adj y * adj2(y * adj y)+    k = norm y+    a = a' `rounddiv` k+    b = b' `rounddiv` k+    c = c' `rounddiv` k+    d = d' `rounddiv` k+    q = Omega a b c d+    r = x - y * q    ++-- ----------------------------------------------------------------------+-- ** Functions++-- | Calculate the remainder for the division of /x/ by /y/.+euclid_mod :: (EuclideanDomain a) => a -> a -> a+euclid_mod x y = r where+  (q,r) = x `divmod` y++infixl 7 `euclid_mod`++-- | Calculate the quotient for the division of /x/ by /y/, ignoring+-- the remainder, if any. This is typically, but not always, used in+-- situations where the remainder is known to be 0 ahead of time.+euclid_div :: (EuclideanDomain a) => a -> a -> a+euclid_div x y = q where+  (q,r) = x `divmod` y++infixl 7 `euclid_div`++-- | Calculate the greatest common divisor in any Euclidean domain.+euclid_gcd :: (EuclideanDomain a) => a -> a -> a+euclid_gcd x y+  | y == 0 = x+  | otherwise = euclid_gcd y r where+    (_,r) = divmod x y++-- | Perform the extended Euclidean algorithm. On inputs /x/ and+-- /y/, this returns (/a/,/b/,/s/,/t/,/d/) such that:+-- +-- * /d/ = gcd(/x/,/y/),+-- +-- * /ax/ + /by/ = /d/,+-- +-- * /sx/ + /ty/ = 0,+-- +-- * /at/ - /bs/ = 1.+extended_euclid :: (EuclideanDomain a) => a -> a -> (a, a, a, a, a)+extended_euclid x y+  | y == 0 = (1, 0, 0, 1, x)+  | otherwise = (b',a'-b'*q,-t',t'*q-s',d) where+    (a',b',s',t',d) = extended_euclid y r+    (q,r) = divmod x y++-- | Find the inverse of a unit in a Euclidean domain. If the given+-- element is not a unit, return 'Nothing'.+euclid_inverse :: (EuclideanDomain a) => a -> Maybe a+euclid_inverse x+  | x == 0    = Nothing+  | r == 0    = Just q+  | otherwise = Nothing+  where+    (q,r) = divmod 1 x++-- | Determine whether an element of a Euclidean domain is a unit.+is_unit :: (EuclideanDomain a) => a -> Bool+is_unit = isJust . euclid_inverse++-- | Compute the inverse of /a/ in /R/\/(p), where /R/ is a Euclidean+-- domain. Note: this works whenever /a/ and /p/ are relatively+-- prime. If /a/ and /p/ are not relatively prime, return 'Nothing'.+inv_mod :: EuclideanDomain a => a -> a -> Maybe a+inv_mod p a = case euclid_inverse d of+  Just d' -> let (q,r) = (b*d') `divmod` p in Just r+  Nothing -> Nothing+  where+    (b,_,_,_,d) = extended_euclid a p++-- ----------------------------------------------------------------------+-- * Auxiliary functions++-- | For /y/ ≠ 0, find the integer /q/ closest to /x/ \/ /y/. This+-- works regardless of whether /x/ and\/or /y/ are positive or+-- negative.  The distance /q/ − /x/ \/ /y/ is guaranteed to be in+-- (-1\/2, 1\/2].+rounddiv :: (Integral a) => a -> a -> a+rounddiv x y = +  -- Note: the use of "quot" and "div" is crucial for the signs to+  -- work out correctly.+  (x + y `quot` 2) `div` y++infixl 7 `rounddiv`
+ Quantum/Synthesis/EulerAngles.hs view
@@ -0,0 +1,51 @@+-- | This module provides functions for converting between matrices in+-- /U/(2) and their Euler angle representation.++module Quantum.Synthesis.EulerAngles where++import Quantum.Synthesis.Ring+import Quantum.Synthesis.Matrix+import Quantum.Synthesis.ArcTan2++-- ----------------------------------------------------------------------+-- * Documentation++-- | Decompose a unitary operator /U/ into Euler angles (α, β, γ, δ).+-- These angles are computed so that+-- +-- * /U/ = [exp /i/α] R[sub /z/](β) R[sub /x/](γ) R[sub /z/](δ).+euler_angles :: (Floating a, ArcTan2 a) => Matrix Two Two (Cplx a) -> (a, a, a, a)+euler_angles op = (alpha, beta, gamma, delta) where+  ((a, b), (c, d)) = from_matrix2x2 op+  beta_plus_delta_over_2 = phase d - alpha+  beta_minus_delta_over_2 = phase c - alpha + pi/2+  alpha = phase (a*d - b*c) / 2+  gamma = 2 * arctan2 (mag b) (mag a)+  delta = phase (b * d * i * adj (a*d - b*c))+  beta = 2 * phase (d * cis (-alpha - delta/2) + c * cis (-alpha+delta/2) * i)++  mag (Cplx a b) = sqrt (a^2 + b^2)+  phase (Cplx a b) = arctan2 b a++  cis x = Cplx (cos x) (sin x)++  adj (Cplx x y) = Cplx x (-y)++-- | Compute the operator+-- +-- * /U/ = [exp /i/α] R[sub /z/](β) R[sub /x/](γ) R[sub /z/](δ).+-- +-- from the given Euler angles.+matrix_of_euler_angles :: (Floating a) => (a, a, a, a) -> Matrix Two Two (Cplx a)  +matrix_of_euler_angles (alpha, beta, gamma, delta) = op where+  op = opa * opb * opc * opd+  opa = cplx_cis alpha `scalarmult` 1+  opb = zrot beta+  opc = hadamard * zrot gamma * hadamard+  opd = zrot delta+  +  cplx_cis theta = Cplx (cos theta) (sin theta)+  hadamard = Cplx (sqrt 0.5) 0 `scalarmult` matrix2x2 (1, 1) (1, -1)+  zrot gamma = matrix2x2 (cplx_cis (-gamma/2), 0) (0, cplx_cis (gamma/2))++
+ Quantum/Synthesis/LaTeX.hs view
@@ -0,0 +1,180 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE OverlappingInstances #-}++-- | This module provides some functionality for pretty-printing+-- certain types to LaTeX format.++module Quantum.Synthesis.LaTeX where++import Quantum.Synthesis.CliffordT+import Quantum.Synthesis.MultiQubitSynthesis+import Quantum.Synthesis.Ring+import Quantum.Synthesis.Matrix+import Quantum.Synthesis.SymReal++import Text.Printf+import Data.Ratio++-- | A type class for things that can be printed to LaTeX format. +-- +-- Minimal complete definition: 'showlatex' or 'showlatex_p'.++-- This is a bit naive at the moment - to do it properly, one should+-- perhaps also supply context information, for example math mode/text+-- mode.+class ShowLaTeX a where+  -- | Print to LaTeX format.+  showlatex :: a -> String+  showlatex x = showlatex_p 0 x ""+  +  -- | Print to LaTeX format, with precedence. Analogous to 'showsPrec'.+  showlatex_p :: Int -> a -> ShowS+  showlatex_p _ x s = showlatex x ++ s++instance ShowLaTeX TwoLevel where+  showlatex (TL_X i j) = printf "X\\level{%d,%d} " (i+1) (j+1)+  showlatex (TL_H i j) = printf "H\\level{%d,%d} " (i+1) (j+1)+  showlatex (TL_T m i j)+    | m' == 0 = ""+    | m' == 1 = printf "T\\level{%d,%d} " (i+1) (j+1)+    | otherwise = printf "T^%d\\level{%d,%d} " m' (i+1) (j+1)+    where m' = m `mod` 8+  showlatex (TL_omega m i)+    | m' == 0 = ""+    | m' == 1 = printf "\\omega\\level{%d} " (i+1)+    | otherwise = printf "\\omega^%d\\level{%d} " m' (i+1)+    where m' = m `mod` 8+  +instance ShowLaTeX [TwoLevel] where+  showlatex = concat . map showlatex++instance ShowLaTeX Integer where+  showlatex = show++instance ShowLaTeX ZOmega where+  showlatex (Omega a b c d) = format_signed_list list2 where+    list = map signedunit [(a,"\\omega^3"),(b,"\\omega^2"),(c,"\\omega"),(d,"")]+    list2 = filter (\(s,a) -> s /= 0) list+    signedunit (a, u) +      | u == ""   = (s, showlatex a')+      | a' == 1   = (s, u)+      | otherwise = (s, showlatex a' ++ u)+      where+        (s,a') = tosigned a+    tosigned a +      | a < 0     = (-1,-a)+      | a == 0    = (0,0)+      | otherwise = (1,a)+    format_signed_list [] = "0"+    format_signed_list ((1,a):t) = a ++ cont t +    format_signed_list ((_,a):t) = "-" ++ a ++ cont t +    cont [] = ""+    cont ((1,a):t) = "+" ++ a ++ cont t+    cont ((0,a):t) = cont t+    cont ((_,a):t) = "-" ++ a ++ cont t++instance (ShowLaTeX a, Nat n) => ShowLaTeX (Matrix n m a) where+  showlatex (Matrix a) = "\\zmatrix{" ++ replicate m 'c' ++ "}{" ++ entries ++ "}" where+    m = length (list_of_vector a)+    entries = concat $ list_of_vector $ vector_map showcolumn (vector_transpose a)+    showcolumn :: ShowLaTeX a => Vector m a -> String+    showcolumn Nil = "\\\\"+    showcolumn (h `Cons` Nil) = showlatex h ++ "\\\\"+    showcolumn (h `Cons` t) = showlatex h ++ " & " ++ showcolumn t++instance ShowLaTeX Rational where+  showlatex r = "\\frac{" ++ showlatex num ++ "}{" ++ showlatex denom ++ "}"+    where+      num = numerator r+      denom = numerator r++instance ShowLaTeX Dyadic where+  showlatex = showlatex . toRational++instance (ShowLaTeX a, Eq a, Ring a) => ShowLaTeX (RootTwo a) where+  showlatex_p d (RootTwo a 0) = showlatex_p d a+  showlatex_p d (RootTwo 0 1) = showString "\\sqrt{2}"+  showlatex_p d (RootTwo 0 (-1)) = showParen (d >= 7) $ showString "-\\sqrt{2}"+  showlatex_p d (RootTwo 0 b) = showParen (d >= 8) $ +    showlatex_p 7 b . showString " \\sqrt{2}"+  showlatex_p d (RootTwo a b) | signum b == 1 = showParen (d >= 7) $+    showlatex_p 6 a . showString " + " . showlatex_p 6 (RootTwo 0 b)+  showlatex_p d (RootTwo a b) | otherwise = showParen (d >= 7) $+    showlatex_p 6 a . showString " - " . showlatex_p 7 (RootTwo 0 (-b))+  ++instance ShowLaTeX (Omega Z2) where+  showlatex (Omega a b c d) = concat $ map show [a,b,c,d]++instance (ShowLaTeX a, Ring a, Eq a) => ShowLaTeX (Cplx a) where+  showlatex_p d (Cplx a 0) = showlatex_p d a+  showlatex_p d (Cplx 0 1) = showString "i"+  showlatex_p d (Cplx 0 (-1)) = showParen (d >= 7) $ showString "-i"+  showlatex_p d (Cplx 0 b) = showParen (d >= 8) $+    showlatex_p 7 b . showString "\\,i"+  showlatex_p d (Cplx a b) | signum b == 1 = showParen (d >= 7) $+    showlatex_p 6 a . showString "+" . showlatex_p 7 (Cplx 0 b)+  showlatex_p d (Cplx a b) | otherwise = showParen (d >= 7) $ +    showlatex_p 6 a . showString "-" . showlatex_p 7 (Cplx 0 (-b))++instance ShowLaTeX Double where+  showlatex x = printf "%0.10f" x++-- This is an overlapping instance+instance Nat n => ShowLaTeX (Matrix n m DOmega) where+  showlatex = showlatex_denomexp++-- This is an overlapping instance+instance Nat n => ShowLaTeX (Matrix n m DRComplex) where+  showlatex = showlatex_denomexp++-- | Generic showlatex-like method that factors out a common+-- denominator exponent.+showlatex_denomexp :: (WholePart a b, ShowLaTeX b, DenomExp a) => a -> String+showlatex_denomexp a+  | k == 0 = showlatex b+  | k == 1 = "\\frac{1}{\\sqrt{2}}" ++ showlatex b+  | otherwise = "\\frac{1}{\\sqrt{2}^{" ++ show k ++ "}}" ++ showlatex b+    where (b, k) = denomexp_decompose a++instance ShowLaTeX [Gate] where+  showlatex [] = "\\epsilon"+  showlatex gates = aux 0 gates where+    aux n (W:t) = aux (n+1) t+    aux 0 []    = ""+    aux 1 []    = "{\\omega}"+    aux n []    = "\\omega^" ++ show n+    aux 0 (h:t) = show h ++ aux 0 t+    aux n t     = aux n [] ++ aux 0 t++instance ShowLaTeX SymReal where+  showlatex_p d (Const x)     = showlatex_p d x+  showlatex_p d (Decimal x s) = showString s+  showlatex_p d (Plus x y)    = showParen (d > 6) $ showlatex_p 6 x . showString "+" . showlatex_p 6 y+  showlatex_p d (Minus x y)   = showParen (d > 6) $ showlatex_p 6 x . showString "-" . showlatex_p 7 y+  showlatex_p d (Times x y)   = showParen (d > 7) $ showlatex_p 7 x . showString "\\cdot" . showlatex_p 7 y+  showlatex_p d (Div x y)     = showParen (d > 7) $ showlatex_p 7 x . showString "/" . showlatex_p 8 y+  showlatex_p d (Power x y)   = showParen (d > 11) $ showlatex_p 12 x . showString "^{" . showlatex_p 0 y . showString "}"+  showlatex_p d (Negate x)    = showParen (d > 5) $ showString "-" . showlatex_p 7 x+  showlatex_p d (Abs x)       = showParen (d > 10) $ showString "|" . showlatex_p 11 x . showString "|"+  showlatex_p d (Signum x)    = showParen (d > 10) $ showString "\\signum " . showlatex_p 11 x+  showlatex_p d (Recip x)     = showParen (d > 7) $ showString "1/" . showlatex_p 8 x+  showlatex_p d Pi            = showString "\\pi"+  showlatex_p d Euler         = showString "e"+  showlatex_p d (Exp x)       = showParen (d > 10) $ showString "e^{" . showlatex_p 0 x . showString "}"+  showlatex_p d (Sqrt x)      = showString "\\sqrt{" . showlatex_p 0 x . showString "}"+  showlatex_p d (Log x)       = showParen (d > 10) $ showString "\\log " . showlatex_p 11 x+  showlatex_p d (Sin x)       = showParen (d > 10) $ showString "\\sin " . showlatex_p 11 x+  showlatex_p d (Tan x)       = showParen (d > 10) $ showString "\\tan " . showlatex_p 11 x+  showlatex_p d (Cos x)       = showParen (d > 10) $ showString "\\cos " . showlatex_p 11 x+  showlatex_p d (ASin x)      = showParen (d > 10) $ showString "\\asin " . showlatex_p 11 x+  showlatex_p d (ATan x)      = showParen (d > 10) $ showString "\\atan " . showlatex_p 11 x+  showlatex_p d (ACos x)      = showParen (d > 10) $ showString "\\acos " . showlatex_p 11 x+  showlatex_p d (Sinh x)      = showParen (d > 10) $ showString "\\sinh " . showlatex_p 11 x+  showlatex_p d (Tanh x)      = showParen (d > 10) $ showString "\\tanh " . showlatex_p 11 x+  showlatex_p d (Cosh x)      = showParen (d > 10) $ showString "\\cosh " . showlatex_p 11 x+  showlatex_p d (ASinh x)     = showParen (d > 10) $ showString "\\asinh " . showlatex_p 11 x+  showlatex_p d (ATanh x)     = showParen (d > 10) $ showString "\\atanh " . showlatex_p 11 x+  showlatex_p d (ACosh x)     = showParen (d > 10) $ showString "\\acosh " . showlatex_p 11 x+  showlatex_p d (ArcTan2 y x) = showParen (d > 10) $ showString "\\arctan2 " . showlatex_p 11 y . showString " " . showlatex_p 11 x
+ Quantum/Synthesis/Matrix.hs view
@@ -0,0 +1,601 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE OverlappingInstances #-}+{-# LANGUAGE IncoherentInstances #-}+{-# LANGUAGE EmptyDataDecls #-}++{-# OPTIONS_GHC -fno-warn-incomplete-patterns #-}++-- | This module provides fixed but arbitrary sized vectors and+-- matrices. The dimensions of the vectors and matrices are determined+-- by the type, for example,+-- +-- > Matrix Two Three Complex+-- +-- for complex 2×3-matrices. The type system ensures that there are no+-- run-time dimension errors.++module Quantum.Synthesis.Matrix where++import Quantum.Synthesis.Ring++-- ----------------------------------------------------------------------+-- * Type-level natural numbers+  +-- $ Note: with Haskell 7.4.2 data-kinds, this could be replaced by a+-- tighter definition; however, the following works just fine in+-- Haskell 7.2.++-- | Type-level representation of zero.+data Zero++-- | Type-level representation of successor.+data Succ a++-- | The natural number 1 as a type.+type One = Succ Zero++-- | The natural number 2 as a type.+type Two = Succ One++-- | The natural number 3 as a type.+type Three = Succ Two++-- | The natural number 4 as a type.+type Four = Succ Three++-- | The natural number 5 as a type.+type Five = Succ Four++-- | The natural number 6 as a type.+type Six = Succ Five++-- | The natural number 7 as a type.+type Seven = Succ Six++-- | The natural number 8 as a type.+type Eight = Succ Seven++-- | The natural number 9 as a type.+type Nine = Succ Eight++-- | The natural number 10 as a type.+type Ten = Succ Nine++-- | The 10th successor of a natural number type. For example, the+-- natural number 18 as a type is+-- +-- > Ten_and Eight+type Ten_and a = Succ (Succ (Succ (Succ (Succ (Succ (Succ (Succ (Succ (Succ a)))))))))++-- | A data type for the natural numbers. Specifically, if /n/ is a+-- type-level natural number, then+-- +-- > NNat n+-- +-- is a singleton type containing only the natural number /n/.+data NNat :: * -> * where+  Zero :: NNat Zero+  Succ :: (Nat n) => NNat n -> NNat (Succ n)++-- | Convert an 'NNat' to an 'Integer'.+fromNNat :: NNat n -> Integer+fromNNat Zero = 0+fromNNat (Succ n) = 1 + fromNNat n++instance Show (NNat n) where+  show = show . fromNNat++-- | A type class for the natural numbers. The members are exactly the+-- type-level natural numbers.+class Nat n where+  -- | Return a term-level natural number corresponding to this+  -- type-level natural number.+  nnat :: NNat n+  +  -- | Return a term-level integer corresponding to this type-level+  -- natural number. The argument is just a dummy argument and is not+  -- evaluated.+  nat :: n -> Integer+  +instance Nat Zero where+  nnat = Zero+  nat n = 0+instance (Nat a) => Nat (Succ a) where+  nnat = Succ nnat+  nat n = 1 + nat (un n) where+    un :: Succ a -> a+    un = undefined++-- | Addition of type-level natural numbers.+type family Plus n m+type instance Zero `Plus` m = m+type instance (Succ n) `Plus` m = Succ (n `Plus` m)++-- | Multiplication of type-level natural numbers.+type family Times n m+type instance Zero `Times` m = Zero+type instance (Succ n) `Times` m = m `Plus` (n `Times` m)++-- ----------------------------------------------------------------------+-- * Fixed-length vectors++-- | @Vector /n/ /a/@ is the type of lists of length /n/ with elements+-- from /a/. We call this a \"vector\" rather than a tuple or list for+-- two reasons: the vectors are homogeneous (all elements have the+-- same type), and they are strict: if any one component is undefined,+-- the whole vector is undefined.+data Vector :: * -> * -> * where+  Nil :: Vector Zero a+  Cons :: !a -> !(Vector n a) -> Vector (Succ n) a++infixr 5 `Cons`++instance (Eq a) => Eq (Vector n a) where+  Nil == Nil = True+  Cons a as == Cons b bs = a == b && as == bs++instance (Show a) => Show (Vector n a) where+  showsPrec d x = showParen (d >= 11) $ showString ("vector " ++ show (list_of_vector x))++instance (ToDyadic a b) => ToDyadic (Vector n a) (Vector n b) where+  maybe_dyadic as = vector_sequence (vector_map maybe_dyadic as)++instance (WholePart a b) => WholePart (Vector n a) (Vector n b) where  +  from_whole = vector_map from_whole+  to_whole = vector_map to_whole+  +instance (DenomExp a) => DenomExp (Vector n a) where+  denomexp as = denomexp (list_of_vector as)+  denomexp_factor as k = vector_map (\a -> denomexp_factor a k) as+  +-- | Construct a vector of length 1.+vector_singleton :: a -> Vector One a+vector_singleton x = x `Cons` Nil++-- | Return the length of a vector. Since this information is+-- contained in the type, the vector argument is never evaluated and+-- can be a dummy (undefined) argument.+vector_length :: (Nat n) => Vector n a -> Integer+vector_length = nat . un where+  un :: Vector n a -> n+  un = undefined++-- | Convert a fixed-length list to an ordinary list.+list_of_vector :: Vector n a -> [a]+list_of_vector Nil = []+list_of_vector (Cons h t) = h : list_of_vector t++-- | Zip two equal length lists.+vector_zipwith :: (a -> b -> c) -> Vector n a -> Vector n b -> Vector n c+vector_zipwith f Nil Nil = Nil+vector_zipwith f (Cons a as) (Cons b bs) = Cons c cs where+  c = f a b+  cs = vector_zipwith f as bs++-- | Map a function over a fixed-length list.+vector_map :: (a -> b) -> Vector n a -> Vector n b+vector_map f Nil = Nil+vector_map f (Cons a as) = Cons (f a) (vector_map f as)++-- | Create the vector (0, 1, …, /n/-1).+vector_enum :: (Num a, Nat n) => Vector n a+vector_enum = aux nnat 0 where+  aux :: (Num a) => NNat n -> a -> Vector n a+  aux Zero a = Nil+  aux (Succ n) a = Cons a (aux n (a+1))++-- | Create the vector (/f/(0), /f/(1), …, /f/(/n/-1)).+vector_of_function :: (Num a, Nat n) => (a -> b) -> Vector n b+vector_of_function f = vector_map f vector_enum++-- | Construct a vector from a list. Note: since the length of the+-- vector is a type-level integer, it cannot be inferred from the+-- length of the input list; instead, it must be specified explicitly+-- in the type. It is an error to apply this function to a list of+-- the wrong length.+vector :: (Nat n) => [a] -> Vector n a+vector = aux nnat where+  aux :: NNat n -> [a] -> Vector n a+  aux Zero [] = Nil+  aux (Succ n) (h:t) = Cons h (aux n t)+  aux _ _ = error "vector: length mismatch"++-- | Return the /i/th element of the vector. Counting starts from 0.+-- Throws an error if the index is out of range.+vector_index :: (Integral i) => Vector n a -> i -> a+vector_index v i = list_of_vector v !! fromIntegral i++-- | Return a fixed-length list consisting of a repetition of the+-- given element. Unlike 'replicate', no count is needed, because this+-- information is already contained in the type. However, the type+-- must of course be inferable from the context.+vector_repeat :: (Nat n) => a -> Vector n a+vector_repeat x = vector_of_function (const x)++-- | Turn a list of columns into a list of rows.+vector_transpose :: (Nat m) => Vector n (Vector m a) -> Vector m (Vector n a)+vector_transpose Nil = vector_repeat Nil+vector_transpose (Cons a as) = vector_zipwith Cons a (vector_transpose as)++-- | Left strict fold over a fixed-length list.+vector_foldl :: (a -> b -> a) -> a -> Vector n b -> a+vector_foldl f x l = foldl f x (list_of_vector l)++-- | Right fold over a fixed-length list.+vector_foldr :: (a -> b -> b) -> b -> Vector n a -> b+vector_foldr f x l = foldr f x (list_of_vector l)++-- | Return the tail of a fixed-length list. Note that the type system+-- ensures that this never fails.+vector_tail :: Vector (Succ n) a -> Vector n a+vector_tail (Cons h t) = t++-- | Return the head of a fixed-length list. Note that the type system+-- ensures that this never fails.+vector_head :: Vector (Succ n) a -> a+vector_head (Cons h t) = h++-- | Append two fixed-length lists.+vector_append :: Vector n a -> Vector m a -> Vector (n `Plus` m) a+vector_append Nil v = v+vector_append (Cons h t) v = Cons h (vector_append t v)++-- | Version of 'sequence' for fixed-length lists.+vector_sequence :: (Monad m) => Vector n (m a) -> m (Vector n a)+vector_sequence Nil = return Nil+vector_sequence (Cons a as) = do+  a' <- a+  as' <- vector_sequence as+  return (Cons a' as')++-- ----------------------------------------------------------------------+-- * Matrices++-- | An /m/×/n/-matrix is a list of /n/ columns, each of which is a+-- list of /m/ scalars.  The type of square matrices of any fixed+-- dimension is an instance of the 'Ring' class, and therefore the+-- usual symbols, such as \"'+'\" and \"'*'\" can be used on+-- them. However, the non-square matrices, the symbols \"'.+.'\" and+-- \"'.*.'\" must be used.+data Matrix m n a = Matrix !(Vector n (Vector m a))+               deriving (Eq)++instance (Nat m, Show a) => Show (Matrix m n a) where+  showsPrec d m = showParen (d >= 11) $ showString ("matrix " ++ show (rows_of_matrix m))+  +-- This is an overlapping instance.+instance (Nat m) => Show (Matrix m n DRootTwo) where+  showsPrec = showsPrec_DenomExp+  +-- This is an overlapping instance.+instance (Nat m) => Show (Matrix m n DRComplex) where+  showsPrec = showsPrec_DenomExp++-- This is an overlapping instance.+instance (Nat m) => Show (Matrix m n DOmega) where+  showsPrec = showsPrec_DenomExp+  +instance (ToDyadic a b) => ToDyadic (Matrix m n a) (Matrix m n b) where+  maybe_dyadic (Matrix a) = do+    b <- maybe_dyadic a+    return (Matrix b)++instance (WholePart a b) => WholePart (Matrix m n a) (Matrix m n b) where+  from_whole (Matrix m) = Matrix (from_whole m)+  to_whole (Matrix m) = Matrix (to_whole m)++instance (DenomExp a) => DenomExp (Matrix m n a) where+  denomexp (Matrix m) = denomexp m+  denomexp_factor (Matrix m) k = Matrix (denomexp_factor m k)++-- | Decompose a matrix into a list of columns.+unMatrix :: Matrix m n a -> (Vector n (Vector m a))+unMatrix (Matrix m) = m++-- | Return the size (/m/, /n/) of a matrix, where /m/ is the number+-- of rows, and /n/ is the number of columns. Since this information+-- is contained in the type, the matrix argument is not evaluated and+-- can be a dummy (undefined) argument.+matrix_size :: (Nat m, Nat n) => Matrix m n a -> (Integer, Integer)+matrix_size op = (nat (m op), nat (n op)) where+  m :: Matrix m n a -> m+  m = undefined+  n :: Matrix m n a -> n+  n = undefined++-- ----------------------------------------------------------------------+-- ** Basic matrix operations++-- | Addition of /m/×/n/-matrices. We use a special symbol because+-- /m/×/n/-matrices do not form a ring; only /n/×/n/-matrices form a+-- ring (in which case the normal symbol \"'+'\" also works).+(.+.) :: (Num a) => Matrix m n a -> Matrix m n a -> Matrix m n a+Matrix a .+. Matrix b = Matrix c where+  c = vector_zipwith (vector_zipwith (+)) a b++infixl 6 .+.++-- | Subtraction of /m/×/n/-matrices. We use a special symbol because+-- /m/×/n/-matrices do not form a ring; only /n/×/n/-matrices form a+-- ring (in which case the normal symbol \"'-'\" also works).+(.-.) :: (Num a) => Matrix m n a -> Matrix m n a -> Matrix m n a+Matrix a .-. Matrix b = Matrix c where+  c = vector_zipwith (vector_zipwith (-)) a b++infixl 6 .-.++-- | Map some function over every element of a matrix.+matrix_map :: (a -> b) -> Matrix m n a -> Matrix m n b+matrix_map f (Matrix a) = Matrix b where+  b = vector_map (vector_map f) a++-- | Create the matrix whose /i/,/j/-entry is (/i/,/j/). Here /i/ and+-- /j/ are 0-based, i.e., the top left entry is (0,0).+matrix_enum :: (Num a, Nat n, Nat m) => Matrix m n (a,a)+matrix_enum = Matrix (vector_of_function f) where+  f i = vector_of_function (\j -> (j,i))++-- | Create the matrix whose /i/,/j/-entry is @f i j@. Here /i/ and+-- /j/ are 0-based, i.e., the top left entry is @f 0 0@.+matrix_of_function :: (Num a, Nat n, Nat m) => (a -> a -> b) -> Matrix m n b+matrix_of_function f = matrix_map (uncurry f) matrix_enum++-- | Multiplication of a scalar and an /m/×/n/-matrix.+scalarmult :: (Num a) => a -> Matrix m n a -> Matrix m n a+scalarmult x m = matrix_map (x *) m++infixl 7 `scalarmult`++-- | Multiplication of /m/×/n/-matrices. We use a special symbol+-- because /m/×/n/-matrices do not form a ring; only /n/×/n/-matrices+-- form a ring (in which case the normal symbol \"'*'\" also works).+(.*.) :: (Num a, Nat m) => Matrix m n a -> Matrix n p a -> Matrix m p a+Matrix a .*. Matrix b = Matrix c where+  c = vector_map (a `mmv`) b+  +  mmv :: (Num a, Nat m) => Vector n (Vector m a) -> Vector n a -> Vector m a+  Nil `mmv` Nil = vector_repeat 0+  (Cons h Nil) `mmv` (Cons k Nil) = k `msv` h+  (Cons h t) `mmv` (Cons k s) = (k `msv` h) `avv` (t `mmv` s)+  +  msv :: (Num b) => b -> Vector n b -> Vector n b+  k `msv` h = vector_map (k*) h+  +  avv :: (Num c) => Vector n c -> Vector n c -> Vector n c+  v `avv` w = vector_zipwith (+) v w++infixl 7 .*.++-- | Return the 0 matrix of the given dimension.+null_matrix :: (Num a, Nat n, Nat m) => Matrix m n a+null_matrix = Matrix (vector_repeat (vector_repeat 0))++-- | Take the transpose of an /m/×/n/-matrix.+matrix_transpose :: (Nat m) => Matrix m n a -> Matrix n m a+matrix_transpose (Matrix a) = Matrix b where+  b = vector_transpose a++-- | Take the adjoint of an /m/×/n/-matrix. Unlike 'adj', this can be+-- applied to non-square matrices.+adjoint :: (Nat m, Adjoint a) => Matrix m n a -> Matrix n m a+adjoint (Matrix a) = Matrix c where+  b = vector_map (vector_map adj) a+  c = vector_transpose b+  +-- | Return the element in the /i/th row and /j/th column of the+-- matrix. Counting of rows and columns starts from 0. Throws an error+-- if the index is out of range.+matrix_index :: (Integral i) => Matrix m n a -> i -> i -> a+matrix_index (Matrix a) i j = a `vector_index` j `vector_index` i++-- | Return a list of all the entries of a matrix, in some fixed but+-- unspecified order.+matrix_entries :: Matrix m n a -> [a]+matrix_entries (Matrix m) = +  concat $ map list_of_vector $ list_of_vector m++-- | Version of 'sequence' for matrices.+matrix_sequence :: (Monad m) => Matrix n p (m a) -> m (Matrix n p a)+matrix_sequence (Matrix m) = do+  m' <- vector_sequence (vector_map vector_sequence m)+  return (Matrix m')++-- | Return the trace of a square matrix.+tr :: (Ring a) => Matrix n n a -> a+tr (Matrix a) = aux a where+  aux :: (Num a) => Vector n (Vector n a) -> a+  aux Nil = 0+  aux ((h `Cons` t) `Cons` s) = h + aux (vector_map vector_tail s)++-- | Return the square of the Hilbert-Schmidt norm of an+-- /m/×/n/-matrix, defined by ‖/M/‖² = tr /M/[sup †]/M/.+hs_sqnorm :: (Ring a, Adjoint a, Nat n) => Matrix n m a -> a+hs_sqnorm m = tr (m .*. adjoint m)++-- ----------------------------------------------------------------------+-- Class instances for the ring of square matrices++instance (Num a, Nat n) => Num (Matrix n n a) where+  (+) = (.+.)+  (*) = (.*.)+  negate = scalarmult (-1)+  (-) = (.-.)+  fromInteger x = matrix_of_function (\i j -> if i == j then fromInteger x else 0)+  abs a = a+  signum a = 1+        +instance (Nat n, Adjoint a) => Adjoint (Matrix n n a) where+  adj (Matrix a) = Matrix c where+    b = vector_map (vector_map adj) a+    c = vector_transpose b++instance (Nat n, Adjoint2 a) => Adjoint2 (Matrix n n a) where+  adj2 (Matrix a) = Matrix b where+    b = vector_map (vector_map adj2) a++instance (HalfRing a, Nat n) => HalfRing (Matrix n n a) where+  half = scalarmult half 1++instance (RootHalfRing a, Nat n) => RootHalfRing (Matrix n n a) where+  roothalf = scalarmult roothalf 1++instance (RootTwoRing a, Nat n) => RootTwoRing (Matrix n n a) where+  roottwo = scalarmult roottwo 1++instance (ComplexRing a, Nat n) => ComplexRing (Matrix n n a) where+  i = scalarmult i 1++-- ----------------------------------------------------------------------+-- ** Operations on block matrices++-- | Stack matrices vertically.+stack_vertical :: Matrix m n a -> Matrix p n a -> Matrix (m `Plus` p) n a+stack_vertical (Matrix a) (Matrix b) = (Matrix c) where+  c = vector_zipwith vector_append a b++-- | Stack matrices horizontally.+stack_horizontal :: Matrix m n a -> Matrix m p a -> Matrix m (n `Plus` p) a+stack_horizontal (Matrix a) (Matrix b) = (Matrix c) where+  c = vector_append a b+  +-- | Repeat a matrix vertically, according to some vector of scalars.+tensor_vertical :: (Num a, Nat n) => Vector p a -> Matrix m n a -> Matrix (p `Times` m) n a+tensor_vertical v m = concat_vertical (vector_map (`scalarmult` m) v)+                               +-- | Vertically concatenate a vector of matrices.+concat_vertical :: (Num a, Nat n) => Vector p (Matrix m n a) -> Matrix (p `Times` m) n a+concat_vertical Nil = null_matrix+concat_vertical (Cons h t) = stack_vertical h (concat_vertical t)++-- | Repeat a matrix horizontally, according to some vector of scalars.+tensor_horizontal :: (Num a, Nat m) => Vector p a -> Matrix m n a -> Matrix m (p `Times` n) a+tensor_horizontal v m = concat_horizontal (vector_map (`scalarmult` m) v)+  +-- | Horizontally concatenate a vector of matrices.+concat_horizontal :: (Num a, Nat m) => Vector p (Matrix m n a) -> Matrix m (p `Times` n) a+concat_horizontal Nil = null_matrix+concat_horizontal (Cons h t) = stack_horizontal h (concat_horizontal t)++-- | Kronecker tensor of two matrices.+tensor :: (Num a, Nat n, Nat (p `Times` m)) => Matrix p q a -> Matrix m n a -> Matrix (p `Times` m) (q `Times` n) a+tensor a b = ab3 where+  Matrix ab1 = matrix_map (`scalarmult` b) a+  ab2 = vector_map concat_vertical ab1+  ab3 = concat_horizontal ab2++-- | Form a diagonal block matrix.+oplus :: (Num a, Nat m, Nat q, Nat n, Nat p) => Matrix p q a -> Matrix m n a -> Matrix (p `Plus` m) (q `Plus` n) a+oplus (a :: Matrix p q a) (b :: Matrix m n a) = +  (a `stack_vertical` (null_matrix :: Matrix m q a)) `stack_horizontal` ((null_matrix :: Matrix p n a) `stack_vertical` b)++-- | Form a controlled gate.+matrix_controlled :: (Eq a, Num a, Nat n) => Matrix n n a -> Matrix (n `Plus` n) (n `Plus` n) a+matrix_controlled (m :: Matrix n n a) = oplus (1 :: Matrix n n a) m++-- ----------------------------------------------------------------------+-- ** Constructors and destructors++-- | A convenient abbreviation for the type of 2×2-matrices.+type U2 a = Matrix Two Two a++-- | A convenient abbreviation for the type of 3×3-matrices.+type SO3 a = Matrix Three Three a++-- | A convenience constructor for matrices: turn a list of columns+-- into a matrix. +-- +-- Note: since the dimensions of the matrix are type-level integers,+-- they cannot be inferred from the dimensions of the input; instead,+-- they must be specified explicitly in the type. It is an error to+-- apply this function to a list of the wrong dimension.+matrix_of_columns :: (Nat n, Nat m) => [[a]] -> Matrix n m a+matrix_of_columns columns = Matrix m where+  m = vector $ map vector columns++-- | A convenience constructor for matrices: turn a list of rows into+-- a matrix.+-- +-- Note: since the dimensions of the matrix are type-level integers,+-- they cannot be inferred from the dimensions of the input; instead,+-- they must be specified explicitly in the type. It is an error to+-- apply this function to a list of the wrong dimension.+matrix_of_rows :: (Nat n, Nat m) => [[a]] -> Matrix n m a+matrix_of_rows = matrix_transpose . matrix_of_columns++-- | A synonym for 'matrix_of_rows'.+matrix :: (Nat n, Nat m) => [[a]] -> Matrix n m a+matrix = matrix_of_rows++-- | Turn a matrix into a list of columns.+columns_of_matrix :: Matrix n m a -> [[a]]+columns_of_matrix (Matrix m) = +  map list_of_vector (list_of_vector m)++-- | Turn a matrix into a list of rows.+rows_of_matrix :: (Nat n) => Matrix n m a -> [[a]]+rows_of_matrix = columns_of_matrix . matrix_transpose++-- | A convenience constructor for 2×2-matrices. The arguments are by+-- rows.+matrix2x2 :: (a, a) -> (a, a) -> Matrix Two Two a+matrix2x2 (a, b) (c, d) = matrix_of_columns [[a,c], [b,d]]++-- | A convenience destructor for 2×2-matrices. The result is by rows.+from_matrix2x2 :: Matrix Two Two a -> ((a, a), (a, a))+from_matrix2x2 (Matrix ((a `Cons` c `Cons` Nil) `Cons` (b `Cons` d `Cons` Nil) `Cons` Nil)) = ((a, b), (c, d))  ++-- | A convenience constructor for 3×3-matrices. The arguments are by+-- rows.+matrix3x3 :: (a, a, a) -> (a, a, a) -> (a, a, a) -> Matrix Three Three a+matrix3x3 (a0, a1, a2) (b0, b1, b2) (c0, c1, c2) = +  matrix_of_columns [[a0, b0, c0], [a1, b1, c1], [a2, b2, c2]]++-- | A convenience constructor for 4×4-matrices. The arguments are by+-- rows.+matrix4x4 :: (a, a, a, a) -> (a, a, a, a) -> (a, a, a, a) -> (a, a, a, a) -> Matrix Four Four a+matrix4x4 (a0, a1, a2, a3) (b0, b1, b2, b3) (c0, c1, c2, c3) (d0, d1, d2, d3) = +  matrix_of_columns [[a0, b0, c0, d0], [a1, b1, c1, d1], [a2, b2, c2, d2], [a3, b3, c3, d3]]++-- | A convenience constructor for 3-dimensional column vectors.+column3 :: (a, a, a) -> Matrix Three One a+column3 (a, b, c) = matrix_of_columns [[a, b, c]]++-- | A convenience destructor for 3-dimensional column vectors. This+-- is the inverse of 'column3'.+from_column3 :: Matrix Three One a -> (a, a, a)+from_column3 (Matrix ((a `Cons` b `Cons` c `Cons` Nil) `Cons` Nil)) = (a, b, c)++-- | A convenience constructor for turning a vector into a column matrix.+column_matrix :: Vector n a -> Matrix n One a+column_matrix v = Matrix (vector_singleton v)++-- ----------------------------------------------------------------------+-- ** Particular matrices++-- | Controlled-not gate.+cnot :: (Num a) => Matrix Four Four a+cnot = matrix4x4 (1,0,0,0)+                 (0,1,0,0)+                 (0,0,0,1)+                 (0,0,1,0)++-- | Swap gate.+swap :: (Num a) => Matrix Four Four a+swap = matrix4x4 (1,0,0,0)+                 (0,0,1,0)+                 (0,1,0,0)+                 (0,0,0,1)++-- | A /z/-rotation gate, /R/[sub /z/](θ) = [exp −/i/θ/Z/\/2].+zrot :: (Eq r, Floating r, Adjoint r) => r -> Matrix Two Two (Cplx r)+zrot theta = matrix2x2 (u, 0)+                       (0, adj u)+  where+    u = Cplx (cos (theta/2)) (-sin (theta/2))
+ Quantum/Synthesis/MultiQubitSynthesis.hs view
@@ -0,0 +1,468 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE GADTs #-}++-- | This module provides functions for the representation and exact+-- synthesis of multi-qubit Clifford+/T/ operators. +-- +-- The multi-qubit Clifford+/T/ exact synthesis algorithm is described+-- in the paper:+-- +-- * Brett Giles, Peter Selinger. Exact synthesis of multiqubit Clifford+T+-- circuits. /Physical Review A/ 87, 032332 (7 pages), 2013. Available+-- from <http://arxiv.org/abs/1212.0506>.+-- +-- It generalizes the single-qubit exact synthesis algorithm of+-- Kliuchnikov, Maslov, and Mosca.++module Quantum.Synthesis.MultiQubitSynthesis where++import Quantum.Synthesis.Matrix+import Quantum.Synthesis.Ring+import Data.List++-- ----------------------------------------------------------------------+-- * Residues++-- | A type class for things that have residues. In a typical+-- instance, /a/ is a ring whose elements are expressed with+-- coefficients in ℤ, and /b/ is a corresponding ring whose elements+-- are expressed with coefficients in ℤ[sub 2].+class Residue a b | a -> b where+  -- | Return the residue of something.+  residue :: a -> b+  +instance Residue Integer Z2 where+  residue = parity+  +instance Residue a b => Residue (Omega a) (Omega b) where+  residue (Omega a b c d) = Omega (residue a) (residue b) (residue c) (residue d)++instance Residue a b => Residue (RootTwo a) (RootTwo b) where+  residue (RootTwo a b) = RootTwo (residue a) (residue b)+  +instance (Residue a a', Residue b b') => Residue (a,b) (a',b') where+  residue (x,y) = (residue x, residue y)+  +instance Residue () () where  +  residue = const ()+  +instance (Residue a b) => Residue [a] [b] where  +  residue = map residue+  +instance (Residue a b) => Residue (Cplx a) (Cplx b) where  +  residue (Cplx a b) = Cplx (residue a) (residue b)+  +instance (Residue a b) => Residue (Vector n a) (Vector n b) where  +  residue = vector_map residue+  +instance (Residue a b) => Residue (Matrix m n a) (Matrix m n b) where+  residue (Matrix m) = Matrix (residue m)+  +-- ----------------------------------------------------------------------+-- * One- and two-level operators+  +-- ----------------------------------------------------------------------  +-- ** Symbolic representation++-- | An index for a row or column of a matrix.+type Index = Int++-- | Symbolic representation of one- and two-level operators. Note+-- that the power /k/ in the 'TL_T' and 'TL_omega' constructors can be+-- positive or negative, and should be regarded modulo 8.+-- +-- Note: when we use a list of 'TwoLevel' operators to express a+-- sequence of operators, the operators are meant to be applied+-- right-to-left, i.e., as in the mathematical notation for matrix+-- multiplication. This is the opposite of the quantum circuit+-- notation.+data TwoLevel = +  TL_X Index Index -- ^ /X/[sub /i/,/j/].+  | TL_H Index Index -- ^ /H/[sub /i/,/j/].+  | TL_T Int Index Index -- ^ (/T/[sub /i/,/j/])[super /k/].+  | TL_omega Int Index -- ^ (ω[sub /i/])[super /k/].+  deriving (Show, Eq)++-- | Invert a 'TwoLevel' operator.+invert_twolevel :: TwoLevel -> TwoLevel+invert_twolevel (TL_X i j) = TL_X i j+invert_twolevel (TL_H i j) = TL_H i j+invert_twolevel (TL_T m i j) = TL_T (-m) i j+invert_twolevel (TL_omega m j) = TL_omega (-m) j++-- | Invert a list of 'TwoLevel' operators.+invert_twolevels :: [TwoLevel] -> [TwoLevel]+invert_twolevels = reverse . map invert_twolevel++-- ----------------------------------------------------------------------+-- ** Constructors for two-level matrices++-- | Construct a two-level matrix with the given entries.+twolevel_matrix :: (Ring a, Nat n) => (a,a) -> (a,a) -> Index -> Index -> Matrix n n a+twolevel_matrix (a,b) (c,d) i j = matrix_of_function f where+  f x y +    | x == i && y == i = a+    | x == i && y == j = b+    | x == j && y == i = c+    | x == j && y == j = d+    | x == y = 1+    | otherwise = 0++-- | Construct a one-level matrix with the given entry.+onelevel_matrix :: (Ring a, Nat n) => a -> Index -> Matrix n n a+onelevel_matrix a i = matrix_of_function f where+  f x y+    | x == i && y == i = a+    | x == y = 1+    | otherwise = 0++-- | Convert a symbolic one- or two-level operator into a matrix.+matrix_of_twolevel :: (ComplexRing a, RootHalfRing a, Nat n) => TwoLevel -> Matrix n n a+matrix_of_twolevel (TL_X i j) = twolevel_matrix (0,1) (1,0) i j+matrix_of_twolevel (TL_H i j) = twolevel_matrix (s,s) (s,-s) i j+  where s = roothalf+matrix_of_twolevel (TL_T k i j) = twolevel_matrix (1,0) (0,omega^(k `mod` 8)) i j+matrix_of_twolevel (TL_omega k i) = onelevel_matrix (omega^(k `mod` 8)) i++-- | Convert a list of symbolic one- or two-level operators into a+-- matrix. Note that the operators are to be applied right-to-left,+-- exactly as in mathematical notation.+matrix_of_twolevels :: (ComplexRing a, RootHalfRing a, Nat n) => [TwoLevel] -> Matrix n n a+matrix_of_twolevels gs = foldl' (*) 1 [ matrix_of_twolevel g | g <- gs ]++-- ----------------------------------------------------------------------+-- * Auxiliary list functions++-- | Replace the /i/th element of a list by /x/.+list_insert :: Index -> a -> [a] -> [a]+list_insert 0 x (h:t) = x:t+list_insert n x (h:t) = h:(list_insert (n-1) x t)+list_insert n x [] = []++-- | Apply a unary operator to element /i/ of a list.+transform_at :: (a -> a) -> Index -> [a] -> [a]+transform_at op i lst = lst' where+  x = lst !! i+  x' = op x+  lst' = list_insert i x' lst++-- | Apply a binary operator to elements /i/ and /j/ of a list.+transform_at2 :: ((a,a) -> (a,a)) -> Index -> Index -> [a] -> [a]+transform_at2 op i j lst = lst' where+  (x,y) = (lst !! i, lst !! j)+  (x',y') = op (x,y)+  lst' = list_insert i x' (list_insert j y' lst)++-- | Split a list into pairs. Return a list of pairs, and a final+-- element if the length of the list was odd.+list_pairs :: [a] -> ([(a,a)], Maybe a)+list_pairs [] = ([], Nothing)+list_pairs [h] = ([], Just h)+list_pairs (h:k:t) = ((h,k):t',r') where (t',r') = list_pairs t++-- ----------------------------------------------------------------------+-- * Functions on ℤ[ω]++-- | Given an element of the form ω[sup /m/], return /m/ ∈ {0,…,7}, or+-- 'Nothing' if not of that form.+log_omega :: ZOmega -> Maybe Int+log_omega (Omega 0 0 0 1) = Just 0+log_omega (Omega 0 0 1 0) = Just 1+log_omega (Omega 0 1 0 0) = Just 2+log_omega (Omega 1 0 0 0) = Just 3+log_omega (Omega 0 0 0 (-1)) = Just 4+log_omega (Omega 0 0 (-1) 0) = Just 5+log_omega (Omega 0 (-1) 0 0) = Just 6+log_omega (Omega (-1) 0 0 0) = Just 7+log_omega _ = Nothing++-- | Multiply a scalar by ω[sup /n/].+omega_power :: (OmegaRing a) => Int -> a -> a+omega_power n x = x * omega^(n `mod` 8)++-- | Divide an element of 'ZOmega' by √2, or throw an error if it is+-- not divisible.+reduce_ZOmega :: ZOmega -> ZOmega+reduce_ZOmega (Omega a b c d) +  | even (a-c) && even (b-d) = Omega a' b' c' d'+  | otherwise = error "reduce_ZOmega: element not reducible"+  where+    a' = (b-d) `div` 2+    b' = (c+a) `div` 2+    c' = (b+d) `div` 2+    d' = (c-a) `div` 2++-- | Apply the /X/ operator to a 2-dimensional vector over 'ZOmega'.+opX_zomega :: (ZOmega, ZOmega) -> (ZOmega, ZOmega)+opX_zomega (x,y) = (y,x)++-- | Apply the /H/ operator to a 2-dimensional vector over+-- 'ZOmega'. This throws an error if the result is not well-defined+-- over 'ZOmega'.+opH_zomega :: (ZOmega, ZOmega) -> (ZOmega, ZOmega)+opH_zomega (x,y) = (reduce_ZOmega (x+y), reduce_ZOmega (x-y))++-- | Apply a 'TwoLevel' operator to a 'ZOmega'-vector, represented as+-- a list. Throws an error if any operation produces a scalar that is+-- not in 'ZOmega'.+apply_twolevel_zomega :: TwoLevel -> [ZOmega] -> [ZOmega]+apply_twolevel_zomega (TL_X i j) w = transform_at2 opX_zomega i j w+apply_twolevel_zomega (TL_H i j) w = transform_at2 opH_zomega i j w+apply_twolevel_zomega (TL_T k i j) w = transform_at (omega_power k) j w+apply_twolevel_zomega (TL_omega k i) w = transform_at (omega_power k) i w++-- | Apply a list of 'TwoLevel' operators to a 'ZOmega'-vector,+-- represented as a list. Throws an error if any operation produces a+-- scalar that is not in 'ZOmega'.+apply_twolevels_zomega :: [TwoLevel] -> [ZOmega] -> [ZOmega]+apply_twolevels_zomega gs w = foldr apply_twolevel_zomega' w gs+  where apply_twolevel_zomega' g w = apply_twolevel_zomega g w++-- ----------------------------------------------------------------------+-- * Functions on residues++-- | The /residue type/ of /t/ ∈ ℤ[ω] is the residue of /t/[sup †]/t/.+-- It is 0000, 0001, or 1010.+data ResidueType = RT_0000 | RT_0001 | RT_1010+                                       deriving (Eq, Ord)++-- | Return the residue's 'ResidueType'.+residue_type :: Omega Z2 -> ResidueType+residue_type r = t where+  (t, _) = residue_type_shift r+  +-- | Return the residue's /shift/.+-- +-- The shift is defined so that: +-- +-- * 0001, 1110, 0011 have shift 0,+-- +-- * 0010, 1101, 0110 have shift 1,+-- +-- * 0100, 1011, 1100 have shift 2, and+-- +-- * 1000, 0111, 1001 have shift 3.+-- +-- Residues of type 'RT_0000' have shift 0.+residue_shift :: Omega Z2 -> Int+residue_shift r = s where+  (_, s) = residue_type_shift r++-- | Return the residue's 'ResidueType' and the shift.+residue_type_shift :: Omega Z2 -> (ResidueType, Int)+residue_type_shift (Omega 0 0 0 0) = (RT_0000, 0)+residue_type_shift (Omega 0 0 0 1) = (RT_0001, 0)+residue_type_shift (Omega 0 0 1 0) = (RT_0001, 1)+residue_type_shift (Omega 0 0 1 1) = (RT_1010, 0)+residue_type_shift (Omega 0 1 0 0) = (RT_0001, 2)+residue_type_shift (Omega 0 1 0 1) = (RT_0000, 0)+residue_type_shift (Omega 0 1 1 0) = (RT_1010, 1)+residue_type_shift (Omega 0 1 1 1) = (RT_0001, 3)+residue_type_shift (Omega 1 0 0 0) = (RT_0001, 3)+residue_type_shift (Omega 1 0 0 1) = (RT_1010, 3)+residue_type_shift (Omega 1 0 1 0) = (RT_0000, 0)+residue_type_shift (Omega 1 0 1 1) = (RT_0001, 2)+residue_type_shift (Omega 1 1 0 0) = (RT_1010, 2)+residue_type_shift (Omega 1 1 0 1) = (RT_0001, 1)+residue_type_shift (Omega 1 1 1 0) = (RT_0001, 0)+residue_type_shift (Omega 1 1 1 1) = (RT_0000, 0)+residue_type_shift _ = undefined  -- to turn off a compiler warning++-- | Given two irreducible residues /a/ and /b/ of the same type, find+-- an index /m/ such that /a/ + ω[sup /m/]/b/ = 0000. If no such index+-- exists, find an index /m/ such that /a/ + ω[sup /m/]/b/ = 1111.+residue_offset :: Omega Z2 -> Omega Z2 -> Int+residue_offset a b = (residue_shift a - residue_shift b) `mod` 4++-- | Check whether a residue is reducible. A residue /r/ is called /reducible/+-- if it is of the form /r/ = √2 ⋅ /r/', i.e., /r/ ∈ {0000, 0101, 1010, 1111}.+reducible :: Omega Z2 -> Bool+reducible (Omega a b c d) = (a==c) && (b==d)++-- ----------------------------------------------------------------------+-- * Exact synthesis++-- | Perform a single row operation as in Lemma 4, applied to rows /i/+-- and /j/.  The entries at rows /i/ and /j/ are /x/ and /y/,+-- respectively, with respective residues /a/ and /b/. A precondition+-- is that /x/ and /y/ are of the same residue type. Returns a list of+-- two-level operations that decreases the denominator exponent.+row_step :: ((Index, Omega Z2, ZOmega), (Index, Omega Z2, ZOmega)) -> [TwoLevel]+row_step ((i,a,x), (j,b,y))+  | reducible a && reducible b = []+  | offs /= 0 = (TL_T offs i j) : row_step ((i,a,x), (j,b',y'))+  | otherwise = (TL_H i j) : row_step ((i,a1,x1), (j,b1,y1))+  where+    offs = residue_offset b a+    y' = omega_power (-offs) y+    b' = residue y'+    (x1,y1) = opH_zomega (x,y)+    (a1,b1) = residue (x1,y1)++-- | Row reduction: Given a unit column vector /v/, generate a+-- sequence of two-level operators that reduces the /i/th standard+-- basis vector /e/[sub /i/] to /v/. Any rows that are already 0 in+-- both vectors are guaranteed not to be touched.+reduce_column :: (Nat n) => Matrix n One (DOmega) -> Index -> [TwoLevel]+reduce_column v i = aux w k where+  vlist = list_of_vector (vector_head (unMatrix v))+  (w, k) = denomexp_decompose vlist+  aux w 0 = m1 ++ m2 where+    j = case findIndices (/= 0) w of+      [j] -> j+      _ -> error "reduce_column: not a unit vector"+    wj = w !! j+    l = case log_omega wj of+      Just l -> l+      Nothing -> error "reduce_column: not a unit vector"+    m1 = if i==j then [] else [TL_X i j]+    m2 = [TL_omega l i]+  aux w k = gates ++ aux w' (k-1) where+    res = residue w+    idx_res = zip3 [0..] res w+    res1010 = [ (i,a,x) | (i,a,x) <- idx_res, residue_type a == RT_1010 ]+    res0001 = [ (i,a,x) | (i,a,x) <- idx_res, residue_type a == RT_0001 ]+    res1010_pairs = case list_pairs res1010 of+      (p, Nothing) -> p+      _ -> error "reduce_column: not a unit vector"+    res0001_pairs = case list_pairs res0001 of+      (p, Nothing) -> p+      _ -> error "reduce_column: not a unit vector"+    m1010 = concat $ map row_step res1010_pairs+    m0001 = concat $ map row_step res0001_pairs+    gates = m1010 ++ m0001+    w' = map (reduce_ZOmega) (apply_twolevels_zomega (invert_twolevels gates) w)++-- | Input an exact /n/×/n/ unitary operator with coefficients in+-- [bold D][ω], and output an equivalent sequence of two-level+-- operators.  This is the algorithm from the Giles-Selinger paper. It+-- has superexponential complexity.+-- +-- Note: the list of 'TwoLevel' operators will be returned in+-- right-to-left order, i.e., as in the mathematical notation for+-- matrix multiplication. This is the opposite of the quantum circuit+-- notation.+synthesis_nqubit :: (Nat n) => Matrix n n DOmega -> [TwoLevel]+synthesis_nqubit m = aux (unMatrix m) 0 where+  aux :: (Nat m) => Vector n (Vector m DOmega) -> Index -> [TwoLevel]+  aux Nil i = []+  aux (c `Cons` cs) i = gates ++ aux (unMatrix m') (i+1)+    where+      gates = reduce_column (column_matrix c) i+      gates_matrix = matrix_of_twolevels (invert_twolevels gates)+      m' = gates_matrix .*. (Matrix cs)++-- ----------------------------------------------------------------------+-- * Alternative algorithm+      +-- $ Section 6 of the Giles-Selinger paper mentions an alternate+-- version of the decomposition algorithm. It requires no ancillas,+-- provided that the determinant of the operator permits this.+      +-- | Symbolic representation of one- and two-level operators, with an+-- alternate set of generators.+-- +-- Note: when we use a list of 'TwoLevel' operators to express a+-- sequence of operators, the operators are meant to be applied+-- right-to-left, i.e., as in the mathematical notation for matrix+-- multiplication. This is the opposite of the quantum circuit+-- notation.+data TwoLevelAlt =+  TL_iX Index Index -- ^ /iX/[sub /i/,/j/].+  | TL_TiHT Int Index Index -- ^ (/T/[super −/m/](iH)T[super /m/])[sub /i/,/j/].+  | TL_W Int Index Index -- ^ /W/[super /m/][sub /i/,/j/].+  | TL_omega_alt Int Index -- ^ (ω[sub /i/])[super /m/].+  deriving (Show, Eq)++-- | Convert from the alternate generators to the original generators.+twolevels_of_twolevelalts :: [TwoLevelAlt] -> [TwoLevel]+twolevels_of_twolevelalts [] = []+twolevels_of_twolevelalts (TL_iX j l : t) = +  TL_X j l : TL_omega 2 j : TL_omega 2 l : twolevels_of_twolevelalts t+twolevels_of_twolevelalts (TL_TiHT m j l : t) =+  TL_T (-m) j l : TL_H j l : TL_omega 2 j : TL_omega 2 l : TL_T m j l : twolevels_of_twolevelalts t+twolevels_of_twolevelalts (TL_W m j l : t) =+  TL_omega m j : TL_omega (-m) l : twolevels_of_twolevelalts t+twolevels_of_twolevelalts (TL_omega_alt m j : t) =+  TL_omega m j : twolevels_of_twolevelalts t++-- | Invert a list of 'TwoLevelAlt' operators, and convert the output+-- to a list of 'TwoLevel' operators.+invert_twolevels_alt :: [TwoLevelAlt] -> [TwoLevel]+invert_twolevels_alt = invert_twolevels . twolevels_of_twolevelalts++-- | Perform a single row operation as in Lemma 4, applied to rows /i/+-- and /j/, using the generators of Section 6.  The entries at rows+-- /i/ and /j/ are /x/ and /y/, respectively, with respective residues+-- /a/ and /b/. A precondition is that /x/ and /y/ are of the same+-- residue type. Returns a list of two-level operations that decreases+-- the denominator exponent.+row_step_alt :: ((Index, Omega Z2, ZOmega), (Index, Omega Z2, ZOmega)) -> [TwoLevelAlt]+row_step_alt ((j,a,x), (l,b,y))+  | reducible a && reducible b = []+  | otherwise = (TL_TiHT m j l) : row_step_alt ((j,a1,x1), (l,b1,y1))+  where+    m = residue_offset a b+    y' = omega_power m y+    (x1,y1') = opH_zomega (-i*x,-i*y')+    y1 = omega_power (-m) y1'+    (a1,b1) = residue (x1,y1)++-- | Row reduction: Given a unit column vector /v/, generate a+-- sequence of two-level operators that reduces the /i/th standard+-- basis vector /e/[sub /i/] to /v/. Any rows that are already 0 in+-- both vectors are guaranteed not to be touched, except possibly row+-- /i/+1 may be multiplied by a scalar.+reduce_column_alt :: (Nat n) => Matrix n One (DOmega) -> Index -> [TwoLevelAlt]+reduce_column_alt v j = aux w k where+  vlist = list_of_vector (vector_head (unMatrix v))+  n = length vlist+  (w, k) = denomexp_decompose vlist+  aux w 0 = m1 ++ m2 where+    l = case findIndices (/= 0) w of+      [l] -> l+      _ -> error "reduce_column: not a unit vector"+    m1 = if j==l then [] else [TL_iX j l]+    wl = if j==l then w !! j else -i*(w !! l)+    m = case log_omega wl of+      Just m -> m+      Nothing -> error "reduce_column: not a unit vector"+    m2 = if j==n-1 then [TL_omega_alt m j] else [TL_W m j (j+1)]+  aux w k = gates ++ aux w' (k-1) where+    res = residue w+    idx_res = zip3 [0..] res w+    res1010 = [ (j,a,x) | (j,a,x) <- idx_res, residue_type a == RT_1010 ]+    res0001 = [ (j,a,x) | (j,a,x) <- idx_res, residue_type a == RT_0001 ]+    res1010_pairs = case list_pairs res1010 of+      (p, Nothing) -> p+      _ -> error "reduce_column: not a unit vector"+    res0001_pairs = case list_pairs res0001 of+      (p, Nothing) -> p+      _ -> error "reduce_column: not a unit vector"+    m1010 = concat $ map row_step_alt res1010_pairs+    m0001 = concat $ map row_step_alt res0001_pairs+    gates = m1010 ++ m0001+    w' = map reduce_ZOmega (apply_twolevels_zomega (invert_twolevels_alt gates) w)++-- | Input an exact /n/×/n/ unitary operator with coefficients in+-- [bold D][ω], and output an equivalent sequence of two-level+-- operators (in the alternative generators, where all but at most one+-- of the generators has determinant 1).  This is the algorithm from+-- the Giles-Selinger paper, Section 6. It has superexponential+-- complexity.+-- +-- Note: the list of 'TwoLevelAlt' operators will be returned in+-- right-to-left order, i.e., as in the mathematical notation for+-- matrix multiplication. This is the opposite of the quantum circuit+-- notation.+synthesis_nqubit_alt :: (Nat n) => Matrix n n DOmega -> [TwoLevelAlt]+synthesis_nqubit_alt m = aux (unMatrix m) 0 where+  aux :: (Nat m) => Vector n (Vector m DOmega) -> Index -> [TwoLevelAlt]+  aux Nil i = []+  aux (c `Cons` cs) i = gates ++ aux (unMatrix m') (i+1)+    where+      gates = reduce_column_alt (column_matrix c) i+      gates_matrix = matrix_of_twolevels (invert_twolevels_alt gates)+      m' = gates_matrix .*. (Matrix cs)
+ Quantum/Synthesis/Newsynth.hs view
@@ -0,0 +1,453 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE ScopedTypeVariables #-}++-- | This module implements an efficient single-qubit Clifford+/T/+-- approximation algorithm. The algorithm is described here:+-- +-- * Peter Selinger. Efficient Clifford+/T/ approximation of+-- single-qubit operators. <http://arxiv.org/abs/1212.6253>.++module Quantum.Synthesis.Newsynth where++import Quantum.Synthesis.Ring+import Quantum.Synthesis.Ring.FixedPrec+import Quantum.Synthesis.Matrix+import Quantum.Synthesis.CliffordT+import Quantum.Synthesis.EuclideanDomain+import Quantum.Synthesis.SymReal++import System.Random+import Data.Number.FixedPrec++-- ----------------------------------------------------------------------+-- * Miscellaneous functions++-- | A useful operation for the 'Maybe' monad, used to ensure that+-- some condition holds (i.e., return 'Nothing' if the condition is+-- false). To be used like this:+-- +-- > do+-- >   x <- something+-- >   y <- something_else+-- >   ensure (x > y)+-- >   ...+ensure :: Bool -> Maybe ()+ensure True = Just ()+ensure False = Nothing++-- | Return the head of a list, if non-empty, or else 'Nothing'.+maybe_head :: [a] -> Maybe a+maybe_head [] = Nothing+maybe_head (h:t) = Just h++-- | Exponentiation via repeated squaring, parameterized by a+-- multiplication function and a unit. Given an associative+-- multiplication function @*@ with unit @e@, the function 'power'+-- @(*)@ /e/ /a/ /n/ efficiently computes /a/[sup /n/] = /a/ @*@ (/a/+-- @*@ (… @*@ (/a/ @*@ /e/)…)).+power :: (a -> a -> a) -> a ->  a -> Integer -> a+power mul unit = aux where+  aux x n+    | n <= 0 = unit+    | n == 1 = x+    | odd n = x `mul` (x `aux` (n-1))+    | otherwise = y `mul` y where y = x `aux` (n `div` 2)+  +-- | Given positive numbers /b/ and /x/, return (/n/, /r/) such that+-- +-- * /x/ = /r/ /b/[sup /n/] and                           +--                                   +-- * 1 ≤ /r/ < /b/.                                  +--                                   +-- In other words, let /n/ = ⌊log[sub /b/] /x/⌋ and +-- /r/ = /x/ /b/[sup −/n/]. This can be more efficient than 'floor'+-- ('logBase' /b/ /x/) depending on the type; moreover, it also works+-- for exact types such as 'Rational' and 'QRootTwo'.+floorlog :: (Fractional b, Ord b) => b -> b -> (Integer, b)+floorlog b x +    | x <= 0            = error "floorlog: argument not positive"    +    | 1 <= x && x < b   = (0, x)+    | 1 <= x*b && x < 1 = (-1, b*x)+    | r < b             = (2*n, r)+    | otherwise         = (2*n+1, r/b)+    where+      (n, r) = floorlog (b^2) x++-- ----------------------------------------------------------------------+-- * Randomized algorithms++-- | A combinator for turning a probabilistic function that succeeds+-- with some small probability into a probabilistic function that+-- always succeeds, by trying again and again.+keeptrying :: (RandomGen g) => (g -> Maybe a) -> (g -> a)+keeptrying f g = case f g1 of+  Just res -> res+  Nothing -> keeptrying f g2+  where+    (g1, g2) = split g++-- | Like 'keeptrying', but also returns a count of the number of attempts.+keeptrying_count :: (RandomGen g) => (g -> Maybe a) -> (g -> (a, Integer))+keeptrying_count f g = aux g 1 where+  aux g n = case f g1 of+    Just res -> (res, n)+    Nothing -> aux g2 n1+    where+      (g1, g2) = split g+      !n1 = n + 1++-- | A combinator for turning a probabilistic function that succeeds+-- with some small probability into a probabilistic function that+-- succeeds with a higher probability, by repeating it /n/ times. +try_for :: (RandomGen g) => Integer -> (g -> Maybe a) -> (g -> Maybe a)+try_for n f g+  | n <= 0 = Nothing+  | otherwise = case f g1 of+      Just res -> Just res+      Nothing -> try_for (n-1) f g2+  where+    (g1, g2) = split g    ++-- ----------------------------------------------------------------------+-- * Square roots in ℤ[√2]++-- | Return a square root of an element of ℤ[√2], if such a square+-- root exists, or else 'Nothing'.+zroottwo_root :: ZRootTwo -> Maybe ZRootTwo+zroottwo_root z@(RootTwo a b) = res where+  d = a^2 - 2*b^2+  r = intsqrt d+  x1 = intsqrt ((a + r) `div` 2)+  x2 = intsqrt ((a - r) `div` 2)+  y1 = intsqrt ((a - r) `div` 4)+  y2 = intsqrt ((a + r) `div` 4)+  w1 = RootTwo x1 y1+  w2 = RootTwo x2 y2+  w3 = RootTwo x1 (-y1)+  w4 = RootTwo x2 (-y2)+  res +    | w1*w1 == z = Just w1+    | w2*w2 == z = Just w2+    | w3*w3 == z = Just w3+    | w4*w4 == z = Just w4+    | otherwise  = Nothing+  +-- ----------------------------------------------------------------------  +-- * Roots of −1 in ℤ[sub /p/]+  +-- | Input an integer /p/, and maybe output a root of −1 modulo /p/.+-- This succeeds with probability at least 1\/2 if /p/ is a positive+-- prime ≡ 1 (mod 4); otherwise, the success probability is+-- unspecified (and may be 0).+root_minus_one_step :: (RandomGen g) => Integer -> g -> Maybe Integer+root_minus_one_step p g = do+  let (b, _) = randomR (1, p-1) g+  let h = power mul 1 b ((p-1) `div` 4)+  ensure $ h `mul` h == p-1  -- succeeds with probability 1/2+  return h+    where+      mul :: Integer -> Integer -> Integer+      mul a b = (a*b) `mod` p+      +-- | Input a positive prime /p/ ≡ 1 (mod 4), and output a root of −1.+root_minus_one :: (RandomGen g) => Integer -> g -> Integer+root_minus_one p = keeptrying (root_minus_one_step p)++-- ----------------------------------------------------------------------+-- * Solving a Diophantine equation++-- | Input ξ ∈ ℤ[√2], and maybe output some /t/ ∈ ℤ[ω] such that +-- /t/[sup †]/t/ = ξ. If ξ ≥ 0, ξ[sup •] ≥ 0 and /p/ = ξ[sup •]ξ is a+-- prime ≡ 1 (mod 4) in ℤ, then this succeeds with probability at least+-- 1\/2.  Otherwise, the success probability is unspecified and may be+-- 0.+dioph_step :: (RandomGen g) => ZRootTwo -> g -> Maybe ZOmega+dioph_step xi g = do+  h <- root_minus_one_step (norm xi) g+  let s = euclid_gcd (fromInteger h+i) (fromZRootTwo xi) :: ZOmega+      ss = zroottwo_of_zomega (adj s * s)+      u = euclid_div xi ss+  v <- zroottwo_root u+  let t = fromZRootTwo v * s+  ensure $ adj t * t == fromZRootTwo xi -- check the answer, just in case+  return t++-- | Input ξ ∈ ℤ[√2] such that ξ ≥ 0, ξ[sup •] ≥ 0, and /p/ = +-- ξ[sup •]ξ is a prime ≡ 1 (mod 4) in ℤ. Output /t/ ∈ ℤ[ω] such that+-- /t/[sup †]/t/ = ξ. If the hypotheses are not satisfied, this will+-- likely loop forever.+dioph :: (RandomGen g) => ZRootTwo -> g -> ZOmega+dioph xi = keeptrying (dioph_step xi)++-- ----------------------------------------------------------------------+-- * Approximations in ℤ[√2]++-- | Input two intervals [/x/₀, /x/₁] ⊆ ℝ and [/y/₀, /y/₁] ⊆ ℝ. Output+-- a list of all points /z/ = /a/ + √2/b/ ∈ ℤ[√2] such that /z/ ∈+-- [/x/₀, /x/₁] and /z/[sup •] ∈ [/y/₀, /y/₁]. The list will be+-- produced lazily, and will be sorted in order of increasing /z/.+-- +-- It is a theorem that there will be at least one solution if ΔxΔy ≥ (1+-- + √2)², and at most one solution if ΔxΔy < 1, where Δx = /x/₁ − /x/₀ ≥ 0+-- and Δy = /y/₁ − /y/₀ ≥ 0. Asymptotically, the expected number of+-- solutions is ΔxΔy/\√8.+-- +-- This function is formulated so that the intervals can be specified+-- exactly (using a type such as 'QRootTwo'), or approximately (using a+-- type such as 'Double' or 'FixedPrec' /e/).+gridpoints :: (RootTwoRing r, Fractional r, Floor r, Ord r) => (r, r) -> (r, r) -> [ZRootTwo]+gridpoints (x0, x1) (y0, y1)+  | dy <= 0 && dx > 0 = +        map adj2 $ gridpoints (y0, y1) (x0, x1)+  | dy >= lambda && even n =+        map (lambdainv^n *) $ gridpoints (lambda^n*x0, lambda^n*x1) (lambda'^n*y0, lambda'^n*y1)+  | dy >= lambda && odd n =+        map (lambdainv^n *) $ gridpoints (lambda^n*x0, lambda^n*x1) (lambda'^n*y1, lambda'^n*y0)+  | dy > 0 && dy < 1 && even n = +        map (lambda^m *) $ gridpoints (lambdainv^m*x0, lambdainv^m*x1) (lambdainv'^m*y0, lambdainv'^m*y1)+  | dy > 0 && dy < 1 && odd n = +        map (lambda^m *) $ gridpoints (lambdainv^m*x0, lambdainv^m*x1) (lambdainv'^m*y1, lambdainv'^m*y0)+  | otherwise =+        [ RootTwo a b | a <- [amin..amax], b <- [bmin a..bmax a], test a b ] +  where+    dx = x1 - x0+    dy = y1 - y0+    (n, _) = floorlog lambda dy+    m = -n+    +    lambda :: (RootTwoRing r) => r+    lambda = 1 + roottwo+    lambda' :: (RootTwoRing r) => r+    lambda' = 1 - roottwo+    lambdainv :: (RootTwoRing r) => r+    lambdainv = -1 + roottwo+    lambdainv' :: (RootTwoRing r) => r+    lambdainv' = -1 - roottwo++    within x (x0, x1) = x0 <= x && x <= x1+    amin = ceiling_of ((x0 + y0) / 2)+    amax = floor_of ((x1 + y1) / 2)+    bmin a = ceiling_of ((fromInteger a - y1) / roottwo)+    bmax a = floor_of ((fromInteger a - y0) / roottwo)+    test a b = fromZRootTwo x `within` (x0, x1) && fromZRootTwo (adj2 x) `within` (y0, y1)+      where x = RootTwo a b++-- | Input two intervals [/x/₀, /x/₁] ⊆ ℝ and [/y/₀, /y/₁] ⊆ ℝ and a+-- source of randomness. Output a random element /z/ = /a/ + √2/b/+-- ∈ ℤ[√2] such that /z/ ∈ [/x/₀, /x/₁] and /z/[sup •] ∈ [/y/₀,+-- /y/₁]. +-- +-- Note: the randomness will not be uniform. To ensure that the set of+-- solutions is non-empty, we must have ΔxΔy ≥ (1 + √2)², where Δx =+-- /x/₁ − /x/₀ ≥ 0 and Δy = /y/₁ − /y/₀ ≥ 0. If there are no solutions+-- at all, the function will return 'Nothing'.+-- +-- This function is formulated so that the intervals can be specified+-- exactly (using a type such as 'QRootTwo'), or approximately (using a+-- type such as 'Double' or 'FixedPrec' /e/).+gridpoint_random :: (RootTwoRing r, Fractional r, Floor r, Ord r, RandomGen g) => (r, r) -> (r, r) -> g -> Maybe ZRootTwo+gridpoint_random (x0, x1) (y0, y1) g = z+  where+    dx = max 0 (x1 - x0)+    dy = max 0 (y1 - y0)+    area = dx * dy+    n = floor_of (area + 1)+    (i,_) = randomR (0, n-1) g+    r = fromInteger i / fromInteger n+    pts = gridpoints (x0 + r * dx, x1) (y0, y1) ++ gridpoints (x0, x1) (y0, y1)+    z = maybe_head pts++-- | Input an integer /e/, two intervals [/x/₀, /x/₁] ⊆ ℝ and [/y/₀,+-- /y/₁] ⊆ ℝ, and a source of randomness. Output random /z/ = /a/ ++-- √2/b/ ∈ ℤ[√2] such that /a/ + √2/b/ ∈ [/x/₀, /x/₁], /a/ - √2/b/ ∈+-- [/y/₀, /y/₁], and /a/-/e/ is even.+-- +-- Note: the randomness will not be uniform. To ensure that the set of+-- solutions is non-empty, we must have ΔxΔy ≥ 2(√2 + 1)², where Δx =+-- /x/₁ − /x/₀ ≥ 0 and Δy = /y/₁ − /y/₀ ≥ 0. If there are no solutions+-- at all, the function will return 'Nothing'.+-- +-- This function is formulated so that the intervals can be specified+-- exactly (using a type such as 'QRootTwo'), or approximately (using a+-- type such as 'Double' or 'FixedPrec' /e/).+gridpoint_random_parity :: (RootTwoRing r, Fractional r, Floor r, Ord r, RandomGen g) => Integer -> (r, r) -> (r, r) -> g -> Maybe ZRootTwo+gridpoint_random_parity e (x0,x1) (y0,y1) g = do+  z' <- gridpoint_random (x0', x1') (-y1', -y0') g+  return (roottwo * z' + fromInteger e2)+  where +    x0' = (x0 - e') / roottwo+    x1' = (x1 - e') / roottwo+    y0' = (y0 - e') / roottwo+    y1' = (y1 - e') / roottwo+    e' = fromInteger e2+    e2 = e `mod` 2++-- ----------------------------------------------------------------------+-- * Approximate synthesis+  +-- ----------------------------------------------------------------------+-- ** The main algorithm++-- | The internal implementation of the approximate synthesis+-- algorithm. The parameters are:+-- +-- * an angle θ, to implement a /R/[sub /z/](θ) = [exp −/i/θ/Z/\/2]+-- gate;+--   +-- * a precision /p/ ≥ 0 in bits, such that ε = 2[sup -/p/];+-- +-- * a source of randomness /g/.+-- +-- With some probability, output a unitary operator in the+-- Clifford+/T/ group that approximates /R/[sub /z/](θ) to within ε in+-- the operator norm. This operator can then be converted to a list of+-- gates with 'to_gates'. Also output log[sub 0.1] of the actual+-- error, or 'Nothing' if the error is 0.+-- +-- This implementation does not use seeding.+-- +-- As a special case, if the /R/[sub /z/](θ) is a Clifford operator+-- (to within the given ε), always return this operator directly.+-- +-- Note: the parameter θ must be of a real number type that has enough+-- precision to perform intermediate calculations; this typically+-- requires precision O(ε[sup 2]).  A more user-friendly function that+-- does this automatically is 'newsynth'.+newsynth_step :: forall r g.(RealFrac r, Floating r, RootHalfRing r, Floor r, Adjoint r, RandomGen g) => r -> r -> g -> Maybe (U2 DOmega, Maybe Double)+newsynth_step prec theta = payload where+  -- We are careful to do all computations that depend only on epsilon+  -- and theta (but not g) outside of aux, to avoid re-computing them+  -- with each attempt.+  +  -- Calculate ε.+  epsilon = 2 ** (-prec)+  +  -- Convert prec to a Double+  dprec = fromRational (toRational prec)+  +  -- Determine k.+  const = 3 + 2 * logBase 2 (1 + sqrt 2) :: Double+  k = ceiling (const + 2 * dprec)+  scale = roottwo^k+  +  -- Normalize θ to be in [-π/4, π/4].+  n = round(theta / (pi/2))+  theta1 = theta - fromInteger n * pi/2+  +  -- Describe the ε-region.+  z @ (x,y) = (cos (theta1 / 2), -sin (theta1 / 2))+  e2 = 1 - epsilon^2/2+  e4 = 1 - epsilon^2/4+  z1 @ (x1,y1) = (e4 * x, e4 * y)+  e' = epsilon / roottwo+  f = e' * sqrt((1+e'/2)*(1-e'/2)) -- == sqrt(1-e4^2)+  w @ (wx,wy) = (-f * y, f * x)+  y_min = y1 - wy+  y_max = y1 + wy+  y'_min = y_min * scale+  y'_max = y_max * scale+  dx = (e4 - e2) * x+  +  find_uU_step = +    -- As a special case, if (1,0) is in the ε-region, return the+    -- identity operator.+    if x >= e2 then \g -> Just 1 else aux++  -- The rest of the computation depends on the random seed g.+  payload g = do+    uU1 <- find_uU_step g  +    let uU = correct uU1 n+    let err = calc_error uU theta+    return (uU, err)+  +  aux g = do+    -- Find a random grid point in the ε-region.+    let (g0,g1) = split g+    beta <- gridpoint_random (y'_min, y'_max) (-roothalf * scale, roothalf * scale) g0+    let  +      beta' = fromZRootTwo beta / scale+      tmp = (beta' - e2 * y) / wy+      x0 = e2 * x + tmp * wx+      x1 = x0 + dx+      x0' = x0 * scale+      x1' = x1 * scale+      (g2,g3) = split g1+      RootTwo c _ = beta+    alpha <- gridpoint_random_parity (c+1) (x0', x1') (-roothalf * scale, roothalf * scale) g2+    +    -- Calculate u, ξ, and solve Diophantine equation to calculate t.+    let  +      u = (fromZRootTwo alpha) + i * (fromZRootTwo beta) :: ZOmega+      xi = zroottwo_of_zomega (2^k - u * adj u)+    t <- dioph_step xi g3+    +    -- If Diophantine equation solved successfully, calculate matrix U.+    let+      u' = fromZOmega u * roothalf^k :: DOmega+      t' = fromZOmega t * roothalf^k :: DOmega+      uU1 = matrix2x2 (u', -(adj t'))+                      (t',  (adj u'))+           +    return uU1+    +  -- Correct for when θ wasn't in [-π/4, π/4].+  correct uU1 n = uU1 * rR^(n `mod` 8) where+    rR = matrix2x2 (omega^7, 0)+                   (0,   omega)+    +  -- Calculate the actual error. Since this is done lazily, this+  -- incurs no overhead in case the error is not actually used.+  calc_error uU theta = log_err where+    uU_fixed = matrix_map fromDOmega uU :: U2 (Cplx r)+    zrot_fixed = zrot theta :: U2 (Cplx r)+    err = sqrt (real (hs_sqnorm (uU_fixed - zrot_fixed)) / 2)+    log_err +      | err <= 0  = Nothing+      | otherwise = Just (log_double err / log 0.1)++-- ----------------------------------------------------------------------+-- ** User-friendly functions++-- | A user-friendly interface to the approximate synthesis+-- algorithm. The parameters are:+-- +-- * an angle θ, to implement a /R/[sub /z/](θ) = [exp −/i/θ/Z/\/2]+-- gate;+--   +-- * a precision /b/ ≥ 0 in bits, such that ε = 2[sup -/b/];+-- +-- * a source of randomness /g/.+-- +-- Output a unitary operator in the Clifford+/T/ group that+-- approximates /R/[sub /z/](θ) to within ε in the operator norm. This+-- operator can then be converted to a list of gates with+-- 'to_gates'.+-- +-- This implementation does not use seeding.+-- +-- Note: the argument /theta/ is given as a symbolic real number. It+-- will automatically be expanded to as many digits as are necessary+-- for the internal calculation. In this way, the caller can specify,+-- e.g., an angle of 'pi'\/128 @::@ 'SymReal', without having to worry+-- about how many digits of π to specify.+newsynth :: (RandomGen g) => Double -> SymReal -> g -> U2 DOmega+newsynth prec theta g = m where+  (m, _, _) = newsynth_stats prec theta g++-- | A version of 'newsynth' that also returns some statistics:+-- log[sub 0.1] of the actual approximation error (or 'Nothing' if the+-- error is 0), and the number of candidates tried.+newsynth_stats :: (RandomGen g) => Double -> SymReal -> g -> (U2 DOmega, Maybe Double, Integer)+newsynth_stats prec theta g = dynamic_fixedprec2 digits f prec theta where+  digits = ceiling (10 + 2 * prec * logBase 10 2)+  f prec theta = (m, err, ct) where+    ((m, err), ct) = keeptrying_count (newsynth_step prec theta) g++-- | A version of 'newsynth' that returns a list of gates instead of a+-- matrix. The inputs are the same as for 'newsynth'.+-- +-- Note: the list of gates will be returned in right-to-left order,+-- i.e., as in the mathematical notation for matrix multiplication.+-- This is the opposite of the quantum circuit notation.+newsynth_gates :: (RandomGen g) => Double -> SymReal -> g -> [Gate]+newsynth_gates prec theta g = synthesis_u2 (newsynth prec theta g)
+ Quantum/Synthesis/Ring.hs view
@@ -0,0 +1,1046 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE IncoherentInstances #-}+{-# LANGUAGE BangPatterns #-}++-- | This module provides type classes for rings. It also provides+-- several specific instances of rings, such as the ring ℤ₂ of+-- integers modulo 2, the ring ℚ of rational numbers, the ring ℤ[½] of+-- dyadic fractions, the ring ℤ[/i/] of Gaussian integers, the ring+-- ℤ[√2] of quadratic integers with radix 2, and the ring ℤ[ω] of+-- cyclotomic integers of degree 8.++module Quantum.Synthesis.Ring where++import Data.Bits+import Data.Complex+import Data.Ratio++-- ----------------------------------------------------------------------+-- * Rings++-- | A type class to denote rings. We make 'Ring' a synonym of+-- Haskell's 'Num' type class, so that we can use the usual notation+-- '+', '-', '*' for the ring operations.  This is not a perfect fit,+-- because Haskell's 'Num' class also contains two non-ring operations+-- 'abs' and 'signum'.  By convention, for rings where these notions+-- don't make sense (or are inconvenient to define), we set 'abs' /x/+-- = /x/ and 'signum' /x/ = 1.++class (Num a) => Ring a+instance (Num a) => Ring a++-- ----------------------------------------------------------------------+-- * Rings with particular elements++-- $ We define several classes of rings with special elements.++-- ----------------------------------------------------------------------+-- ** Rings with ½++-- | A type class for rings that contain ½.+class (Ring a) => HalfRing a where+  -- | The value ½.+  half :: a++instance HalfRing Double where+  half = 0.5++instance HalfRing Float where+  half = 0.5++instance HalfRing Rational where+  half = 1%2++instance (HalfRing a, RealFloat a) => HalfRing (Complex a) where+  half = half :+ 0++-- ----------------------------------------------------------------------+-- ** Rings with √2++-- | A type class for rings that contain √2.+class (Ring a) => RootTwoRing a where+  -- | The square root of 2.+  roottwo :: a+  +instance RootTwoRing Double where+  roottwo = sqrt 2++instance RootTwoRing Float where+  roottwo = sqrt 2++instance (RootTwoRing a, RealFloat a) => RootTwoRing (Complex a) where+  roottwo = roottwo :+ 0++-- ----------------------------------------------------------------------+-- ** Rings with 1\/√2++-- | A type class for rings that contain 1\/√2.+class (HalfRing a, RootTwoRing a) => RootHalfRing a where+  -- | The square root of ½.+  roothalf :: a+  +instance RootHalfRing Double where+  roothalf = sqrt 0.5++instance RootHalfRing Float where+  roothalf = sqrt 0.5++instance (RootHalfRing a, RealFloat a) => RootHalfRing (Complex a) where+  roothalf = roothalf :+ 0+++-- ----------------------------------------------------------------------+-- ** Rings with /i/++-- | A type class for rings that contain a square root of -1.+class (Ring a) => ComplexRing a where+  -- | The complex unit.+  i :: a+       +instance (Ring a, RealFloat a) => ComplexRing (Complex a) where+  i = 0 :+ 1++-- ----------------------------------------------------------------------+-- ** Rings with ω++-- | A type class for rings that contain a square root of /i/, or+-- equivalently, a fourth root of −1.+class (Ring a) => OmegaRing a where+  -- | The square root of /i/.+  omega :: a+  +instance (ComplexRing a, RootHalfRing a) => OmegaRing a where+  omega = roothalf * (1 + i)++-- ----------------------------------------------------------------------+-- * Rings with particular automorphisms++-- ----------------------------------------------------------------------+-- ** Rings with complex conjugation++-- | A type class for rings with complex conjugation, i.e., an+-- automorphism mapping /i/ to −/i/. +-- +-- When instances of this type class are vectors or matrices, the+-- conjugation also exchanges the roles of rows and columns (in other+-- words, it is the adjoint).+-- +-- For rings that are not complex, the conjugation can be defined to+-- be the identity function.+class Adjoint a where+  -- | Compute the adjoint (complex conjugate transpose).+  adj :: a -> a++instance Adjoint Integer where+  adj x = x+  +instance Adjoint Int where+  adj x = x+  +instance Adjoint Double where+  adj x = x+  +instance Adjoint Float where+  adj x = x+  +instance Adjoint Rational where  +  adj x = x+  +instance (Adjoint a, Ring a) => Adjoint (Complex a) where+  adj (a :+ b) = adj a :+ (-adj b)++-- ----------------------------------------------------------------------+-- ** Rings with √2-conjugation++-- | A type class for rings with a √2-conjugation, i.e., an+-- automorphism mapping √2 to −√2. +-- +-- When instances of this type class are vectors or matrices, the+-- √2-conjugation does /not/ exchange the roles of rows and columns.+-- +-- For rings that have no √2, the conjugation can be defined to be the+-- identity function.+class Adjoint2 a where+  -- | Compute the adjoint, mapping /a/ + /b/√2 to /a/ −/b/√2.+  adj2 :: a -> a++instance Adjoint2 Integer where+  adj2 x = x++instance Adjoint2 Int where+  adj2 x = x+  +instance Adjoint2 Rational where  +  adj2 x = x+  +-- ----------------------------------------------------------------------+-- * Normed rings++-- | A (number-theoretic) /norm/ on a ring /R/ is a function /N/ : /R/+-- → ℤ such that /N/(/rs/) = /N/(/r/)/N/(/s/), for all /r/, /s/ ∈ /R/.+-- The norm also satisfies /N/(/r/) = 0 iff /r/ = 0, and /N/(/r/) = ±1+-- iff /r/ is a unit of the ring.+class (Ring r) => NormedRing r where+  norm :: r -> Integer+  +instance NormedRing Integer where+  norm x = x+  +-- ----------------------------------------------------------------------+-- * Floor and ceiling+  +-- | The 'floor' and 'ceiling' functions provided by the standard+-- Haskell libraries are predicated on many unnecessary assumptions.+-- This type class provides an alternative.+-- +-- Minimal complete definition: 'floor_of' or 'ceiling_of'.+class (Ring r) => Floor r where+  -- | Compute the floor of /x/, i.e., the greatest integer /n/ such+  -- that /n/ ≤ /x/.+  floor_of :: r -> Integer+  floor_of x = -(ceiling_of (-x))+  -- | Compute the ceiling of /x/, i.e., the least integer /n/ such+  -- that /x/ ≤ /n/.+  ceiling_of :: r -> Integer+  ceiling_of x = -(floor_of (-x))++instance Floor Integer where+  floor_of = id+  ceiling_of = id++instance Floor Rational where+  floor_of = floor+  ceiling_of = ceiling++instance Floor Float where+  floor_of = floor+  ceiling_of = ceiling++instance Floor Double where+  floor_of = floor+  ceiling_of = ceiling++-- ----------------------------------------------------------------------+-- * Particular rings++-- ----------------------------------------------------------------------+-- ** The ring ℤ₂ of integers modulo 2++-- | The ring ℤ₂ of integers modulo 2. +data Z2 = Even | Odd+        deriving (Eq)+                     +instance Show Z2 where+  show Even = "0"+  show Odd = "1"++instance Num Z2 where+  Even + x = x+  x + Even = x+  Odd + Odd = Even+  Even * x = Even+  x * Even = Even+  Odd * Odd = Odd+  negate x = x+  fromInteger n = if even n then Even else Odd+  abs x = x+  signum x = 1++instance Adjoint Z2 where+  adj x = x++instance Adjoint2 Z2 where+  adj2 x = x++-- ----------------------------------------------------------------------+-- ** The ring [bold D] of dyadic fractions++-- | A dyadic fraction is a rational number whose denominator is a+-- power of 2. We denote the dyadic fractions by [bold D] = ℤ[½].+-- +-- We internally represent a dyadic fraction /a/\/2[sup /n/] as a pair+-- (/a/,/n/). Note that this representation is not unique. When it is+-- necessary to choose a canonical representative, we choose the least+-- possible /n/≥0.+data Dyadic = Dyadic !Integer !Integer++-- | Given a dyadic fraction /r/, return (/a/,/n/) such that /r/ =+-- /a/\/2[sup /n/], where /n/≥0 is chosen as small as possible.+decompose_dyadic :: Dyadic -> (Integer, Integer)+decompose_dyadic (Dyadic a n) +  | a == 0 = (0, 0)+  | n >= k = (b, n-k)+  | otherwise = (shiftL b (fromInteger (k-n)), 0)+  where+    (b,k) = lobit a++-- | Given a dyadic fraction /r/ and an integer /k/≥0, such that /a/ =+-- /r/2[sup /k/] is an integer, return /a/. If /a/ is not an integer,+-- the behavior is undefined.+integer_of_dyadic :: Dyadic -> Integer -> Integer+integer_of_dyadic (Dyadic a n) k =+  shift a (fromInteger (k-n))++instance Real Dyadic where+  toRational (Dyadic a n) +    | n >= 0    = a % 2^n+    | otherwise = a * 2^(-n) % 1++instance Show Dyadic where+  showsPrec d a = showsPrec_rational d (toRational a)++instance Eq Dyadic where+  Dyadic a n == Dyadic b m = a * 2^(k-n) == b * 2^(k-m) where+    k = max n m++instance Ord Dyadic where+  compare (Dyadic a n) (Dyadic b m) = compare (a * 2^(k-n)) (b * 2^(k-m)) where+    k = max n m++instance Num Dyadic where+  Dyadic a n + Dyadic b m +    | n < m     = Dyadic c m+    | otherwise = Dyadic d n+    where+      c = shiftL a (fromInteger (m-n)) + b+      d = a + shiftL b (fromInteger (n-m))+  Dyadic a n * Dyadic b m = Dyadic (a*b) (n+m)+  negate (Dyadic a n) = Dyadic (-a) n+  abs x = if x >= 0 then x else -x+  signum x = case compare 0 x of { LT -> 1; EQ -> 0; GT -> -1 }+  fromInteger n = Dyadic n 0++instance HalfRing Dyadic where+  half = Dyadic 1 1++instance Adjoint Dyadic where+  adj x = x++instance Adjoint2 Dyadic where+  adj2 x = x++-- | The unique ring homomorphism from ℤ[½] to any ring containing+-- ½. This exists because ℤ[½] is the free such ring.++-- Implementation note: we can't just use fromInteger a * half^n,+-- because this can give potentially bad round-off errors in case+-- half^n underflows in the target type. Moreover, this does not work+-- if n is negative.+fromDyadic :: (HalfRing a) => Dyadic -> a+fromDyadic x = aux (fromInteger a) n where+  (a,n) = decompose_dyadic x+  aux !a !n+    | n>0       = aux (half*a) (n-1)+    | n==0      = a+    | otherwise = aux (2*a) (n+1)++-- ----------------------------------------------------------------------+-- ** The ring ℚ of rational numbers++-- | We define our own variant of the rational numbers, which is an+-- identical copy of the type 'Rational' from the standard Haskell+-- library, except that it has a more sensible 'Show' instance.+newtype Rationals = ToRationals { unRationals :: Rational }+                  deriving (Num, Eq, Ord, Fractional, Real, RealFrac, HalfRing, Adjoint, Adjoint2, ToQOmega, Floor)++-- | An auxiliary function for printing rational numbers, using+-- correct precedences, and omitting denominators of 1.+showsPrec_rational :: (Show a, Integral a) => Int -> Ratio a -> ShowS+showsPrec_rational d a+  | denom == 1 = showsPrec d numer+  | numer < 0  = showParen (d >= 7) $ showString "-" . showsPrec_rational 7 (-a)+  | otherwise  = showParen (d >= 8) $+                 showsPrec 7 numer . showString "/" . showsPrec 8 denom+    where+      numer = numerator a+      denom = denominator a++instance Show Rationals where+  showsPrec d (ToRationals a) = showsPrec_rational d a++-- | Conversion from 'Rationals' to any 'Fractional' type.+fromRationals :: (Fractional a) => Rationals -> a+fromRationals = fromRational . unRationals++-- ----------------------------------------------------------------------+-- ** The ring /R/[√2]+  +-- | The ring /R/[√2], where /R/ is any ring. The value 'RootTwo' /a/+-- /b/ represents /a/ + /b/ √2.+data RootTwo a = RootTwo !a !a+                deriving (Eq)++instance (Eq a, Num a) => Num (RootTwo a) where+  RootTwo a b + RootTwo a' b' = RootTwo a'' b'' where+    a'' = a + a'+    b'' = b + b'+  RootTwo a b * RootTwo a' b' = RootTwo a'' b'' where+    a'' = a * a' + 2 * b * b'+    b'' = a * b' + a' * b+  negate (RootTwo a b) = RootTwo a' b' where+    a' = -a+    b' = -b+  fromInteger n = RootTwo n' 0 where+    n' = fromInteger n+  abs f = f * signum f+  signum f@(RootTwo a b)+    | sa == 0 && sb == 0 = 0+    | sa /= -1 && sb /= -1 = 1+    | sa /= 1 && sb /= 1 = -1+    | sa /= -1 && sb /= 1 && signum (a*a - 2*b*b) /= -1 = 1+    | sa /= 1 && sb /= -1 && signum (a*a - 2*b*b) /= 1  = 1+    | otherwise = -1+    where+      sa = signum a+      sb = signum b++instance (Eq a, Ring a) => Ord (RootTwo a) where+  x <= y  =  signum (y-x) /= (-1)+  +instance (Show a, Eq a, Ring a) => Show (RootTwo a) where+  showsPrec d (RootTwo a 0) = showsPrec d a+  showsPrec d (RootTwo 0 1) = showString "roottwo"+  showsPrec d (RootTwo 0 (-1)) = showParen (d >= 7) $ showString "-roottwo"+  showsPrec d (RootTwo 0 b) = showParen (d >= 8) $ +    showsPrec 7 b . showString "*roottwo"+  showsPrec d (RootTwo a b) | signum b == 1 = showParen (d >= 7) $+    showsPrec 6 a . showString " + " . showsPrec 6 (RootTwo 0 b)+  showsPrec d (RootTwo a b) | otherwise = showParen (d >= 7) $+    showsPrec 6 a . showString " - " . showsPrec 7 (RootTwo 0 (-b))++instance (Eq a, Fractional a) => Fractional (RootTwo a) where+  recip (RootTwo a b) = RootTwo (a/k) (-b/k) where+    k = a^2 - 2*b^2+  fromRational r = RootTwo (fromRational r) 0++instance (Eq a, Ring a) => RootTwoRing (RootTwo a) where+  roottwo = RootTwo 0 1++instance (Eq a, HalfRing a) => HalfRing (RootTwo a) where+  half = RootTwo half 0+  +instance (Eq a, HalfRing a) => RootHalfRing (RootTwo a) where+  roothalf = RootTwo 0 half+  +instance (Eq a, ComplexRing a) => ComplexRing (RootTwo a) where+  i = RootTwo i 0++instance (Adjoint a) => Adjoint (RootTwo a) where  +  adj (RootTwo a b) = RootTwo (adj a) (adj b)++instance (Adjoint2 a, Num a) => Adjoint2 (RootTwo a) where  +  adj2 (RootTwo a b) = RootTwo (adj2 a) (-adj2 b)++instance (Eq a, NormedRing a) => NormedRing (RootTwo a) where+  norm (RootTwo a b) = (norm a)^2 - 2 * (norm b)^2++-- ----------------------------------------------------------------------+-- ** The ring ℤ[√2]++-- | The ring ℤ[√2].+type ZRootTwo = RootTwo Integer++-- | The unique ring homomorphism from ℤ[√2] to any ring containing+-- √2. This exists because ℤ[√2] is the free such ring.+fromZRootTwo :: (RootTwoRing a) => ZRootTwo -> a+fromZRootTwo (RootTwo x y) = fromInteger x + roottwo * fromInteger y++-- ----------------------------------------------------------------------+-- ** The ring [bold D][√2]++-- | The ring [bold D][√2] = ℤ[1\/√2]. +type DRootTwo = RootTwo Dyadic++-- | The unique ring homomorphism from [bold D][√2] to any ring containing+-- 1\/√2. This exists because [bold D][√2] = ℤ[1\/√2] is the free such ring.+fromDRootTwo :: (RootHalfRing a) => DRootTwo -> a+fromDRootTwo (RootTwo x y) = fromDyadic x + roottwo * fromDyadic y++-- ----------------------------------------------------------------------+-- ** The field ℚ[√2]++-- | The field ℚ[√2].+type QRootTwo = RootTwo Rationals++instance Floor QRootTwo where+  floor_of x@(RootTwo a b)+    | r'+1 <= x  = r+1+    | r' <= x    = r+    | otherwise = r-1 +   where +    a' = floor a+    b' = intsqrt (floor (2 * b^2))+    r | b >= 0    = a' + b'+      | otherwise = a' - b'+    r' = fromInteger r++-- | The unique ring homomorphism from ℚ[√2] to any ring containing+-- the rational numbers and √2. This exists because ℚ[√2] is the free+-- such ring.+fromQRootTwo :: (RootTwoRing a, Fractional a) => QRootTwo -> a+fromQRootTwo (RootTwo x y) = fromRationals x + roottwo * fromRationals y++-- ----------------------------------------------------------------------+-- ** The ring /R/[/i/]++-- | The ring /R/[/i/], where /R/ is any ring. The reason we do not+-- use the 'Complex' /a/ type from the standard Haskell libraries is+-- that it assumes too much, for example, it assumes /a/ is a member+-- of the 'RealFloat' class. Also, this allows us to define a more+-- sensible 'Show' instance.+data Cplx a = Cplx !a !a+            deriving (Eq)++instance (Eq a, Show a, Num a) => Show (Cplx a) where+  showsPrec d (Cplx a 0) = showsPrec d a+  showsPrec d (Cplx 0 1) = showString "i"+  showsPrec d (Cplx 0 (-1)) = showParen (d >= 7) $ showString "-i"+  showsPrec d (Cplx 0 b) = showParen (d >= 8) $ +    showsPrec 7 b . showString "*i"+  showsPrec d (Cplx a b) | signum b == 1 = showParen (d >= 7) $+    showsPrec 6 a . showString " + " . showsPrec 6 (Cplx 0 b)+  showsPrec d (Cplx a b) | otherwise = showParen (d >= 7) $+    showsPrec 6 a . showString " - " . showsPrec 7 (Cplx 0 (-b))++instance (Num a) => Num (Cplx a) where+  Cplx a b + Cplx a' b' = Cplx a'' b'' where+    a'' = a + a'+    b'' = b + b'+  Cplx a b * Cplx a' b' = Cplx a'' b'' where+    a'' = a * a' - b * b'+    b'' = a * b' + a' * b+  negate (Cplx a b) = Cplx a' b' where+    a' = -a+    b' = -b+  fromInteger n = Cplx n' 0 where+    n' = fromInteger n+  abs x = x+  signum x = 1++instance (Fractional a) => Fractional (Cplx a) where+  recip (Cplx a b) = Cplx (a/d) (-b/d) where+    d = a^2 + b^2+  fromRational a = Cplx (fromRational a) 0++instance (Ring a) => ComplexRing (Cplx a) where+  i = Cplx 0 1++instance (HalfRing a) => HalfRing (Cplx a) where+  half = Cplx half 0++instance (RootHalfRing a) => RootHalfRing (Cplx a) where+  roothalf = Cplx roothalf 0++instance (RootTwoRing a) => RootTwoRing (Cplx a) where+  roottwo = Cplx roottwo 0++instance (Adjoint a, Ring a) => Adjoint (Cplx a) where+  adj (Cplx a b) = (Cplx (adj a) (-(adj b)))++instance (Adjoint2 a, Ring a) => Adjoint2 (Cplx a) where+  adj2 (Cplx a b) = (Cplx (adj2 a) (adj2 b))++instance (NormedRing a) => NormedRing (Cplx a) where+  norm (Cplx a b) = (norm a)^2 + (norm b)^2++-- ----------------------------------------------------------------------+-- ** The ring ℤ[/i/] of Gaussian integers++-- | The ring ℤ[/i/] of Gaussian integers.+type ZComplex = Cplx Integer++-- | The unique ring homomorphism from ℤ[/i/] to any ring containing+-- /i/. This exists because ℤ[/i/] is the free such ring.+fromZComplex :: (ComplexRing a) => ZComplex -> a+fromZComplex (Cplx a b) = fromInteger a + i * fromInteger b++-- ----------------------------------------------------------------------+-- ** The ring [bold D][/i/]++-- | The ring [bold D][/i/] = ℤ[½, /i/] of Gaussian dyadic fractions.+type DComplex = Cplx Dyadic++-- | The unique ring homomorphism from [bold D][/i/] to any ring containing+-- ½ and /i/. This exists because [bold D][/i/] is the free such ring.+fromDComplex :: (ComplexRing a, HalfRing a) => DComplex -> a+fromDComplex (Cplx a b) = fromDyadic a + i * fromDyadic b++-- ----------------------------------------------------------------------+-- ** The ring ℚ[/i/] of Gaussian rationals++-- | The ring ℚ[/i/] of Gaussian rationals.+type QComplex = Cplx Rationals++-- | The unique ring homomorphism from ℚ[/i/] to any ring containing+-- the rational numbers and /i/. This exists because ℚ[/i/] is the+-- free such ring.+fromQComplex :: (ComplexRing a, Fractional a) => QComplex -> a+fromQComplex (Cplx a b) = fromRationals a + i * fromRationals b++-- ----------------------------------------------------------------------+-- ** The ring [bold D][√2, /i/]++-- | The ring [bold D][√2, /i/] = ℤ[1\/√2, /i/].+type DRComplex = Cplx DRootTwo++-- | The unique ring homomorphism from [bold D][√2, /i/] to any ring+-- containing 1\/√2 and /i/. This exists because [bold D][√2, /i/] =+-- ℤ[1\/√2, /i/] is the free such ring.+fromDRComplex :: (RootHalfRing a, ComplexRing a) => DRComplex -> a+fromDRComplex (Cplx a b) = fromDRootTwo a + i * fromDRootTwo b++-- ----------------------------------------------------------------------+-- ** The ring ℚ[√2, /i/]++-- | The field ℚ[√2, /i/].+type QRComplex = Cplx QRootTwo++-- | The unique ring homomorphism from ℚ[√2, /i/] to any ring+-- containing the rational numbers, √2, and /i/. This exists because+-- ℚ[√2, /i/] is the free such ring.+fromQRComplex :: (RootTwoRing a, ComplexRing a, Fractional a) => QRComplex -> a+fromQRComplex (Cplx a b) = fromQRootTwo a + i * fromQRootTwo b++-- ----------------------------------------------------------------------+-- ** The ring ℂ of complex numbers++-- $ We provide two versions of the complex numbers using floating+-- point arithmetic.++-- | Double precision complex floating point numbers.+type CDouble = Cplx Double++-- | Single precision complex floating point numbers.+type CFloat = Cplx Float++-- ----------------------------------------------------------------------+-- ** The ring /R/[ω]++-- | The ring /R/[ω], where /R/ is any ring, and ω = [exp iπ/4] is an+-- 8th root of unity. The value 'Omega' /a/ /b/ /c/ /d/ represents+-- /a/ω[sup 3]+/b/ω[sup 2]+/c/ω+/d/.+data Omega a = Omega !a !a !a !a+            deriving (Eq)++-- | An inverse to the embedding /R/ ↦ /R/[ω]: return the \"real+-- rational\" part. +-- In other words, map /a/ω[sup 3]+/b/ω[sup 2]+/c/ω+/d/ to /d/.+omega_real :: Omega a -> a+omega_real (Omega a b c d) = d++instance (Show a, Ring a) => Show (Omega a) where+  showsPrec p (Omega a b c d) = +    showParen (p >= 11) $ showString "Omega " . +                         showsPrec 11 a . showString " " . +                         showsPrec 11 b . showString " " . +                         showsPrec 11 c . showString " " . +                         showsPrec 11 d++instance (Num a) => Num (Omega a) where+  Omega a b c d + Omega a' b' c' d' = Omega a'' b'' c'' d'' where+    a'' = a + a'+    b'' = b + b'+    c'' = c + c'+    d'' = d + d'+  Omega a b c d * Omega a' b' c' d' = Omega a'' b'' c'' d'' where  +    a'' = a*d' + b*c' + c*b' + d*a'+    b'' = b*d' + c*c' + d*b' - a*a'+    c'' = c*d' + d*c' - a*b' - b*a'+    d'' = d*d' - a*c' - b*b' - c*a'+  negate (Omega a b c d) = Omega (-a) (-b) (-c) (-d) where+  fromInteger n = Omega 0 0 0 n' where+    n' = fromInteger n+  abs x = x+  signum x = 1++instance (Fractional a) => Fractional (Omega a) where+  recip (Omega a b c d) = x1 * x2 * x3 * Omega 0 0 0 (1/denom)+    where+      x1 = Omega (-c) (-b) (-a) d+      x2 = Omega (-a) b (-c) d+      x3 = Omega c (-b) a d+      denom = (a^2+b^2+c^2+d^2)^2-2*(a*b+b*c+c*d-d*a)^2+  fromRational r = fromInteger a / fromInteger b where+    a = numerator r+    b = denominator r++instance (HalfRing a) => HalfRing (Omega a) where+  half = Omega 0 0 0 half++instance (HalfRing a) => RootHalfRing (Omega a) where+  roothalf = Omega (-half) 0 half 0++instance (Ring a) => RootTwoRing (Omega a) where+  roottwo = Omega (-1) 0 1 0++instance (Ring a) => ComplexRing (Omega a) where+  i = Omega 0 1 0 0++instance (Adjoint a, Ring a) => Adjoint (Omega a) where+  adj (Omega a b c d) = Omega (-(adj c)) (-(adj b)) (-(adj a)) (adj d)++instance (Adjoint2 a, Ring a) => Adjoint2 (Omega a) where+  adj2 (Omega a b c d) = Omega (-adj2 a) (adj2 b) (-adj2 c) (adj2 d)++instance (NormedRing a) => NormedRing (Omega a) where+  norm (Omega x y z w) = (a^2+b^2+c^2+d^2)^2-2*(a*b+b*c+c*d-d*a)^2+    where+      a = norm x+      b = norm y+      c = norm z+      d = norm w++-- This is an undecidable instance, but is not overlapping. Note: we+-- do not declare OmegaRing (Omega a), because this usually follows+-- from (ComplexRing a, RootHalfRing a). But there are some+-- exceptions, such as OmegaRing (Omega Z2) and OmegaRing (Omega+-- Integer).+instance OmegaRing (Omega Z2) where+  omega = Omega 0 0 1 0++-- ----------------------------------------------------------------------+-- ** The ring ℤ[ω]++-- | The ring ℤ[ω] of /cyclotomic integers/ of degree 8. Such rings+-- were first studied by Kummer around 1840, and used in his proof of+-- special cases of Fermat's Last Theorem.  See also:+-- +-- * <http://fermatslasttheorem.blogspot.com/2006/05/basic-properties-of-cyclotomic.html>+-- +-- * <http://fermatslasttheorem.blogspot.com/2006/02/cyclotomic-integers.html>+-- +-- * Harold M. Edwards, \"Fermat's Last Theorem: A Genetic+-- Introduction to Algebraic Number Theory\".+type ZOmega = Omega Integer++-- | The unique ring homomorphism from ℤ[ω] to any ring containing+-- ω. This exists because ℤ[ω] is the free such ring.+fromZOmega :: (OmegaRing a) => ZOmega -> a+fromZOmega (Omega a b c d) = fromInteger a * omega^3 + fromInteger b * omega^2 + fromInteger c * omega + fromInteger d++-- This is an undecidable instance, but is not overlapping.+instance OmegaRing ZOmega where+  omega = Omega 0 0 1 0++-- | Inverse of the embedding ℤ[√2] → ℤ[ω]. Note that ℤ[√2] = ℤ[ω] ∩+-- ℝ. This function takes an element of ℤ[ω] that is real, and+-- converts it to an element of ℤ[√2]. It throws an error if the input+-- is not real.+zroottwo_of_zomega :: ZOmega -> ZRootTwo+zroottwo_of_zomega (Omega a b c d)+  | a == -c && b == 0  = RootTwo d c+  | otherwise = error "zroottwo_of_zomega: non-real value"+  +-- ----------------------------------------------------------------------+-- ** The ring [bold D][ω]++-- | The ring [bold D][ω]. Here [bold D]=ℤ[½] is the ring of dyadic+-- fractions. In fact, [bold D][ω] is isomorphic to the ring [bold D][√2,+-- i], but they have different 'Show' instances.+type DOmega = Omega Dyadic++-- | The unique ring homomorphism from [bold D][ω] to any ring containing+-- 1\/√2 and /i/. This exists because [bold D][ω] is the free such ring.+fromDOmega :: (RootHalfRing a, ComplexRing a) => DOmega -> a+fromDOmega (Omega a b c d) = fromDyadic a * omega^3 + fromDyadic b * omega^2 + fromDyadic c * omega + fromDyadic d++-- This is an overlapping instance. It is nicer to show an element of+-- D[ω] by pulling out a common denominator exponent. But in cases+-- where the typechecker cannot infer this, then it will just fall+-- back to the more general method.+instance Show DOmega where+  showsPrec = showsPrec_DenomExp+  +-- This is an overlapping instance. See previous comment.+instance Show DRComplex where+  showsPrec = showsPrec_DenomExp++-- ----------------------------------------------------------------------+-- ** The field ℚ[ω]++-- | The field ℚ[ω] of /cyclotomic rationals/ of degree 8.+type QOmega = Omega Rationals++-- | The unique ring homomorphism from ℚ[ω] to any ring containing the+-- rational numbers, √2, and /i/. This exists because ℚ[ω] is the free+-- such ring.+fromQOmega :: (RootHalfRing a, ComplexRing a, Fractional a) => QOmega -> a+fromQOmega (Omega a b c d) = fromRationals a * omega^3 + fromRationals b * omega^2 + fromRationals c * omega + fromRationals d++-- ----------------------------------------------------------------------+-- * Conversion to dyadic++-- | A type class relating \"rational\" types to their dyadic+-- counterparts.+class ToDyadic a b | a -> b where+  -- | Convert a \"rational\" value to a \"dyadic\" value, if the+  -- denominator is a power of 2. Otherwise, return 'Nothing'.+  maybe_dyadic :: a -> Maybe b++-- | Convert a \"rational\" value to a \"dyadic\" value, if the+-- denominator is a power of 2. Otherwise, throw an error.+to_dyadic :: (ToDyadic a b) => a -> b+to_dyadic a = case maybe_dyadic a of+  Just b -> b+  Nothing -> error "to_dyadic: denominator not a power of 2"++instance ToDyadic Dyadic Dyadic where+  maybe_dyadic = return++instance ToDyadic Rational Dyadic where+  maybe_dyadic x = do+    k <- log2 denom+    return (Dyadic numer k)+    where denom = denominator x+          numer = numerator x++instance ToDyadic Rationals Dyadic where+  maybe_dyadic = maybe_dyadic . unRationals++instance (ToDyadic a b) => ToDyadic (RootTwo a) (RootTwo b) where+  maybe_dyadic (RootTwo x y) = do+    x' <- maybe_dyadic x+    y' <- maybe_dyadic y+    return (RootTwo x' y')++instance (ToDyadic a b) => ToDyadic (Cplx a) (Cplx b) where+  maybe_dyadic (Cplx x y) = do+    x' <- maybe_dyadic x+    y' <- maybe_dyadic y+    return (Cplx x' y')++instance (ToDyadic a b) => ToDyadic (Omega a) (Omega b) where+  maybe_dyadic (Omega x y z w) = do+    x' <- maybe_dyadic x+    y' <- maybe_dyadic y+    z' <- maybe_dyadic z+    w' <- maybe_dyadic w+    return (Omega x' y' z' w')++-- ----------------------------------------------------------------------+-- * Real part+    +-- | A type class for rings that have a \"real\" component. A typical+-- instance is /a/ = 'DRComplex' with /b/ = 'DRootTwo'.+class RealPart a b | a -> b where+  -- | Take the real part.+  real :: a -> b++instance RealPart (Cplx a) a where+  real (Cplx a b) = a++instance (HalfRing a) => RealPart (Omega a) (RootTwo a) where+  real (Omega a b c d) = RootTwo d (half * (c - a))++-- ----------------------------------------------------------------------+-- * Rings of integers+  +-- | A type class for rings that have a distinguished subring \"of+-- integers\". A typical instance is /a/ = 'DRootTwo', which has /b/ =+-- 'ZRootTwo' as its ring of integers.+class WholePart a b | a -> b where  +  -- | The embedding of the ring of integers into the larger ring.+  from_whole :: b -> a+  -- | The inverse of 'from_whole'. Throws an error if the given+  -- element is not actually an integer in the ring.+  to_whole :: a -> b+  +instance WholePart Dyadic Integer where+  from_whole = fromInteger+  to_whole d +    | n == 0 = a+    | otherwise = error "to_whole: non-integral value"+    where+      (a,n) = decompose_dyadic d++instance WholePart DRootTwo ZRootTwo where+  from_whole = fromZRootTwo+  to_whole (RootTwo x y) = RootTwo (to_whole x) (to_whole y)+  +instance WholePart DOmega ZOmega where+  from_whole = fromZOmega+  to_whole (Omega x y z w) = Omega (to_whole x) (to_whole y) (to_whole z) (to_whole w)+  +instance (WholePart a a', WholePart b b') => WholePart (a,b) (a',b') where+  from_whole (x,y) = (from_whole x, from_whole y)+  to_whole (x,y) = (to_whole x, to_whole y)+  +instance WholePart () () where  +  from_whole = const ()+  to_whole = const ()+  +instance (WholePart a b) => WholePart [a] [b] where  +  from_whole = map from_whole+  to_whole = map to_whole+  +instance (WholePart a b) => WholePart (Cplx a) (Cplx b) where  +  from_whole (Cplx a b) = Cplx (from_whole a) (from_whole b)+  to_whole (Cplx a b) = Cplx (to_whole a) (to_whole b)+  +-- ----------------------------------------------------------------------+-- * Common denominators+  +-- | A type class for things from which a common power of 1\/√2 (a+-- least denominator exponent) can be factored out. Typical instances+-- are 'DRootTwo', 'DRComplex', as well as tuples, lists, vectors, and+-- matrices thereof.+class DenomExp a where+  -- | Calculate the least denominator exponent /k/ of /a/. Returns+  -- the smallest /k/≥0 such that /a/ = /b/\/√2[sup /k/] for some+  -- integral /b/.+  denomexp :: a -> Integer+  +  -- | Factor out a /k/th power of 1\/√2 from /a/. In other words,+  -- calculate /a/√2[sup /k/].+  denomexp_factor :: a -> Integer -> a++-- | Calculate and factor out the least denominator exponent /k/ of+-- /a/. Return (/b/,/k/), where /a/ = /b/\/(√2)[sup /k/] and /k/≥0.+denomexp_decompose :: (WholePart a b, DenomExp a) => a -> (b, Integer)+denomexp_decompose a = (b, k) where+  k = denomexp a+  b = to_whole (denomexp_factor a k)++-- | Generic 'show'-like method that factors out a common denominator+-- exponent.+showsPrec_DenomExp :: (WholePart a b, Show b, DenomExp a) => Int -> a -> ShowS+showsPrec_DenomExp d a +  | k == 0 = showsPrec d b+  | k == 1 = showParen (d >= 8) $ +             showString "roothalf * " . showsPrec 7 b+  | otherwise = showParen (d >= 8) $+                showString ("roothalf^" ++ show k ++ " * ") . showsPrec 7 b+  where (b, k) = denomexp_decompose a++instance DenomExp DRootTwo where+  denomexp (RootTwo x y) = k'+    where+      (a,k) = decompose_dyadic x+      (b,l) = decompose_dyadic y+      k' = maximum [2*k, 2*l-1]+  denomexp_factor a k = a * roottwo^k++instance DenomExp DOmega where+  denomexp (Omega x y z w) = k'+      where+        (a,ak) = decompose_dyadic x+        (b,bk) = decompose_dyadic y+        (c,ck) = decompose_dyadic z+        (d,dk) = decompose_dyadic w+        k = maximum [ak, bk, ck, dk]+        a' = if k == ak then a else 0+        b' = if k == bk then b else 0+        c' = if k == ck then c else 0+        d' = if k == dk then d else 0+        k' | k>0 && even (a'-c') && even (b'-d') = 2*k-1+           | otherwise = 2*k+  denomexp_factor a k = a * roottwo^k+        +instance (DenomExp a, DenomExp b) => DenomExp (a,b) where+  denomexp (a,b) = denomexp a `max` denomexp b+  denomexp_factor (a,b) k = (denomexp_factor a k, denomexp_factor b k)++instance DenomExp () where+  denomexp _ = 0+  denomexp_factor _ k = ()++instance (DenomExp a) => DenomExp [a] where+  denomexp as = maximum (0 : [ denomexp a | a <- as ])+  denomexp_factor as k = [ denomexp_factor a k | a <- as ]++instance (DenomExp a) => DenomExp (Cplx a) where+  denomexp (Cplx a b) = denomexp a `max` denomexp b+  denomexp_factor (Cplx a b) k = Cplx (denomexp_factor a k) (denomexp_factor b k)++-- ----------------------------------------------------------------------+-- * Conversion to ℚ[ω]++-- $ 'QOmega' is the largest one of our \"exact\" arithmetic types. We+-- define a 'toQOmega' family of functions for converting just about+-- anything to 'QOmega'.++-- | A type class for things that can be exactly converted to ℚ[ω].+class ToQOmega a where+  -- | Conversion to 'QOmega'.+  toQOmega :: a -> QOmega++instance ToQOmega Integer where+  toQOmega = fromInteger++instance ToQOmega Rational where+  toQOmega = fromRational++instance (ToQOmega a) => ToQOmega (RootTwo a) where+  toQOmega (RootTwo a b) = toQOmega a + roottwo * toQOmega b+  +instance ToQOmega Dyadic where+  toQOmega (Dyadic a n)+    | n >= 0    = toQOmega a * half^n+    | otherwise = toQOmega a * 2^(-n)++instance (ToQOmega a) => ToQOmega (Cplx a) where+  toQOmega (Cplx a b) = toQOmega a + i * toQOmega b++instance (ToQOmega a) => ToQOmega (Omega a) where+  toQOmega (Omega a b c d) = omega^3 * a' + omega^2 * b' + omega * c' + d'+    where+      a' = toQOmega a+      b' = toQOmega b+      c' = toQOmega c+      d' = toQOmega d++-- ----------------------------------------------------------------------+-- * Parity+    +-- | A type class for things that have parity.+class Parity a where+  -- | Return the parity of something.+  parity :: a -> Z2++instance Integral a => Parity a where+  parity x = if even x then 0 else 1+  +instance Parity ZRootTwo where+  parity (RootTwo a b) = parity a++-- ----------------------------------------------------------------------+-- * Auxiliary functions++-- | If /n/≠0, return (/a/,/k/) such that /a/ is odd and /n/ =+-- /a/⋅2[sup /k/]. If /n/=0, return (/0/,/0/).+lobit :: Integer -> (Integer, Integer)+lobit 0 = (0,0)+lobit n = aux n 0 where+  aux n !k+    | n .&. 0xffffffff == 0  = aux (shiftR n 32) (k+32)+    | n .&. 0xff == 0        = aux (shiftR n 8) (k+8)+    | even n                 = aux (shiftR n 1) (k+1)+    | otherwise              = (n,k)+        +-- | If /n/ is of the form 2[sup /k/], return /k/. Otherwise, return+-- 'Nothing'.+log2 :: Integer -> Maybe Integer+log2 n+  | a == 1 = Just k+  | otherwise = Nothing+    where+      (a,k) = lobit n++-- | For /n/ ≥ 0, return the floor of the square root of /n/. This is+-- done using integer arithmetic, so there are no rounding errors.+intsqrt :: (Integral n) => n -> n+intsqrt n +  | n <= 0 = 0+  | otherwise = iterate 1 +    where+      iterate m+        | m_sq <= n && m_sq + 2*m + 1 > n = m+        | otherwise = iterate ((m + n `div` m) `div` 2)+          where+            m_sq = m*m+
+ Quantum/Synthesis/Ring/FixedPrec.hs view
@@ -0,0 +1,26 @@+-- | This module provides ring instances for "Data.Number.FixedPrec".++module Quantum.Synthesis.Ring.FixedPrec where++import Quantum.Synthesis.Ring++import Data.Number.FixedPrec++instance Precision e => RootHalfRing (FixedPrec e) where+  roothalf = sqrt 0.5++instance Precision e => RootTwoRing (FixedPrec e) where+  roottwo = sqrt 2++instance Precision e => HalfRing (FixedPrec e) where+  half = 0.5++instance Precision e => Adjoint (FixedPrec e) where+  adj x = x+  +instance Precision e => Adjoint2 (FixedPrec e) where+  adj2 x = x++instance Precision e => Floor (FixedPrec e) where+  floor_of = floor+  ceiling_of = ceiling
+ Quantum/Synthesis/Ring/SymReal.hs view
@@ -0,0 +1,21 @@+-- | This module provides ring instances for "Quantum.Synthesis.SymReal".++module Quantum.Synthesis.Ring.SymReal where++import Quantum.Synthesis.Ring+import Quantum.Synthesis.SymReal++instance RootHalfRing SymReal where+  roothalf = sqrt 0.5++instance RootTwoRing SymReal where+  roottwo = sqrt 2++instance HalfRing SymReal where+  half = 0.5++instance Adjoint SymReal where+  adj x = x+  +instance Adjoint2 SymReal where+  adj2 x = x
+ Quantum/Synthesis/RotationDecomposition.hs view
@@ -0,0 +1,160 @@+-- | This module provides functions for decomposing a unitary /n/×/n/+-- operator into one- and two-level unitaries. +-- +-- The algorithm is adapted from Section 4.5.1 of Nielsen and+-- Chuang. In addition to what is described in Nielsen and Chuang, our+-- algorithm produces two-level operators that can be decomposed using+-- only two Euler angles. The algorithm produces at most /n/(/n/−1)\/2+-- two-level operators of type /R/[sub /z/](δ)/R/[sub /x/](γ), as well+-- as /n/ one-level operators of type [exp /i/θ]. Therefore, the+-- decomposition of a unitary /n/×/n/ operator yields /n/[sup 2] real+-- parameters, which is optimal.++module Quantum.Synthesis.RotationDecomposition where++import Quantum.Synthesis.Matrix+import Quantum.Synthesis.MultiQubitSynthesis+import Quantum.Synthesis.Ring+import Quantum.Synthesis.EulerAngles+import Quantum.Synthesis.ArcTan2++import Data.List+import System.Random++-- ----------------------------------------------------------------------+-- * Elementary rotations++-- | An elementary rotation is either a combined /x/- and+-- /z/-rotation, applied at indices /j/ and /k/, or a phase change+-- applied at index /j/.+-- +-- * 'ERot_zx' δ γ /j/ /k/ represents the operator +-- /R/[sub /z/](δ)/R/[sub /x/](γ), applied to levels /j/ and /k/.+-- +-- \[image ERot_zx.png]+-- +-- * 'ERot_phase' θ /j/ represents the operator [exp /i/θ] applied to level+-- /j/.+-- +-- \[image ERot_phase.png]+-- +-- Note: when we use a list of 'ElementaryRot's to express a sequence of+-- operators, the operators are meant to be applied right-to-left,+-- i.e., as in the mathematical notation for matrix multiplication.+-- This is the opposite of the quantum circuit notation.+data ElementaryRot a = +  ERot_zx a a Index Index+  | ERot_phase a Index+    deriving (Show)++-- | Convert a symbolic elementary rotation to a concrete matrix.+matrix_of_elementary :: (Ring a, Floating a, Nat n) => ElementaryRot a -> Matrix n n (Cplx a)+matrix_of_elementary (ERot_zx delta gamma j k) = +  twolevel_matrix (a, b) (c, d) j k where+  a = ed' * cg+  b = -i * ed' * sg+  c = -i * ed * sg+  d = ed * cg+  cg = Cplx (cos (gamma/2)) 0+  sg = Cplx (sin (gamma/2)) 0+  ed = Cplx cd sd+  ed' = Cplx cd (-sd)+  cd = cos (delta/2) +  sd = sin (delta/2)+matrix_of_elementary (ERot_phase theta j) = +  onelevel_matrix (Cplx c s) j where+    c = cos theta+    s = sin theta++-- | Convert a sequence of elementary rotations to an /n/×/n/-matrix.+matrix_of_elementaries :: (Ring a, Floating a, Nat n) => [ElementaryRot a] -> Matrix n n (Cplx a)+matrix_of_elementaries ops =+  foldl' (*) 1 [ matrix_of_elementary op | op <- ops ]++-- ----------------------------------------------------------------------+-- * Decomposition into elementary rotations++-- | Convert an /n/×/n/-matrix to a sequence of elementary rotations.+-- +-- Note: the list of elementary rotations will be returned in+-- right-to-left order, i.e., as in the mathematical notation for+-- matrix multiplication.  This is the opposite of the quantum circuit+-- notation.+rotation_decomposition :: (Eq a, Fractional a, Floating a, Adjoint a, ArcTan2 a, Nat n) => Matrix n n (Cplx a) -> [ElementaryRot a]+rotation_decomposition op = concat gates ++ reverse gates' where+  (op', gates) = mapAccumL rowop op [ (i,j) | j <- [0..n-2], i <- [j+1..n-1] ]+  gates' = [ get_phase op' i | i <- [0..n-1] ]+  (n', _) = matrix_size op+  n = fromInteger n'++-- ----------------------------------------------------------------------+-- * Auxiliary functions++-- | Construct a two-level /n/×/n/-matrix from a given 2×2-matrix and+-- indices /j/ and /k/.+twolevel_matrix_of_matrix :: (Ring a, Nat n) => Matrix Two Two a -> Index -> Index -> Matrix n n a+twolevel_matrix_of_matrix u j k = op where+  op = twolevel_matrix (a,b) (c,d) j k+  ((a,b), (c,d)) = from_matrix2x2 u+  +-- | Extract the phase of the /j/th diagonal entry of the given+-- matrix.+get_phase :: (ArcTan2 a) => Matrix n n (Cplx a) -> Index -> ElementaryRot a+get_phase op j = ERot_phase theta j where+  a = matrix_index op j j+  theta = arctan2 y x+  Cplx x y = a+             +-- | Perform a two-level operation on rows /j/ and /k/ of a matrix /U/,+-- such that the resulting matrix has a 0 in the (/j/,/k/)-position.+-- Return the inverse of the two-level operation used, as well as the+-- updated matrix.+rowop :: (Eq a, Fractional a, Floating a, Adjoint a, ArcTan2 a, Nat n) => Matrix n n (Cplx a) -> (Index, Index) -> (Matrix n n (Cplx a), [ElementaryRot a])+rowop op (j,k) +  | b == 0 = (op, [])+  | otherwise = (op', gates) +  where+    a = matrix_index op k k+    b = matrix_index op j k+    matrix = 1/Cplx (sqrt(real (a * adj a + b * adj b))) 0 `scalarmult` matrix2x2 (adj a, adj b) (b, -a)+    (alpha, beta, gamma, delta) = euler_angles matrix+    matrix2 = matrix_of_euler_angles (0, 0, gamma, delta)+    op' = twolevel_matrix_of_matrix matrix2 k j .*. op+    gates = [ ERot_zx (-delta) (-gamma) k j ]++-- ----------------------------------------------------------------------+-- * Testing++-- | Return a \"random\" unitary /n/×/n/-matrix. These matrices will+-- not quite be uniformly distributed; this function is primarily+-- meant to generate test cases. +random_unitary :: (RandomGen g, Nat n, Floating a, Random a) => g -> Matrix n n (Cplx a)+random_unitary g = op where+  op = matrix_of_elementaries gates+  gates = random_gates g (20*n^2)+  random_gates g 0 = []+  random_gates g m = h:t where+    (gamma, g1) = randomR (0, 2*pi) g+    (delta, g1') = randomR (0, 2*pi) g1+    (c, g2) = randomR (0, 1) g1'+    (j, g3) = randomR (0, n-2) g2+    (k, g4) = randomR (j+1, n-1) g3+    h = case c :: Int of+      0 -> ERot_zx delta gamma j k+      _ -> ERot_phase delta j+    t = random_gates g4 (m-1)+  (n', _) = matrix_size op+  n = fromInteger n'++-- | Generate a random matrix, decompose it, and then re-calculate the+-- matrix from the decomposition.+test :: IO ()+test = do+  g <- newStdGen+  let m = random_unitary g :: Matrix Four Four CDouble+  let gates = rotation_decomposition m+  let m' = matrix_of_elementaries gates :: Matrix Four Four CDouble+  putStrLn $ "m = " ++ show m+  putStrLn $ "gates = " ++ show gates+  putStrLn $ "m' = " ++ show m'+
+ Quantum/Synthesis/SymReal.hs view
@@ -0,0 +1,493 @@+{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE Rank2Types #-}++-- | This module provides a symbolic representation of real number+-- expressions, as well as a type class of things that can be+-- converted to arbitrary precision real numbers.+module Quantum.Synthesis.SymReal where++import Quantum.Synthesis.ArcTan2++import Control.Monad+import Data.Char (isAlpha, isAlphaNum, isDigit)+import Data.Number.FixedPrec+import Text.ParserCombinators.ReadP+import Data.Ratio++-- ----------------------------------------------------------------------+-- * Symbolic real number expressions++-- | A type to represent symbolic expressions for real numbers.+-- +-- Caution: equality '==' at this type denotes symbolic equality of+-- expressions, not equality of the defined real numbers.+data SymReal =+  Const Integer             -- ^ An integer constant.+  | Decimal Rational String -- ^ A decimal constant. This has a rational value and a string representation.+  | Plus SymReal SymReal    -- ^ /x/ @+@ /y/.+  | Minus SymReal SymReal   -- ^ /x/ @−@ /y/.+  | Times SymReal SymReal   -- ^ /x/ @*@ /y/.+  | Div SymReal SymReal     -- ^ /x/ @\/@ /y/.+  | Negate SymReal          -- ^ \−/x/.+  | Abs SymReal             -- ^ |/x/|.+  | Signum SymReal          -- ^ signum(/x/).+  | Recip SymReal           -- ^ 1\//x/.+  | Pi                      -- ^ π.+  | Euler                   -- ^ /e/.+  | Exp SymReal             -- ^ \[exp /x/].+  | Sqrt SymReal            -- ^ √/x/.+  | Log SymReal             -- ^ log /x/.+  | Power SymReal SymReal   -- ^ /x/[sup /y/].+  | Sin SymReal             -- ^ sin /x/.+  | Tan SymReal             -- ^ cos /x/.+  | Cos SymReal             -- ^ cos /x/.+  | ASin SymReal            -- ^ asin /x/.+  | ATan SymReal            -- ^ atan /x/.+  | ACos SymReal            -- ^ acos /x/.+  | Sinh SymReal            -- ^ sinh /x/.+  | Tanh SymReal            -- ^ tanh /x/.+  | Cosh SymReal            -- ^ cosh /x/.+  | ASinh SymReal           -- ^ asinh /x/.+  | ATanh SymReal           -- ^ atanh /x/.+  | ACosh SymReal           -- ^ acosh /x/.+  | ArcTan2 SymReal SymReal -- ^ arctan2 /x/ /y/.+    deriving (Eq)++instance Show SymReal where+  showsPrec d (Const x)     = showsPrec d x+  showsPrec d (Decimal x s) = showString s+  showsPrec d (Plus x y)    = showParen (d > 6) $ showsPrec 6 x . showString "+" . showsPrec 6 y+  showsPrec d (Minus x y)   = showParen (d > 6) $ showsPrec 6 x . showString "-" . showsPrec 7 y+  showsPrec d (Times x y)   = showParen (d > 7) $ showsPrec 7 x . showString "*" . showsPrec 7 y+  showsPrec d (Div x y)     = showParen (d > 7) $ showsPrec 7 x . showString "/" . showsPrec 8 y+  showsPrec d (Power x y)   = showParen (d > 8) $ showsPrec 9 x . showString "**" . showsPrec 9 y+  showsPrec d (Negate x)    = showParen (d > 5) $ showString "-" . showsPrec 7 x+  showsPrec d (Abs x)       = showParen (d > 10) $ showString "abs " . showsPrec 11 x+  showsPrec d (Signum x)    = showParen (d > 10) $ showString "signum " . showsPrec 11 x+  showsPrec d (Recip x)     = showParen (d > 7) $ showString "1/" . showsPrec 8 x+  showsPrec d Pi            = showString "pi"+  showsPrec d Euler         = showString "e"+  showsPrec d (Exp x)       = showParen (d > 10) $ showString "exp " . showsPrec 11 x+  showsPrec d (Sqrt x)      = showParen (d > 10) $ showString "sqrt " . showsPrec 11 x+  showsPrec d (Log x)       = showParen (d > 10) $ showString "log " . showsPrec 11 x+  showsPrec d (Sin x)       = showParen (d > 10) $ showString "sin " . showsPrec 11 x+  showsPrec d (Tan x)       = showParen (d > 10) $ showString "tan " . showsPrec 11 x+  showsPrec d (Cos x)       = showParen (d > 10) $ showString "cos " . showsPrec 11 x+  showsPrec d (ASin x)      = showParen (d > 10) $ showString "asin " . showsPrec 11 x+  showsPrec d (ATan x)      = showParen (d > 10) $ showString "atan " . showsPrec 11 x+  showsPrec d (ACos x)      = showParen (d > 10) $ showString "acos " . showsPrec 11 x+  showsPrec d (Sinh x)      = showParen (d > 10) $ showString "sinh " . showsPrec 11 x+  showsPrec d (Tanh x)      = showParen (d > 10) $ showString "tanh " . showsPrec 11 x+  showsPrec d (Cosh x)      = showParen (d > 10) $ showString "cosh " . showsPrec 11 x+  showsPrec d (ASinh x)     = showParen (d > 10) $ showString "asinh " . showsPrec 11 x+  showsPrec d (ATanh x)     = showParen (d > 10) $ showString "atanh " . showsPrec 11 x+  showsPrec d (ACosh x)     = showParen (d > 10) $ showString "acosh " . showsPrec 11 x+  showsPrec d (ArcTan2 y x) = showParen (d > 10) $ showString "arctan2 " . showsPrec 11 y . showString " " . showsPrec 11 x++instance Num SymReal where+  (+) = Plus+  (*) = Times+  (-) = Minus+  negate = Negate+  abs = Abs+  signum = Signum+  fromInteger = Const+  +instance Fractional SymReal where+  (/) = Div+  recip = Recip+  fromRational x = Const (numerator x) `Div` Const (denominator x)+  +instance Floating SymReal where+  pi = Pi+  exp = Exp+  sqrt = Sqrt+  log = Log+  (**) = Power+  logBase x y = log y / log x+  sin = Sin+  tan = Tan+  cos = Cos+  asin = ASin+  atan = ATan+  acos = ACos+  sinh = Sinh+  tanh = Tanh+  cosh = Cosh+  asinh = ASinh+  atanh = ATanh+  acosh = ACosh++instance ArcTan2 SymReal where+  arctan2 y x = ArcTan2 y x++-- ----------------------------------------------------------------------+-- * Conversion to real number types++-- | A type class for things that can be converted to a real number at+-- arbitrary precision.+class ToReal a where+  to_real :: (Floating r, ArcTan2 r) => a -> r++instance ToReal SymReal where+  to_real (Const x) = fromInteger x+  to_real (Decimal x s) = fromRational x+  to_real (Plus x y) = to_real x + to_real y+  to_real (Minus x y) = to_real x - to_real y+  to_real (Times x y) = to_real x * to_real y+  to_real (Negate x) = -(to_real x)+  to_real (Abs x) = abs (to_real x)+  to_real (Signum x) = signum (to_real x)+  to_real (Div x y) =  to_real x / to_real y+  to_real (Recip x) = recip (to_real x)+  to_real Pi = pi+  to_real Euler = exp 1+  to_real (Exp x) = exp (to_real x)+  to_real (Sqrt x) = sqrt (to_real x)+  to_real (Log x) = log (to_real x)+  to_real (Power x y) = to_real x ** to_real y+  to_real (Sin x) = sin (to_real x)+  to_real (Tan x) = tan (to_real x)+  to_real (Cos x) = cos (to_real x)+  to_real (ASin x) = asin (to_real x)+  to_real (ATan x) = atan (to_real x)+  to_real (ACos x) = acos (to_real x)+  to_real (Sinh x) = sinh (to_real x)+  to_real (Tanh x) = tanh (to_real x)+  to_real (Cosh x) = cosh (to_real x)+  to_real (ASinh x) = asinh (to_real x)+  to_real (ATanh x) = atanh (to_real x)+  to_real (ACosh x) = acosh (to_real x)+  to_real (ArcTan2 y x) = arctan2 (to_real y) (to_real x)+  +instance ToReal Rational where+  to_real = fromRational+  +instance ToReal Integer where+  to_real = fromInteger+  +instance ToReal Int where+  to_real = fromIntegral+  +instance ToReal Double where+  to_real = fromRational . toRational++instance ToReal Float where+  to_real = fromRational . toRational++instance (Precision e) => ToReal (FixedPrec e) where+  to_real = fromRational . toRational++instance ToReal String where+  to_real x = case parse_SymReal x of+    Just n -> to_real n+    Nothing -> error "ToReal String: string does not parse"++-- ----------------------------------------------------------------------+-- ** Dynamic conversion to FixedPrec++-- | It would be useful to have a function for converting a symbolic+-- real number to a fixed-precision real number with a chosen+-- precision, such that the precision /e/ depends on a parameter /d/:+-- +-- > to_fixedprec :: (ToReal r) => Integer -> r -> FixedPrec e+-- > to_fixedprec d x = ...+--  +-- However, since /e/ is a type, /d/ is a term, and Haskell is not+-- dependently typed, this cannot be done directly.+-- +-- The function 'dynamic_fixedprec' is the closest thing we have to a+-- workaround. The call @dynamic_fixedprec@ /d/ /f/ /x/ calls+-- /f/(/x/'), where /x/' is the value /x/ converted to /d/ digits of+-- precision.  In other words, we have+-- +-- > dynamic_fixedprec d f x = f (to_fixedprec d x),+-- +-- with the restriction that the precision /e/ cannot occur freely in+-- the result type of /f/.+dynamic_fixedprec :: forall a r.(ToReal r) => Integer -> (forall e.(Precision e) => FixedPrec e -> a) -> r -> a+dynamic_fixedprec d f x = loop d (undefined :: P0)+  where +    loop :: forall e.(Precision e) => Integer -> e -> a+    loop d e+      | d >= 1000 = loop (d-1000) (undefined :: PPlus1000 e)+      | d >= 100  = loop (d-100)  (undefined :: PPlus100 e)+      | d >= 10   = loop (d-10) (undefined :: PPlus10 e)+      | d > 0     = loop (d-1) (undefined :: PPlus1 e)+      | otherwise = f (to_real x :: FixedPrec e)++-- | Like 'dynamic_fixedprec', but take two real number arguments. In+-- terms of the fictitious function @to_fixedprec@, we have:+-- +-- > dynamic_fixedprec2 d f x y = f (to_fixedprec d x) (to_fixedprec d y).+dynamic_fixedprec2 :: forall a r s.(ToReal r, ToReal s) => Integer -> (forall e.(Precision e) => FixedPrec e -> FixedPrec e -> a) -> r -> s -> a+dynamic_fixedprec2 d f x y = loop d (undefined :: P0)+  where +    loop :: forall e.(Precision e) => Integer -> e -> a+    loop d e+      | d >= 1000 = loop (d-1000) (undefined :: PPlus1000 e)+      | d >= 100  = loop (d-100)  (undefined :: PPlus100 e)+      | d >= 10   = loop (d-10) (undefined :: PPlus10 e)+      | d > 0     = loop (d-1) (undefined :: PPlus1 e)+      | otherwise = f (to_real x :: FixedPrec e) (to_real y :: FixedPrec e)++-- ----------------------------------------------------------------------+-- * A parser for real number expressions+  +-- ----------------------------------------------------------------------+-- ** Grammar specification++-- $ Each function in this section corresponds to a production rule+-- for a context-free grammar. The type of each function is 'ReadP'+-- /a/, where /a/ is the type of the semantic value produced by the+-- grammar for that expression.+-- +-- The parser uses simple precedences. +-- +-- * Unary \"+\" and \"−\" have precedence 6. +-- +-- * Binary \"+\" and \"−\" have precedence 6 and are left+-- associative.+-- +-- * Binary \"*\" and \"\/\" have precedence 7 and are left+-- associative.+-- +-- * Binary \"**\" and \"^\" have precedence 8 and are right+-- associative.+-- +-- * All unary operators other than \"+\" and \"−\" have precedence+-- 10.+-- +-- We use /exp6/ to denote an expression whose+-- top-level operator has precedence 6 or higher, /exp7/ to denote an+-- expression whose top-level operator has precedence 7 or higher, and+-- so on.+-- +-- We also allow whitespace between lexicographic entities. For+-- simplicity, whitespace is not shown in the production rules,+-- although it appears in the code.++-- | /integer/ ::= /digit/ /digit/*.+integer :: ReadP SymReal+integer = do+  s <- munch1 isDigit+  let n = read s+  return (Const (fromInteger n))++-- | /float/ ::= /digit/* \".\" /digit/*.+-- +-- There must be at least one digit, either before or after the decimal point.+float :: ReadP SymReal+float = do+  (s1, _) <- gather $ do+    munch isDigit+  char '.'+  (s2, _) <- gather $ do+    munch isDigit+  when (length s1 == 0 && length s2 == 0) $ do+    pfail+  let num = read (s1++s2) :: Integer+  let denom = 10^(length s2)+  let s1' = if s1 == [] then "0" else s1+  let s2' = if s2 == [] then "0" else s2+  return (Decimal (num % denom) (s1' ++ "." ++ s2'))++-- | /const_pi/ ::= \"pi\".+const_pi :: ReadP SymReal+const_pi = do+  string "pi"+  return Pi++-- | /const_e/ ::= \"e\".+const_e :: ReadP SymReal+const_e = do+  string "e"+  return Euler++-- | /negative/ ::= \"−\".+negative :: ReadP (SymReal -> SymReal)+negative = do+  string "-"+  skipSpaces+  return Negate++-- | /positive/ ::= \"+\".+positive :: ReadP (SymReal -> SymReal)+positive = do+  string "+"+  skipSpaces+  return id++-- | /plus_term/ ::= \"+\" /exp7/.+plus_term :: ReadP (SymReal -> SymReal)+plus_term = do+  skipSpaces+  string "+"+  skipSpaces+  n2 <- exp7+  return (\n1 -> Plus n1 n2)++-- | /minus_term/ ::= \"−\" /exp7/.+minus_term :: ReadP (SymReal -> SymReal)+minus_term = do+  skipSpaces+  string "-"+  skipSpaces+  n2 <- exp7+  return (\n1 -> Minus n1 n2)++-- | /times_term/ ::= \"*\" /exp8/.+times_term :: ReadP (SymReal -> SymReal)+times_term = do+  skipSpaces+  string "*"+  skipSpaces+  n2 <- exp8+  return (\n1 -> Times n1 n2)++-- | /div_term/ ::= \"\/\" /exp8/.+div_term :: ReadP (SymReal -> SymReal)+div_term = do+  skipSpaces+  string "/"+  skipSpaces+  n2 <- exp8+  return (\n1 -> Div n1 n2)++-- | /power_term/ ::= /exp10/ \"**\" | /exp10/ \"^\".+power_term :: ReadP (SymReal -> SymReal)+power_term = do+  n1 <- exp10+  skipSpaces+  string "**" +++ string "^"+  skipSpaces+  return (\n2 -> Power n1 n2)++-- | /unary_fun/ ::= /unary_op/ /exp10/.+unary_fun :: ReadP SymReal+unary_fun = do+  skipSpaces+  op <- unary_op+  skipSpaces+  n <- exp10+  return (op n)++-- | /unary_op/ ::= \"abs\" | \"signum\" | ...+unary_op :: ReadP (SymReal -> SymReal)+unary_op = +  choice [ do { string s; return op } | (s, op) <- ops ]+  where +    ops = [ ("abs", Abs),+            ("signum", Signum),+            ("recip", Recip),+            ("exp", Exp),+            ("sqrt", Sqrt),+            ("log", Log),+            ("sin", Sin),+            ("tan", Tan),+            ("cos", Cos),+            ("asin", ASin),+            ("atan", ATan),+            ("acos", ACos),+            ("sinh", Sinh),+            ("tanh", Tanh),+            ("cosh", Cosh),+            ("asinh", ASinh),+            ("atanh", ATanh),+            ("acosh", ACosh) ]++-- | /binary_fun/ ::= /binary_op/ /exp10/ /exp10/.+binary_fun :: ReadP SymReal+binary_fun = do+  skipSpaces+  op <- binary_op+  skipSpaces+  n <- exp10+  skipSpaces+  m <- exp10+  return (op n m)++-- | /binary_op/ ::= \"abs\" | \"signum\" | ...+binary_op :: ReadP (SymReal -> SymReal -> SymReal)+binary_op = +  choice [ do { string s; return op } | (s, op) <- ops ]+  where +    ops = [ ("arctan2", ArcTan2) ]++-- | /exp6/ ::= (/negative/ | /positive/)? /exp7/ ( /plus_term/ | /minus_term/ )*.+-- +-- An expression whose top-level operator has precedence 6 or+-- above. The operators of precedence 6 are \"+\" and \"−\".+exp6 :: ReadP SymReal+exp6 = do+  sign <- option id (negative +++ positive)+  n1 <- exp7+  ops <- many $ do+    plus_term +++ minus_term+  return (foldl (\x f -> f x) (sign n1) ops)++-- | /exp7/ ::= /exp8/ ( /times_term/ | /div_term/ )*.+-- +-- An expression whose top-level operator has precedence 7 or+-- above. The operators of precedence 6 are \"*\" and \"\/\".+exp7 :: ReadP SymReal+exp7 = do+  n1 <- exp8+  ops <- many $ do+    times_term +++ div_term+  return (foldl (\x f -> f x) n1 ops)++-- | /exp8/ ::= ( /power_term/ )* /exp10/+-- +-- An expression whose top-level operator has precedence 8 or+-- above. The operators of precedence 6 are \"**\" and \"^\".+exp8 :: ReadP SymReal+exp8 = do+  ops <- many $ do+    power_term+  n2 <- exp10+  return (foldr (\f x -> f x) n2 ops)++-- | /exp10/ ::= /parenthesized/ | /const_pi/ | /const_e/ | /integer/ | /float/ | /unary_fun/ | /binary_fun/.+-- +-- An expression whose top-level operator has precedence 10 or+-- above. Such expressions are constants, applications of unary+-- operators (except unary \"−\" and \"+\"), and parenthesized+-- expressions.+exp10 :: ReadP SymReal+exp10 = parenthesized +++ const_pi +++ const_e +++ integer +++ float +++ unary_fun +++ binary_fun++-- | /parenthesized/ ::= \"(\" /exp6/ \")\".+parenthesized :: ReadP SymReal+parenthesized = do+  string "("+  skipSpaces+  n <- exp6+  skipSpaces+  string ")"+  return n++-- | /expression/ ::= /exp6/ /end-of-line/.+--   +-- This is a top-level expression.+expression :: ReadP SymReal+expression = do+  skipSpaces+  s <- exp6+  skipSpaces+  eof+  return s++-- ----------------------------------------------------------------------+-- ** Top-level parser++-- | Parse a symbolic real number expression. Typical strings that can+-- be parsed are @\"1.0\"@, @\"pi\/128\"@, @\"(1+sin(pi\/3))^2\"@, etc. If+-- the expression cannot be parsed, return 'Nothing'.+parse_SymReal :: String -> Maybe SymReal+parse_SymReal s =+  case readP_to_S expression s of+    (h, ""):_ -> Just h+    _ -> Nothing
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Superdoc+main = superdocMain
+ images/E.png view

binary file changed (absent → 1496 bytes)

+ images/ERot_phase.png view

binary file changed (absent → 1572 bytes)

+ images/ERot_zx.png view

binary file changed (absent → 3957 bytes)

+ newsynth.cabal view
@@ -0,0 +1,102 @@+-- The name of the package.+name:                newsynth++-- The package version.  See the Haskell package versioning policy (PVP) +-- for standards guiding when and how versions should be incremented.+-- http://www.haskell.org/haskellwiki/Package_versioning_policy+-- PVP summary:      +-+------- breaking API changes+--                   | | +----- non-breaking API additions+--                   | | | +--- code changes with no API change+version:             0.1.0.0++-- A short (one-line) description of the package.+synopsis:            Exact and approximate synthesis of quantum circuits++-- A longer description of the package.+description:         ++  A library of algorithms for exact and approximate synthesis of+  quantum circuits over the Clifford+T gate set. This includes, among+  other things:+  .+  * "Quantum.Synthesis.Newsynth": an efficient single-qubit+    approximate synthesis algorithm. From P. Selinger, \"Efficient+    Clifford+T approximation of single-qubit operators\",+    <http://arxiv.org/abs/1212.6253>.+  .+  * "Quantum.Synthesis.MultiQubitSynthesis": multi-qubit exact+    synthesis algorithms. From B. Giles and P. Selinger, \"Exact+    synthesis of multiqubit Clifford+/T/ circuits\", Physical Review A+    87, 032332, 2013, <http://arxiv.org/abs/1212.0506>.+  .+  * "Quantum.Synthesis.CliffordT": the computation of+    Matsumoto-Amano normal forms. From K. Matsumoto and K. Amano,+    \"Representation of Quantum Circuits with Clifford and π\/8+    Gates\", <http://arxiv.org/abs/0806.3834>.+  .+  * "Quantum.Synthesis.RotationDecomposition": an algorithm for+    decomposing multi-qubit unitary operators into one- and two-level+    unitaries. See e.g. Section 4.5.1 of M. A. Nielsen and+    I. L. Chuang, \"Quantum Computation and Quantum Information\",+    Cambridge University Press, 2002.+  .+  This package also provides an easy-to-use command line tool for+  single-qubit approximate synthesis.++-- URL for the project homepage or repository.+homepage:            http://www.mathstat.dal.ca/~selinger/newsynth/++-- The license under which the package is released.+license:             GPL-3++-- The file containing the license text.+license-file:        LICENSE++-- The package author(s).+author:              Peter Selinger++-- An email address to which users can send suggestions, bug reports, and +-- patches.+maintainer:          selinger@mathstat.dal.ca++-- A copyright notice.+copyright:           Copyright (c) 2012-2013 Peter Selinger++-- A classification category for future use by the package catalogue+-- Hackage. These categories have not yet been specified, but the+-- upper levels of the module hierarchy make a good start.+category:            Quantum++-- The type of build used by this package.+build-type:          Custom++-- Constraint on the version of Cabal needed to build this package.+cabal-version:       >=1.8++-- A list of additional files to be included in source distributions+-- built with setup sdist.+extra-source-files:  images/*.png ChangeLog++library+  -- Modules exported by the library.+  exposed-modules:     Quantum.Synthesis.Newsynth, Quantum.Synthesis.Matrix, Quantum.Synthesis.LaTeX, Quantum.Synthesis.RotationDecomposition, Quantum.Synthesis.ArcTan2, Quantum.Synthesis.EulerAngles, Quantum.Synthesis.EuclideanDomain, Quantum.Synthesis.SymReal, Quantum.Synthesis.Ring, Quantum.Synthesis.Clifford, Quantum.Synthesis.MultiQubitSynthesis, Quantum.Synthesis.CliffordT, Quantum.Synthesis.Ring.FixedPrec, Quantum.Synthesis.Ring.SymReal+  +  -- Modules included in this library but not exported.+  -- other-modules:       +  +  -- Other library packages from which modules are imported.+  build-depends:       base ==4.6.*, random ==1.0.*, fixedprec ==0.2.*, superdoc ==0.1.*+++executable newsynth+  -- .hs or .lhs file containing the Main module.+  main-is:             newsynth.hs++  -- Root directories for the module hierarchy.+  hs-source-dirs:      programs++  -- Modules included in this executable, other than Main.+  other-modules:       CommandLine+  +  -- Other library packages from which modules are imported.+  build-depends:       base ==4.6.*, random ==1.0.*, time ==1.4.*, superdoc ==0.1.*, newsynth
+ programs/CommandLine.hs view
@@ -0,0 +1,81 @@+-- | This module provides some functions that are useful in the+-- processing of command line options, and that are shared between+-- several algorithms.++module CommandLine where++-- import other stuff+import System.Exit+import System.IO+import Data.List+import Data.Char++-- ----------------------------------------------------------------------+-- * Formatting of lists and strings++-- | A general list-to-string function. Example:+-- +-- > string_of_list "{" ", " "}" "{}" show [1,2,3] = "{1, 2, 3}"+string_of_list :: String -> String -> String -> String -> (t -> String) -> [t] -> String+string_of_list lpar comma rpar nil string_of_elt lst =+  let string_of_tail lst =+        case lst of+          [] -> ""+          h:t -> comma ++ string_of_elt h ++ string_of_tail t+  in+  case lst of+    [] -> nil+    h:t -> lpar ++ string_of_elt h ++ string_of_tail t ++ rpar++-- ----------------------------------------------------------------------+-- * Option processing+      +-- | Exit with an error message after a command line error. This also+-- outputs information on where to find command line help.+optfail :: String -> IO a+optfail msg = do+  hPutStr stderr msg+  hPutStrLn stderr "Try --help for more info."+  exitFailure++-- | Parse a string to an integer, or return 'Nothing' on failure.+parse_int :: (Integral r) => String -> Maybe r+parse_int s = case reads s of+  [(n, "")] -> Just (fromInteger n)+  _ -> Nothing++-- | Parse a string to a list of integers, or return 'Nothing' on failure.+parse_list_int :: String -> Maybe [Int]      +parse_list_int s = case reads s of+  [(ns, "")] -> Just ns+  _ -> Nothing++-- | Parse a string to a 'Double', or return 'Nothing' on failure.+parse_double :: String -> Maybe Double+parse_double s = case reads s of+  [(n, "")] -> Just n+  _ -> Nothing++-- | In an association list, find the key that best matches the given+-- string. If one key matches exactly, return the corresponding+-- key-value pair. Otherwise, return a list of all key-value pairs+-- whose key have the given string as a prefix. This list could be of+-- length 0 (no match), 1 (unique match), or greater (ambiguous key).+-- Note: the keys in the association list must be lower case. The+-- input string is converted to lower case as well, resulting in+-- case-insensitive matching.+match_enum :: [(String, a)] -> String -> [(String, a)]+match_enum list key =+  case lookup s list of+    Just v -> [(s,v)]+    Nothing -> filter (\(k,v) -> isPrefixOf s k) list+  where+    s = map toLower key+    +-- | Pretty-print a list of possible values for a parameter. The+-- first argument is the name of the parameter, and the second+-- argument is its enumeration.+show_enum :: String -> [(String, a)] -> String    +show_enum param list =+  "Possible values for " ++ param ++ " are: " +++  string_of_list "" ", " "" "no possible values" fst list ++ ".\n"
+ programs/newsynth.hs view
@@ -0,0 +1,270 @@+-- | This module provides a command line interface to the+-- decomposition library.++module Main where++import Quantum.Synthesis.Newsynth+import Quantum.Synthesis.SymReal+import Quantum.Synthesis.CliffordT+import Quantum.Synthesis.Ring+import Quantum.Synthesis.Matrix+import Quantum.Synthesis.LaTeX++import CommandLine++-- import other stuff+import Control.Monad+import Data.Time+import System.Console.GetOpt+import System.Environment    +import System.Exit+import System.Random+import Text.Printf++-- ----------------------------------------------------------------------+-- * Option processing++-- | A data type to hold values set by command line options.+data Options = Options {+  opt_digits :: Double,       -- ^ Requested precision in digits (default: 10).+  opt_theta :: SymReal,       -- ^ Angle to approximate.+  opt_hex   :: Bool,          -- ^ Output operator in hex coding? (default: ASCII).+  opt_stats :: Bool,          -- ^ Output statistics?+  opt_latex :: Bool,          -- ^ Additional LaTeX output?+  opt_count :: Integer,       -- ^ Repetition count for stats (default: 1).+  opt_rseed :: Maybe StdGen   -- ^ An optional random seed.+} deriving Show++-- | The default options.+defaultOptions :: Options+defaultOptions = Options+  { opt_digits = 10,+    opt_theta = 0.0,+    opt_hex   = False,+    opt_stats = False,+    opt_latex = False,+    opt_count = 1,+    opt_rseed = Nothing+  }++-- | The list of command line options, in the format required by 'getOpt'.+options :: [OptDescr (Options -> IO Options)]+options =+  [ Option ['h'] ["help"]    (NoArg help)           "print usage info and exit",+    Option ['d'] ["digits"]  (ReqArg digits "<n>")  "set precision in decimal digits (default: 10)",+    Option ['b'] ["bits"]    (ReqArg bits "<n>")    "set precision in bits",+    Option ['e'] ["epsilon"] (ReqArg epsilon "<n>") "set precision as epsilon (default: 1e-10)",+    Option ['x'] ["hex"]     (NoArg hex)            "output hexadecimal coding (default: ASCII)",+    Option ['s'] ["stats"]   (NoArg stats)          "output statistics",+    Option ['l'] ["latex"]   (NoArg latex)          "additional output in LaTeX format",+    Option ['c'] ["count"]   (ReqArg count "<n>")   "average statistics over <n> runs (default: 1)",+    Option ['r'] ["rseed"]   (ReqArg rseed "\"<s>\"") "set optional random seed (default: random)"+  ]+    where+      help :: Options -> IO Options+      help o = do+        usage+        exitSuccess++      digits :: String -> Options -> IO Options+      digits string o =+        case parse_double string of+          Just n | n >= 0 -> return o { opt_digits = n }+          Just n -> optfail ("Number of digits must not be negative -- " ++ string ++ "\n")+          _ -> optfail ("Invalid digits -- " ++ string ++ "\n")++      bits :: String -> Options -> IO Options+      bits string o =+        case parse_double string of+          Just n | n >= 0 -> return o { opt_digits = n * logBase 10 2 }+          Just n -> optfail ("Number of bits must not be negative -- " ++ string ++ "\n")+          _ -> optfail ("Invalid bits -- " ++ string ++ "\n")++      epsilon :: String -> Options -> IO Options+      epsilon string o =+        case parse_double string of+          Just eps | eps < 1 && eps > 0 -> return o { opt_digits = -logBase 10 eps }+          Just n -> optfail ("Epsilon must be between 0 and 1 -- " ++ string ++ "\n")+          _ -> optfail ("Invalid epsilon -- " ++ string ++ "\n")++      hex :: Options -> IO Options+      hex o = return o { opt_hex = True }++      stats :: Options -> IO Options+      stats o = return o { opt_stats = True }++      latex :: Options -> IO Options+      latex o = return o { opt_latex = True }++      count :: String -> Options -> IO Options+      count string o =+        case parse_int string of+          Just n | n >= 1 -> return o { opt_count = n }+          Just n -> optfail ("Invalid count, must be positive -- " ++ string ++ "\n")+          _ -> optfail ("Invalid count -- " ++ string ++ "\n")++      rseed :: String -> Options -> IO Options+      rseed string o =+        case reads string of+          [(g, "")] -> return o { opt_rseed = Just g }+          _ -> optfail ("Invalid random seed -- " ++ string ++ "\n")++-- | Process /argv/-style command line options into an 'Options' structure.+dopts :: [String] -> IO Options+dopts argv = do+  let (o, args, errs) = getOpt Permute options argv+  opts <- foldM (flip id) defaultOptions o+  when (errs /= []) $ do+    optfail (concat errs)+  case args of+    [] -> optfail "Missing argument: theta.\n"+    [string] -> do+      case parse_SymReal string of+        Just theta -> return opts { opt_theta = theta }+        _ -> optfail ("Invalid theta -- " ++ string ++ "\n")+    h1:h2:[] -> optfail ("Too many non-option arguments -- " ++ h1 ++ ", " ++ h2 ++ "\n")+    h1:h2:_ -> optfail ("Too many non-option arguments -- " ++ h1 ++ ", " ++ h2 ++ "...\n")++-- | Print usage message to 'stdout'.+usage :: IO ()+usage = do+  putStr (usageInfo header options) +    where header = +            "Usage: newsynth [OPTION...] <theta>\n"+            ++ "Arguments:\n"+            ++ " <theta>                   z-rotation angle. May be symbolic, e.g. pi/128\n"+            ++ "Options:"++-- ----------------------------------------------------------------------+-- * The main function++-- | Main function: read options, then execute the appropriate tasks.+main :: IO()+main = do+  -- Read options.+  argv <- getArgs+  options <- dopts argv+  let digits = opt_digits options+  let prec = digits * logBase 2 10+  let theta = opt_theta options+  let count = opt_count options+  let exponent = ceiling digits+  +  -- Set random seed.+  g <- case opt_rseed options of+    Nothing -> newStdGen+    Just g -> return g+  +  -- Expand random seed opt_count times.+  let gs = expand_seed count g++  -- Do it for each element of gs.+  stats <- sequence $ flip map (zip gs [1..]) $ \(g,n) -> do+    +    when (count > 1 && (opt_stats options || opt_latex options)) $ do+      putStrLn ("Solution " ++ show n ++ ":")+    +    -- Payload.+    t0 <- getCurrentTime+    let (m,err,ct) = newsynth_stats prec theta g+        gates = to_gates m+    if opt_hex options then+      printf "%x\n" (convert gates :: Integer)+      else+      putStrLn (if gates == [] then "I" else convert gates)+    t1 <- getCurrentTime++    -- Print optional statistics+    let tcount = length $ filter (==T) gates+    let secs = diffUTCTime t1 t0+  +    when (opt_stats options || opt_latex options) $ do+      putStrLn ("Random seed: " ++ show g)+      putStrLn ("T-count: " ++ show tcount)+    +    when (opt_stats options) $ do+      putStrLn ("Theta: " ++ show theta)+      putStrLn ("Epsilon: " ++ show_exp 10 exponent (Just digits))+      putStrLn ("Matrix: " ++ show m)+      putStrLn ("Actual error: " ++ show_exp 10 exponent err)+      putStrLn ("Runtime: " ++ show secs)+      putStrLn ("Candidates tried: " ++ show ct)+      putStrLn ("Time/candidate: " ++ show (secs / fromInteger ct))++    -- Optional LaTeX output+    when (opt_latex options) $ do+      putStrLn ("LaTeX Gates: " ++ showlatex gates)+      putStrLn ("LaTeX Theta: " ++ showlatex theta)+      putStrLn ("LaTeX Epsilon: " ++ showlatex_exp 5 exponent (Just digits))+      putStrLn ("LaTeX Matrix: " ++ showlatex (convert gates :: U2 DOmega))+      putStrLn ("LaTeX Actual error: " ++ showlatex_exp 5 exponent err)+      putStrLn ("LaTeX Runtime: " ++ show (round_to 2 secs))+      putStrLn ("LaTeX Candidates tried: " ++ show ct)+      putStrLn ("LaTeX Time/candidate: " ++ show (round_to 4 (secs / fromInteger ct)))+      +    when (count > 1 && (opt_stats options || opt_latex options)) $ do+      putStrLn ""++    return (ct, secs)++  -- If count > 1, show summary stats.+  when (count > 1) $ do+    let (cts, secss) = unzip stats+    let ct_total = sum cts+    let secs_total = sum secss+    +    when (opt_stats options || opt_latex options) $ do+      putStrLn "Summary:"+      putStrLn ("Number of runs: " ++ show count)+      putStrLn ("Total runtime: " ++ show secs_total)+    +    when (opt_stats options) $ do+      putStrLn ("Average runtime: " ++ show (secs_total / fromInteger count))+      putStrLn ("Average candidates tried: " ++ show (fromInteger ct_total / fromInteger count :: Double))+      putStrLn ("Average time/candidate: " ++ show (secs_total / fromInteger ct_total))++    when (opt_latex options) $ do+      putStrLn ("LaTeX Average runtime: " ++ show (round_to 2 (secs_total / fromInteger count)))+      putStrLn ("LaTeX Average candidates tried: " ++ show (fromInteger ct_total / fromInteger count :: Double))+      putStrLn ("LaTeX Average time/candidate: " ++ show (round_to 4 (secs_total / fromInteger ct_total)))++-- ----------------------------------------------------------------------+-- * Miscellaneous++-- | Round a 'RealFrac' value to the given number of decimals.                +round_to :: (RealFrac r) => Integer -> r -> r               +round_to n x = fromInteger (round (x * 10^n)) / 10^n++-- | Show the number 10[sup -/x/] in the format 10^(-n) or+-- 1.23*10^(-n), with precision /d/ and exponent -/n/. A value of+-- 'Nothing' is treated as 0.+-- +-- For example, display @0.316*10^(-13)@ instead of @10^(-13.5)@.+show_exp :: (Show r, RealFrac r, Floating r, PrintfArg r) => Integer -> Integer -> Maybe r -> String+show_exp d n x+  | y == 1    = "10^(" ++ show (-n) ++ ")"+  | otherwise = printf ("%." ++ show d ++ "f") y ++ "*10^(" ++ show (-n) ++ ")"+  where+    y = case x of+      Nothing -> 0+      Just x -> round_to d (10 ** (fromInteger n - x))+  +-- | Show the number 10[sup -/x/] in the format @10^{-n}@ or+-- @1.23\\cdot 10^{-n}@, with precision /d/ and exponent -/n/. A value+-- of 'Nothing' is treated as 0.+showlatex_exp :: (Show r, RealFrac r, Floating r, PrintfArg r) => Integer -> Integer -> Maybe r -> String+showlatex_exp d n x +  | y == 1    = "10^{" ++ show (-n) ++ "}"+  | otherwise = printf ("%." ++ show d ++ "f") y ++ "\\cdot 10^{" ++ show (-n) ++ "}"+  where+    y = case x of+      Nothing -> 0+      Just x -> round_to d (10 ** (fromInteger n - x))++-- | Expand a random seed /g/ into a list [/g/[sub 1], …, +-- /g/[sub /n/]] of /n/ random seeds. This is done in such a way that+-- /g/[sub 1] = /g/.+expand_seed :: (RandomGen g) => Integer -> g -> [g]+expand_seed 0 g = []+expand_seed n g = g:expand_seed (n-1) g' where+  (g', _) = split g