canon (empty) → 0.1.0.0
raw patch · 12 files changed
+3213/−0 lines, 12 filesdep +arithmoidep +arraydep +basesetup-changedbinary-added
Dependencies added: arithmoi, array, base, containers, polynomial
Files
- Changes +2/−0
- LICENSE +675/−0
- Math/NumberTheory/Canon.hs +879/−0
- Math/NumberTheory/Canon/Additive.hs +109/−0
- Math/NumberTheory/Canon/AurifCyclo.hs +490/−0
- Math/NumberTheory/Canon/Internals.hs +590/−0
- Math/NumberTheory/Canon/Simple.hs +292/−0
- README.md +1/−0
- Setup.hs +2/−0
- canon.cabal +38/−0
- goBigOrGoHome.odp binary
- test-suite/CanonManualTests.hs +135/−0
+ Changes view
@@ -0,0 +1,2 @@+0.1.0.0:+ First release
+ LICENSE view
@@ -0,0 +1,675 @@+ GNU GENERAL PUBLIC LICENSE+ Version 3, 29 June 2007++ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>+ Everyone is permitted to copy and distribute verbatim copies+ of this license document, but changing it is not allowed.++ Preamble++ The GNU General Public License is a free, copyleft license for+software and other kinds of works.++ The licenses for most software and other practical works are designed+to take away your freedom to share and change the works. By contrast,+the GNU General Public License is intended to guarantee your freedom to+share and change all versions of a program--to make sure it remains free+software for all its users. We, the Free Software Foundation, use the+GNU General Public License for most of our software; it applies also to+any other work released this way by its authors. You can apply it to+your programs, too.++ When we speak of free software, we are referring to freedom, not+price. Our General Public Licenses are designed to make sure that you+have the freedom to distribute copies of free software (and charge for+them if you wish), that you receive source code or can get it if you+want it, that you can change the software or use pieces of it in new+free programs, and that you know you can do these things.++ To protect your rights, we need to prevent others from denying you+these rights or asking you to surrender the rights. Therefore, you have+certain responsibilities if you distribute copies of the software, or if+you modify it: responsibilities to respect the freedom of others.++ For example, if you distribute copies of such a program, whether+gratis or for a fee, you must pass on to the recipients the same+freedoms that you received. You must make sure that they, too, receive+or can get the source code. And you must show them these terms so they+know their rights.++ Developers that use the GNU GPL protect your rights with two steps:+(1) assert copyright on the software, and (2) offer you this License+giving you legal permission to copy, distribute and/or modify it.++ For the developers' and authors' protection, the GPL clearly explains+that there is no warranty for this free software. For both users' and+authors' sake, the GPL requires that modified versions be marked as+changed, so that their problems will not be attributed erroneously to+authors of previous versions.++ Some devices are designed to deny users access to install or run+modified versions of the software inside them, although the manufacturer+can do so. This is fundamentally incompatible with the aim of+protecting users' freedom to change the software. The systematic+pattern of such abuse occurs in the area of products for individuals to+use, which is precisely where it is most unacceptable. Therefore, we+have designed this version of the GPL to prohibit the practice for those+products. If such problems arise substantially in other domains, we+stand ready to extend this provision to those domains in future versions+of the GPL, as needed to protect the freedom of users.++ Finally, every program is threatened constantly by software patents.+States should not allow patents to restrict development and use of+software on general-purpose computers, but in those that do, we wish to+avoid the special danger that patents applied to a free program could+make it effectively proprietary. To prevent this, the GPL assures that+patents cannot be used to render the program non-free.++ The precise terms and conditions for copying, distribution and+modification follow.++ TERMS AND CONDITIONS++ 0. Definitions.++ "This License" refers to version 3 of the GNU General Public License.++ "Copyright" also means copyright-like laws that apply to other kinds of+works, such as semiconductor masks.++ "The Program" refers to any copyrightable work licensed under this+License. Each licensee is addressed as "you". "Licensees" and+"recipients" may be individuals or organizations.++ To "modify" a work means to copy from or adapt all or part of the work+in a fashion requiring copyright permission, other than the making of an+exact copy. The resulting work is called a "modified version" of the+earlier work or a work "based on" the earlier work.++ A "covered work" means either the unmodified Program or a work based+on the Program.++ To "propagate" a work means to do anything with it that, without+permission, would make you directly or secondarily liable for+infringement under applicable copyright law, except executing it on a+computer or modifying a private copy. Propagation includes copying,+distribution (with or without modification), making available to the+public, and in some countries other activities as well.++ To "convey" a work means any kind of propagation that enables other+parties to make or receive copies. Mere interaction with a user through+a computer network, with no transfer of a copy, is not conveying.++ An interactive user interface displays "Appropriate Legal Notices"+to the extent that it includes a convenient and prominently visible+feature that (1) displays an appropriate copyright notice, and (2)+tells the user that there is no warranty for the work (except to the+extent that warranties are provided), that licensees may convey the+work under this License, and how to view a copy of this License. If+the interface presents a list of user commands or options, such as a+menu, a prominent item in the list meets this criterion.++ 1. Source Code.++ The "source code" for a work means the preferred form of the work+for making modifications to it. "Object code" means any non-source+form of a work.++ A "Standard Interface" means an interface that either is an official+standard defined by a recognized standards body, or, in the case of+interfaces specified for a particular programming language, one that+is widely used among developers working in that language.++ The "System Libraries" of an executable work include anything, other+than the work as a whole, that (a) is included in the normal form of+packaging a Major Component, but which is not part of that Major+Component, and (b) serves only to enable use of the work with that+Major Component, or to implement a Standard Interface for which an+implementation is available to the public in source code form. A+"Major Component", in this context, means a major essential component+(kernel, window system, and so on) of the specific operating system+(if any) on which the executable work runs, or a compiler used to+produce the work, or an object code interpreter used to run it.++ The "Corresponding Source" for a work in object code form means all+the source code needed to generate, install, and (for an executable+work) run the object code and to modify the work, including scripts to+control those activities. However, it does not include the work's+System Libraries, or general-purpose tools or generally available free+programs which are used unmodified in performing those activities but+which are not part of the work. For example, Corresponding Source+includes interface definition files associated with source files for+the work, and the source code for shared libraries and dynamically+linked subprograms that the work is specifically designed to require,+such as by intimate data communication or control flow between those+subprograms and other parts of the work.++ The Corresponding Source need not include anything that users+can regenerate automatically from other parts of the Corresponding+Source.++ The Corresponding Source for a work in source code form is that+same work.++ 2. Basic Permissions.++ All rights granted under this License are granted for the term of+copyright on the Program, and are irrevocable provided the stated+conditions are met. This License explicitly affirms your unlimited+permission to run the unmodified Program. The output from running a+covered work is covered by this License only if the output, given its+content, constitutes a covered work. This License acknowledges your+rights of fair use or other equivalent, as provided by copyright law.++ You may make, run and propagate covered works that you do not+convey, without conditions so long as your license otherwise remains+in force. You may convey covered works to others for the sole purpose+of having them make modifications exclusively for you, or provide you+with facilities for running those works, provided that you comply with+the terms of this License in conveying all material for which you do+not control copyright. Those thus making or running the covered works+for you must do so exclusively on your behalf, under your direction+and control, on terms that prohibit them from making any copies of+your copyrighted material outside their relationship with you.++ Conveying under any other circumstances is permitted solely under+the conditions stated below. Sublicensing is not allowed; section 10+makes it unnecessary.++ 3. Protecting Users' Legal Rights From Anti-Circumvention Law.++ No covered work shall be deemed part of an effective technological+measure under any applicable law fulfilling obligations under article+11 of the WIPO copyright treaty adopted on 20 December 1996, or+similar laws prohibiting or restricting circumvention of such+measures.++ When you convey a covered work, you waive any legal power to forbid+circumvention of technological measures to the extent such circumvention+is effected by exercising rights under this License with respect to+the covered work, and you disclaim any intention to limit operation or+modification of the work as a means of enforcing, against the work's+users, your or third parties' legal rights to forbid circumvention of+technological measures.++ 4. Conveying Verbatim Copies.++ You may convey verbatim copies of the Program's source code as you+receive it, in any medium, provided that you conspicuously and+appropriately publish on each copy an appropriate copyright notice;+keep intact all notices stating that this License and any+non-permissive terms added in accord with section 7 apply to the code;+keep intact all notices of the absence of any warranty; and give all+recipients a copy of this License along with the Program.++ You may charge any price or no price for each copy that you convey,+and you may offer support or warranty protection for a fee.++ 5. Conveying Modified Source Versions.++ You may convey a work based on the Program, or the modifications to+produce it from the Program, in the form of source code under the+terms of section 4, provided that you also meet all of these conditions:++ a) The work must carry prominent notices stating that you modified+ it, and giving a relevant date.++ b) The work must carry prominent notices stating that it is+ released under this License and any conditions added under section+ 7. This requirement modifies the requirement in section 4 to+ "keep intact all notices".++ c) You must license the entire work, as a whole, under this+ License to anyone who comes into possession of a copy. This+ License will therefore apply, along with any applicable section 7+ additional terms, to the whole of the work, and all its parts,+ regardless of how they are packaged. This License gives no+ permission to license the work in any other way, but it does not+ invalidate such permission if you have separately received it.++ d) If the work has interactive user interfaces, each must display+ Appropriate Legal Notices; however, if the Program has interactive+ interfaces that do not display Appropriate Legal Notices, your+ work need not make them do so.++ A compilation of a covered work with other separate and independent+works, which are not by their nature extensions of the covered work,+and which are not combined with it such as to form a larger program,+in or on a volume of a storage or distribution medium, is called an+"aggregate" if the compilation and its resulting copyright are not+used to limit the access or legal rights of the compilation's users+beyond what the individual works permit. Inclusion of a covered work+in an aggregate does not cause this License to apply to the other+parts of the aggregate.++ 6. Conveying Non-Source Forms.++ You may convey a covered work in object code form under the terms+of sections 4 and 5, provided that you also convey the+machine-readable Corresponding Source under the terms of this License,+in one of these ways:++ a) Convey the object code in, or embodied in, a physical product+ (including a physical distribution medium), accompanied by the+ Corresponding Source fixed on a durable physical medium+ customarily used for software interchange.++ b) Convey the object code in, or embodied in, a physical product+ (including a physical distribution medium), accompanied by a+ written offer, valid for at least three years and valid for as+ long as you offer spare parts or customer support for that product+ model, to give anyone who possesses the object code either (1) a+ copy of the Corresponding Source for all the software in the+ product that is covered by this License, on a durable physical+ medium customarily used for software interchange, for a price no+ more than your reasonable cost of physically performing this+ conveying of source, or (2) access to copy the+ Corresponding Source from a network server at no charge.++ c) Convey individual copies of the object code with a copy of the+ written offer to provide the Corresponding Source. This+ alternative is allowed only occasionally and noncommercially, and+ only if you received the object code with such an offer, in accord+ with subsection 6b.++ d) Convey the object code by offering access from a designated+ place (gratis or for a charge), and offer equivalent access to the+ Corresponding Source in the same way through the same place at no+ further charge. You need not require recipients to copy the+ Corresponding Source along with the object code. If the place to+ copy the object code is a network server, the Corresponding Source+ may be on a different server (operated by you or a third party)+ that supports equivalent copying facilities, provided you maintain+ clear directions next to the object code saying where to find the+ Corresponding Source. Regardless of what server hosts the+ Corresponding Source, you remain obligated to ensure that it is+ available for as long as needed to satisfy these requirements.++ e) Convey the object code using peer-to-peer transmission, provided+ you inform other peers where the object code and Corresponding+ Source of the work are being offered to the general public at no+ charge under subsection 6d.++ A separable portion of the object code, whose source code is excluded+from the Corresponding Source as a System Library, need not be+included in conveying the object code work.++ A "User Product" is either (1) a "consumer product", which means any+tangible personal property which is normally used for personal, family,+or household purposes, or (2) anything designed or sold for incorporation+into a dwelling. In determining whether a product is a consumer product,+doubtful cases shall be resolved in favor of coverage. For a particular+product received by a particular user, "normally used" refers to a+typical or common use of that class of product, regardless of the status+of the particular user or of the way in which the particular user+actually uses, or expects or is expected to use, the product. A product+is a consumer product regardless of whether the product has substantial+commercial, industrial or non-consumer uses, unless such uses represent+the only significant mode of use of the product.++ "Installation Information" for a User Product means any methods,+procedures, authorization keys, or other information required to install+and execute modified versions of a covered work in that User Product from+a modified version of its Corresponding Source. The information must+suffice to ensure that the continued functioning of the modified object+code is in no case prevented or interfered with solely because+modification has been made.++ If you convey an object code work under this section in, or with, or+specifically for use in, a User Product, and the conveying occurs as+part of a transaction in which the right of possession and use of the+User Product is transferred to the recipient in perpetuity or for a+fixed term (regardless of how the transaction is characterized), the+Corresponding Source conveyed under this section must be accompanied+by the Installation Information. But this requirement does not apply+if neither you nor any third party retains the ability to install+modified object code on the User Product (for example, the work has+been installed in ROM).++ The requirement to provide Installation Information does not include a+requirement to continue to provide support service, warranty, or updates+for a work that has been modified or installed by the recipient, or for+the User Product in which it has been modified or installed. Access to a+network may be denied when the modification itself materially and+adversely affects the operation of the network or violates the rules and+protocols for communication across the network.++ Corresponding Source conveyed, and Installation Information provided,+in accord with this section must be in a format that is publicly+documented (and with an implementation available to the public in+source code form), and must require no special password or key for+unpacking, reading or copying.++ 7. Additional Terms.++ "Additional permissions" are terms that supplement the terms of this+License by making exceptions from one or more of its conditions.+Additional permissions that are applicable to the entire Program shall+be treated as though they were included in this License, to the extent+that they are valid under applicable law. If additional permissions+apply only to part of the Program, that part may be used separately+under those permissions, but the entire Program remains governed by+this License without regard to the additional permissions.++ When you convey a copy of a covered work, you may at your option+remove any additional permissions from that copy, or from any part of+it. (Additional permissions may be written to require their own+removal in certain cases when you modify the work.) You may place+additional permissions on material, added by you to a covered work,+for which you have or can give appropriate copyright permission.++ Notwithstanding any other provision of this License, for material you+add to a covered work, you may (if authorized by the copyright holders of+that material) supplement the terms of this License with terms:++ a) Disclaiming warranty or limiting liability differently from the+ terms of sections 15 and 16 of this License; or++ b) Requiring preservation of specified reasonable legal notices or+ author attributions in that material or in the Appropriate Legal+ Notices displayed by works containing it; or++ c) Prohibiting misrepresentation of the origin of that material, or+ requiring that modified versions of such material be marked in+ reasonable ways as different from the original version; or++ d) Limiting the use for publicity purposes of names of licensors or+ authors of the material; or++ e) Declining to grant rights under trademark law for use of some+ trade names, trademarks, or service marks; or++ f) Requiring indemnification of licensors and authors of that+ material by anyone who conveys the material (or modified versions of+ it) with contractual assumptions of liability to the recipient, for+ any liability that these contractual assumptions directly impose on+ those licensors and authors.++ All other non-permissive additional terms are considered "further+restrictions" within the meaning of section 10. If the Program as you+received it, or any part of it, contains a notice stating that it is+governed by this License along with a term that is a further+restriction, you may remove that term. If a license document contains+a further restriction but permits relicensing or conveying under this+License, you may add to a covered work material governed by the terms+of that license document, provided that the further restriction does+not survive such relicensing or conveying.++ If you add terms to a covered work in accord with this section, you+must place, in the relevant source files, a statement of the+additional terms that apply to those files, or a notice indicating+where to find the applicable terms.++ Additional terms, permissive or non-permissive, may be stated in the+form of a separately written license, or stated as exceptions;+the above requirements apply either way.++ 8. Termination.++ You may not propagate or modify a covered work except as expressly+provided under this License. Any attempt otherwise to propagate or+modify it is void, and will automatically terminate your rights under+this License (including any patent licenses granted under the third+paragraph of section 11).++ However, if you cease all violation of this License, then your+license from a particular copyright holder is reinstated (a)+provisionally, unless and until the copyright holder explicitly and+finally terminates your license, and (b) permanently, if the copyright+holder fails to notify you of the violation by some reasonable means+prior to 60 days after the cessation.++ Moreover, your license from a particular copyright holder is+reinstated permanently if the copyright holder notifies you of the+violation by some reasonable means, this is the first time you have+received notice of violation of this License (for any work) from that+copyright holder, and you cure the violation prior to 30 days after+your receipt of the notice.++ Termination of your rights under this section does not terminate the+licenses of parties who have received copies or rights from you under+this License. If your rights have been terminated and not permanently+reinstated, you do not qualify to receive new licenses for the same+material under section 10.++ 9. Acceptance Not Required for Having Copies.++ You are not required to accept this License in order to receive or+run a copy of the Program. Ancillary propagation of a covered work+occurring solely as a consequence of using peer-to-peer transmission+to receive a copy likewise does not require acceptance. However,+nothing other than this License grants you permission to propagate or+modify any covered work. These actions infringe copyright if you do+not accept this License. Therefore, by modifying or propagating a+covered work, you indicate your acceptance of this License to do so.++ 10. Automatic Licensing of Downstream Recipients.++ Each time you convey a covered work, the recipient automatically+receives a license from the original licensors, to run, modify and+propagate that work, subject to this License. You are not responsible+for enforcing compliance by third parties with this License.++ An "entity transaction" is a transaction transferring control of an+organization, or substantially all assets of one, or subdividing an+organization, or merging organizations. If propagation of a covered+work results from an entity transaction, each party to that+transaction who receives a copy of the work also receives whatever+licenses to the work the party's predecessor in interest had or could+give under the previous paragraph, plus a right to possession of the+Corresponding Source of the work from the predecessor in interest, if+the predecessor has it or can get it with reasonable efforts.++ You may not impose any further restrictions on the exercise of the+rights granted or affirmed under this License. For example, you may+not impose a license fee, royalty, or other charge for exercise of+rights granted under this License, and you may not initiate litigation+(including a cross-claim or counterclaim in a lawsuit) alleging that+any patent claim is infringed by making, using, selling, offering for+sale, or importing the Program or any portion of it.++ 11. Patents.++ A "contributor" is a copyright holder who authorizes use under this+License of the Program or a work on which the Program is based. The+work thus licensed is called the contributor's "contributor version".++ A contributor's "essential patent claims" are all patent claims+owned or controlled by the contributor, whether already acquired or+hereafter acquired, that would be infringed by some manner, permitted+by this License, of making, using, or selling its contributor version,+but do not include claims that would be infringed only as a+consequence of further modification of the contributor version. For+purposes of this definition, "control" includes the right to grant+patent sublicenses in a manner consistent with the requirements of+this License.++ Each contributor grants you a non-exclusive, worldwide, royalty-free+patent license under the contributor's essential patent claims, to+make, use, sell, offer for sale, import and otherwise run, modify and+propagate the contents of its contributor version.++ In the following three paragraphs, a "patent license" is any express+agreement or commitment, however denominated, not to enforce a patent+(such as an express permission to practice a patent or covenant not to+sue for patent infringement). To "grant" such a patent license to a+party means to make such an agreement or commitment not to enforce a+patent against the party.++ If you convey a covered work, knowingly relying on a patent license,+and the Corresponding Source of the work is not available for anyone+to copy, free of charge and under the terms of this License, through a+publicly available network server or other readily accessible means,+then you must either (1) cause the Corresponding Source to be so+available, or (2) arrange to deprive yourself of the benefit of the+patent license for this particular work, or (3) arrange, in a manner+consistent with the requirements of this License, to extend the patent+license to downstream recipients. "Knowingly relying" means you have+actual knowledge that, but for the patent license, your conveying the+covered work in a country, or your recipient's use of the covered work+in a country, would infringe one or more identifiable patents in that+country that you have reason to believe are valid.++ If, pursuant to or in connection with a single transaction or+arrangement, you convey, or propagate by procuring conveyance of, a+covered work, and grant a patent license to some of the parties+receiving the covered work authorizing them to use, propagate, modify+or convey a specific copy of the covered work, then the patent license+you grant is automatically extended to all recipients of the covered+work and works based on it.++ A patent license is "discriminatory" if it does not include within+the scope of its coverage, prohibits the exercise of, or is+conditioned on the non-exercise of one or more of the rights that are+specifically granted under this License. You may not convey a covered+work if you are a party to an arrangement with a third party that is+in the business of distributing software, under which you make payment+to the third party based on the extent of your activity of conveying+the work, and under which the third party grants, to any of the+parties who would receive the covered work from you, a discriminatory+patent license (a) in connection with copies of the covered work+conveyed by you (or copies made from those copies), or (b) primarily+for and in connection with specific products or compilations that+contain the covered work, unless you entered into that arrangement,+or that patent license was granted, prior to 28 March 2007.++ Nothing in this License shall be construed as excluding or limiting+any implied license or other defenses to infringement that may+otherwise be available to you under applicable patent law.++ 12. No Surrender of Others' Freedom.++ If conditions are imposed on you (whether by court order, agreement or+otherwise) that contradict the conditions of this License, they do not+excuse you from the conditions of this License. If you cannot convey a+covered work so as to satisfy simultaneously your obligations under this+License and any other pertinent obligations, then as a consequence you may+not convey it at all. For example, if you agree to terms that obligate you+to collect a royalty for further conveying from those to whom you convey+the Program, the only way you could satisfy both those terms and this+License would be to refrain entirely from conveying the Program.++ 13. Use with the GNU Affero General Public License.++ Notwithstanding any other provision of this License, you have+permission to link or combine any covered work with a work licensed+under version 3 of the GNU Affero General Public License into a single+combined work, and to convey the resulting work. The terms of this+License will continue to apply to the part which is the covered work,+but the special requirements of the GNU Affero General Public License,+section 13, concerning interaction through a network will apply to the+combination as such.++ 14. Revised Versions of this License.++ The Free Software Foundation may publish revised and/or new versions of+the GNU General Public License from time to time. Such new versions will+be similar in spirit to the present version, but may differ in detail to+address new problems or concerns.++ Each version is given a distinguishing version number. If the+Program specifies that a certain numbered version of the GNU General+Public License "or any later version" applies to it, you have the+option of following the terms and conditions either of that numbered+version or of any later version published by the Free Software+Foundation. If the Program does not specify a version number of the+GNU General Public License, you may choose any version ever published+by the Free Software Foundation.++ If the Program specifies that a proxy can decide which future+versions of the GNU General Public License can be used, that proxy's+public statement of acceptance of a version permanently authorizes you+to choose that version for the Program.++ Later license versions may give you additional or different+permissions. However, no additional obligations are imposed on any+author or copyright holder as a result of your choosing to follow a+later version.++ 15. Disclaimer of Warranty.++ THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY+APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR+PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM+IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.++ 16. Limitation of Liability.++ IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF+SUCH DAMAGES.++ 17. Interpretation of Sections 15 and 16.++ If the disclaimer of warranty and limitation of liability provided+above cannot be given local legal effect according to their terms,+reviewing courts shall apply local law that most closely approximates+an absolute waiver of all civil liability in connection with the+Program, unless a warranty or assumption of liability accompanies a+copy of the Program in return for a fee.++ END OF TERMS AND CONDITIONS++ How to Apply These Terms to Your New Programs++ If you develop a new program, and you want it to be of the greatest+possible use to the public, the best way to achieve this is to make it+free software which everyone can redistribute and change under these terms.++ To do so, attach the following notices to the program. It is safest+to attach them to the start of each source file to most effectively+state the exclusion of warranty; and each file should have at least+the "copyright" line and a pointer to where the full notice is found.++ {one line to give the program's name and a brief idea of what it does.}+ Copyright (C) {year} {name of author}++ This program is free software: you can redistribute it and/or modify+ it under the terms of the GNU General Public License as published by+ the Free Software Foundation, either version 3 of the License, or+ (at your option) any later version.++ This program is distributed in the hope that it will be useful,+ but WITHOUT ANY WARRANTY; without even the implied warranty of+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the+ GNU General Public License for more details.++ You should have received a copy of the GNU General Public License+ along with this program. If not, see <http://www.gnu.org/licenses/>.++Also add information on how to contact you by electronic and paper mail.++ If the program does terminal interaction, make it output a short+notice like this when it starts in an interactive mode:++ {project} Copyright (C) {year} {fullname}+ 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>.+
+ Math/NumberTheory/Canon.hs view
@@ -0,0 +1,879 @@+-- |+-- Module: Math.NumberTheory.Canon+-- Copyright: (c) 2015-2018 Frederick Schneider+-- Licence: MIT+-- Maintainer: Frederick Schneider <frederick.schneider2011@gmail.com> +-- Stability: Provisional+--+-- A Canon is an exponentation-based representation for arbitrarily massive numbers, including prime towers.++{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, PatternSynonyms, ViewPatterns, RankNTypes #-}++module Math.NumberTheory.Canon ( + makeCanon, makeC,+ canonToGCR, cToGCR,++ cMult, cDiv, cAdd, cSubtract, cExp,+ cReciprocal,+ cGCD, cLCM, cMod, cOdd, cTotient, cPhi,+ cLog, cLogDouble,+ cNegative, cPositive,+ cIntegral, cRational, cIrrational,+ cSimplify, cSimplified, + cDepth, + cSplit, cNumerator, cDenominator,+ cCanonical, cBare, cBareStatus, cValueType,+ cIsPrimeTower, cPrimeTowerLevel,+ + cTetration, cPentation, cHexation, cHyperOp,+ (>^), (<^), (%), (<^>), (<<^>>), (<<<^>>>) +)+where++import Math.NumberTheory.Primes.Testing (isPrime)+import Data.List (intersperse)+import GHC.Real (Ratio(..))+import Math.NumberTheory.Canon.Internals+import Math.NumberTheory.Canon.Additive+import Math.NumberTheory.Canon.AurifCyclo (CycloMap, crCycloInitMap)+import Math.NumberTheory.Canon.Simple (CanonConv(..))++-- | CanonValueType: 3 possibilities for this GADT. Imaginary/complex numbers are not supported+data CanonValueType = IntegralC | NonIntRationalC | IrrationalC deriving (Eq, Ord, Show)++-- | GCR_ stands for Generalized Canonical Representation+type GCR_ = [GCRE_]+type GCRE_ = (Integer, Canon)++-- | Canon: GADT for either Bare or some variation of a canonical form.+data Canon = Bare Integer BareStatus | Canonical GCR_ CanonValueType ++-- | BareStatus: A "Bare Simplified" number means a prime number, +/-1 or 0. The code must set the flag properly+-- A "Bare NotSimplified" number is an integer that has not been checked (to see if it can be factored).+data BareStatus = Simplified | NotSimplified deriving (Eq, Ord, Show)++makeCanon, makeC, makeCanonFull, makeDefCanonForExpnt :: Integer -> Canon++-- | Create a Canon from an integer. This may involve expensive factorization.+makeCanon n = makeCanonI n False++-- | Shorthand for makeCanon+makeC = makeCanon++-- | Make a Canon and attempt a full factorization+makeCanonFull n = makeCanonI n True++makeCanonI :: Integer -> Bool -> Canon+makeCanonI n b = crToC (crFromI n) b +-- TODO: next step: enhance this once we can partially factor numbers++cCutoff :: Integer+cCutoff = 1000000++-- | Create type of Canon based on whether it exceeds a cutoff+makeDefCanonForExpnt e | e > cCutoff = Bare e (getBareStatus e) + | otherwise = makeCanonFull e++-- | Convert from underlying canonical rep. to Canon. The 2nd param indicates whether or not to force factorization/simplification.+crToC :: CR_ -> Bool -> Canon+crToC POne _ = Bare 1 Simplified+crToC c b | crSimplified c = Bare (fst $ head c) Simplified -- a little ugly+ | otherwise = Canonical g (gcrCVT g)+ where g = map (\(p,e) -> (p, convPred e)) c+ convPred e | b = makeCanonFull e -- do complete factorization+ | otherwise = makeDefCanonForExpnt e+ -- Leave exponents "Bare" with flag based on if whether it's "simplified"+ -- (can't be reduced any further)++-- | Instances for Canon+instance Eq Canon where+ x == y = cEq x y++instance Show Canon where + show (Bare n NotSimplified) = "(" ++ (show n) ++ ")" -- Note the extra characters. This does not mean the figure is negative.+ show (Bare n Simplified) = show n+ show c | denom == c1 = s numer False + | otherwise = s numer True ++ " / " ++ s denom True+ where (numer, denom) = cSplit c + s (Bare n f) _ = show $ Bare n f+ s v w | w = "(" ++ catList ++ ")" + | otherwise = catList -- w = with(out) parens+ where catList = concat $ intersperse " * " $ map sE $ cToGCR v -- sE means showElem+ sE (p, e) | ptLevel > 2 = sp ++ " <^> " ++ s (makeCanonFull $ ptLevel) True+ | otherwise = case e of+ Bare 1 _ -> sp + Bare _ _ -> sp ++ "^" ++ se+ _ -> sp ++ " ^ (" ++ se ++ ")"+ where ptLevel = cPrimeTowerLevelI e p 1 + sp = show p+ se = show e+ -- TODO: try for add'l collapse into <<^>>+ +instance Enum Canon where+ toEnum n = makeCanon $ fromIntegral n+ fromEnum c = fromIntegral $ cToI c++instance Ord Canon where+ compare x y = cCmp x y ++instance Real Canon where+ toRational c | cIrrational c = toRational $ cToD c+ | otherwise = (cToI $ cNumerator c) :% (cToI $ cDenominator c)++instance Integral Canon where+ toInteger c | cIntegral c = cToI c + | otherwise = floor $ cToD c+ quotRem n m = fst $ cQuotRem n m crCycloInitMap -- tries to use map but ultimately throws it away + -- ToDo: mod n m = fst $ cModBAD n m crCycloInitMap -- fix this "bad" logic and use this instead of the original function+ mod n m = cMod n m+ +instance Fractional Canon where+ fromRational (n :% d) = makeCanon n / makeCanon d + (/) x y = fst $ cDiv x y crCycloInitMap -- tries to use map but ultimately throws it away++instance Num Canon where -- tries to use the map but ultimately throws it away when using +, - and * operators+ fromInteger n = makeCanon n+ x + y = fst $ cAdd x y crCycloInitMap+ x - y = fst $ cSubtract x y crCycloInitMap+ x * y = fst $ cMult x y crCycloInitMap+ + negate x = cNegate x+ abs x = cAbs x+ signum x = cSignum x++-- | Is the Canon a more complex expression? +cCanonical :: Canon -> Bool+cCanonical (Canonical _ _ ) = True+cCanonical _ = False++-- | Checks if the Canon just a "Bare" Integer.+cBare :: Canon -> Bool+cBare (Bare _ _ ) = True+cBare _ = False++-- | Returns the status for "Bare" numbers.+cBareStatus :: Canon -> BareStatus+cBareStatus (Bare _ b) = b+cBareStatus _ = error "cBareStatus: Can only checked for 'Bare' Canons"++-- | Return the CanonValueType (Integral, etc).+cValueType :: Canon -> CanonValueType+cValueType (Bare _ _ ) = IntegralC+cValueType (Canonical _ v ) = v++-- | Split a Canon into the numerator and denominator.+cSplit :: Canon -> (Canon, Canon)+cSplit c = (cNumerator c, cDenominator c)++-- | Check for equality.+cEq:: Canon -> Canon -> Bool +cEq (Bare x _ ) (Bare y _ ) = x == y+cEq (Bare _ Simplified) (Canonical _ _ ) = False+cEq (Canonical _ _ ) (Bare _ Simplified) = False++cEq (Bare x NotSimplified) y | cValueType y /= IntegralC = False+ | otherwise = cEq (makeCanon x) y++cEq x (Bare y NotSimplified) | cValueType x /= IntegralC = False+ | otherwise = cEq x (makeCanon y)++cEq (Canonical x a ) (Canonical y b) = if a /= b then False else gcrEqCheck x y++-- | Check if a Canon is an odd integer. Note: If the Canon is not integral, return False +cOdd :: Canon -> Bool+cOdd (Bare x _) = mod x 2 == 1+cOdd (Canonical c IntegralC ) = gcrOdd c+cOdd (Canonical _ _ ) = False++-- | GCD and LCM functions for Canon+cGCD, cLCM :: Canon -> Canon -> Canon+cGCD x y = cLGApply gcrGCD x y+cLCM x y = cLGApply gcrLCM x y++-- | Compute log as a Rational number.+cLog :: Canon -> Rational+cLog x = gcrLog $ cToGCR x ++-- | Compute log as a Double.+cLogDouble :: Canon -> Double+cLogDouble x = gcrLogDouble $ cToGCR x ++-- | Compare Function+cCmp :: Canon -> Canon -> Ordering+cCmp (Bare x _) (Bare y _) = compare x y+cCmp x y = gcrCmp (cToGCR x) (cToGCR y)++-- | QuotRem Function+cQuotRem :: Canon -> Canon -> CycloMap -> ((Canon, Canon), CycloMap)+cQuotRem x y m | cIntegral x && cIntegral y = ((gcrToC q', md'), m'')+ | otherwise = error "cQuotRem: Must both parameters must be integral"+ where (q', md', m'') = case gcrDiv (cToGCR x) gy of+ -- ToDo: Left _ -> (q, md, mm) -- fix "cModBAD" and stop pointing to orig fcn+ Left _ -> (q, md, m')+ Right quotient -> (quotient, c0, m)+ where gy = cToGCR y+ -- ToDo: fix (md, mm) = cModBAD x y m' -- Better to compute quotient this way .. to take adv. of alg. forms+ md = cMod x y+ q = gcrDivStrict (cToGCR d) gy -- equivalent to: (x - x%y) / y.+ (d, m') = cSubtract x md m++-- | Mod function+cMod :: Canon -> Canon -> Canon+cMod c m = if (cIntegral c) && (cIntegral m) then (makeCanon $ cModI c (cToI m))+ else error "cMod: Must both parameters must be integral"++cModI :: Canon -> Integer -> Integer+cModI _ 0 = error "cModI: Divide by zero error when computing n mod 0"+cModI _ 1 = 0+cModI _ (-1) = 0+cModI Pc1 PIntPos = 1+cModI Pc0 _ = 0+cModI c m | cn && mn = -1 * cModI (cAbs c) (abs m)+ | (cn && not mn) ||+ (mn && not cn) = (signum m) * ( (abs m) - (cModI' (cAbs c) (abs m)) )+ | otherwise = cModI' c m+ where cn = cNegative c+ mn = m < 0++cModI' :: Canon -> Integer -> Integer+cModI' (Bare n _ ) m = mod n m+cModI' (Canonical c IntegralC ) m = mod (product $ map (\x -> pmI (fst x) (mmt $ snd x) m) c) m+ where tm = totient m+ mmt e = cModI e tm -- optimization+cModI' (Canonical _ _ ) _ = error "cModI': Logic error: Canonical var has to be integral at this point" ++-- | Totient functions+cTotient, cPhi :: Canon -> CycloMap -> (Canon, CycloMap)+cTotient c m | (not $ cIntegral c) || cNegative c = error "Not defined for non-integral or negative numbers"+ | c == c0 = (c0, m)+ | otherwise = f (cToGCR c) c1 m+ where f [] prd m' = (prd, m') + f ((p,e):gs) prd m' = f gs wp mw + -- f is equivalent to the crTotient function but with threading of CycloMap + -- => product $ map (\(p,e) -> (p-1) * p^(e-1)) cr+ where cp = makeC p -- "Canon-ize" this. Generally, this should be a prime already+ (pM1, mp) = cSubtract cp c1 m'+ (eM1, me) = cSubtract e c1 mp + (pxeM1, mpm) = cExp cp eM1 False me+ (nprd, mprd) = cMult pM1 pxeM1 mpm + (wp, mw) = cMult prd nprd mprd++cPhi = cTotient++-- | Hyperoperations (including tetration and beyond): https://en.wikipedia.org/wiki/Hyperoperation+-- | The thinking around the operators is that they should look progressively scarier :)+infixr 9 <^>, <<^>>, <<<^>>>+(<^>), (<<^>>), (<<<^>>>) :: Canon -> Integer -> Canon+a <^> b = fst $ cTetration a b crCycloInitMap+a <<^>> b = fst $ cPentation a b crCycloInitMap+a <<<^>>> b = fst $ cHexation a b crCycloInitMap++cTetration, cPentation, cHexation :: Canon -> Integer -> CycloMap -> (Canon, CycloMap)++-- | Tetration function+cTetration = cHyperOp 4 ++-- | Pentation Function+cPentation = cHyperOp 5++-- | Hexation Function+cHexation = cHyperOp 6++-- | Generalized Hyperoperation Function+cHyperOp :: Integer -> Canon -> Integer -> CycloMap -> (Canon, CycloMap)+cHyperOp n a b m | b < -1 = error "Hyperoperations not defined when b < -1"+ | n < 0 = error "Hyperoperations require the level n >= 0"+ | a /= c0 && a /= c1 && + b > 1 && (a /= c2 && b == 2) = c n cb m+ | otherwise = (sp n a b, m)+ where cb = makeCanon b+ -- Function for regular cases+ c 1 b' m' = cAdd a b' m' -- addition+ c 2 b' m' = cMult a b' m' -- multiplication+ c 3 b' m' = (a <^ b', m') -- exponentiation: ToDo: Plug in the CycloMap logic for expo.+ -- Tetration and beyond+ c _ Pc1 m' = (a, m')+ c n' b' m' = c (n'-1) r m'' -- TODO: Find a way to optimize this+ where (r, m'') = c n' (b'-1) m'++ -- Function for special cases: + -- Note: When n (first param) is zero, that is known as "succession"+ -- Cases when a is zero ...+ sp 0 Pc0 b' = makeCanon (b' + 1)+ sp 1 Pc0 b' = makeCanon b'+ sp 2 Pc0 _ = c0+ sp 3 Pc0 b' = if b' == 0 then c1 else c0+ sp _ Pc0 b' = if (mod b' 2) == 1 then c0 else c1+ -- Cases when b is zero ...+ sp 0 _ 0 = c1 + sp 1 a' 0 = a'+ sp 2 _ 0 = c0 + sp _ _ 0 = c1 + -- Cases when b is -1 ...+ sp 0 _ (-1) = c0+ sp 1 a' (-1) = a' - 1+ sp 2 a' (-1) = cNegate a'+ sp 3 a' (-1) = cReciprocal a'+ sp _ _ (-1) = c0+ -- Other Cases ...+ sp h Pc2 2 | h == 0 = makeCanon 3+ | otherwise = makeCanon 4 -- recursive identity+ sp _ Pc1 _ = c1+ sp _ a' 1 = a'+ sp _ _ _ = error "Can't compute this hyperoperation. b must be >= -1" ++infixl 7 %+-- | Mod operator+(%) :: (Integral a) => a -> a -> a+n % m = mod n m ++-- | Exponentation operator declaration+infixr 9 <^ +-- Note: Even with Flexible Contexts switched on, it doesn't infer a bare number to be an Integer++-- | Dedicated multi-param typeclass for exponentiation operator.+class CanonExpnt a b c | a b -> c where + -- | Exponentiation operator+ (<^) :: a -> b -> c++instance CanonExpnt Canon Canon Canon where+ p <^ e = fst $ cExp p e True crCycloInitMap + +instance CanonExpnt Integer Integer Canon where+ p <^ e = fst $ cExp (makeCanon p) (makeDefCanonForExpnt e) True crCycloInitMap++instance CanonExpnt Canon Integer Canon where+ p <^ e = fst $ cExp p (makeDefCanonForExpnt e) True crCycloInitMap++instance CanonExpnt Integer Canon Canon where+ p <^ e = fst $ cExp (makeCanon p) e True crCycloInitMap++-- | Operator declaration: r >^ n means: attempt to take the rth root of n +infixr 9 >^ ++-- | Dedicated multi-param typeclass for radical or root operator.+class CanonRoot a b c | a b -> c where + -- | Root operator+ (>^) :: a -> b -> c++instance CanonRoot Canon Canon Canon where+ r >^ n = cRoot n r+ +instance CanonRoot Integer Integer Canon where+ r >^ n = cRoot (makeCanon n) (makeCanon r)+ +instance CanonRoot Integer Canon Canon where+ r >^ n = cRoot n (makeCanon r) ++instance CanonRoot Canon Integer Canon where+ r >^ n = cRoot (makeCanon n) r ++-- | Check if underlying rep is simplified+crSimplified :: CR_ -> Bool+crSimplified POne = True+crSimplified PZero = True +crSimplified PN1 = True +crSimplified c = crPrime c++-- | Convert a Canon back to its underlying rep (if possible).+cToCR :: Canon -> CR_+cToCR (Canonical c v) | v /= IrrationalC = gcrToCR c + | otherwise = error "cToCR: Cannot convert irrational canons to underlying data structure"+cToCR (Bare 1 _ ) = cr1+cToCR (Bare n NotSimplified) = crFromI n+cToCR (Bare n Simplified) = [(n,1)] -- not ideal++-- | Convert generalized canon rep to Canon.+gcrToC :: GCR_ -> Canon+gcrToC g | gcrBare g = Bare (gcrToI g) Simplified+ | otherwise = Canonical g (gcrCVT g)++-- | For generalized canon rep, determine the CanonValueType. +gcrCVT :: GCR_ -> CanonValueType +gcrCVT POne = IntegralC+gcrCVT g = g' g IntegralC -- start Integral, can only get "worse"+ where g' _ IrrationalC = IrrationalC -- short-circuits once IrrationalC is found+ g' POne v = v+ g' ((_,ce):cs) v = g' cs (dcv v ce) -- check the exponents for expr's value type+ g' _ _ = error "gcrCVT : Logic error. Patterns should have been exhaustive"++ -- checking exponents+ dcv IrrationalC _ = IrrationalC+ dcv _ (Canonical _ IrrationalC) = IrrationalC+ dcv _ (Canonical _ NonIntRationalC) = IrrationalC+ dcv IntegralC (Bare n _ ) = if n < 0 then NonIntRationalC else IntegralC+ dcv v (Bare _ _ ) = v+ dcv v c = if cNegative c then NonIntRationalC else v++c1, c0, cN1, c2 :: Canon+c1 = makeCanon 1+c0 = makeCanon 0+cN1 = makeCanon (-1)+c2 = makeCanon 2++-- | Convert Canon to Integer if possible.+cToI :: Canon -> Integer+cToI (Bare i _ ) = i+cToI (Canonical c v) | v == IntegralC = gcrToI c + | otherwise = error "Can't convert non-integral canon to an integer"++-- | Convert Canon To Double.+cToD :: Canon -> Double+cToD (Bare i _ ) = fromIntegral i+cToD (Canonical c _ ) = gcrToD c ++-- | Multiply Function: Generally speaking this will be much cheaper.+cMult :: Canon -> Canon -> CycloMap -> (Canon, CycloMap) +cMult Pc0 _ m = (c0, m)+cMult _ Pc0 m = (c0, m)+cMult Pc1 y m = (y, m)+cMult x Pc1 m = (x, m)+cMult x y m = (gcrToC g, m') + where (g, m') = gcrMult (cToGCR x) (cToGCR y) m++-- | Addition and subtraction is generally much more expensive because it requires refactorization.+-- There is logic to look for algebraic forms which can greatly reduce simplify factorization.+cAdd, cSubtract :: Canon -> Canon -> CycloMap -> (Canon, CycloMap)+cAdd = cApplyAdtvOp True +cSubtract = cApplyAdtvOp False ++-- | Internal Function to compute sum or difference based on first param. Much heavy lifting under the hood here.+cApplyAdtvOp :: Bool -> Canon -> Canon -> CycloMap -> (Canon, CycloMap)+cApplyAdtvOp _ x Pc0 m = (x, m)+cApplyAdtvOp True Pc0 y m = (y, m) -- True -> (+)+cApplyAdtvOp False Pc0 y m = (negate y, m) -- False -> (-) +cApplyAdtvOp b x y m = (gcd' * r, m')+ where gcd' = cGCD x y + x' = x / gcd'+ y' = y / gcd'+ r = crToC c False+ (c, m') = crApplyAdtvOptConv b (cToCR x') (cToCR y') m -- costly bit++-- | Exponentiation: This does allow for negative exponentiation if the Bool flag is True.+cExp :: Canon -> Canon -> Bool -> CycloMap -> (Canon, CycloMap)+cExp c e b m | cNegative e && (not b) + = error "Per param flag, negative exponentiation is not allowed here."+ | cIrrational c && cIrrational e + = error "cExp: Raising an irrational number to an irrational power is not currently supported."+ | otherwise = cExp' c e m+ where cExp' Pc0 e' m' | cPositive e' = (c0, m')+ | otherwise = error "0^e where e <= 0 is either undefined or illegal"+ cExp' Pc1 _ m' = (c1, m')+ cExp' _ Pc0 m' = (c1, m')+ cExp' c' e' m' = (gcrToC g, mg)+ where (g, mg) = gE (cToGCR c') e' m' + gE :: GCR_ -> Canon -> CycloMap -> (GCR_, CycloMap)+ gE g' e' m' | gcrNegative g' + = case cValueType e' of -- gcr exponentiation+ IntegralC -> if cOdd e' then (gcreN1:absTail, m'')+ else (absTail, m'')+ NonIntRationalC -> if cOdd d then (gcreN1:absTail, m'')+ else error "gE: Imaginary numbers not supported"+ IrrationalC -> error "gE: Raising neg numbers to irr. powers not supported" + | otherwise + = f g' m' -- equivalent to multiplying each exp by e' (with CycloMap threaded)+ where (absTail, m'') = gE (gcrAbs g') e' m'+ (_, d) = cSplit e' -- even denominator means you will have an imag. number+ f [] mf = ([], mf) + f ((p,x):gs) mf = (fp, mf')+ where (prd, mx) = cMult e' x mf+ (t, mn) = f gs mx+ (fp, mf') = gcrMult [(p, prd)] t mn++-- | Functions to check if a canon is negative/positive+cNegative, cPositive :: Canon -> Bool++cNegative (Bare n _ ) = n < 0+cNegative (Canonical c _ ) = gcrNegative c++cPositive (Bare n _ ) = n > 0+cPositive (Canonical c _ ) = gcrPositive c++-- | Functions for negation, absolute value and signum+cNegate, cAbs, cSignum :: Canon -> Canon ++cNegate (Bare 0 _) = c0+cNegate (Bare 1 _) = cN1+cNegate (Bare x Simplified) = Canonical (gcreN1 : [(x, c1)]) IntegralC -- prepend a "-1", not ideal+cNegate (Bare x NotSimplified) = Bare (-1 * x) NotSimplified +cNegate (Canonical x v) = gcrNegateCanonical x v+ +cAbs x | cNegative x = cNegate x+ | otherwise = x++cSignum (Bare 0 _) = c0+cSignum g | cNegative g = cN1+ | otherwise = c1++-- This internal function works for either gcrGCD or gcrLCM.+cLGApply :: (GCR_ -> GCR_ -> GCR_) -> Canon -> Canon -> Canon+cLGApply _ Pc0 y = y+cLGApply _ x Pc0 = x+cLGApply f x y | cNegative x || + cNegative y = gcrToC $ f (cToGCR $ cAbs x) (cToGCR $ cAbs y)+ | otherwise = gcrToC $ f (cToGCR x) (cToGCR y)++-- | Div function : Multiply by the reciprocal.+cDiv :: Canon -> Canon -> CycloMap -> (Canon, CycloMap)+cDiv _ Pc0 _ = error "cDiv: Division by zero error"+cDiv x y m = cMult (cReciprocal y) x m -- multiply by the reciprocal++-- | Compute reciprocal (by negating exponents).+cReciprocal :: Canon -> Canon+cReciprocal x = fst $ cExp x cN1 True crCycloInitMap -- raise number to (-1)st power++-- | Functions to check if a Canon is integral, (ir)rational, "simplified" or a prime tower+cIntegral, cIrrational, cRational, cSimplified, cIsPrimeTower :: Canon -> Bool++cIntegral (Bare _ _ ) = True+cIntegral (Canonical _ v ) = v == IntegralC++cIrrational (Canonical _ IrrationalC ) = True+cIrrational _ = False++cRational c = not $ cIrrational c++cSimplified (Bare _ Simplified) = True+cSimplified (Bare _ NotSimplified) = True+cSimplified (Canonical c _ ) = gcrSimplified c++cIsPrimeTower c = cPrimeTowerLevel c > 0 -- x^x would not be, but x^x^x would be++-- | cNumerator and cDenominator are for processing "rational" canon reps.+cNumerator, cDenominator :: Canon -> Canon++cNumerator (Canonical c _ ) = gcrToC $ filter (\x -> cPositive $ snd x) c -- filter positive exponents+cNumerator b = b ++cDenominator (Canonical c _ ) = gcrToC $ map (\(p,e) -> (p, cN1 * e)) $ filter (\(_,e) -> cNegative e) c -- negate negative expnts+cDenominator _ = c1 ++-- | Determines the depth/height of maximum prime tower in the Canon.+cDepth :: Canon-> Integer+cDepth (Bare _ _ ) = 1+cDepth (Canonical c _ ) = 1 + gcrDepth c++-- | Force the expression to be simplified. This could potentially be very expensive.+cSimplify :: Canon -> Canon+cSimplify (Bare n NotSimplified) = makeCanonFull n+cSimplify (Canonical c _ ) = gcrToC $ gcrSimplify c+cSimplify g = g -- Bare number already simplified : Fix when expr come into play++-- | Compute the rth-root of a Canon.+cRoot :: Canon -> Canon -> Canon +cRoot c r | cNegative r = error "r-th roots are not allowed when r <= 0" + | otherwise = gcrToC $ gcrRootI (cToGCR c) (cToGCR r)++-- | This is used for tetration, etc. It defaults to zero for non-integral reps.+cPrimeTowerLevel :: Canon -> Integer +cPrimeTowerLevel (Bare _ Simplified) = 1+cPrimeTowerLevel (Canonical g IntegralC) = case gcrPrimePower g of+ False -> 0+ True -> cPrimeTowerLevelI (snd $ head g) (fst $ head g) (1 :: Integer)+cPrimeTowerLevel _ = 0++-- | Internal workhorse function+cPrimeTowerLevelI :: Canon -> Integer -> Integer -> Integer+cPrimeTowerLevelI (Bare b _ ) n l | b == n = l + 1 + | otherwise = 0+cPrimeTowerLevelI (Canonical g IntegralC) n l | gcrPrimePower g == False = 0 + | n /= (fst $ head g) = 0+ | otherwise = cPrimeTowerLevelI (snd $ head g) n (l+1)+cPrimeTowerLevelI _ _ _ = 0++-- | Functions to Convert Canon to Generalized Canon Rep+canonToGCR, cToGCR :: Canon -> GCR_+canonToGCR (Canonical x _) = x+canonToGCR (Bare x NotSimplified) = canonToGCR $ makeCanon x -- ToDo: Thread in CycloMap?+canonToGCR (Bare x Simplified) | x == 1 = gcr1 + | otherwise = [(x, c1)]+cToGCR = canonToGCR++-- Warning: Don't call this for 0 or +/- 1. The value type will not change by negating the value +gcrNegateCanonical :: GCR_ -> CanonValueType -> Canon +gcrNegateCanonical g v | gcrNegative g = case gcrPrime (tail g) of+ True -> Bare (fst $ head $ tail g) Simplified+ False -> Canonical (tail g) v + | otherwise = Canonical (gcreN1 : g) v -- just prepend++gcrNegate :: GCR_ -> GCR_+gcrNegate Pg0 = gcr0+gcrNegate x | gcrNegative x = tail x + | otherwise = gcreN1 : x ++gcrNegative :: GCR_ -> Bool+gcrNegative PgNeg = True+gcrNegative _ = False++gcrPositive :: GCR_ -> Bool+gcrPositive PNeg = False+gcrPositive PZero = False+gcrPositive _ = True++gcrMult :: GCR_ -> GCR_ -> CycloMap -> (GCR_, CycloMap)+gcrMult x POne m = (x, m)+gcrMult POne y m = (y, m)+gcrMult _ Pg0 m = (gcr0, m)+gcrMult Pg0 _ m = (gcr0, m)+gcrMult x@(xh@(xp,xe):xs) y@(yh@(yp,ye):ys) m = case compare xp yp of+ LT -> (xh:g, m') + where (g, m') = gcrMult xs y m+ EQ -> if gcrNegative x || expSum == c0 + then gcrMult xs ys m -- cancel double negs/exponents adding to zero+ else ((xp, expSum):gf, mf) + where (expSum, m') = cAdd xe ye m + (gf, mf) = gcrMult xs ys m'+ GT -> (yh:g, m') + where (g, m') = gcrMult x ys m+gcrMult x y _ = error e + where e = "Non-exhaustive pattern error in gcrMult. Params: " ++ (show x) ++ "*" ++ (show y)++gcr1, gcr0 :: GCR_+gcr1 = []+gcr0 = [(0, c1)] ++gcreN1 :: GCRE_+gcreN1 = (-1, c1)++gcrToI :: GCR_ -> Integer+gcrToI g = product $ map f g+ where f (p, e) | ce > 0 = p ^ ce + | otherwise = error negExpErr+ where ce = cToI e + negExpErr = "gcrToI: Negative exponent found trying to convert " ++ (show g)++gcrToD :: GCR_ -> Double+gcrToD g = product $ map (\(p,e) -> (fromIntegral p) ** cToD e) g+ +gcrCmp :: GCR_ -> GCR_ -> Ordering+gcrCmp POne y = gcrCmpTo1 y True+gcrCmp x POne = gcrCmpTo1 x False+gcrCmp x y | x == y = EQ + | xN && yN = compare (gcrToC $ tail y) (gcrToC $ tail x)+ | xN = LT+ | yN = GT + | gcrIsZero x = LT+ | gcrIsZero y = GT+ | otherwise = case compare (gcrLogDouble x) (gcrLogDouble y) of+ -- If equal: we have to break out the big guns, both evaluate to infinity+ EQ -> compare (gcrLog'' x) (gcrLog'' y) + cmp -> cmp++ where xN = gcrNegative x+ yN = gcrNegative y ++ -- This is much more expensive but accurate. You have an "infinity" result issue potentially with gcrLogDouble+ gcrLog'' g = sum $ map f g+ f (p,e) = (toRational $ logDouble $ fromIntegral p) * (toRational e)+ logDouble :: Double -> Double+ logDouble n = log n+ +gcrCmpTo1 :: GCR_ -> Bool -> Ordering+gcrCmpTo1 PNeg b = if b then GT else LT+gcrCmpTo1 Pg0 b = if b then GT else LT+gcrCmpTo1 _ b = if b then LT else GT ++gcrLog :: GCR_ -> Rational+gcrLog g = crLog $ gcrToCR g ++-- | These internal functions should not be called directly. The definition of GCD and LCM are extended to handle non-integers.+gcrGCD, gcrLCM :: GCR_ -> GCR_ -> GCR_+gcrGCD POne _ = gcr1+gcrGCD _ POne = gcr1+gcrGCD x y = case compare xp yp of+ LT -> gcrGCD xs y+ EQ -> (xp, min xe ye) : gcrGCD xs ys + GT -> gcrGCD x ys+ where ((xp,xe):xs) = x+ ((yp,ye):ys) = y +gcrLCM POne y = y+gcrLCM x POne = x +gcrLCM x y = case compare xp yp of+ LT -> xh : gcrLCM xs y+ EQ -> (xp, max xe ye) : gcrLCM xs ys+ GT -> yh : gcrLCM x ys+ where (xh@(xp,xe) : xs) = x+ (yh@(yp,ye) : ys) = y ++gcrLogDouble :: GCR_ -> Double+gcrLogDouble g = sum $ map (\(p,e) -> (log $ fromIntegral p) * (cToD e)) g++divisionError :: String+divisionError = "gcrDiv: As requested per param, the dividend must be a multiple of the divisor." ++divByZeroError :: String+divByZeroError = "gcrDiv: Division by zero error!"++zeroDivZeroError :: String+zeroDivZeroError = "gcrDiv: Zero divided by zero is undefined!"++gcrDivStrict :: GCR_ -> GCR_ -> GCR_+gcrDivStrict x y = case (gcrDiv x y) of+ Left errorMsg -> error errorMsg+ Right results -> results++gcrDiv :: GCR_ -> GCR_ -> Either String GCR_+gcrDiv Pg0 Pg0 = Left zeroDivZeroError +gcrDiv Pg0 _ = Right gcr0+gcrDiv _ Pg0 = Left divByZeroError+gcrDiv n d = g' n d + where g' x POne = Right x+ g' POne _ = Left divisionError+ g' x y + | gcrNegative y = g' (gcrNegate x) (gcrAbs y)+ | otherwise = case compare xp yp of + LT -> case (g' xs y) of+ Left _ -> Left divisionError+ Right res -> Right ((xp, xe) : res)+ EQ | xe > ye -> case (g' xs ys) of+ Left _ -> Left divisionError+ Right res -> Right ((xp, xe - ye) : res)+ EQ | xe == ye -> gcrDiv xs ys+ _ -> Left divisionError + where ((xp,xe):xs) = x+ ((yp,ye):ys) = y ++-- GCR functions (GCR is an acronym for generalized canonical representation)+gcrAbs :: GCR_ -> GCR_+gcrAbs x | gcrNegative x = tail x+ | otherwise = x++gcrToCR :: GCR_ -> CR_+gcrToCR c = map (\(p,e) -> (p, cToI e)) c++gcrBare :: GCR_ -> Bool+gcrBare PBare = True+gcrBare POne = True+gcrBare _ = False++gcrPrime :: GCR_ -> Bool+gcrPrime PgPrime = True+gcrPrime _ = False ++gcrPrimePower :: GCR_ -> Bool+gcrPrimePower PgPPower = True+gcrPrimePower _ = False ++gcrIsZero :: GCR_ -> Bool+gcrIsZero Pg0 = True;+gcrIsZero _ = False ++gcrOdd, gcrEven :: GCR_ -> Bool+gcrOdd Pg0 = False+gcrOdd POne = True+gcrOdd c | gcrNegative c = gcrOdd (gcrAbs c)+ | otherwise = cp /= 2 + where (cp,_):_ = c++gcrEven g = not (gcrOdd g)++gcrEqCheck :: GCR_ -> GCR_ -> Bool+gcrEqCheck POne POne = True+gcrEqCheck POne _ = False+gcrEqCheck _ POne = False +gcrEqCheck ((xp,xe):xs) ((yp,ye):ys) | xp /= yp || xe /= ye = False + | otherwise = gcrEqCheck xs ys+gcrEqCheck x y = error e+ where e = "Non-exhaustive patterns in gcrEqCheck comparing " ++ (show x) ++ " to " ++ (show y)++gcrDepth :: GCR_ -> Integer+gcrDepth g = maximum $ map (\(_,e) -> cDepth e) g++gcrSimplified :: GCR_ -> Bool+gcrSimplified g = all (\(_,e) -> cSimplified e) g ++gcrSimplify :: GCR_ -> GCR_+gcrSimplify g = map (\(p,e) -> (p, cSimplify e)) g++gcrRootI :: GCR_ -> GCR_ -> GCR_+gcrRootI POne _ = gcr1 +gcrRootI c r | not $ gcrNegative c = case gcrDiv (cToGCR ce) r of+ Left _ -> error e + Right quotient -> (cp, gcrToC quotient) : gcrRootI cs r+ | gcrEven r = error "Imaginary numbers not allowed: Even root of negative number requested."+ | otherwise = gcreN1 : gcrRootI (gcrAbs c) r+ where ((cp,ce):cs) = c + e = "gcrRootI: All expnts must be multiples of " ++ (show r) ++ ". Not so with " ++ (show c)++-- | Check if the number is simplified rather than factoring it. Simplified is equivalent to having one term in the list.+getBareStatus :: Integer -> BareStatus+getBareStatus n | n < -1 = NotSimplified + | n <= 1 || isPrime n = Simplified + | otherwise = NotSimplified++-- | Instance of CanonConv class +instance CanonConv Canon where+ toSC c = toSC $ cToCR c+ toRC c = toRC $ cToCR c+ +-- | Canon of form x^1. (Does not match on 1 itself)+pattern PBare :: forall t. [(t, Canon)]+pattern PBare <- [(_, Bare 1 _)] ++-- | Canon of form p^e where e >= 1. p has already been verified to be prime.+pattern PgPPower :: forall t a. (Num a, Ord a) => [(a, t)]+pattern PgPPower <- [(compare 1 -> LT, _ )]++-- | Canon of form p^1 where p is prime+pattern PgPrime :: forall a. (Num a, Ord a) => [(a, Canon)]+pattern PgPrime <- [(compare 1 -> LT, Bare 1 _)] ++-- | Pattern looks for Canons beginning with negative 1+pattern PgNeg :: forall a. (Num a, Eq a) => [(a, Canon)]+pattern PgNeg <- ((-1, Bare 1 _):_) ++-- | Pattern for "generalized" zero+pattern Pg0 :: forall a. (Num a, Eq a) => [(a, Canon)]+pattern Pg0 <- [(0, Bare 1 _)] -- internal pattern for zero++-- | Patterns for 0 and 1+pattern Pc0 :: Canon+pattern Pc0 <- Bare 0 _++pattern Pc1 :: Canon+pattern Pc1 <- Bare 1 _ ++pattern Pc2 :: Canon+pattern Pc2 <- Bare 2 _++-- ToDo: Fix this Mod function. "Proper" rewrite has terrible performance+{-+pattern PcN1 :: Canon -- this pattern is only used in the "bad" function+pattern PcN1 <- Canonical [(-1, Bare 1 _)] _++cModBAD :: Canon -> Canon -> CycloMap -> (Canon, CycloMap)+cModBAD c m cm | cIntegral c && cIntegral m = f c m cm+ | otherwise = error "cModBAD: Must both parameters must be integral"+ where f _ Pc0 _ = error "cModBAD: Divide by zero error when computing n mod 0"+ f _ Pc1 cm' = (0, cm')+ f _ PcN1 cm' = (0, cm')+ f Pc0 _ cm' = (0, cm')+ f c' m' cm' | m' == c0 = error "cModBAD: Divide by zero error when computing n mod 0"+ | ma == c1 = (c0, cm')+ | ca == ma = (c0, cm')+ | cn && mn = (cNegate mrn, cmn) -- both (n)egative+ | (not cn) && (not mn) &&+ ca < ma = (ca, cm')+ | (cn && not mn) ||+ (mn && not cn) = ((cSignum m') * (makeC $ maI - mrm), cmm) -- (m)ixed sign: TODO: CycloMap threading+ | otherwise = (makeC io, mo)+ where (cn, mn) = (cNegative c', cNegative m')+ (ca, ma) = (cAbs c', cAbs m')+ (mrn, cmn) = f ca ma cm'+ (mrm, cmm) = f' ca maI cm'+ maI = cToI ma+ (io, mo) = f' c' (cToI m') cm'+ f' (Bare n _ ) mI cm' = (mod n mI, cm')+ f' ic@(Canonical c' IntegralC) mI cm' | cNegative ic = error "The canon must be positive here"+ | otherwise = (mod ip mI, cmf)+ where (ip, cmf) = i c' cm'' (1 :: Integer) -- performs fold-like product+ i [] cmi pri = (pri, cmi) -- with CycloMap threading+ i (l:ls) cmi pri | pri == 0 = (pri, cmi)+ | otherwise = i ls cmv (pri * v)+ where (v, cmv) = pf l cmi+ pf (p,e) mp = (pmI p (cToI v) mI, cmv)+ where (v, cmv) = f e tm mp+ (tm, cm'') = cTotient ic cm'+ f' (Canonical _ _ ) _ _ = error "cModBAD: Logic error: Canonical var has to be integral at this point"+-}+
+ Math/NumberTheory/Canon/Additive.hs view
@@ -0,0 +1,109 @@+-- |+-- Module: Math.NumberTheory.Canon.Additive+-- Copyright: (c) 2015-2018 Frederick Schneider+-- Licence: MIT+-- Maintainer: Frederick Schneider <frederick.schneider2011@gmail.com>+-- Stability: Provisional+--+-- Mostly functions for the addition and subtraction of CRs (Canonical Representations of numbers)++module Math.NumberTheory.Canon.Additive (+ crAdd,+ crSubtract,+ crAddR,+ crSubtractR, + crApplyAdtvOpt,+ crApplyAdtvOptConv, + crApplyAdtvOptR,+ crQuotRem +)+where++import Math.NumberTheory.Canon.Internals+import Math.NumberTheory.Canon.AurifCyclo (crCycloAurifApply, CycloMap)++-- | Functions for computing sums and differences. +crAdd, crSubtract, crAddR, crSubtractR :: CR_ -> CR_ -> CycloMap -> (CR_, CycloMap)+crAdd = crApplyAdtvOpt True +crSubtract = crApplyAdtvOpt False+crAddR = crApplyAdtvOptR True+crSubtractR = crApplyAdtvOptR False++-- | crApplyAdtvOptR performs addition/subtraction on two rational canons. +{-+ Like the nonR version, we take the GCD to try to simplify the expression we need to + convert to an integer and back. Here's a breakdown of the steps ...+ + nx ny nx*dy op ny*dx nf1 op nf2+ x op y => -- op -- => -------------- = ---------- => + dx dy dx * dy dx * dy++ ngcd * (nf1r op nf2r) ngcd * nf n+ -------------------- => --------- => -+ dx * dy dx * dy d++-} +crApplyAdtvOptR :: Bool -> CR_ -> CR_ -> CycloMap -> (CR_, CycloMap)+crApplyAdtvOptR _ x PZero m = (x, m)+crApplyAdtvOptR True PZero y m = (y, m) -- True -> (+)+crApplyAdtvOptR False PZero y m = (crNegate y, m) -- False -> (-) +crApplyAdtvOptR b x y m = (crDivRational n d, m')+ where (nx, dx) = crSplit x+ (ny, dy) = crSplit y + nf1 = crMult nx dy+ nf2 = crMult ny dx + ngcd = crGCD nf1 nf2+ nf1r = crDivStrict nf1 ngcd+ nf2r = crDivStrict nf2 ngcd + (nf, m') = crApplyAdtvOpt b nf1r nf2r m -- costly bit+ n = crMult ngcd nf+ d = crMult dx dy ++-- | crApplyAdtvOpt: Simplifies/Factorizes expressions x +/- y.+crApplyAdtvOpt :: Bool -> CR_ -> CR_ -> CycloMap -> (CR_, CycloMap)+crApplyAdtvOpt _ x PZero m = (x, m)+crApplyAdtvOpt True PZero y m = (y, m) -- True -> (+)+crApplyAdtvOpt False PZero y m = (crNegate y, m) -- False -> (-) +crApplyAdtvOpt b x y m = (crMult gcd' r, m')+ where gcd' = crGCD x y + xres = crDivStrict x gcd'+ yres = crDivStrict y gcd'+ (r, m') = crApplyAdtvOptConv b xres yres m -- costly bit + +logThreshold :: Double+logThreshold = 10 * (log 10) -- 'n' digit number++-- | crApplyAdtvOptConv is setup to convert different cases in a standard manner. All 8 combinations of signs and operators are covered here.+{-+ p1 + p2 => p1 + p2, p1 - p2 => p1 - p2+ p1 + n2 => p1 - p2, p1 - n2 => p1 + p2++ n1 + n2 => -(p1 + p2), n1 - n2 => (p2 - p1) + n1 + p2 => (p2 - p1), n1 - p2 => -(p2 + p1) +-} +crApplyAdtvOptConv :: Bool -> CR_ -> CR_ -> CycloMap -> (CR_, CycloMap)+crApplyAdtvOptConv b x y m + | gi < 2 || mL <= logThreshold + = (crSimpleApply op x y, m) -- no algebraic optimization we can perform+ | crPositive x = if crPositive y then crCycloAurifApply b ax ay g gi m+ else crCycloAurifApply (not b) ax ay g gi m+ | (crNegative y) && b = (crNegate c1, m1)+ | (crNegative y) && not b = crCycloAurifApply b ay ax g gi m+ | b = crCycloAurifApply (not b) ay ax g gi m+ | otherwise = (crNegate c2, m2) + where op = if b then (+) else (-)+ (ax, ay) = (crAbs x, crAbs y)+ gi = gcd (crMaxRoot ax) (crMaxRoot ay)+ g = crFromInteger $ fromIntegral gi+ mL = max (crLogDouble ax) (crLogDouble ay)+ (c1, m1) = crCycloAurifApply b ax ay g gi m+ (c2, m2) = crCycloAurifApply (not b) ax ay g gi m -- corresponds to "otherwise"++-- | Quot Rem function for Canon Rep. Optimization: Check first if q is a multiple of r. If so, we avoid the potentially expensive conversion.+crQuotRem :: CR_ -> CR_ -> CycloMap -> ((CR_, CR_), CycloMap)+crQuotRem x y m = case (crDiv x y) of+ Left _ -> ((q, md), m') + Right quotient -> ((quotient, cr0), m)+ where md = crMod x y -- Better to compute quotient this way .. to take adv. of alg. forms+ q = crDivStrict d y -- (x - x%y) / y. + (d, m') = crSubtract x md m
+ Math/NumberTheory/Canon/AurifCyclo.hs view
@@ -0,0 +1,490 @@+-- |+-- Module: Math.NumberTheory.Canon.AurifCyclo+-- Copyright: (c) 2015-2018 Frederick Schneider+-- Licence: MIT+-- Maintainer: Frederick Schneider <frederick.schneider2011@gmail.com>+-- Stability: Provisional+--+-- Aurifeullian and Cyclotomic factorization method functions.++{-# LANGUAGE PatternSynonyms, ViewPatterns #-}++module Math.NumberTheory.Canon.AurifCyclo (+ aurCandDec, aurCandDecCr, + aurDec, aurDecCr,+ applyCycloPair, applyCycloPairWithMap,+ cyclo, cycloWithMap,+ cycloDivSet, cycloDivSetWithMap,+ chineseAurif, chineseAurifWithMap, chineseAurifCr,++ crCycloAurifApply, applyCrCycloPair, divvy,+ CycloMap, fromCycloMap, fromCM, showCyclo, crCycloInitMap+)+where++import Math.NumberTheory.Canon.Internals+import Math.NumberTheory.Moduli.Jacobi (JacobiSymbol(..), jacobi)+import Data.Array (array, (!), Array(), elems) -- to do: convert to unboxed? https://wiki.haskell.org/Arrays+import GHC.Real (numerator, denominator)+import Math.Polynomial( Poly(), poly, multPoly, quotPoly, Endianness(..), polyCoeffs)+import Data.List (sort, sortBy, (\\))+import qualified Data.Map as M++-- CR_ Rep of 2+cr2 :: CR_+cr2 = crFromI 2++-- | This function checks if the inputs along with operator flag have a cyclotomic or Aurifeuillian form to greatly simplify factoring.+-- If they do not, potentially much more expesive simple factorization is used via crSimpleApply.+-- Note: The cyclotomic map is threaded into the functions+crCycloAurifApply :: Bool -> CR_ -> CR_ -> CR_ -> Integer -> CycloMap -> (CR_, CycloMap)+crCycloAurifApply b x y g gi m+ -- Optimization for prime g: If g is a prime (and exp not of from x^2 + y^2) but not Aurifeullian (verify + | (crPrime g) && not (g == cr2 && b) + = eA ([term1, termNM1], m) -- split into (x +/- y) and (x^(n-1) ... -/+ y^(n-1)) ++ -- Factorize: grtx^g - grty^g via cyclotomic polynomials + | not b = eA (cycA grtx grty g) ++ -- Factorize x^n + y^n using cyclotomic polynomials (if n = 2^x*m where x >= 0 and m > 2)+ | b && not gpwr2 = eA (cycA (oddRoot x) (-1 * oddRoot y) odd') ++ | otherwise = (crSimpleApply op x y, m)+ where op = if b then (+) else (-)+ ((gp, _):gs) = g+ gpwr2 = gp == 2 && gs == [] + gth_root v = crToI $ crRoot v gi+ grtx = gth_root x+ grty = gth_root y+ + -- used when factoring x^p +/- 1 where p is prime+ term1 = integerApply op (crRoot x gi) (crRoot y gi) -- a +/- b+ termNM1 = div (integerApply op x y) term1 -- divide a^g +/- b^g by the term above++ cycA x' y' n = (sort ia, m') -- sort the integers returned from low to high, should help if there are larger terms+ where (ia, m') = applyCrCycloPair x' y' n m+ eA (a,mp) = (foldr1 crMult $ map crFromI v, m') -- eA stands for "enriched apply"+ where (v, m') = case aurCandDecU x y gi g b of+ Nothing -> auL a mp -- can't do anything Brent Aurif-wise, try Chinese method+ Just (a1, a2) -> auL (divvy a a1 a2) mp -- meld in the 2 Aurif factors with input array+ auL al ma = case c of -- aL stands for "augmented list)+ Nothing -> (al, mp') -- just return what was passed in+ Just (a3, a4) -> (divvy al a3 a4, mp') -- additional "Chinese" factors+ where (c, mp') = chineseAurifCr x y b ma ++ odd' | gp == 2 = tail g -- grabs number sans any power of 2+ | otherwise = g+ oddRoot v = crToI $ crRoot v (crToI odd')+ +{-+The following functions implement Richard Brent's algorithm for computing Aurifeullian factors.+His logic is used in conjuction with cyclotomic polynomial decomposition.++http://maths-people.anu.edu.au/~brent/pd/rpb127.pdf+http://maths-people.anu.edu.au/~brent/pd/rpb135.pdf+-}++-- | Integer wrapper for aurCandDecCr+aurCandDec :: Integer -> Integer -> Bool -> Maybe (Integer, Integer)+aurCandDec xi yi b = aurCandDecCr (crFromI xi) (crFromI yi) b++-- | This function checks if the input is a candidate for Aurifeuillian decomposition.+-- If so, split it into two and evaluate it. Otherwise, return nothing. +-- The code will "prep" the input params so they will be relatively prime.+aurCandDecCr :: CR_ -> CR_ -> Bool -> Maybe (Integer, Integer)+aurCandDecCr xp yp b = aurCandDecU x y n (crFromI n) b + where n = gcd (crMaxRoot $ crAbs x) (crMaxRoot $ crAbs y)+ gxy = crGCD xp yp + (x, y) = (crDivStrict xp gxy, crDivStrict yp gxy) -- this will fix the input to be relatively prime++-- U belows means unsafe. Don't call this directly. The function assumes that x and y are relatively prime. Currently uses Brent logic only+aurCandDecU :: CR_ -> CR_ -> Integer -> CR_ -> Bool -> Maybe (Integer, Integer)+aurCandDecU x y n cr b| (nm4 == 1 && b) || (nm4 /= 1 && not b) ||+ (xdg == x && ydg == y) || (m /= 0)+ = Nothing -- + | otherwise = case aurDecI n' cr' of+ Nothing -> Nothing+ Just (gamma, delta) -> apply gamma delta+ where + -- override of n, to attempt decomp for g = gcd when number of form: g^gd +/-1, + -- this will only work when either x or y = 1 and not for any other divisor of g. + -- If both terms are not 1, we just attempt an Aurif. decomp for n+ + -- need to integrate chineseAurif, it does something different + (n', cr') | x /= cr1 && y /= cr1 = (n, cr) + | otherwise = (gcd1i, gcd1)+ where x1 = if y /= cr1 then y else x+ gcd1 = crRadical $ crGCD x1 cr + gcd1i = crToI gcd1 + nm4 = mod n' 4 + divTry a = case crDiv a (crExp cr' n' False) of -- check to divide by n^n, if not return original+ Left _ -> a + Right quotient -> quotient+ xdg = divTry x+ ydg = divTry y+ mrGCD = gcd (crMaxRoot $ crAbs xdg) (crMaxRoot $ crAbs ydg)+ m = mod mrGCD (2*n')++ -- need to consider cyclotomic translations, order the terms + (x', ml) | ydg /= y = ( (crDivRational ydg x), if (not b) then (-1) else 1)+ | otherwise = ( (crDivRational xdg y), 1); + + {- The more familar form of the below is (C(x)^2 - nxD(x)^2):+ gm(x)^2 -nx*dt(x)^2 => + gamma +/- sqrt(nx) * delta -}+ xrtn = crMult cr' (crRoot x' n')+ xrtnr = crToRational xrtn+ sqrtnxr = crToRational $ crRoot (crMult cr' xrtn) 2++ apply gm dt = Just (ml * numerator f1, numerator f2)+ where f1 = c - sqrtnxr * d+ f2 = c + sqrtnxr * d+ c = aA gm xrtnr+ d = aA dt xrtnr + -- aA means applyArray. array/lists are treated like polynomials (zero-base assumed) + aA a x'' = f (elems a) 1 0 + where f [] _ a' = a'+ f (c:cs) m' a' = f cs (m'*x'') (a' + (toRational c)*m') + +-- Example Aurif. decomp: C5(x) = x^2 + 3x + 1, D5(x) = x + 1 => Cyclotomic5(x) = C5(x)^2 − 5x*D5(x)^2 ++-- | This function returns a pair of polynomials (in array form) or Nothing (if it's squareful). +-- An illogical n (n <= 1) will generate an error.+aurDec :: Integer -> Maybe (Array Integer Integer, Array Integer Integer)+aurDec n | n <= 1 = error "aurifDecomp: n must be greater than 1"+ | otherwise = aurDecI n (crFromI n)++-- | CR_ wrapper for aurDec +aurDecCr :: CR_ -> Maybe (Array Integer Integer, Array Integer Integer) +aurDecCr cr = aurDec (crToI cr)++-- | Internal Aurifeullian Decomposition Workhorse Function+aurDecI :: Integer -> CR_ -> Maybe (Array Integer Integer, Array Integer Integer) +aurDecI n cr | crHasSquare cr || n < 2 || n' < 2+ = Nothing+ | otherwise = Just (gamma, delta)+ where nm4 = mod n 4 + n' = if (nm4 == 1) then n else (2*n)+ d = div (totient n') 2+ dm2 = mod d 2+ + -- max gamma and delta subscripts to explicitly compute (add'l terms come from symmetry)+ mg | dm2 == 1 = div (d-1) 2+ | otherwise = div d 2+ md | dm2 == 1 = div (d-1) 2+ | otherwise = div (d-2) 2+ + -- create q array of size td2: q(2k+1) = jacobi n (2k+1), q(2k+2 = mn * totient (2k+2)+ q = array (1, d) ([(i, f i) | i <- [1..d]]) + where f i | mod i 2 == 1 = convJacobi $ jacobi n i+ | otherwise = eQ+ where eQ = moeb (crFromI $ div n' g) * (totient g) * (cos' $ (n-1)*i) + g = gcd n' i ++ -- moebius fcn: 0 if has square, otherwise based on number of distinct prime factors+ moeb cr' | crHasSquare cr' = 0 + | mod (length cr') 2 == 1 = -1 + | otherwise = 1++ cos' c | m8 == 2 || m8 == 6 = 0 -- "cosine" function+ | m8 == 4 = -1+ | m8 == 0 = 1+ | otherwise = error "Logic error: bad/odd value passed to cos'" + where m8 = mod c 8++ -- These two arrarys have mutually recursive definitions+ gamma = array (0, d) ([(0,1)] ++ [(i, gf i) | i <- [1..d]]) + where gf k | k > mg = gamma!(d-k) + | otherwise = div gTerm (2*k)+ where gTerm = sum $ map f [0..k-1]+ where f j = n * q!(2*k-2*j-1) * delta!j - q!(2*k-2*j) * gamma!j++ delta = array (0, d-1) ([(0,1)] ++ [(i, df i) | i <- [1..d-1]]) + where df k | k > md = delta!(d-k-1) + | otherwise = div dTerm (2*k+1)+ where dTerm = gamma!k + sum (map f [0..k-1]) + where f j = q!(2*k-2*j+1) * gamma!j - q!(2*k-2*j) * delta!j++ {- Pseudocode for computing gammas G and deltas D + G(0) = 1+ D(0) = 1++ Evaluate G(k) for 1 .. floor(d/2)+ G(k) = (1/2k) * sum(n*q(2k-2j-1)*D(j) - q(2k-2j)G(j)) [for j= 0 to k-1)++ Evaluate D(k) for 1 .. floor(d-1/2)+ D(k) = (1/2k+1) * ( G(k) + sum(q(2k+1-2j)*D(j) - q(2k-2j)D(j)) )++ Evaluate G(k) for (floor(d/2)+1) to d => G(k) = G(d-k)+ Evaluate D(k) for (floor(d+1/2)) to d-1 => D(k) = D(d-k) + + Cyc(n) = C(x)^2 -nxD(x)^2 where gamma and delta are the coeffs for C(x) and D(x) respectively+ -} ++-- | Internal function requires two integers (computed via Aurif. methods) along with a list of Integers. The product of +-- the Integers must be a divisor of the list's product otherwise an error will be generated.+-- It's called divvy because it splits the 2 integers across the array using the gcd.+-- This will help factoring because the larger term(s) will be broken up into smaller pieces.+divvy :: [Integer] -> Integer -> Integer -> [Integer]+divvy a x y = d (sortBy rev a) (abs x) (abs y) + where rev a' b' = if (a' > b') then LT else GT + d [] x' y' | x' == 1 && y' == 1 = []+ | abs x' == 1 && abs y' == 1 = [x' * y']+ | otherwise = error "Empty list passed as first param but x' and y' weren't both 1"+ d (c:cs) x' y' | x' == 1 && y' == 1 = c:cs + | otherwise = v ++ d cs (div x' gnx) (div y' gny)+ where v = filter (>1) $ [div q gny, gnx, gny]+ gnx = gcd c x'+ q = div c gnx+ gny = gcd q y' +{- Test: Run divvy on this:+ let a = [4,7,8401,62324358089100907319521,682969,61,374857981681,547]+ let x = 50031545486420197+ let y = 50031544711579219 -}++{- Cyclotomic factorizations for numbers of the form x^n +/- y^n+ Example:+ x^15-y^15 = (x-y) (x^2+x y+y^2) (x^4+x^3 y+x^2 y^2+x y^3+y^4) (x^8-x^7*y +x^5*y^3 - x^4*y^4 +x^3*y^5 -x*y^7+y^8)++ C(15) is a form of the last term where y = 1+ It's possible in some cases to do an additional Aurifeullian factorization (of the last term). -}++-- | CycloPair: Pair of an Integer and its Corresponding Cyclotomic Polynomial+type CycloPair = (Integer, Poly Float)++-- | CycloMapInternal: Map internal to CycloMap newtype+type CycloMapInternal = M.Map CR_ CycloPair++-- | CycloMap is a newtype hiding the details of a map of CR_ to pairs of integers and corresponding cyclotomic polynomials.+newtype CycloMap = MakeCM CycloMapInternal deriving (Eq, Show)++-- | Unwrap the CycloMap newtype.+fromCM, fromCycloMap :: CycloMap -> CycloMapInternal+fromCM (MakeCM cm) = cm+fromCycloMap = fromCM++-- | This is an initial map with the cyclotomic polynomials for 1 and 2.+crCycloInitMap :: CycloMap+crCycloInitMap = MakeCM $ M.insert cr1 (1, poly LE ([-1.0, 1.0] :: [Float])) M.empty++-- Two internal functions for the map internals+cmLookup :: CR_ -> CycloMap -> Maybe CycloPair +cmLookup c m = M.lookup c (fromCM m)++cmInsert :: CR_ -> CycloPair -> CycloMap -> CycloMap+cmInsert c p m = MakeCM $ M.insert c p (fromCM m)++-- Computing the cyclotomic polynomials for the divisor set of a number:+-- Begin with with a map of 2 elements to the cylotomic polynomial for 1 and 2+-- Check the radical and return the map: crCycloRad+-- If the cr doesn't equal the radical, then open it up to the other factors and the square-free factors must be in the map+-- Identity: x^n -1 is the product of cyclotomic polynomial for each d where d | n. ++-- | Integer wrapper for crCyclo with default CycloMap parameter+cyclo :: Integer -> (CycloPair, CycloMap)+cyclo n = crCyclo (crFromI n) crCycloInitMap++-- | Integer wrapper for crCyclo +cycloWithMap :: Integer -> CycloMap -> (CycloPair, CycloMap)+cycloWithMap n m = crCyclo (crFromI n) m++-- | Integer wrapper for crCycloDivSet with default CycloMap parameter+cycloDivSet :: Integer -> CycloMap+cycloDivSet n = fst $ crCycloDivSet (crFromI n) crCycloInitMap++-- | Integer wrapper for crCycloDivSet+cycloDivSetWithMap :: Integer -> CycloMap -> (CycloMap, CycloMap)+cycloDivSetWithMap n m = crCycloDivSet (crFromI n) m++-- | Return pair of expon. multiplier and radical's polynomial along with working cyclotomic map.+crCyclo :: CR_ -> CycloMap -> (CycloPair, CycloMap)+crCyclo cr m | crPositive cr = ((crToI $ crDivStrict cr r, p), m')+ | otherwise = error "crCyclo: Positive integer needed"+ where r = crRadical cr+ ((_,p), m') = crCycloRad r m ++-- | Return a pair of cyclo maps just for the divisors and then a master map.+crCycloDivSet :: CR_ -> CycloMap -> (CycloMap, CycloMap)+crCycloDivSet cr m | crPositive cr = m2+ | otherwise = error "crCycloDivSet: Positive integer needed" + where (_,m2) = c+ crd = crDivisors cr+ c = case cmLookup cr m of + Nothing -> c' -- need to compute it+ Just p -> (p, (mf, m)) -- found it. add a filtered version. ToDo: optimize this+ where mf = MakeCM $ M.fromList $ filter (\(n,_) -> elem n crd) $ M.toList $ fromCM m++ {- Performance note: The filter version above tended to be somewhat faster than the lookup+ version below so I used that+ where mf = MakeCM $ M.fromList $ map l crd -- "lookup version"+ where l d = case cmLookup d m of+ Nothing -> error $ e d+ Just p -> (d, p) + e d = error "crCycloDivSet: Logic error: Divisor = '" ++ (show d) ++ "' not found!" + -} ++ c' | r == cr = (pr, (pm, pm))+ | otherwise = (cm, (mm, mm)) + where (cm, mm) = mfn sqFulDivs pm+ r = crRadical cr+ (pr, pm) = crCycloRad r m+ sqFulDivs = crd \\ crDivisors r -- "squareful" divisors+ mfn [] _ = error "Logic Error in mfn: Empty list is forbidden"+ mfn (n:ns) mp | ns == [] = cp+ | otherwise = mfn ns mp' + where cp@(_, mp') = crCycloAll n mp++-- | Compute the "radical" divisors first and then the non-square free entries.+crCycloRad :: CR_ -> CycloMap -> (CycloPair, CycloMap)+crCycloRad cr m = case cmLookup cr m of + Nothing -> c' -- need to compute it+ Just p -> (p, m) -- found it+ where c' | cs == [] = (cycpr, cmInsert cr cycpr m) + | otherwise = (cyc_n, cmInsert cr cyc_n mp)+ where (_ : cs) = cr+ -- for primes, because the tail of the cr is [] meaning only one prime factor+ r = fromInteger $ crToI $ crRadical cr + -- ToDo: Optimize cycpr to be quotient of (r^n -1)/(r-1) when r is a big prime + cycpr = (1, poly LE (replicate r 1.0)) --prime : ToDo: Optimize this to be quotient when a+ -- for composites+ -- Create polynomial of the form : x^n -1+ xNm1 = poly LE ( (-1.0:(replicate (r-1) 0.0) ++ [1.0]) :: [Float] )+ (cPrd, mp) = mf (init $ crDivisors cr)+ cyc_n = (1, quotPoly xNm1 cPrd)++ -- mf (Memo Fold) takes a list of divisors and returns the pair: (cyclotomic product, memoized map)+ mf (n:ns) = f ns p m' + where ((_,p), m') = crCycloRad n m + f (n':ns') p' mp = f ns' (multPoly p' p'') m'' -- cycloMap-threaded mult. fold+ where ((_,p''), m'') = crCycloRad n' mp+ f _ p' mp = (p', mp)+ mf [] = error "Logic error: Blank list can't be passed to mf aka crCycloMemoFold" ++-- Return a pair of Integer and its Cyclotomic Polynomial while efficiently building up a cyc. poly. map+crCycloAll :: CR_ -> CycloMap -> (CycloPair, CycloMap)+crCycloAll cr m | p == cr1 = case cmLookup cr m of + Nothing -> error "Logic error: Radical value not found for crCycloAll"+ Just cb -> (cb, m) -- found it + | otherwise = (crp, cmInsert cr crp md)+ where (p, d) = crPullSq+ ((i,y), md) = case cmLookup d m of + Nothing -> crCycloAll d m -- need to compute it+ Just c -> (c, m) -- found it + -- Optimization: (1, x + 1) => (4, x + 1). Note: Cyc2(x) = x + 1+ -- The first value is the exponential multiplier Cyc8(x) = C2(x^4) = (x^4) + 1+ crp = ((fst $ head p) * i, y) + crPullSq = f [] cr+ where f h [] = (cr1, h)+ f h (c@(cp,ce):cs) | ce > 1 = ([(cp, 1)], h ++ (cp, ce-1):cs) + | otherwise = f (h ++ [c]) cs++-- | These "apply cyclo" functions will use cyclotomic polynomial methods to factor x^e - b^e.+applyCrCycloPair :: Integer -> Integer -> CR_ -> CycloMap -> ([Integer], CycloMap)+applyCrCycloPair l r cr m = (applyCrCycloPairI l r cr (M.elems $ fromCM md), mn)+ where (md, mn) = crCycloDivSet cr m++applyCrCycloPairI :: Integer -> Integer -> CR_ -> [CycloPair] -> [Integer]+applyCrCycloPairI l r cr cds = map applyPoly cds+ where nd = crTotient cr+ pA v = a where a = array (0,nd) ([(0,1)] ++ [(i, v*a!(i-1)) | i <- [1..nd]]) -- array of powers+ lpa = pA l+ rpa = pA r+ applyPoly (m,p) = foldr1 (+) (map f $ zip fmtdCy [0..])+ where f (a, b) | a == 0 = 0+ | otherwise = a * lpa!(m*b) * rpa!(m*(maxExp - b))+ fmtdCy = map ceiling $ polyCoeffs LE p -- format poly from mult poly pair+ maxExp = toInteger $ length fmtdCy - 1++-- | Wraps applyCycloPairWithMap with default CycloMap argument.+applyCycloPair :: Integer -> Integer -> Integer -> [Integer]+applyCycloPair x y e = fst $ applyCycloPairWithMap x y e crCycloInitMap++-- | This will use cyclotomic polynomial methods to factor x^e - b^e.+applyCycloPairWithMap :: Integer -> Integer -> Integer -> CycloMap -> ([Integer], CycloMap) +applyCycloPairWithMap x y e m = applyCrCycloPair x y (crFromI e) m++-- | This will display the cyclotomic polynomials for a CR.+showCyclo :: CR_ -> CycloMap -> [Char]+showCyclo n m = p $ map (\x -> (ceiling x) :: Integer) $ polyCoeffs LE (snd $ fst $ crCyclo n m)+ where p (c:cs) = show c ++ (p' cs (1 :: Int)) -- "LE" endianness is assumed here+ p _ = []+ p' (c:cs) s | c == 0 = r+ | otherwise = (if c > 0 then " + " else " - ") ++ (if ac == 1 then "" else show ac) +++ "x" ++ (if s == 1 then "" else "^" ++ show s) ++ r+ where r = p' cs (s+1) + ac = abs c+ p' _ _ = []++-- | All the exponents must be even to return True.+crSquareFlag :: CR_ -> Bool+crSquareFlag = all (\(_, ce) -> mod ce 2 == 0) ++-- | Wrapper for chineseAurifWithMap with default CycloMap parameter+chineseAurif :: Integer -> Integer -> Bool -> Maybe (Integer, Integer)+chineseAurif x y b = fst $ chineseAurifWithMap x y b crCycloInitMap++-- | Integer wrapper for chineseAurifCr+chineseAurifWithMap :: Integer -> Integer -> Bool -> CycloMap -> (Maybe (Integer, Integer), CycloMap)+chineseAurifWithMap x y b m = chineseAurifCr (crFromI x) (crFromI y) b m++-- The source for this algorithm is the paper by Sun Qi, Ren Debin, Hong Shaofang, Yuan Pingzhi and Han Qing+-- http://www.jams.or.jp/scm/contents/Vol-2-3/2-3-16.pdf (The formula at 2.7 there is implemented below)+-- This will handle a subset of the cases that the main Aurif. routines handle++-- | This function reduces the two CR parameters by gcd before calling an internal function to find a "Chinese" Aurifeullian factorization.+chineseAurifCr :: CR_ -> CR_ -> Bool -> CycloMap -> (Maybe (Integer, Integer), CycloMap)+chineseAurifCr xp yp b m = case c of+ Nothing -> chineseAurifI mbyx n myx (crToI myx) b m' -- if first try fails, try the reverse+ r -> (r, m') + where (c, m') = chineseAurifI mbxy n mxy (crToI mxy) b m+ gcdxy = crGCD xp yp+ (x, y) = (crDivStrict xp gcdxy, crDivStrict yp gcdxy) -- strip out any commonality+ n = gcd (crMaxRoot $ crAbs x) (crMaxRoot $ crAbs y) + ncr = crFromI n + mbxy = crRoot (crDivRational x y) n+ mxy = crGCD (crNumer mbxy) ncr + mbyx = crRecip mbxy+ myx = crGCD (crNumer mbyx) ncr ++-- | Internal function to find factor of mb^n +/- 1 (mb would be M from paper, mb meaning m "big".+-- Solution forms: Any non-zero multiple of (q^2m * p) ^ (p * k) op (r^2n)^(p * k) where k is an odd, postive number, m, n > 0.+-- This will work if the op is "+" when mod p 4 = 3 OR when op is "-" for when mod p 4 = 1.+chineseAurifI :: CR_ -> Integer -> CR_ -> Integer -> Bool -> CycloMap -> (Maybe (Integer, Integer), CycloMap)+chineseAurifI mbcr n mcr m b mp | mod n 2 == 0 || mod m 2 == 0 || -- n and m must both be odd+ m < 3 || km /= 0 || -- m must be odd and > 1 and m | n+ (mm4 == 1 && b) || -- sign and modulus+ (mm4 == 3 && not b) || -- must be in synch + mbdm == cr0 || not (crSquareFlag mbdm) -- mb/m must be a square and integral+ = (Nothing, mp)+ | otherwise = case cv - (gd1 * gd2) of+ 0 -> (Just (gd1, gd2), mp')+ _ -> (Nothing, mp) -- addl check since paper doesn't indicate rationals are supported+ where mm4 = mod m 4+ e = toRational $ if (mm4 == 3) then (-1) else mm4+ (k, km) = quotRem n m+ mbdm = case crDiv mbcr mcr of+ Left _ -> cr0 -- error condition if not a multiple+ Right q -> q+ r = crToRational $ crRoot (crMult mbcr mcr) 2 -- sqrt (m*M) from paper+ mb = crToRational mbcr+ jR c = toRational $ convJacobi $ jacobi c m + eM = e * mb+ v1 = (toRational m) * mb^(div (k * (m + 1)) 2)+ v2 = t * s+ where t = (jR 2) * r * (mb ^ (div (k-1) 2))+ s = sum $ map (\c -> (jR c) * eM^(k*c)) + $ filter (\c -> gcd c m == 1) [1..m] -- rel. prime+ ncr = crFromI n+ -- get cyclotomic value+ cv = head $ applyCrCycloPairI (numerator eM) (denominator eM) ncr [cp]+ (cp, mp') = crCyclo ncr mp + gd1 = gcd cv (numerator $ v1 - v2) -- delta1 + gd2 = gcd cv (numerator $ v1 + v2) -- delta2++-- workaround after arithmoi changes+convJacobi :: JacobiSymbol -> Integer+convJacobi j = case j of+ MinusOne -> -1+ Zero -> 0+ One -> 1+
+ Math/NumberTheory/Canon/Internals.hs view
@@ -0,0 +1,590 @@+-- |+-- Module: Math.NumberTheory.Canon.Internals+-- Copyright: (c) 2015-2018 Frederick Schneider+-- Licence: MIT+-- Maintainer: Frederick Schneider <frederick.schneider2011@gmail.com>+-- Stability: Provisional+--+-- This module defines the internal canonical representation of numbers (CR_), a product of pairs (prime and exponent). +-- It's not meant to be called directly.++{-# LANGUAGE PatternSynonyms, ViewPatterns, ScopedTypeVariables, DataKinds, RankNTypes #-}++module Math.NumberTheory.Canon.Internals (+ CanonRep_, CR_,+ crValidIntegral, crValidIntegralViaUserFunc,+ crValidRational, crValidRationalViaUserFunc,+ crExp,+ crRoot, + crMaxRoot,+ crShow,+ ceShow,+ crFromInteger, crFromI,+ crToInteger, crToI,+ crCmp, + crMult,+ crNegate,+ crAbs,+ crDivStrict,+ crSignum,+ crNumer,+ crDenom,+ crSplit,+ crDivRational,+ crIntegral,+ crShowRational,+ crToRational,+ crGCD,+ crLCM,+ crNegative,+ crPositive,+ crLog,+ crLogDouble,+ crDiv,+ crRadical,+ integerApply,+ crSimpleApply,+ crPrime,+ crHasSquare,+ crRecip,+ crMin,+ crMax,+ crValid,+ crMod, crModI,+ + crDivisors,+ crNumDivisors,+ crWhichDivisor,+ crNthDivisor,+ crDivsPlus,+ crTau,+ crTotient,+ crPhi,+ + crN1,+ cr0,+ cr1,+ creN1,+ + pattern PZero,+ pattern PZeroBad,+ pattern POne,+ pattern PNeg,+ pattern PNotPos,+ pattern PN1,+ + pattern PIntNeg,+ pattern PIntNPos,+ pattern PIntPos,++ -- functions deprecated from arithmoi that needed to be included here+ totient,+ pmI -- stands for powerModInteger+) +where++{- Canon or canon rep is short for canonical representation.++ In this context, it refers to the prime factorization of a number.+ Example: 250 = 2 * 5^3. It would be represented internally here as: [(2,1),(5,3)]++ So, this library along with Canon.hs can be used as shorthand for extremely large numbers.+ Multiplicative actions are cheap. But, addition and subtraction is potentially very expensive + as any "raw" sums or differences are factorized. There are optimizations for sums/differences+ including those with special forms (algebraic, Aurifeuillean).++ Here are the possibilities:+ + Zero: [(0,1)]+ + One [] + + Other Positive + Numbers: A set of (p,e) pairs where the p's are prime and in ascending order. + For "integers", the e(xponents) must be positive+ For "rationals", the e(xponents) must not be zero+ All "integers" are "rationals" but not vice-versa.+ + Negative Numbers: (-1,1):P where P is a canon rep for any positive number+++ Note: Much of the code is spent on handling these special cases (Zero, One, Negatives). + + Each integer and rational will have a unique "canon rep".+ + Caveats: The behavior is undefined when directly working with "canon reps" which are not valid.+ The behavior of using rational CRs directly (when integral CRs are specified) is also undefined.+ + The Canon library should be used as it hides these internal details.+-}++import Data.List (intersperse)+import Math.NumberTheory.Primes.Factorisation (factorise, factorise')+import Data.List (sortBy)+import Math.NumberTheory.Primes.Testing (isPrime)+import GHC.Real (Ratio(..))++-- | Canon element: prime and exponent pair+type CanonElement_ = (Integer,Integer)++-- | Canonical representation: list of canon elements+type CanonRep_ = [CanonElement_]++-- | Shorthand for canonical representation+type CR_ = CanonRep_++-- | Pattern to match the CR_ equivalent of 1+pattern POne :: forall t. [t]+pattern POne = []++-- | Pattern to match the CR_ equivalent of zero+pattern PZero :: forall a a1. (Num a, Num a1, Eq a, Eq a1) => [(a1, a)]+pattern PZero = [(0,1)]++-- | Pattern to match the CR_ equivalent of -1+pattern PN1 :: forall a a1. (Num a, Num a1, Eq a, Eq a1) => [(a1, a)]+pattern PN1 = [(0,1)]++-- | Pattern to match a badly formed zero, meaning it's an invalid CR_+pattern PZeroBad :: forall t a. (Num a, Eq a) => [(a, t)]+pattern PZeroBad <- ((0,_):_) -- MUST check after PZero++-- No longer necessary+-- pattern PNotIntegral :: forall a t. (Num a, Ord a) => [(t, a)]+-- pattern PNotIntegral <- ( (_, compare 0 -> GT):_ ) -- negative exponent in the 2nd member of pair++-- | Pattern to match a non-positive CR_+pattern PNotPos :: forall t a. (Num a, Ord a) => [(a, t)]+pattern PNotPos <- ( (compare 1 -> GT, _):_ ) -- first term is 0, -1 and so not positive++-- | Pattern to match a negative number+pattern PIntNeg :: forall a. (Num a, Ord a) => a+pattern PIntNeg <- (compare 0 -> GT)++-- | Pattern to match a positive number+pattern PIntPos :: forall a. (Num a, Ord a) => a+pattern PIntPos <- (compare 0 -> LT)++-- | Pattern to match a non-positive number+pattern PIntNPos :: forall a. (Num a, Ord a) => a+pattern PIntNPos <- (compare 1 -> GT)++-- | Canonical values for a few special numbers+creN1, cre0 :: CanonElement_+creN1 = (-1,1)+cre0 = (0,1)++crN1, cr0, cr1 :: CanonRep_++-- | Canon rep for -1 +crN1 = [creN1]++-- | Canon rep for 0+cr0 = [cre0]++-- | Canon rep for 1 +cr1 = [] -- Yes, a canonical "1" is just an empty list.++-- | Pattern for a negative CR_+pattern PNeg :: forall a a1. (Num a, Num a1, Eq a, Eq a1) => [(a1, a)]+pattern PNeg <- ((-1,1):_) ++crNegative, crPositive :: CR_ -> Bool++-- | Check if a CR_ is negative.+crNegative PNeg = True+crNegative _ = False++-- | Check if a CR_ is positive.+crPositive PZero = False+crPositive x = not $ crNegative x++-- | Canon rep validity check: +-- The 2nd param checks the validity of the base, the 3rd of the exponent.+-- The base pred should be some kind of prime or psuedo-prime test unless you knew for +-- certain the bases are prime. There are two choices for the exp pred: +-- positiveOnly (True) or nonZero (False) (which allows for "rationals"). +crValid :: CR_ -> (Integer -> Bool) -> Bool -> Bool+crValid POne _ _ = True+crValid PZero _ _ = True+crValid PZeroBad _ _ = False+crValid c bp ef + | crNegative c = f (tail c) 1+ | otherwise = f c 1+ where f POne _ = True+ f ((cp,ce):cs) n | cp <= n || not (expPred ef ce) || not (bp cp) = False+ | otherwise = f cs cp+ f _ _ = error "Logic error in crValid'. Issue with pattern matching?"+ expPred b e = if b then (e > 0) else (e /= 0)++crValidIntegral, crValidRational :: CR_ -> Bool+crValidIntegralViaUserFunc, crValidRationalViaUserFunc :: CR_ -> (Integer -> Bool) -> Bool++-- | Checks if a CR_ represents an integral number.+crValidIntegral n = crValid n isPrime True++-- | Checks if a CR_ is Integral and valid per user-supplied criterion.+crValidIntegralViaUserFunc n f = crValid n f True++-- | Checks if a CR_ is represents a rational number (inclusive of integral numbers).+crValidRational n = crValid n isPrime False++-- | Checks if a CR_ is Rational and valid per user-supplied criterion.+crValidRationalViaUserFunc n f = crValid n f False++crFromInteger, crFromI :: Integer -> CR_++-- | Factor the number to convert it to a canonical rep. This is of course can be extremely expensive.+crFromInteger 0 = cr0+crFromInteger n = map (\(p, e) -> (p, toInteger e)) $ sortBy sf $ factorise n+ -- the prime factors must be in ascending order+ where sf (p1, _) (p2, _) | p1 < p2 = LT+ | otherwise = GT++-- | Shorthand for crFromInteger function+crFromI n = crFromInteger n ++crToInteger, crToI :: CR_ -> Integer++-- | Converts a canon rep back to an Integer.+crToInteger POne = 1+crToInteger PZero = 0+crToInteger c | (head c) == creN1 = -1 * (crToInteger $ tail c) -- negative number+ | otherwise = product $ map (\(x,y) -> x ^ y) c++-- | Alias to crToInteger.+crToI = crToInteger++-- | Compute the modulus between a CR_ and Integer and return an Integer.+crModI :: CR_ -> Integer -> Integer+crModI _ 0 = error "Divide by zero error when computing n mod 0"+crModI _ 1 = 0+crModI _ (-1) = 0+crModI POne PIntPos = 1+crModI PZero _ = 0+crModI c m | cn && mn = -1 * crModI (crAbs c) am+ | (cn && not mn) ||+ (mn && not cn) = (signum m) * (am - f (crAbs c) am)+ | otherwise = f c m+ where cn = crNegative c+ mn = m < 0+ am = abs m+ f c' m' = mod (product $ map (\(x,y) -> pmI x (mmt y) m') c') m'+ mmt e | e >= 1 = mod e $ totient m -- optimization+ | otherwise = error "Negative exponents are not allowed in crModI" ++-- | Compute modulus with all CR_ parameters. This wraps crModI.+crMod :: CR_ -> CR_ -> CR_+crMod c m = crFromI $ crModI c (crToI m)+ +-- | Display a Canon Element (as either p^e or p).+ceShow :: CanonElement_ -> String+ceShow (p,e) = show p ++ if e == 1 then "" + else "^" ++ (if e < 0 then "(" ++ se ++ ")" else se)+ where se = show e++crShow, crShowRational :: CR_ -> String++-- | Display a canonical representation.+crShow POne = show (1 :: Integer)+crShow x | null (tail x) = ceShow $ head x+ | otherwise = concat $ intersperse " * " $ map ceShow x ++-- | Display a Canonical Rep rationally, as a quotient of its numerator and denominator.+crShowRational c | d == cr1 = crShow n+ | otherwise = crShow n ++ "\n/\n" ++ crShow d+ where (n, d) = crSplit c ++crNegate, crAbs, crSignum :: CR_ -> CR_++-- | Negate a CR_.+crNegate PZero = cr0+crNegate x | crNegative x = tail x + | otherwise = creN1 : x ++-- | Take the Absolute Value of a CR_.+crAbs x | crNegative x = tail x+ | otherwise = x++-- | Compute the signum and return as CR_.+crSignum PZero = cr0;+crSignum x | crNegative x = crN1+ | otherwise = cr1++-- | CR_ Compare Function +crCmp :: CR_ -> CR_ -> Ordering+crCmp POne y = crCmp1 y True+crCmp x POne = crCmp1 x False+crCmp x y | x == y = EQ + | xN && yN = crCmp (tail y) (tail x)+ | xN = LT+ | yN = GT + | eqZero x = LT+ | eqZero y = GT+ | otherwise = case compare (crLogDouble x) (crLogDouble y) of+ EQ -> compare (crLog x) (crLog y) -- We have to break out the big guns, both evaluate to infinity+ cmp -> cmp+ where xN = crNegative x+ yN = crNegative y + eqZero PZero = True;+ eqZero _ = False++-- Internal: Compare when either term is 1.+crCmp1 :: CR_ -> Bool -> Ordering+crCmp1 PNeg b = if b then GT else LT+crCmp1 PZero b = if b then GT else LT+crCmp1 _ b = if b then LT else GT ++crMin, crMax :: CR_ -> CR_ -> CR_++-- | Min function+crMin x y = case crCmp x y of+ GT -> y+ _ -> x++-- | Max functon +crMax x y = case crCmp x y of+ LT -> y+ _ -> x + +divisionError, divByZeroError, zeroDivZeroError, negativeLogError :: String+divisionError = "For this function, the dividend must be a multiple of the divisor." +divByZeroError = "Division by zero error!"+zeroDivZeroError = "Zero divided by zero is undefined!"+negativeLogError = "The log of a negative number is undefined!"++-- | Strict division: Generates error if exact division is not possible.+crDivStrict :: CR_ -> CR_ -> CR_+crDivStrict x y = case crDiv x y of+ Left errorMsg -> error errorMsg+ Right quotient -> quotient++-- | Attempt to take the quotient.+crDiv :: CR_ -> CR_ -> Either String CR_+crDiv PZero PZero = Left zeroDivZeroError +crDiv PZero _ = Right cr0+crDiv _ PZero = Left divByZeroError+crDiv x' y' = f x' y'+ where -- call this after handling zeroes above, then division just occurs within here+ f x POne = Right x+ f POne _ = Left divisionError+ f x y | crNegative y = f (crNegate x) (crAbs y)+ | otherwise = case compare xp yp of + LT -> case f xs y of+ Left _ -> Left divisionError+ Right r -> Right ((xp, xe):r)+ EQ| (xe > ye) -> case f xs ys of+ Left _ -> Left divisionError+ Right r -> Right ((xp,xe-ye):r)+ EQ| (xe == ye) -> f xs ys+ _ -> Left divisionError+ where ((xp,xe):xs) = x+ ((yp,ye):ys) = y ++crMult, crDivRational, crGCD, crLCM :: CR_ -> CR_ -> CR_++-- | Multiply two crs by summing the exponents for each prime.+crMult PZero _ = cr0+crMult _ PZero = cr0+crMult POne y = y+crMult x POne = x+crMult x y = case compare xp yp of+ LT -> xh : crMult xs y+ -- cancel double negs or expnts adding to zero+ EQ -> if crNegative x || expSum == 0 then r+ else (xp, expSum) : r+ where r = crMult xs ys+ GT -> yh : crMult x ys+ where (xh@(xp,xe):xs) = x+ (yh@(yp,ye):ys) = y+ expSum = xe + ye++-- | Division of rationals is equivalent to multiplying with negated exponents.+crDivRational x y = crMult (crRecip y) x -- multiply by the reciprocal++-- | For the GCD (Greatest Common Divisor), take the lesser of two exponents for each prime encountered.+crGCD PZero y = y+crGCD x PZero = x+crGCD x y | crNegative x || crNegative y = f (crAbs x) (crAbs y)+ | otherwise = f x y+ where f POne _ = cr1+ f _ POne = cr1+ f x' y' = case compare xp yp of+ LT -> f xs y'+ EQ -> (xp, min xe ye) : f xs ys + GT -> f x' ys+ where ((xp,xe):xs) = x'+ ((yp,ye):ys) = y' ++-- | For the LCM (Least Common Multiple), take the max of two exponents for each prime encountered.+crLCM PZero y = y+crLCM x PZero = x+crLCM x y | crNegative x || crNegative y = f (crAbs x) (crAbs y)+ | otherwise = f x y+ where f POne y' = y'+ f x' POne = x'+ f x' y' = case compare xp yp of+ LT -> xh : f xs y'+ EQ -> (xp, max xe ye) : f xs ys+ GT -> yh : f x' ys+ where (xh@(xp,xe):xs) = x'+ (yh@(yp,ye):ys) = y' ++-- | Take the reciprocal by raising a CR to the -1 power (equivalent to multiplying exponents by -1).+crRecip :: CR_ -> CR_+crRecip x = crExp x (-1) True++rootError :: CR_ -> Integer -> String+rootError c r = "crRoot: All exponents must be multiples of " ++ (show r) ++ ". Not so with " ++ (show c)++-- | Attempt to compute a particular root of a CR_.+crRoot :: CR_ -> Integer -> CR_ +crRoot _ PIntNeg = error "r-th roots are not allowed when r <= 0" +crRoot POne _ = cr1 +crRoot c r+ | crNegative c = if mod r 2 == 1 then creN1 : crRoot (crAbs c) r + else error "Imaginary numbers not allowed: Even root of negative number requested"+ | otherwise = if mod ce r == 0 then (cp, div ce r) : crRoot cs r+ else error $ rootError c r+ where ((cp,ce):cs) = c ++-- | Takes the maximum root of the number. Generally, the abs value would be passed to the function.+crMaxRoot :: CR_ -> Integer+crMaxRoot c = foldr (\x -> flip gcd $ snd x) 0 c++-- | Exponentiation. Note: this does allow for negative exponentiation if bool flag is True.+crExp :: CR_ -> Integer -> Bool -> CR_+crExp _ PIntNeg False = error "Per param flag, negative exponentiation is not allowed here."+crExp PZero PIntNPos _ = error "0^e where e <= 0 is either undefined or illegal"+crExp PZero _ _ = cr0+crExp POne _ _ = cr1+crExp _ 0 _ = cr1+crExp c em _ = ce c + where ce c' | crNegative c' = if mod em 2 == 1 then creN1 : absTail+ else absTail+ | otherwise = map (\(p,e) -> (p, e * em)) c'+ where absTail = ce $ crAbs c'++-- | This log function is much more expensive but accurate. You have an "infinity" problem potentially with crLogDouble.+crLog :: CR_ -> Rational+crLog PNeg = error negativeLogError+crLog c = sum $ map (\(p,e) -> (toRational $ logDouble $ fromIntegral p) * (fromIntegral e)) c+ where logDouble :: Double -> Double+ logDouble n = log n++-- | Returns log of CR_ as a Double.+crLogDouble :: CR_ -> Double+crLogDouble PNeg = error negativeLogError+crLogDouble c = sum $ map (\(x,y) -> log (fromIntegral x) * fromIntegral y) c+ +crNumer, crDenom, crRadical :: CR_ -> CR_+ +-- | Compute numerator (by filtering on positive exponents).+crNumer c = filter (\(_,e) -> e > 0) c++-- | Compute denominator. (Grab the primes with negative exponents and then flip the sign of the exponents.)+crDenom c = map (\(p,e) -> (p, (-1) * e)) $ filter (\(_,e) -> e < 0) c++-- | Check if a CR_ represents an integer.+crIntegral :: CR_ -> Bool+crIntegral x = all (\(_,e) -> e > 0) x -- all exponents must be positive++-- | Split a CR_ into its Numerator and Denominator.+crSplit :: CR_ -> (CR_, CR_)+crSplit c = (crNumer c, crDenom c)++-- | Convert a CR_ to a Rational number.+crToRational :: CR_ -> Rational+crToRational c = (crToI $ crNumer c) :% (crToI $ crDenom c)++-- | Compute the Radical of a CR_ (http://en.wikipedia.org/wiki/Radical_of_an_integer).+-- Its the product of the unique primes in its factorization.+crRadical n = map (\(p,_) -> (p, 1)) n ++-- | The Op(eration) is intended to be plus or minus.+integerApply :: (Integer -> Integer -> Integer) -> CR_ -> CR_ -> Integer+integerApply op x y = op (crToI x) (crToI y)++-- | Calls integerApply and returns a CR_.+crSimpleApply :: (Integer -> Integer -> Integer) -> CR_ -> CR_ -> CR_+crSimpleApply op x y = crFromI $ integerApply op x y++pattern PPrime :: forall a a1. (Eq a, Num a, Num a1, Ord a1) => [(a1, a)]+pattern PPrime <- [(compare 1 -> LT, 1)] -- of form x^1 where x > 1 -- prime (assumption p has already been verified to be prime)++crPrime, crHasSquare :: CR_ -> Bool++-- | Check if a number is a prime.+crPrime PPrime = True+crPrime _ = False++-- | Checks if a number has a squared (or higher) factor.+crHasSquare = any (\(_,e) -> e > 1) + +-- | Divisor functions -- should be called with integral CRs (no negative exponents).+crNumDivisors, crTau, crTotient, crPhi :: CR_ -> Integer++crNumDivisors cr = product $ map (\(_,e) -> 1 + e) cr -- does return 1 for cr1+crTau = crNumDivisors+crTotient cr = product $ map (\(p,e) -> (p-1) * p^(e-1)) cr+crPhi = crTotient++-- | Computes the nth divisor. This is zero based. +-- Note: This is deterministic but it's not ordered by the value of the divisor.+crNthDivisor :: Integer -> CR_ -> CR_+crNthDivisor 0 _ = cr1+crNthDivisor _ POne = error "Bad div num requested"+crNthDivisor n c | m == 0 = r+ | otherwise = (cb,m) : r+ where (cb,ce):cs = c+ r = crNthDivisor (div n (ce + 1)) cs -- first param is the next n+ m = mod n (ce + 1) ++-- | Consider this to be an inverse of the crNthDivisor function.+crWhichDivisor :: CR_ -> CR_ -> Integer+crWhichDivisor d c | crPositive d == False ||+ crPositive c == False = error "crWhichDivisor: Both params must be positive"+ | otherwise = f d c + where f POne _ = 0+ f _ POne = error "Not a valid divisor" + f d' c' | dp < cp || + (dp == cp && de > ce) = error "Not a valid divisor"+ | dp == cp = de + (ce + 1) * (f ds cs)+ | otherwise = (ce + 1) * (f d cs)+ where ((dp, de):ds) = d'+ ((cp, ce):cs) = c' ++-- | Efficiently computes all of the divisors based on the canonical representation.+crDivisors :: CR_ -> [CR_]+crDivisors c = foldr1 cartProd $ map pwrDivList c+ where cartProd xs ys = [x ++ y | y <- ys, x <- xs]+ pwrDivList (n,e) = [if y == 0 then cr1 else [(n,y)] | y <- [0..(fromInteger e)]]++-- | Like the crDivisors function, except that it returns pairs of the CR_ and resp. numeric value, instead of just the CR_.+crDivsPlus :: CR_ -> [(CR_, Integer)]+crDivsPlus c = foldr1 cartProd (map pwrDivList c)+ where cartProd xs ys = [(xl ++ yl, xv * yv) | (yl, yv) <- ys, (xl, xv) <- xs] + pwrDivList (e,n) = map tr $ pwrList e n+ powers x = 1 : map (* x) (powers x)+ pwrList n e = [(n,y) | y <- zip [0..e] (take (e'+1) $ powers n)] + where e' = fromInteger e+ tr (a,b) = (if fb == 0 then cr1 else [(a, fb)], sb) -- this just transforms the data structure+ where (fb, sb) = b++-- | Compute totient: Logic from deprecated arithmoi function used here.+totient :: Integer -> Integer+totient n+ | n < 1 = error "Totient only defined for positive numbers"+ | n == 1 = 1+ | otherwise = product $ map (\(p,e) -> (p-1) * p ^ (e-1)) $ factorise' n ++-- | powerModInteger adapted here from deprecated arithmoi function.+pmI :: Integer -> Integer -> Integer -> Integer+pmI x p m | x < 1 || p < 0 || m < 1 = error "pmI (powerModInteger) requires: x >= 1 &&, p >= 0, m >= 1"+ | otherwise = f p 1 (mod x m) -- last is the running exp of mod initially+ where f q w e | w == 0 || q == 0 = w+ | q == 1 = mod (w*e) m+ | otherwise = f (div q 2) nw (mod (e*e) m) + where nw | mod q 2 == 1 = mod (w*e) m + | otherwise = w+
+ Math/NumberTheory/Canon/Simple.hs view
@@ -0,0 +1,292 @@+-- |+-- Module: Math.NumberTheory.Canon.Simple+-- Copyright: (c) 2015-2018 Frederick Schneider+-- Licence: MIT+-- Maintainer: Frederick Schneider <frederick.schneider2011@gmail.com>+-- Stability: Provisional+--+-- This a wrapper for the Canonical Representation type found in the Internals module. +-- If you want to work with arbitrarily nested prime towers, you can use the Math.NumberTheory.Canon module.++{-# LANGUAGE MultiParamTypeClasses,FunctionalDependencies, FlexibleInstances, PatternSynonyms, ViewPatterns #-}++module Math.NumberTheory.Canon.Simple ( + SimpleCanon(..), SC,+ toSimpleCanon, toSC, toSimpleCanonViaUserFunc,+ fromSimpleCanon, fromSC,+ CanonConv,++ scGCD, scLCM,+ scLog, scLogDouble,+ scNegative, scPositive,+ scToInteger, scToI,+ + RationalSimpleCanon(..), RC,+ toRationalSimpleCanon, toRC, toRationalSimpleCanonViaUserFunc,+ fromRationalSimpleCanon, fromRC, + rcNegative, rcPositive, + + getNumer, getDenom, getNumerDenom,+ getNumerAsRC, getDenomAsRC, getNumerDenomAsRCs, + rcLog, rcLogDouble,+ + (>^), (<^), (%) +)+where++import GHC.Real (Ratio(..))+import Math.NumberTheory.Canon.Internals+import Math.NumberTheory.Canon.Additive+import Math.NumberTheory.Canon.AurifCyclo (crCycloInitMap)++-- | SimpleCanon is a new type wrapping a canonical representation.+newtype SimpleCanon = MakeSC CR_ deriving (Eq)++-- | This function allow you to specify a user function when converting a canon rep to an SC.+toSimpleCanonViaUserFunc :: CR_ -> (Integer -> Bool) -> SimpleCanon+toSimpleCanonViaUserFunc c f | crValidIntegralViaUserFunc c f == False = error $ invalidError + | otherwise = MakeSC c+ where invalidError = "toSimpleCanonViaUserFunc: Invalid integral canonical rep passed to constructor: " ++ (show c) ++-- | Grab the canon rep from a SimpleCanon.+fromSimpleCanon, fromSC :: SimpleCanon -> CR_+fromSimpleCanon (MakeSC i) = i+fromSC = fromSimpleCanon++-- | Shorthand type declaration+type SC = SimpleCanon++-- | Define various instances+instance Show SimpleCanon where + show c = crShow $ fromSC c+ +instance Enum SimpleCanon where+ toEnum n = toSimpleCanon $ crFromI $ fromIntegral n+ fromEnum c = fromIntegral $ crToI $ fromSC c++instance Ord SimpleCanon where+ compare x y = crCmp (fromSC x) (fromSC y)++instance Real SimpleCanon where+ toRational c = scToI c :% 1++instance Integral SimpleCanon where+ toInteger c = scToI c+ quotRem n m = (MakeSC n', MakeSC m') + where (n', m') = fst $ crQuotRem (fromSC n) (fromSC m) crCycloInitMap+ mod n m = MakeSC $ crMod (fromSC n) (fromSC m)+ +instance Fractional SimpleCanon where+ fromRational (n :% d) | m == 0 = MakeSC $ crFromI q+ | otherwise = error "Modulus not zero. Use Rational SimpleCanons for non-integers."+ where (q, m) = quotRem n d+ (/) x y = MakeSC $ crDivStrict (fromSC x) (fromSC y)++instance Num SimpleCanon where+ fromInteger n = MakeSC $ crFromI n -- to do: check where called?+ x + y = MakeSC $ fst $ crAdd (fromSC x) (fromSC y) crCycloInitMap -- discard the map info+ x - y = MakeSC $ fst $ crSubtract (fromSC x) (fromSC y) crCycloInitMap -- discard the map info+ x * y = MakeSC $ crMult (fromSC x) (fromSC y)+ + negate x = MakeSC $ crNegate $ fromSC x+ abs x = MakeSC $ crAbs $ fromSC x+ signum x = MakeSC $ crSignum $ fromSC x++-- | Convert a SimpleCanon back to an Integer.+scToInteger, scToI :: SimpleCanon -> Integer+scToI c = crToI $ fromSC c+scToInteger = scToI++-- | SimpleCanon GCD and LCM functions+scGCD, scLCM :: SimpleCanon -> SimpleCanon -> SimpleCanon+scGCD x y = MakeSC $ crGCD (fromSC x) (fromSC y)+scLCM x y = MakeSC $ crLCM (fromSC x) (fromSC y)++-- | Wrappers for underlying canon rep functions+scNegative, scPositive :: SimpleCanon -> Bool+scNegative c = crNegative $ fromSC c+scPositive c = crPositive $ fromSC c++-- | Wrapper for underlying CR function+scLog :: SimpleCanon -> Rational+scLog x = crLog $ fromSC x ++-- | Wrapper for underlying CR function+scLogDouble :: SimpleCanon -> Double+scLogDouble x = crLogDouble $ fromSC x ++-- | Newtype for RationalSimpleCanon. The underlying canon rep is the same.+newtype RationalSimpleCanon = MakeRC CR_ deriving (Eq)++-- | Convert canon rep to RationalSimpleCanon via a user-supplied criterion function.+toRationalSimpleCanonViaUserFunc :: CR_ -> (Integer -> Bool) -> RationalSimpleCanon+toRationalSimpleCanonViaUserFunc c f | crValidRationalViaUserFunc c f == False = error $ invalidError + | otherwise = MakeRC c+ where invalidError = + "toRationalSimpleCanonViaUserFunc: Invalid rational canonical rep passed to constructor: " + ++ (show c) ++ " (user predicate supplied)" ++-- | Convert RC back to underlying canon rep.+fromRationalSimpleCanon, fromRC :: RationalSimpleCanon -> CR_+fromRC (MakeRC i) = i+fromRationalSimpleCanon = fromRC++-- | Shorthand type name +type RC = RationalSimpleCanon++-- | Define various instances for RationSimpleCanon.+instance Show RationalSimpleCanon where + show rc = crShowRational $ fromRC rc+ +instance Enum RationalSimpleCanon where+ toEnum n = toRC $ crFromI $ fromIntegral n+ fromEnum c = fromIntegral $ toInteger c -- if not integral, this will fail++-- | Caveat: These functions will error out (in)directly if there are any negative exponents.+instance Integral RationalSimpleCanon where+ toInteger rc = crToI $ fromRC rc+ quotRem n m | crIntegral $ fromRC n = (MakeRC n', MakeRC m') + | otherwise = error "Can't perform 'quotRem' on non-integral RationalSimpleCanon"+ where (n', m') = fst $ crQuotRem (fromRC n) (fromRC m) crCycloInitMap+ mod n m | crIntegral $ fromRC n = MakeRC $ crMod (fromRC n) (fromRC m) + | otherwise = error "Can't perform 'mod' on non-integral RationalSimpleCanon"++instance Fractional RationalSimpleCanon where+ fromRational (n :% d) = MakeRC $ crDivRational (crFromI n) (crFromI d)+ (/) x y = MakeRC $ crDivRational (fromRC x) (fromRC y)++instance Ord RationalSimpleCanon where+ compare x y = crCmp (fromRC x) (fromRC y)+ +instance Real RationalSimpleCanon where+ toRational rc = crToRational $ fromRC rc+ +instance Num RationalSimpleCanon where+ fromInteger n = MakeRC $ crFromI n+ x + y = MakeRC $ fst $ crAddR (fromRC x) (fromRC y) crCycloInitMap+ x - y = MakeRC $ fst $ crSubtractR (fromRC x) (fromRC y) crCycloInitMap+ x * y = MakeRC $ crMult (fromRC x) (fromRC y) + + negate x = MakeRC $ crNegate $ fromRC x+ abs x = MakeRC $ crAbs $ fromRC x + signum x = MakeRC $ crSignum $ fromRC x++-- | Calls underlying canon rep function.+rcLog :: RationalSimpleCanon -> Rational+rcLog c = crLog $ fromRC c ++-- | Calls underlying canon rep function. +rcLogDouble :: RationalSimpleCanon -> Double+rcLogDouble c = crLogDouble $ fromRC c++-- | Calls underlying canon rep function. +getNumerAsRC :: RationalSimpleCanon -> RationalSimpleCanon+getNumerAsRC c = MakeRC $ crNumer $ fromRC c+ +-- | Calls underlying canon rep function. +getDenomAsRC :: RationalSimpleCanon -> RationalSimpleCanon+getDenomAsRC c = MakeRC $ crDenom $ fromRC c++-- | Pulls numerator or denominator from RC and converts it to a SimpleCanon.+getNumer, getDenom :: RationalSimpleCanon -> SimpleCanon+getNumer c = MakeSC $ crNumer $ fromRC c+getDenom c = MakeSC $ crDenom $ fromRC c ++-- | Wraps crSplit function and returns a pair of SimpleCanons.+getNumerDenom :: RationalSimpleCanon -> (SimpleCanon, SimpleCanon)+getNumerDenom c = (MakeSC n, MakeSC d) + where (n, d) = crSplit $ fromRC c+ +-- | Wraps crSplit function and returns a pair of RationalSimpleCanons. +getNumerDenomAsRCs :: RationalSimpleCanon -> (RationalSimpleCanon, RationalSimpleCanon)+getNumerDenomAsRCs c = (MakeRC n, MakeRC d) + where (n, d) = crSplit $ fromRC c ++-- | Wraps underlying canon rep functions.+rcNegative, rcPositive :: RationalSimpleCanon -> Bool+rcNegative x = crNegative $ fromRC x+rcPositive x = crPositive $ fromRC x ++-- | Modulus operator+infixl 7 %+(%) :: (Integral a) => a -> a -> a+n % m = mod n m +++-- | Multi-param typeclass for exponentiation+infixr 9 <^++class SimpleCanonExpnt a b c | a b -> c where + -- | Exponentiation Operator+ (<^) :: a -> b -> c++instance SimpleCanonExpnt Integer Integer SimpleCanon where+ p <^ e = MakeSC $ crExp (crFromI p) e False++instance SimpleCanonExpnt SimpleCanon Integer SimpleCanon where+ p <^ e = MakeSC $ crExp (fromSC p) e False++instance SimpleCanonExpnt RationalSimpleCanon Integer RationalSimpleCanon where+ p <^ e = MakeRC $ crExp (fromRC p) e True++instance SimpleCanonExpnt RationalSimpleCanon SimpleCanon RationalSimpleCanon where+ p <^ e = MakeRC $ crExp (fromRC p) (crToI $ fromSC e) True++-- | Multi-param typeclass for radical/root operator+infixr 9 >^ -- r >^ n means attempt to take the rth root of n++class SimpleCanonRoot a b c | a b -> c where+ -- | Root Operator+ (>^) :: a -> b -> c++instance SimpleCanonRoot SimpleCanon SimpleCanon SimpleCanon where+ r >^ n = MakeSC $ crRoot (fromSC n) (toInteger r)+ +instance SimpleCanonRoot Integer Integer SimpleCanon where+ r >^ n = MakeSC $ crRoot (crFromI n) r+ +instance SimpleCanonRoot Integer SimpleCanon SimpleCanon where+ r >^ n = MakeSC $ crRoot (fromSC n) r++instance SimpleCanonRoot SimpleCanon Integer SimpleCanon where+ r >^ n = MakeSC $ crRoot (crFromI n) (toInteger r) + +instance SimpleCanonRoot Integer RationalSimpleCanon RationalSimpleCanon where+ r >^ n = MakeRC $ crRoot (fromRC n) r++-- | Convert CanonConv types to SimpleCanon.+toSimpleCanon :: (CanonConv a) => a -> SimpleCanon+toSimpleCanon = toSC++-- | Convert CanonConv types to RationalSimpleCanon.+toRationalSimpleCanon :: (CanonConv a) => a -> RationalSimpleCanon+toRationalSimpleCanon = toRC++-- | Typeclass for converting to SimpleCanon and RationalSimpleCanon +class CanonConv c where+ -- | Convert a type to a SimpleCanon.+ toSC :: c -> SimpleCanon ++ -- | Convert a type to a RationalSimpleCanon.+ toRC :: c -> RationalSimpleCanon + +instance CanonConv SimpleCanon where + toSC c = c+ toRC c = MakeRC $ fromSC c++instance CanonConv CR_ where + toSC cr | crValidIntegral cr = MakeSC cr+ | otherwise = error invalidError+ where invalidError = "Invalid integral canonical rep passed to constructor: " ++ (show cr) + + toRC cr | crValidRational cr = MakeRC cr+ | otherwise = error invalidRepRatioError + where invalidRepRatioError = "toRC: Invalid canonical rep passed to constructor: " ++ (show cr) + +instance CanonConv RationalSimpleCanon where + toSC rc | crValidIntegral frc = MakeSC frc+ | otherwise = error invalidError+ where frc = fromRC rc+ invalidError = "Invalid integral canonical rep passed to constructor: " ++ (show rc) + toRC rc = rc+
+ README.md view
@@ -0,0 +1,1 @@+# canon
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ canon.cabal view
@@ -0,0 +1,38 @@+-- Initial canon.cabal generated by cabal init. For further documentation,+-- see http://haskell.org/cabal/users-guide/++name: canon+version: 0.1.0.0+synopsis: Massive Number Arithmetic+description: This library allows one to manipulate of practically unlimited by keeping them in factored "canonical" form.+ For manipulating sums and differences, there is additional code to factor expressions of special forms.++homepage: https://github.com/grandpascorpion/canon+license: MIT+license-file: LICENSE+author: Frederick Schneider+maintainer: frederick dot schneider2011 at gmail dot com+-- copyright: +category: Math+build-type: Simple+extra-source-files: Changes, README.md+cabal-version: >=1.10++library+ build-depends : base >= 4.9.1.0 && < 5+ , arithmoi >= 0.6.0.1 && < 0.7+ , polynomial >= 0.7.3 && < 0.8+ , array >= 0.5.1.1 && < 0.6+ , containers >= 0.5.7.1 && < 0.6++ exposed-modules : Math.NumberTheory.Canon+ Math.NumberTheory.Canon.Simple+ Math.NumberTheory.Canon.AurifCyclo+ Math.NumberTheory.Canon.Internals+-- other-modules : Math.NumberTheory.Canon.Internals+ Math.NumberTheory.Canon.Additive++ ghc-options : -O2 -Wall+ ghc-prof-options : -O2 -auto++ default-language: Haskell2010
+ goBigOrGoHome.odp view
binary file changed (absent → 55747 bytes)
+ test-suite/CanonManualTests.hs view
@@ -0,0 +1,135 @@+-- |+-- Module: Math.NumberTheory.CanonTests+-- Copyright: (c) 2018 Andrew Lelechenko+-- Licence: MIT+-- Maintainer: Frederick Schneider <frederick.schneider2011@gmail.com> +-- Stability: Provisional+--+-- Tests for Math.NumberTheory.Canon, etc+--++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++import Math.NumberTheory.Canon+import Math.NumberTheory.Canon.AurifCyclo+import Data.Array (array)+import Control.Monad (forever)++trueS, falseS :: String+trueS = "true"+falseS = "FALSE"++divvyTest, aCaDTest, aDTest, cATest, cATestM, aCPTest :: Bool+divvyTest = ans == divvy a x y+ where a = [4,7,8401,62324358089100907319521,682969,61,374857981681,547]+ x = 50031545486420197+ y = 50031544711579219+ ans = [765980456641,81365467681,374857981681,301,2269,31,271,547,61,7,4]++aCaDTest = r == aurCandDec (2^58) 1 True+ where r = Just (536838145,536903681)++aDTest = r == aurDec 5+ where r = Just (array (0,2) [(0,1),(1,3),(2,1)],array (0,1) [(0,1),(1,1)])++--chineseAurif: Any non-zero multiple of (q^2m * p) ^ (p * k) + (r^2n)^(p * k) where k is odd pos, m, n > 0+-- q and r do not both equal 1++cATest = r == chineseAurif (44^253) 1 True -- Equivalent to 44^253 + 1. Should these two factors below. Example from Sichuan 5 paper+ where r = Just (+ 3708082114051284931014527275382936962949050019900548504948093002539948192457694962513241254988377338102340862648630965276420678480576906389289483833735873261700512602622143146599971,+ 10004590597907985573943582945748620748239251502916976018978239877682278432398712396908419662400063010898795149152694380672517014008143612228221361453877714927361019333917217066917231+ )++cATestM = r == chineseAurif ((q*p)^(p*m)) (s^(p*m)) False -- Equivalent to 117^221 - 4^221.+ where (p, q, m, s) = (13, 9, 17, 4)+ r = Just (+ 6777177566148891825866597484460647604677595599339380316570867390255387995775996275965289196556216301121524986265286043390240164764352880415910164128699689628228970391585087480017942725930708815238591,+ 1760066887061200649747400254189929691148780525097352459006167529137160481395362942899568367387185016629815611235574356459997927884738618313122568841325847458481738026543495018488267067010986061866931+ ) ++aCPTest = boolFlag + where (_, _, boolFlag) = verifyApplyCycloPair 5 3 2310++verifyApplyCycloPair :: Integer -> Integer -> Integer -> (Integer, [Integer], Bool)+verifyApplyCycloPair x y e = (v, factors, v == product factors) -- should be True+ where v = x^e - y^e+ factors = applyCycloPair x y e++{-+-- fix this+f p e = aurCandDec (p^(e*p)) 1 True +f 11 15+Just (1271895306126722839332303077175663680408337203898205913979279374031646228168717,1271895436663409913619808282471754544629810369181390985798959949855832620283803)++oddExpAur p em = aurCandDec (p^(em*p)) 1 True+evenExpAur p em = aurCandDec (p^(em*p)) 1 False+-}++canonBasicProperties :: Int -> [(Int, String)]+canonBasicProperties m | m <= 0 = error "Positive ints only"+ | otherwise = formatTestOutput tests+ where tests = [mc * rc == c1,+ mc / mc == c1,+ (toInteger $ mc - mc) == 0,+ mc + mc == mc * 2,+ mc ^ 2 == mc * mc,++ mc + c1 > mc,+ mc - 1 < mc,+ mc == negate nmc, -- double negative+ mc == abs nmc,+ mc == cReciprocal rc, -- double reciprocal++ signum mc == c1,+ signum nmc == negate c1,+ toInteger mc == toInteger m, -- undo the conversion+ cNegative nmc && cPositive mc,+ cIntegral mc && cIntegral nmc,++ --XXXX cIrrational sr && + --XXXX sr <^ mp1 == c1 &&+ oc /= ec && (oc || ec) + ] + mc = makeCanon $ toInteger $ m+ rc = cReciprocal mc + nmc = negate mc+ -- mp1 = mc + c1+ -- sr = mp1 >^ mp1 -- take the xth root of x (x is > 1). Must be irrational+ oc = cOdd mc+ ec = cOdd $ mc + c1 -- even check+ c1 = makeCanon 1++canonBasicProperties2 :: Int -> Int -> [(Int, String)]+canonBasicProperties2 m n | m <= 0 || n <= 0 = error "m and n must be positive"+ | otherwise = formatTestOutput tests+ where + tests = [mc * nc / g == cLCM mc nc,+ mc == q * nc + r,+ mc == nc >^ mxn, -- test the root operator+ r == mod mc nc+ ]+ mc = makeCanon $ toInteger m+ nc = makeCanon $ toInteger n+ g = makeCanon $ toInteger $ gcd m n+ (q,r) = quotRem mc nc+ mxn = mc <^ nc -- exponentiation++formatTestOutput :: [Bool] -> [(Int, String)]+formatTestOutput tests = zip [1..] $ map (\b -> if b then trueS else falseS) tests++main :: IO ()+main = forever $ do+ print "Canon Basic Properties (Enter 1 param): " + p <- getLine+ print $ canonBasicProperties (read p :: Int) + print ""+ print "Canon Basic Properties (Enter 2 params, one each line): "+ p1 <- getLine+ p2 <- getLine+ print $ canonBasicProperties2 (read p1 :: Int) (read p2 :: Int)+ print ""+ print "Canon Specific Tests (0 params): "+ print $ formatTestOutput [divvyTest, aCaDTest, aDTest, cATest, cATestM, aCPTest]+ print ""+