diff --git a/ChangeLog b/ChangeLog
new file mode 100644
--- /dev/null
+++ b/ChangeLog
@@ -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.
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -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>.
+
diff --git a/Quantum/Synthesis/ArcTan2.hs b/Quantum/Synthesis/ArcTan2.hs
new file mode 100644
--- /dev/null
+++ b/Quantum/Synthesis/ArcTan2.hs
@@ -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)
diff --git a/Quantum/Synthesis/Clifford.hs b/Quantum/Synthesis/Clifford.hs
new file mode 100644
--- /dev/null
+++ b/Quantum/Synthesis/Clifford.hs
@@ -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)
diff --git a/Quantum/Synthesis/CliffordT.hs b/Quantum/Synthesis/CliffordT.hs
new file mode 100644
--- /dev/null
+++ b/Quantum/Synthesis/CliffordT.hs
@@ -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
diff --git a/Quantum/Synthesis/EuclideanDomain.hs b/Quantum/Synthesis/EuclideanDomain.hs
new file mode 100644
--- /dev/null
+++ b/Quantum/Synthesis/EuclideanDomain.hs
@@ -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`
diff --git a/Quantum/Synthesis/EulerAngles.hs b/Quantum/Synthesis/EulerAngles.hs
new file mode 100644
--- /dev/null
+++ b/Quantum/Synthesis/EulerAngles.hs
@@ -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))
+
+  
diff --git a/Quantum/Synthesis/LaTeX.hs b/Quantum/Synthesis/LaTeX.hs
new file mode 100644
--- /dev/null
+++ b/Quantum/Synthesis/LaTeX.hs
@@ -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
diff --git a/Quantum/Synthesis/Matrix.hs b/Quantum/Synthesis/Matrix.hs
new file mode 100644
--- /dev/null
+++ b/Quantum/Synthesis/Matrix.hs
@@ -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))
diff --git a/Quantum/Synthesis/MultiQubitSynthesis.hs b/Quantum/Synthesis/MultiQubitSynthesis.hs
new file mode 100644
--- /dev/null
+++ b/Quantum/Synthesis/MultiQubitSynthesis.hs
@@ -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)
diff --git a/Quantum/Synthesis/Newsynth.hs b/Quantum/Synthesis/Newsynth.hs
new file mode 100644
--- /dev/null
+++ b/Quantum/Synthesis/Newsynth.hs
@@ -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)
diff --git a/Quantum/Synthesis/Ring.hs b/Quantum/Synthesis/Ring.hs
new file mode 100644
--- /dev/null
+++ b/Quantum/Synthesis/Ring.hs
@@ -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
+
diff --git a/Quantum/Synthesis/Ring/FixedPrec.hs b/Quantum/Synthesis/Ring/FixedPrec.hs
new file mode 100644
--- /dev/null
+++ b/Quantum/Synthesis/Ring/FixedPrec.hs
@@ -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
diff --git a/Quantum/Synthesis/Ring/SymReal.hs b/Quantum/Synthesis/Ring/SymReal.hs
new file mode 100644
--- /dev/null
+++ b/Quantum/Synthesis/Ring/SymReal.hs
@@ -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
diff --git a/Quantum/Synthesis/RotationDecomposition.hs b/Quantum/Synthesis/RotationDecomposition.hs
new file mode 100644
--- /dev/null
+++ b/Quantum/Synthesis/RotationDecomposition.hs
@@ -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'
+  
diff --git a/Quantum/Synthesis/SymReal.hs b/Quantum/Synthesis/SymReal.hs
new file mode 100644
--- /dev/null
+++ b/Quantum/Synthesis/SymReal.hs
@@ -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
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Superdoc
+main = superdocMain
diff --git a/images/E.png b/images/E.png
new file mode 100644
Binary files /dev/null and b/images/E.png differ
diff --git a/images/ERot_phase.png b/images/ERot_phase.png
new file mode 100644
Binary files /dev/null and b/images/ERot_phase.png differ
diff --git a/images/ERot_zx.png b/images/ERot_zx.png
new file mode 100644
Binary files /dev/null and b/images/ERot_zx.png differ
diff --git a/newsynth.cabal b/newsynth.cabal
new file mode 100644
--- /dev/null
+++ b/newsynth.cabal
@@ -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
diff --git a/programs/CommandLine.hs b/programs/CommandLine.hs
new file mode 100644
--- /dev/null
+++ b/programs/CommandLine.hs
@@ -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"
diff --git a/programs/newsynth.hs b/programs/newsynth.hs
new file mode 100644
--- /dev/null
+++ b/programs/newsynth.hs
@@ -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