diff --git a/Changes b/Changes new file mode 100644 --- /dev/null +++ b/Changes @@ -0,0 +1,2 @@ +0.1.0.0: + First release diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,675 @@ + GNU GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + 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 . + +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 +. + + 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 +. + diff --git a/Math/NumberTheory/Canon.hs b/Math/NumberTheory/Canon.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Canon.hs @@ -0,0 +1,879 @@ +-- | +-- Module: Math.NumberTheory.Canon +-- Copyright: (c) 2015-2018 Frederick Schneider +-- Licence: MIT +-- Maintainer: Frederick Schneider +-- 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" +-} + diff --git a/Math/NumberTheory/Canon/Additive.hs b/Math/NumberTheory/Canon/Additive.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Canon/Additive.hs @@ -0,0 +1,109 @@ +-- | +-- Module: Math.NumberTheory.Canon.Additive +-- Copyright: (c) 2015-2018 Frederick Schneider +-- Licence: MIT +-- Maintainer: Frederick Schneider +-- 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 diff --git a/Math/NumberTheory/Canon/AurifCyclo.hs b/Math/NumberTheory/Canon/AurifCyclo.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Canon/AurifCyclo.hs @@ -0,0 +1,490 @@ +-- | +-- Module: Math.NumberTheory.Canon.AurifCyclo +-- Copyright: (c) 2015-2018 Frederick Schneider +-- Licence: MIT +-- Maintainer: Frederick Schneider +-- 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 + diff --git a/Math/NumberTheory/Canon/Internals.hs b/Math/NumberTheory/Canon/Internals.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Canon/Internals.hs @@ -0,0 +1,590 @@ +-- | +-- Module: Math.NumberTheory.Canon.Internals +-- Copyright: (c) 2015-2018 Frederick Schneider +-- Licence: MIT +-- Maintainer: Frederick Schneider +-- 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 + diff --git a/Math/NumberTheory/Canon/Simple.hs b/Math/NumberTheory/Canon/Simple.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Canon/Simple.hs @@ -0,0 +1,292 @@ +-- | +-- Module: Math.NumberTheory.Canon.Simple +-- Copyright: (c) 2015-2018 Frederick Schneider +-- Licence: MIT +-- Maintainer: Frederick Schneider +-- 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 + diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,1 @@ +# canon diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/canon.cabal b/canon.cabal new file mode 100644 --- /dev/null +++ b/canon.cabal @@ -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 diff --git a/goBigOrGoHome.odp b/goBigOrGoHome.odp new file mode 100644 Binary files /dev/null and b/goBigOrGoHome.odp differ diff --git a/test-suite/CanonManualTests.hs b/test-suite/CanonManualTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/CanonManualTests.hs @@ -0,0 +1,135 @@ +-- | +-- Module: Math.NumberTheory.CanonTests +-- Copyright: (c) 2018 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Frederick Schneider +-- 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 "" +