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. + + + Copyright (C) + + 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: + + Copyright (C) + 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/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,5 @@ +#!/usr/bin/env runhaskell + +import qualified Distribution.Simple + +main = Distribution.Simple.defaultMain diff --git a/changelog b/changelog new file mode 100644 --- /dev/null +++ b/changelog @@ -0,0 +1,16 @@ +2011-03-01 Dr. Alistair Ward + +0.0.0.1 + * First version of the package. +0.0.0.2 + * Created modules; "Factory.Test.QuickCheck.Bounds", "Factory.Math.Implementations.Pi.Borwein.*" and "Factory.Test.Performance.Statistics". + * Created a new module "Factory.Data.PrimeFactors", and migrated definitions from both "Factory.Math.PrimeFactorisation" and "Factory.Math.Implementations.PrimeFactorisation". + * Created class 'Factory.Math.Factorial.Factorial' and new module "Factory.Math.Implementations.Factorial". + Moved existing implementation (Bisection) into new module, with new implementation (PrimeFactorisation). + * Added function 'Factory.Math.Summation.sumR'. + * Added a parameter to functions 'Factory.Math.DivideAndConquer.divideAndConquer' and 'Factory.Data.Bounds.divideAndConquer' to permit asymmetric bisection. + * Added methods to class "Factory.Math.Pi.Algorithm" to permit the retrieval of /Pi/ as a 'Rational' or a 'String'. + * Renamed 'Factory.Math.Precision.capPrecision' to 'Factory.Math.Precision.simplify'. + * Removed module "Factory.Test.Performance.Exponential". + * Removed function 'Factory.Math.Power.raise', which was no more efficient than ghc's implementation of '(^)'. +0.0.0.3 diff --git a/copyright b/copyright new file mode 100644 --- /dev/null +++ b/copyright @@ -0,0 +1,11 @@ +Author: + Dr. Alistair Ward . + +Copyright: + Copyright (C) 2011 Dr. Alistair Ward. All Rights Reserved. + +Home-page: + http://functionalley.eu + +License: + GNU GENERAL PUBLIC LICENSE Version 3; see '/usr/share/common-licenses/GPL-3' or '/usr/share/doc/licenses/gpl-3.0.txt' where available, or the local packaged file 'LICENSE'. diff --git a/factory.cabal b/factory.cabal new file mode 100644 --- /dev/null +++ b/factory.cabal @@ -0,0 +1,144 @@ +--Package-properties +Name: factory +Version: 0.0.0.2 +Cabal-Version: >= 1.6 +Copyright: (C) 2011 Dr. Alistair Ward +License: GPL +License-file: LICENSE +Author: Dr. Alistair Ward +Stability: Unstable interface, incomplete features. +Synopsis: Rational arithmetic in an irrational world. +Build-Type: Simple +Description: A library of number-theory functions, for; factorials, square-roots, Pi, primality-testing, prime-factorisation ... +Category: Math, Number Theory +Tested-With: GHC == 6.12, GHC == 7.0 +Homepage: http://functionalley.eu +Maintainer: factory functionalley eu +Bug-reports: factory functionalley eu +Extra-Source-Files: changelog, copyright, makefile + +flag llvm + Description: Whether the 'llvm' compiler-backend has been installed and is required for code-generation. + manual: True + default: False + +flag threaded + Description: Enable parallelized code. + default: True + +Library + hs-source-dirs: src + + Exposed-modules: + Factory.Data.Bounds + Factory.Data.Exponential + Factory.Data.MonicPolynomial + Factory.Data.Monomial + Factory.Data.Polynomial + Factory.Data.PrimeFactors + Factory.Data.QuotientRing + Factory.Data.Ring + Factory.Math.ArithmeticGeometricMean + Factory.Math.DivideAndConquer + Factory.Math.Factorial + Factory.Math.Fibonacci + Factory.Math.Implementations.Factorial + Factory.Math.Implementations.Primality + Factory.Math.Implementations.PrimeFactorisation + Factory.Math.Implementations.SquareRoot + Factory.Math.MultiplicativeOrder + Factory.Math.Pi + Factory.Math.Implementations.Pi.AGM.Algorithm + Factory.Math.Implementations.Pi.AGM.BrentSalamin + Factory.Math.Implementations.Pi.BBP.Algorithm + Factory.Math.Implementations.Pi.BBP.Base65536 + Factory.Math.Implementations.Pi.BBP.Bellard + Factory.Math.Implementations.Pi.BBP.Implementation + Factory.Math.Implementations.Pi.BBP.Series + Factory.Math.Implementations.Pi.Borwein.Algorithm + Factory.Math.Implementations.Pi.Borwein.Borwein1993 + Factory.Math.Implementations.Pi.Borwein.Implementation + Factory.Math.Implementations.Pi.Borwein.Series + Factory.Math.Implementations.Pi.Ramanujan.Algorithm + Factory.Math.Implementations.Pi.Ramanujan.Chudnovsky + Factory.Math.Implementations.Pi.Ramanujan.Classic + Factory.Math.Implementations.Pi.Ramanujan.Implementation + Factory.Math.Implementations.Pi.Ramanujan.Series + Factory.Math.Implementations.Pi.Spigot.Algorithm + Factory.Math.Implementations.Pi.Spigot.Gosper + Factory.Math.Implementations.Pi.Spigot.RabinowitzWagon + Factory.Math.Implementations.Pi.Spigot.Series + Factory.Math.Implementations.Pi.Spigot.Spigot + Factory.Math.Power + Factory.Math.Precision + Factory.Math.Primality + Factory.Math.PrimeFactorisation + Factory.Math.Radix + Factory.Math.SquareRoot + Factory.Math.Statistics + Factory.Math.Summation + + Build-depends: + array, + base == 4.*, + deepseq >= 1.1, + containers, + primes >= 0.1, + toolshed == 0.11.* + + if flag(threaded) + Build-depends: parallel >= 3.0 + else + Build-depends: parallel + + GHC-options: -Wall -O2 + GHC-prof-options: -prof -auto-all -caf-all + + if impl(ghc >= 7.0) && flag(llvm) + GHC-options: -fllvm + +Executable factory + hs-source-dirs: src + + Main-Is: Main.hs + + Other-modules: + Factory.Test.CommandOptions + Factory.Test.Performance.Factorial + Factory.Test.Performance.Pi + Factory.Test.Performance.Primality + Factory.Test.Performance.PrimeFactorisation + Factory.Test.Performance.SquareRoot + Factory.Test.Performance.Statistics + Factory.Test.QuickCheck.ArithmeticGeometricMean + Factory.Test.QuickCheck.Bounds + Factory.Test.QuickCheck.Factorial + Factory.Test.QuickCheck.MonicPolynomial + Factory.Test.QuickCheck.Pi + Factory.Test.QuickCheck.Polynomial + Factory.Test.QuickCheck.Power + Factory.Test.QuickCheck.Primality + Factory.Test.QuickCheck.PrimeFactorisation + Factory.Test.QuickCheck.Radix + Factory.Test.QuickCheck.QuickChecks + Factory.Test.QuickCheck.SquareRoot + Factory.Test.QuickCheck.Statistics + Factory.Test.QuickCheck.Summation + + Build-depends: + Cabal >= 1.6 && < 2, + haskell98, + QuickCheck >= 2.2 + + GHC-options: -Wall -O2 + GHC-prof-options: -prof -auto-all -caf-all + + if flag(threaded) + GHC-options: -threaded -feager-blackholing + + if impl(ghc >= 7.0) + GHC-options: -rtsopts + + if flag(llvm) + GHC-options: -fllvm + diff --git a/makefile b/makefile new file mode 100644 --- /dev/null +++ b/makefile @@ -0,0 +1,53 @@ +# Copyright (C) 2011 Dr. Alistair Ward +# +# 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 . + +.PHONY: all build clean configure copy haddock help hlint install prof sdist + +all: install + +install: build haddock + @[ -z "$$CABAL_INSTALL_OPTIONS" ] || echo "INFO: CABAL_INSTALL_OPTIONS='$$CABAL_INSTALL_OPTIONS'" + runhaskell Setup.hs $@ $$CABAL_INSTALL_OPTIONS + +prof: + CABAL_CONFIGURE_OPTIONS="--enable-library-profiling --enable-executable-profiling $$CABAL_CONFIGURE_OPTIONS" make install + +copy: build + @[ -z "$$CABAL_COPY_OPTIONS" ] || echo "INFO: CABAL_COPY_OPTIONS='$$CABAL_COPY_OPTIONS'" + runhaskell Setup.hs $@ $$CABAL_COPY_OPTIONS + +build: configure + @[ -z "$$CABAL_BUILD_OPTIONS" ] || echo "INFO: CABAL_BUILD_OPTIONS='$$CABAL_BUILD_OPTIONS'" + runhaskell Setup.hs $@ $$CABAL_BUILD_OPTIONS + +configure: factory.cabal Setup.hs + @[ -z "$$CABAL_CONFIGURE_OPTIONS" ] || echo "INFO: CABAL_CONFIGURE_OPTIONS='$$CABAL_CONFIGURE_OPTIONS'" + runhaskell Setup.hs $@ $$CABAL_CONFIGURE_OPTIONS #--user + +haddock: configure + PATH=~/.cabal/bin:$$PATH runhaskell Setup.hs $@ --hyperlink-source #Amend path to find 'HsColour', as required for 'hyperlink-source'. + +hlint: + @$@ src/ + +sdist: configure + runhaskell Setup.hs $@ + +clean: + runhaskell Setup.hs $@ + +help: + @grep '^[a-zA-Z].*:' makefile | sed -e 's/:.*//' + diff --git a/src/Factory/Data/Bounds.hs b/src/Factory/Data/Bounds.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Data/Bounds.hs @@ -0,0 +1,170 @@ +{-# LANGUAGE CPP #-} +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] + + * Describes a /bounded/ range of, typically integral, quantities. + + * Operations have been defined, on the list of /consecutive/ quantities delimited by these two bounds. + + * The point is that if the list is composed from /consecutive/ quantities, the intermediate values can be inferred, rather than physically represented. + + [@CAVEATS@] + + * The API was driven top-down by its caller's requirements, rather than a bottom-up attempt to provide a complete interface. + consequently there may be omissions from the view point of future callers. +-} + +module Factory.Data.Bounds ( +-- * Types +-- ** Type-synonyms + Bounds, +-- * Functions +-- divideAndConquer, + elem', + length', + normalise, + product', + splitAt', + toList, +-- ** Accessors + minBound', + maxBound' +-- ** Predicates +-- isReversed +) where + +import Control.Arrow((***)) +import qualified Data.Monoid +import qualified Data.Ratio + +#if MIN_VERSION_parallel(3,0,0) +import qualified Control.Parallel.Strategies +#endif + +#if MIN_VERSION_base(4,3,0) +import Data.Tuple(swap) +#else +-- | Swap the components of a pair. +swap :: (a, b) -> (b, a) +swap (a, b) = (b, a) +#endif + +-- | Defines a range of consecutive values, bracketed by /inclusive/ bounds. +type Bounds limit = (limit, limit) + +-- | Accessor. +{-# INLINE minBound' #-} +minBound' :: Bounds a -> a +minBound' = fst + +-- | Accessor. +{-# INLINE maxBound' #-} +maxBound' :: Bounds a -> a +maxBound' = snd + +-- | 'True' if the specified value is within the inclusive 'Bounds'. +elem' :: Ord limit => limit -> Bounds limit -> Bool +elem' x = uncurry (&&) . ((<= x) *** (x <=)) + +-- | 'True' if /minBound'/ exceeds /maxBound'/ extent. +isReversed :: Ord limit => Bounds limit -> Bool +isReversed = uncurry (>) + +-- | Swap the limits where they were originally reversed, but otherwise do nothing. +normalise :: Ord limit => Bounds limit -> Bounds limit +normalise b + | isReversed b = swap b + | otherwise = b + +-- | Bisect the bounds at the specified limit; which should be between the two existing limits. +splitAt' :: (Num limit, Ord limit) => limit -> Bounds limit -> (Bounds limit, Bounds limit) +splitAt' i bounds@(l, r) + | any ($ i) [(< l), (>= r)] = error $ "Factory.Data.Bounds.splitAt':\tunsuitable index=" ++ show i ++ " for bounds=" ++ show bounds ++ "." + | otherwise = ((l, i), (i + 1, r)) + +-- | The length of 'toList'. +{-# INLINE length' #-} +length' :: (Num limit, Ord limit) => Bounds limit -> limit +length' (l, r) = r + 1 - l + +-- | Converts 'Bounds' to a list by enumerating the values. +{-# INLINE toList #-} +toList :: Enum limit => Bounds limit -> [limit] +toList = uncurry enumFromTo + +{- | + * Reduces 'Bounds' to a single integral value encapsulated in a 'Data.Monoid.Monoid', + using a /divide-and-conquer/ strategy, + bisecting the /bounds/ and recursively evaluating each part; . + + * By choosing a 'ratio' other than @(1 % 2)@, the bisection can be made asymmetrical. + The specified ratio represents the length of the left-hand portion over the original list-length; + eg. @(1 % 3)@ results in the first part, half the length of the second. + + * This process of recursive bisection, is terminated beneath the specified minimum length, + after which the 'Bounds' are expanded into the corresponding list, and the /monoid/'s binary operator is directly /folded/ over it. + + * One can view this as a , + in which 'Bounds' is exploded into a binary tree-structure + (each leaf of which contains a list of up to 'minLength' integers, and each node of which contains an associative binary operator), + and then collapsed to a scalar, by application of the operators. +-} +divideAndConquer :: (Integral i, Data.Monoid.Monoid monoid) + => (i -> monoid) -- ^ The monoid's constructor. + -> Data.Ratio.Ratio i -- ^ The ratio of the original span, at which to bisect the 'Bounds'. + -> i -- ^ For efficiency, the bounds will not be bisected, when it's length has been reduced to this value. + -> Bounds i + -> monoid -- ^ The resulting scalar. +divideAndConquer monoidConstructor ratio minLength + | any ($ ratio) [ + (< 0), + (>= 1) + ] = error $ "Factory.Data.Bounds.divideAndConquer:\tunsuitable ratio='" ++ show ratio ++ "'." + | minLength < 1 = error $ "Factory.Data.Bounds.divideAndConquer:\tunsuitable minLength=" ++ show minLength ++ "." + | otherwise = slave + where + slave bounds@(l, r) + | length' bounds <= minLength = Data.Monoid.mconcat . map monoidConstructor $ toList bounds --Fold the monoid's binary operator over the delimited list. + | otherwise = uncurry Data.Monoid.mappend . +#if MIN_VERSION_parallel(3,0,0) + Control.Parallel.Strategies.withStrategy ( + Control.Parallel.Strategies.parTuple2 Control.Parallel.Strategies.rseq Control.Parallel.Strategies.rseq + ) . +#endif + (slave *** slave) $ splitAt' ( + l + (r - l) * Data.Ratio.numerator ratio `div` Data.Ratio.denominator ratio --Use the ratio to generate the split-index. + ) bounds --Apply the monoid's binary operator to the two operands resulting from bisection. + +{- | + * Multiplies the consecutive sequence of integers within 'Bounds'. + + * Since the result can be large, 'divideAndConquer' is used to form operands of a similar order of magnitude, + thus improving the efficiency of the big-number multiplication. +-} +product' :: Integral i + => Data.Ratio.Ratio i -- ^ The ratio at which to bisect the 'Bounds'. + -> i -- ^ For efficiency, the bounds will not be bisected, when it's length has been reduced to this value. + -> Bounds i + -> i -- ^ The resulting product. +product' ratio minLength bounds + | elem' 0 bounds = 0 + | otherwise = Data.Monoid.getProduct $ divideAndConquer Data.Monoid.Product ratio minLength bounds + diff --git a/src/Factory/Data/Exponential.hs b/src/Factory/Data/Exponential.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Data/Exponential.hs @@ -0,0 +1,89 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] + + * Describes a simple numeric type, designed to contain an /exponential/ number. + + * . +-} + +module Factory.Data.Exponential( +-- * Types +-- ** Type-synonyms + Exponential, +-- * Functions + evaluate, + invert, +-- ** Accessors + getBase, + getExponent, +-- ** Constructors + rightIdentity, +-- ** Operators + (<^), + (=~) +) where + +import qualified Control.Arrow + +infix 4 =~ --Same as (==). +infixr 8 <^ --Same as (^). + +-- | Describes an /exponential/, in terms of its /base/ and /exponent/. +type Exponential base exponent = (base, exponent) + +-- | Accessor. +{-# INLINE getBase #-} +getBase :: Exponential base exponent -> base +getBase = fst + +-- | Accessor. +{-# INLINE getExponent #-} +getExponent :: Exponential base exponent -> exponent +getExponent = snd + +{- | + * Construct an 'Exponential' merely raised to the 1st power. + + * The value of the resulting exponential is the same as specified 'base'; . +-} +rightIdentity :: Num exponent => base -> Exponential base exponent +rightIdentity x = (x, 1) + +-- | Evaluate the specified 'Exponential', returning the resulting number. +{-# INLINE evaluate #-} +evaluate :: (Num base, Integral exponent) => Exponential base exponent -> base +evaluate = uncurry (^) + +-- | 'True' if the /bases/ are equal. +(=~) :: Eq base => Exponential base exponent -> Exponential base exponent -> Bool +(l, _) =~ (r, _) = l == r + +-- | Raise the specified 'Exponential' to a power. +(<^) :: Num exponent + => Exponential base exponent -- ^ The operand. + -> exponent -- ^ The power to which the exponential is to be raised. + -> Exponential base exponent -- ^ The result. +(b, e) <^ power = (b, e * power) + +-- | Invert the value, by negating the exponent. +invert :: Num exponent => Exponential base exponent -> Exponential base exponent +invert = Control.Arrow.second negate + diff --git a/src/Factory/Data/MonicPolynomial.hs b/src/Factory/Data/MonicPolynomial.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Data/MonicPolynomial.hs @@ -0,0 +1,83 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] + + * Describes a /monic polynomial; ; + ie. in which the /coefficient/ of the /leading term/ is one. +-} + +module Factory.Data.MonicPolynomial( +-- * Types +-- ** Data-types, + MonicPolynomial(getPolynomial), --Hide the data-constructor. +-- * Functions +-- ** Constructors + mkMonicPolynomial +) where + +import Control.Arrow((***)) +import qualified Control.Arrow +import qualified Factory.Data.Monomial as Data.Monomial +import Factory.Data.Polynomial((*=)) +import qualified Factory.Data.Polynomial as Data.Polynomial +import qualified Factory.Data.QuotientRing as Data.QuotientRing +import Factory.Data.Ring((=*=), (=+=), (=-=)) +import qualified Factory.Data.Ring as Data.Ring + +-- | A type of 'Data.Polynomial.Polynomial', in which the /leading term/ is required to have a /coefficient/ of one. +newtype MonicPolynomial c e = MkMonicPolynomial { + getPolynomial :: Data.Polynomial.Polynomial c e +} deriving (Eq, Show) + +-- | Constructs an arbitrary /monic polynomial/. +mkMonicPolynomial :: (Num c, Ord e, Show e) => Data.Polynomial.Polynomial c e -> MonicPolynomial c e +mkMonicPolynomial polynomial + | not $ Data.Polynomial.isMonic polynomial = error $ "Factory.Data.MonicPolynomial.mkMonicPolynomial:\tnot monic; " ++ show polynomial + | otherwise = MkMonicPolynomial polynomial + +{- + * This instance-declaration merely delegates to the 'Data.Polynomial.Polynomial' payload. + + * CAVEAT: it's not strictly an instance of this class, since the result of some methods isn't /monic/. +-} +instance ( + Num c, + Num e, + Ord e, + Show e + ) => Data.Ring.Ring (MonicPolynomial c e) where + MkMonicPolynomial l =*= MkMonicPolynomial r = MkMonicPolynomial $ l =*= r + MkMonicPolynomial l =+= MkMonicPolynomial r = mkMonicPolynomial $ l =+= r --CAVEAT: potentially non-monic. +-- additiveInverse (MkMonicPolynomial p) = MkMonicPolynomial $ Data.Ring.additiveInverse p --CAVEAT: not monic ! + additiveInverse _ = error "Factory.Data.MonicPolynomial.additiveInverse:\tresult isn't monic" + multiplicativeIdentity = MkMonicPolynomial Data.Ring.multiplicativeIdentity + additiveIdentity = MkMonicPolynomial Data.Ring.additiveIdentity --CAVEAT: not monic ! + +-- Since the /leading term/ of the /denominator/ is one, the /coefficient/ isn't required to implement 'Fractional'. +instance (Num c, Num e, Ord e) => Data.QuotientRing.QuotientRing (MonicPolynomial c e) where + MkMonicPolynomial polynomialN `quotRem'` MkMonicPolynomial polynomialD = (MkMonicPolynomial *** MkMonicPolynomial) $ longDivide polynomialN where +-- longDivide :: (Num c, Num e, Ord e) => Polynomial c e -> (Polynomial c e, Polynomial c e) + longDivide numerator + | Data.Polynomial.isZero numerator || Data.Monomial.getExponent quotient < 0 = (Data.Polynomial.zero, numerator) + | otherwise = Control.Arrow.first (Data.Polynomial.lift (quotient :)) $ longDivide (numerator =-= polynomialD *= quotient) + where +-- quotient :: Num e => Data.Monomial.Monomial c e + quotient = Data.Polynomial.getLeadingTerm numerator `Data.Monomial.shiftExponent` negate (Data.Monomial.getExponent $ Data.Polynomial.getLeadingTerm polynomialD) + diff --git a/src/Factory/Data/Monomial.hs b/src/Factory/Data/Monomial.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Data/Monomial.hs @@ -0,0 +1,148 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] + + * Describes a and operations on it. + + * A /monomial/ is merely a /polynomial/ with a single non-zero term; cf. /Binomial/. +-} + +module Factory.Data.Monomial( +-- * Types +-- ** Type-synonyms + Monomial, +-- * Functions + double, + mod', + negateCoefficient, + realCoefficientToFrac, + shiftCoefficient, + shiftExponent, + square, +-- ** Accessors + getExponent, + getCoefficient, +-- ** Operators + (<=>), + (), + (<*>), + (=~), +-- ** Predicates + isMonomial +) where + +import qualified Control.Arrow +import qualified Factory.Math.Power as Math.Power + +infix 4 <=> --Same as (==). +infix 4 =~ --Same as (==). +infixl 7 --Same as (/). +infixl 7 <*> --Same as (*). + +{- | + * The type of an arbitrary monomial. + + * CAVEAT: though a /monomial/ has an integral power, this contraint is only imposed at the function-level. +-} +type Monomial coefficient exponent = (coefficient, exponent) + +-- | Accessor. +{-# INLINE getCoefficient #-} +getCoefficient :: Monomial c e -> c +getCoefficient = fst + +-- | Accessor. +{-# INLINE getExponent #-} +getExponent :: Monomial c e -> e +getExponent = snd + +{- | + * 'True' if the /exponent/ is both integral and non-/negative/. + + * CAVEAT: one can't even call this function unless the /exponent/ is integral. +-} +isMonomial :: Integral e => Monomial c e -> Bool +isMonomial = (>= 0) . getExponent + +-- | Compares the /exponents/ of the specified 'Monomial's. +{-# INLINE (<=>) #-} +(<=>) :: Ord e => Monomial c e -> Monomial c e -> Ordering +(_, l) <=> (_, r) = l `compare` r + +-- | 'True' if the /exponents/ are equal. +(=~) :: Eq e => Monomial c e -> Monomial c e -> Bool +(_, l) =~ (_, r) = l == r + +-- | Multiply the two specified 'Monomial's. +{-# INLINE (<*>) #-} +(<*>) :: (Num c, Num e) => Monomial c e -> Monomial c e -> Monomial c e +(cL, eL) <*> (cR, eR) = (cL * cR, eL + eR) + +-- | Divide the two specified 'Monomial's. +() :: (Fractional c, Num e) + => Monomial c e -- ^ Numerator. + -> Monomial c e -- ^ Denominator. + -> Monomial c e +(cN, eN) (1, eD) = (cN, eN - eD) +(cN, eN) (cD, eD) = (cN / cD, eN - eD) + +-- | Square the specified 'Monomial'. +square :: (Num c, Num e) => Monomial c e -> Monomial c e +square (c, e) = (Math.Power.square c, 2 * e) + +-- | Double the specified 'Monomial'. +{-# INLINE double #-} +double :: Num c => Monomial c e -> Monomial c e +double (c, e) = (2 * c, e) + +-- | Shift the /coefficient/, by the specified amount. +{-# INLINE shiftCoefficient #-} +shiftCoefficient :: Num c + => Monomial c e + -> c -- ^ The magnitude of the shift. + -> Monomial c e +--m `shiftCoefficient` i = Control.Arrow.first (+ i) m --CAVEAT: Too slow. +(c, e) `shiftCoefficient` i = (c + i, e) + +-- | Shift the /exponent/, by the specified amount. +{-# INLINE shiftExponent #-} +shiftExponent :: Num e + => Monomial c e + -> e -- ^ The magnitude of the shift. + -> Monomial c e +--m `shiftExponent` i = Control.Arrow.second (+ i) m --CAVEAT: Too slow. +(c, e) `shiftExponent` i = (c, e + i) + +-- | Negate the coefficient. +negateCoefficient :: Num c => Monomial c e -> Monomial c e +negateCoefficient = Control.Arrow.first negate + +-- | Reduce the coefficient using /modular/ arithmetic. +{-# INLINE mod' #-} +mod' :: Integral c + => Monomial c e + -> c -- ^ Modulus. + -> Monomial c e +monomial `mod'` modulus = Control.Arrow.first (`mod` modulus) monomial + +-- | Convert the type of the /coefficient/. +realCoefficientToFrac :: (Real r, Fractional f) => Monomial r e -> Monomial f e +realCoefficientToFrac = Control.Arrow.first realToFrac + diff --git a/src/Factory/Data/Polynomial.hs b/src/Factory/Data/Polynomial.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Data/Polynomial.hs @@ -0,0 +1,364 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] + + * Describes a polynomial and operations on it. + + * . + + * . +-} + +module Factory.Data.Polynomial( +-- * Types +-- ** Type-synonyms +-- MonomialList, +-- ** Data-types, + Polynomial, +-- * Constants + zero, + one, +-- * Functions + evaluate, + getDegree, + getLeadingTerm, + lift, + mod', + normalise, +-- pruneCoefficients, + raiseModulo, + realCoefficientsToFrac, + terms, +-- ** Constructors + mkConstant, + mkLinear, + mkPolynomial, +-- ** Operators + (*=), +-- ** Predicates + areCongruentModulo, + inAscendingOrder, + inDescendingOrder, +-- inOrder, + isMonic, + isMonomial, + isNormalised, + isPolynomial, +-- isReduced, + isZero +) where + +import Control.Arrow((&&&)) +import qualified Control.Arrow +import qualified Data.List +import Factory.Data.Monomial((<*>), (), (<=>), (=~)) +import qualified Factory.Data.Monomial as Data.Monomial +import qualified Factory.Data.QuotientRing as Data.QuotientRing +import Factory.Data.Ring((=*=), (=+=), (=-=)) +import qualified Factory.Data.Ring as Data.Ring + +infixl 7 *= --Same as (*). + +-- | The guts of a 'Polynomial'. +type MonomialList coefficient exponent = [Data.Monomial.Monomial coefficient exponent] + +{- | + * The type of an arbitrary /univariate/ polynomial; + actually it's more general, since it permits negative powers (s). + It can't describe /multivariate/ polynomials, which would require a list of /exponents/. + Rather than requiring the /exponent/ to implement the /type-class/ 'Integral', this is implemented at the function-level, as required. + + * The structure permits gaps between /exponents/, + in which /coefficients/ are inferred to be zero, thus enabling efficient representation of sparse polynomials. + + * CAVEAT: the 'MonomialList' is required to; + be ordered by /descending/ exponent (ie. reverse ); + have had zero coefficients removed; + and to have had /like/ terms merged; + so the raw data-constructor isn't exported. +-} +newtype {- Integral exponent => -} Polynomial coefficient exponent = MkPolynomial { + getMonomialList :: MonomialList coefficient exponent -- ^ Accessor. +} deriving (Eq, Show) + +-- | Makes /Polynomial/ a 'Data.Ring.Ring', over the /field/ composed from all possible /coefficients/; . +instance (Num c, Num e, Ord e) => Data.Ring.Ring (Polynomial c e) where + MkPolynomial [] =*= _ = zero + _ =*= MkPolynomial [] = zero + polynomialL =*= polynomialR +-- | polynomialL == one = polynomialR --Counterproductive. +-- | polynomialR == one = polynomialL --Counterproductive. + | terms polynomialL > terms polynomialR = polynomialL `times` polynomialR + | otherwise = polynomialR `times` polynomialL + where + l `times` r = {-# SCC "times" #-} Data.Ring.sum' (recip 2) {-TODO-} 10 {-empirical-} . map (l *=) $ getMonomialList r + + MkPolynomial [] =+= p = p + p =+= MkPolynomial [] = p + MkPolynomial listL =+= MkPolynomial listR = {-# SCC "merge" #-} MkPolynomial $ merge listL listR where + merge [] r = r + merge l [] = l + merge l@(lh : ls) r@(rh : rs) = case lh <=> rh of + GT -> lh : merge ls r + LT -> rh : merge l rs + _ -> case lh `Data.Monomial.shiftCoefficient` Data.Monomial.getCoefficient rh of + (0, _) -> merge ls rs + monomial -> monomial : merge ls rs + + additiveInverse = lift (Data.Monomial.negateCoefficient `map`) + multiplicativeIdentity = one + additiveIdentity = zero + +{- + Override the default implementation, + in order to take advantage of the symmetry under reflection about the main diagonal, + in the square matrix of products formed from the multiplication of each term by each term. + Eg: + (ax^3 + bx^2 + cx + d)^2 = [ + (a^2x^6 + abx^5 + acx^4 + adx^3) + + (bax^5 + b^2x^4 + bcx^3 + bdx^2) + + (cax^4 + cbx^3 + c^2x^2 + cdx) + + (dax^3 + dbx^2 + dcx + d^2) + ] + + = (a^2x^6 + b^2x^4 + c^2x^2 + d^2) + 2 * [ax^3 * (bx^2 + cx + d) + bx^2 * (cx + d) + cx * (d)] +-} + square (MkPolynomial []) = zero + square p = Data.Ring.sum' (recip 2) {-TODO-} 10 {-empirical-} $ diagonal : corners where + diagonal = {-# SCC "diagonal" #-} map Data.Monomial.square `lift` p + corners = {-# SCC "corners" #-} uncurry ( + zipWith (*=) + ) $ map MkPolynomial . init {-remove terminal null-} . Data.List.tails . tail &&& map Data.Monomial.double $ getMonomialList p + +-- | Defines the ability to divide /polynomials/. +instance (Fractional c, Num e, Ord e) => Data.QuotientRing.QuotientRing (Polynomial c e) where +{- + Uses /Euclidian division/. + . + . +-} + _ `quotRem'` MkPolynomial [] = error "Factory.Data.Polynomial.quotRem':\tzero denominator." + polynomialN `quotRem'` polynomialD = longDivide polynomialN where +-- longDivide :: (Fractional c, Num e, Ord e) => Polynomial c e -> (Polynomial c e, Polynomial c e) + longDivide (MkPolynomial []) = (zero, zero) --Exactly divides. + longDivide numerator + | Data.Monomial.getExponent quotient < 0 = (zero, numerator) --Indivisible remainder. + | otherwise = Control.Arrow.first (lift (quotient :)) $ longDivide (numerator =-= polynomialD *= quotient ) + where +-- quotient :: (Fractional c, Num e) => Data.Monomial.Monomial c e + quotient = getLeadingTerm numerator getLeadingTerm polynomialD + +{- | + * Transforms the data behind the constructor. + + * CAVEAT: similar to 'Data.Functor.fmap', but 'Polynomial' isn't an instance of 'Data.Functor.Functor' since we may want to operate on both /type-parameters/. + + * CAVEAT: the caller is required to re-'normalise' the resulting polynomial depending on the nature of the transformation of the data. +-} +lift :: (MonomialList c1 e1 -> MonomialList c2 e2) -> Polynomial c1 e1 -> Polynomial c2 e2 +lift transform = MkPolynomial . transform . getMonomialList + +-- | Returns the number of non-zero terms in the polynomial. +terms :: Polynomial c e -> Int +terms (MkPolynomial l) = length l + +-- | Return the highest-degree monomial. +getLeadingTerm :: Polynomial c e -> Data.Monomial.Monomial c e +getLeadingTerm (MkPolynomial []) = error "Factory.Data.Polynomial.getLeadingTerm:\tzero polynomial has no leading term." +getLeadingTerm (MkPolynomial (m : _)) = m + +-- | Removes terms with a /coefficient/ of zero. +pruneCoefficients :: Num c => Polynomial c e -> Polynomial c e +pruneCoefficients (MkPolynomial []) = zero +pruneCoefficients p = filter ((/= 0) . Data.Monomial.getCoefficient) `lift` p + +-- | Sorts into /descending order/ of exponents, groups /like/ exponents, and calls 'pruneCoefficients'. +normalise :: (Num c, Ord e) => Polynomial c e -> Polynomial c e +normalise = pruneCoefficients . lift ( + map ( + foldr ((+) . Data.Monomial.getCoefficient) 0 &&& Data.Monomial.getExponent . head + ) . Data.List.groupBy (=~) . Data.List.sortBy (flip (<=>)) + ) + +-- | Constructs an arbitrary /zeroeth-degree polynomial/, ie. independent of the /indeterminate/. +mkConstant :: (Num c, Num e) => c -> Polynomial c e +mkConstant 0 = zero +mkConstant c = MkPolynomial [(c, 0)] + +-- | Constructs an arbitrary /first-degree polynomial/. +mkLinear :: (Num c, Num e) + => c -- ^ Gradient. + -> c -- ^ Constant. + -> Polynomial c e +mkLinear m c = pruneCoefficients $ MkPolynomial [(m, 1), (c, 0)] + +-- | Constructs an arbitrary /polynomial/. +mkPolynomial :: (Num c, Ord e) => MonomialList c e -> Polynomial c e +mkPolynomial [] = zero +mkPolynomial l = normalise $ MkPolynomial l + +-- | Constructs a /polynomial/ with zero terms. +zero :: Polynomial c e +zero = MkPolynomial [] + +-- | Constructs a constant /monomial/, independent of the /indeterminate/. +one :: (Num c, Num e) => Polynomial c e +one = mkConstant 1 + +-- | 'True' if all /exponents/ are in the order defined by the specified comparator. +inOrder :: (e -> e -> Bool) -> Polynomial c e -> Bool +inOrder comparator p + | any ($ p) [isZero, isMonomial] = True + | otherwise = and . uncurry (zipWith comparator) . (init &&& tail) . map Data.Monomial.getExponent $ getMonomialList p + +-- | 'True' if the /exponents/ of successive terms are in /ascending/ order. +inAscendingOrder :: Ord e => Polynomial c e -> Bool +inAscendingOrder = inOrder (<=) + +-- | 'True' if the /exponents/ of successive terms are in /descending/ order. +inDescendingOrder :: Ord e => Polynomial c e -> Bool +inDescendingOrder = inOrder (>=) + +-- | 'True' if no term has a /coefficient/ of zero. +isReduced :: Num c => Polynomial c e -> Bool +isReduced = all ((/= 0) . Data.Monomial.getCoefficient) . getMonomialList + +-- | 'True' if no term has a /coefficient/ of zero and the /exponents/ of successive terms are in /descending/ order. +isNormalised :: (Num c, Ord e) => Polynomial c e -> Bool +isNormalised polynomial = all ($ polynomial) [isReduced, inDescendingOrder] + +{- | + * 'True' if the /leading coefficient/ is one. + + * . +-} +isMonic :: Num c => Polynomial c e -> Bool +isMonic (MkPolynomial []) = False --All coefficients are zero, and have therefore been removed. +isMonic p = (== 1) . Data.Monomial.getCoefficient $ getLeadingTerm p + +-- | 'True' if there are zero terms. +isZero :: Polynomial c e -> Bool +isZero (MkPolynomial []) = True +isZero _ = False + +-- | 'True' if there's exactly one term. +isMonomial :: Polynomial c e -> Bool +isMonomial (MkPolynomial []) = True +isMonomial _ = False + +-- | 'True' if all /exponents/ are /positive/ integers as required. +isPolynomial :: Integral e => Polynomial c e -> Bool +isPolynomial = all Data.Monomial.isMonomial . getMonomialList + +{- | + * 'True' if the two specified /polynomials/ are /congruent/ in /modulo/-arithmetic. + + * . +-} +areCongruentModulo :: (Integral c, Num e, Ord e) + => Polynomial c e -- ^ LHS. + -> Polynomial c e -- ^ RHS. + -> c -- ^ Modulus. + -> Bool +areCongruentModulo _ _ 0 = error "Factory.Data.Polynomial.areCongruentModulo:\tzero modulus." +areCongruentModulo _ _ 1 = True +areCongruentModulo l r modulus + | l == r = True + | otherwise = all ((== 0) . (`mod` modulus) . Data.Monomial.getCoefficient) . getMonomialList $ l =-= r + +{- | + * Return the /degree/ (AKA /order/) of the /polynomial/. + + * . + + * . +-} +getDegree :: Num e => Polynomial c e -> e +getDegree (MkPolynomial []) = -1 --CAVEAT: debatable, but makes some operations more robust and consistent. +getDegree p = Data.Monomial.getExponent $ getLeadingTerm p + +{- | + * Scale-up the specified /polynomial/ by a constant /monomial/ factor. + + * . +-} +(*=) :: (Num c, Num e) => Polynomial c e -> Data.Monomial.Monomial c e -> Polynomial c e +polynomial *= monomial + | Data.Monomial.getCoefficient monomial == 1 = map (`Data.Monomial.shiftExponent` Data.Monomial.getExponent monomial) `lift` polynomial + | otherwise = map (monomial <*>) `lift` polynomial + +{- | + * Raise a /polynomial/ to the specified positive integral power, but using /modulo/-arithmetic. + + * Whilst one could naively implement this as @(x Data.Ring.=^ n) `mod` m@, this will result in arithmetic operatons on unnecessarily big integers. +-} +raiseModulo :: (Integral c, Integral power, Num e, Ord e) + => Polynomial c e -- ^ The base. + -> power -- ^ The exponent to which the base should be raised. + -> c -- ^ The modulus. + -> Polynomial c e -- ^ The result. +raiseModulo _ _ 0 = error "Factory.Data.Polynomial.raiseModulo:\tzero modulus." +raiseModulo _ _ 1 = zero +raiseModulo _ 0 modulus = mkConstant $ 1 `mod` modulus +raiseModulo polynomial power modulus + | power < 0 = error $ "Factory.Data.Polynomial.raiseModulo:\tthe result isn't guaranteed to be a polynomial, for power=" ++ show power + | first `elem` [zero, one] = first --Eg 'raiseModulo (mkPolynomial [(3,1)]) 100 3' or 'raiseModulo (mkPolynomial [(3,1),(1,0)]) 100 3'. + | otherwise = slave power + where +-- first :: Integral c => Polynomial c e + first = polynomial `mod'` modulus + +-- slave :: (Integral c, Integral power, Num e, Ord e) => power -> Polynomial c e + slave 1 = first + slave n = (`mod'` modulus) . (if r == 0 {-even-} then id else (polynomial =*=)) . Data.Ring.square $ slave q {-recurse-} where + (q, r) = n `quotRem` 2 + +-- | Reduces all the coefficients using /modular/ arithmetic. +mod' :: Integral c + => Polynomial c e + -> c -- ^ Modulus. + -> Polynomial c e +mod' p modulus = pruneCoefficients $ map (`Data.Monomial.mod'` modulus) `lift` p + +{- | + * Evaluate the /polynomial/ at a specific /indeterminate/. + + * CAVEAT: requires positive exponents; but it wouldn't really be a /polynomial/ otherwise. + + * If the /polynomial/ is very sparse, this may be inefficient, + since it /memoizes/ the complete sequence of powers up to the polynomial's /degree/. +-} +evaluate :: (Num n, Integral e) + => n -- ^ The /indeterminate/. + -> Polynomial n e + -> n -- ^ The Result. +evaluate x = foldr ((+) . raise) 0 . getMonomialList where + powers = iterate (* x) 1 + + raise monomial + | exponent' < 0 = error $ "Factory.Data.Polynomial.evaluate.raise:\tnegative exponent; " ++ show exponent' + | otherwise = Data.Monomial.getCoefficient monomial * Data.List.genericIndex powers exponent' + where + exponent' = Data.Monomial.getExponent monomial + +-- | Convert the type of the /coefficient/s. +realCoefficientsToFrac :: (Real r, Fractional f) => Polynomial r e -> Polynomial f e +realCoefficientsToFrac = lift (Data.Monomial.realCoefficientToFrac `map`) + diff --git a/src/Factory/Data/PrimeFactors.hs b/src/Factory/Data/PrimeFactors.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Data/PrimeFactors.hs @@ -0,0 +1,153 @@ +{-# LANGUAGE CPP #-} +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] + + * Describes a list of /prime factors/. + + * The product of this list of prime-factors represents the /composite/ integer from which they were originally extracted. +-} + +module Factory.Data.PrimeFactors( +-- * Types +-- ** Type-synonyms + Factors, +-- * Functions + insert', +-- invert, +-- merge, + product', + reduce, +-- reduceSorted, +-- sumExponents, +-- ** Operators + (>*<), + (>/<), + (>^) +) where + +import qualified Control.Arrow +import Control.Arrow((&&&)) +import qualified Data.List +import qualified Data.Ord +import qualified Factory.Math.DivideAndConquer as Math.DivideAndConquer +import qualified Factory.Data.Exponential as Data.Exponential +import Factory.Data.Exponential((<^), (=~)) + +#if MIN_VERSION_toolshed(11,1,1) +import qualified ToolShed.ListPlus as ListPlus +#endif + +infixl 7 >/<, >*< --Same as (/). +infixr 8 >^ --Same as (^). + +{- | + * Each element of this list represents one /prime-factor/, expressed as an /exponential/ with a /prime/ base, of the original integer. + + * Whilst it only makes sense for both the /base/ and /exponent/ to be integral, these constrains are applied at the function-level as required. +-} +type Factors base exponent = [Data.Exponential.Exponential base exponent] + +{- | + * Sorts a list representing a product of /prime factors/ by increasing /base/. + + * Multiplies 'Data.Exponential.Exponential's of similar /base/. +-} +reduce :: (Ord base, Num exponent, Ord exponent) => Factors base exponent -> Factors base exponent +reduce = reduceSorted . Data.List.sort {-primarily by base-} + +-- | Multiplies 'Data.Exponential.Exponential's of similar /base/. +reduceSorted :: (Eq base, Num exponent) => Factors base exponent -> Factors base exponent +--reduceSorted = map (Data.Exponential.getBase . head &&& sumExponents) . Data.List.groupBy (=~) --Slow +reduceSorted [] = [] +reduceSorted (x : xs) + | null matched = x : reduceSorted remainder + | otherwise = Control.Arrow.second (+ sumExponents matched) x : reduceSorted remainder + where + (matched, remainder) = span (=~ x) xs + +{- | + * Insert a 'Data.Exponential.Exponential', into a list representing a product of /prime factors/, multiplying with any incumbent of like /base/. + + * The list should be sorted by increasing /base/. + + * Preserves the sort-order. + + * CAVEAT: this is tolerably efficient for the odd insertion; to insert a list, use '>*<'. +-} +insert' :: (Ord base, Num exponent) => Data.Exponential.Exponential base exponent -> Factors base exponent -> Factors base exponent +insert' e [] = [e] +insert' e l@(x : xs) = case Data.Ord.comparing Data.Exponential.getBase e x of + LT -> e : l + GT -> x : insert' e xs --Recurse. + _ -> Control.Arrow.second (+ Data.Exponential.getExponent e) x : xs --Multiply by adding exponents. + +{- | + * Multiplies two lists each representing a product of /prime factors/, and sorted by increasing /base/. + + * Preserves the sort-order. +-} +(>*<) :: (Ord base, Num exponent, Ord exponent) => Factors base exponent -> Factors base exponent -> Factors base exponent +l >*< r = +#if MIN_VERSION_toolshed(11,1,1) + reduceSorted $ ListPlus.merge l r +#else + reduce $ l ++ r --CAVEAT: concatenation disorders the list, necessitating a re-sort. +#endif + +-- | Invert the product of a list /prime factors/, by negating each of the /exponents/. +invert :: Num exponent => Factors base exponent -> Factors base exponent +invert = map Data.Exponential.invert + +{- | + * Divides two lists, each representing a product of /prime factors/, and sorted by increasing /base/. + + * Preserves the sort-order. +-} +(>/<) :: (Integral base, Integral exponent) + => Factors base exponent -- ^ The list of /prime factors/ in the /numerator/. + -> Factors base exponent -- ^ The list of /prime factors/ in the /denominator/. + -> (Factors base exponent, Factors base exponent) -- ^ The ratio of /numerator/ and /denominator/, after like /prime factors/ are cancelled. +numerator >/< denominator = filter ( + (> 0) . Data.Exponential.getExponent + ) &&& invert . filter ( + (< 0) . Data.Exponential.getExponent + ) $ numerator >*< invert denominator + +{- | + * Raise the product of a list /prime factors/ to the specified power. + + * CAVEAT: this merely involves raising each element to the specified power; cf. raising a /polynomial/ to a power. +-} +(>^) :: Num exponent => Factors base exponent -> exponent -> Factors base exponent +factors >^ power = map (<^ power) factors + +-- | Sum the /exponents/ of the specified list; as required to multiply exponentials with identical /base/. +sumExponents :: Num exponent => Factors base exponent -> exponent +sumExponents = foldr ((+) . Data.Exponential.getExponent) 0 + +-- | Multiply a list of /prime factors/. +product' :: (Num base, Integral exponent) + => Math.DivideAndConquer.BisectionRatio + -> Math.DivideAndConquer.MinLength + -> Factors base exponent -- ^ The list on which to operate. + -> base -- ^ The result. +product' bisectionRatio minLength = Math.DivideAndConquer.product' bisectionRatio minLength . map Data.Exponential.evaluate + diff --git a/src/Factory/Data/QuotientRing.hs b/src/Factory/Data/QuotientRing.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Data/QuotientRing.hs @@ -0,0 +1,79 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] + + * Describes a /Quotient Ring/; . + + * This is a /ring/ composed from a residue-class resulting from /modular/ division. +-} + +module Factory.Data.QuotientRing ( +-- * Type-classes + QuotientRing(..), +-- * Functions + quot', + rem', +-- ** Predicates + areCongruentModulo, + isDivisibleBy +) where + +import Factory.Data.Ring((=-=)) +import qualified Factory.Data.Ring as Data.Ring + +-- | Defines a sub-class of 'Data.Ring.Ring', in which division is implemented. +class Data.Ring.Ring q => QuotientRing q where + quotRem' :: q -> q -> (q, q) -- ^ Divides the first operand by the second, to yield a pair composed from the /quotient/ and the /remainder/. + +-- | Returns the /quotient/, after division of the two specified 'QuotientRing's. +quot' :: QuotientRing q + => q -- ^ Numerator. + -> q -- ^ Denominator. + -> q +quot' numerator = fst . quotRem' numerator + +-- | Returns the /remainder/, after division of the two specified 'QuotientRing's. +rem' :: QuotientRing q + => q -- ^ Numerator. + -> q -- ^ Denominator. + -> q +rem' numerator = snd . quotRem' numerator + +{- | + * 'True' if the two specified 'QuotientRing's are /congruent/ in /modulo/-arithmetic, where the /modulus/ is a third 'QuotientRing'. + + * . +-} +areCongruentModulo :: (Eq q, QuotientRing q) + => q -- ^ LHS. + -> q -- ^ RHS. + -> q -- ^ Modulus. + -> Bool +areCongruentModulo l r modulus + | l == r = True --Only required for efficiency. + | otherwise = (l =-= r) `isDivisibleBy` modulus + +-- | 'True' if the second operand /divides/ the first. +isDivisibleBy :: (Eq q, QuotientRing q) + => q -- ^ Numerator. + -> q -- ^ Denominator. + -> Bool +numerator `isDivisibleBy` denominator = rem' numerator denominator == Data.Ring.additiveIdentity + diff --git a/src/Factory/Data/Ring.hs b/src/Factory/Data/Ring.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Data/Ring.hs @@ -0,0 +1,113 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] + + * Describes a /ring/ and operations on its members. + + * . + + * . +-} + +module Factory.Data.Ring( +-- * Type-classes + Ring(..), +-- * Types +-- ** Data.types +-- Product, +-- Sum, +-- * Functions + product', + sum', +-- ** Operators + (=^) +) where + +import qualified Data.Monoid +import qualified Factory.Math.DivideAndConquer as Math.DivideAndConquer + +infixl 6 =+= --Same as (+). +infixl 6 =-= --Same as (-). +infixl 7 =*= --Same as (*). +infixr 8 =^ --Same as (^). + +{- | + * Define both the operations applicable to all members of the /ring/, and its mandatory members. + + * Minimal definition; '=+=', '=*=', 'additiveInverse', 'multiplicativeIdentity', 'additiveIdentity'. +-} +class Ring r where + (=+=) :: r -> r -> r -- ^ Addition of two members; required to be /commutative/; . + (=*=) :: r -> r -> r -- ^ Multiplication of two members. + additiveInverse :: r -> r -- ^ The operand required to yield /zero/ under addition; . + multiplicativeIdentity :: r -- ^ The /identity/-member under multiplication; . + additiveIdentity :: r -- ^ The /identity/-member under addition (AKA /zero/); . + + (=-=) :: r -> r -> r -- ^ Subtract the two specified /ring/-members. + l =-= r = l =+= additiveInverse r --Default implementation. + + square :: r -> r -- ^ Square the ring. + square r = r =*= r --Default implementation; there may be a more efficient one. + +{- | + * Raise a /ring/-member to the specified positive integral power. + + * Exponentiation is implemented as a sequence of either squares of, or multiplications by, the /ring/-member; + . +-} +(=^) :: (Ring r, Eq r, Integral power) => r -> power -> r +_ =^ 0 = multiplicativeIdentity +ring =^ power + | power < 0 = error $ "Factory.Data.Ring.(=^):\tthe result isn't guaranteed to be a ring-member, for power=" ++ show power + | ring `elem` [additiveIdentity, multiplicativeIdentity] = ring + | otherwise = slave power + where + slave 1 = ring + slave n = (if r == 0 {-even-} then id else (=*= ring)) . square $ slave q where + (q, r) = n `quotRem` 2 + +-- | Does for 'Ring', what 'Data.Monoid.Product' does for type 'Num', in that it makes it an instance of 'Data.Monoid.Monoid' under multiplication. +newtype Product p = MkProduct { + getProduct :: p -- ^ Access the polymorphic payload. +} deriving (Read, Show) + +instance Ring r => Data.Monoid.Monoid (Product r) where + mempty = MkProduct multiplicativeIdentity + MkProduct x `mappend` MkProduct y = MkProduct $ x =*= y + +-- | Returns the /product/ of the list of /ring/-members. +product' :: Ring r => Math.DivideAndConquer.BisectionRatio -> Math.DivideAndConquer.MinLength -> [r] -> r +--product' _ _ = getProduct . Data.Monoid.mconcat . map MkProduct +product' ratio minLength = getProduct . Math.DivideAndConquer.divideAndConquer ratio minLength . map MkProduct + +-- | Does for 'Ring', what 'Data.Monoid.Sum' does for type 'Num', in that it makes it an instance of 'Data.Monoid.Monoid' under addition. +newtype Sum s = MkSum { + getSum :: s -- ^ Access the polymorphic payload. +} deriving (Read, Show) + +instance Ring r => Data.Monoid.Monoid (Sum r) where + mempty = MkSum additiveIdentity + MkSum x `mappend` MkSum y = MkSum $ x =+= y + +-- | Returns the /sum/ of the list of /ring/-members. +sum' :: Ring r => Math.DivideAndConquer.BisectionRatio -> Math.DivideAndConquer.MinLength -> [r] -> r +--sum' _ _ = getSum . Data.Monoid.mconcat . map MkSum +sum' ratio minLength = getSum . Math.DivideAndConquer.divideAndConquer ratio minLength . map MkSum + diff --git a/src/Factory/Math/ArithmeticGeometricMean.hs b/src/Factory/Math/ArithmeticGeometricMean.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/ArithmeticGeometricMean.hs @@ -0,0 +1,99 @@ +{-# LANGUAGE CPP #-} +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Determines the /Arithmetic-geometric mean/; . +-} + +module Factory.Math.ArithmeticGeometricMean( +-- * Types +-- ** Type-synonyms + ArithmeticMean, + GeometricMean, + AGM, +-- * Functions + convergeToAGM, + spread, +-- ** Accessors + getArithmeticMean, + getGeometricMean, +-- ** Predicates + isValid +) where + +import Control.Arrow((&&&)) +import qualified Data.Ratio +import qualified Factory.Math.Precision as Math.Precision +import qualified Factory.Math.SquareRoot as Math.SquareRoot + +#if MIN_VERSION_parallel(3,0,0) +import qualified Control.Parallel.Strategies +#endif + +-- | The type of the /arithmetic mean/; . +type ArithmeticMean = Data.Ratio.Rational + +-- | The type of the /geometric mean/; . +type GeometricMean = Data.Ratio.Rational + +-- | Encapsulates both /arithmetic/ and /geometric/ means. +type AGM = (ArithmeticMean, GeometricMean) + +-- | Accessor. +{-# INLINE getArithmeticMean #-} +getArithmeticMean :: AGM -> ArithmeticMean +getArithmeticMean = fst + +-- | Accessor. +{-# INLINE getGeometricMean #-} +getGeometricMean :: AGM -> GeometricMean +getGeometricMean = snd + +-- | Returns an infinite list which converges on the /Arithmetic-geometric mean/. +convergeToAGM :: Math.SquareRoot.Algorithm squareRootAlgorithm => squareRootAlgorithm -> Math.Precision.DecimalDigits -> AGM -> [AGM] +convergeToAGM squareRootAlgorithm decimalDigits agm + | decimalDigits <= 0 = error $ "Factory.Math.ArithmeticGeometricMean.convergeToAGM:\tinvalid number of decimal digits; " ++ show decimalDigits + | not $ isValid agm = error $ "Factory.Math.ArithmeticGeometricMean.convergeToAGM:\tboth means must be positive for a real geometric mean; " ++ show agm + | spread agm == 0 = repeat agm + | otherwise = let + simplify :: Data.Ratio.Rational -> Data.Ratio.Rational + simplify = Math.Precision.simplify (decimalDigits - 1 {-ignore single integral digit-}) --This makes a gigantic difference to performance. + + findArithmeticMean :: AGM -> ArithmeticMean + findArithmeticMean = (/ 2) . uncurry (+) + + findGeometricMean :: AGM -> GeometricMean + findGeometricMean = Math.SquareRoot.squareRoot squareRootAlgorithm decimalDigits . uncurry (*) + in iterate ( +#if MIN_VERSION_parallel(3,0,0) + Control.Parallel.Strategies.withStrategy ( + Control.Parallel.Strategies.parTuple2 Control.Parallel.Strategies.rdeepseq Control.Parallel.Strategies.rdeepseq + ) . +#endif + (simplify . findArithmeticMean &&& simplify . findGeometricMean) + ) agm + +-- | Returns the bounds within which the 'AGM' has been constrained. +spread :: AGM -> Data.Ratio.Rational +spread = uncurry (-) + +-- | Checks that both /means/ are positive, as required for the /geometric mean/ to be consistently /real/. +isValid :: AGM -> Bool +isValid (a, g) = all (>= 0) [a, g] + diff --git a/src/Factory/Math/DivideAndConquer.hs b/src/Factory/Math/DivideAndConquer.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/DivideAndConquer.hs @@ -0,0 +1,130 @@ +{-# LANGUAGE CPP #-} +{- + Copyright (C) 2010 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] + + * Provides a polymorphic algorithm, to /unfold/ a list into a tree, to which an /associative binary operator/ is then applied to re-/fold/ the tree to a /scalar/. + + * Implementations of this strategy have been provided for /addition/ and /multiplication/, + though other associative binary operators, like 'gcd' or 'lcm' could also be used. + + * Where the contents of the list are consecutive, a more efficient implementation is available in /Factory.Data.Bounds/. +-} + +module Factory.Math.DivideAndConquer( +-- * Types +-- ** Type-synonyms + BisectionRatio, + MinLength, +-- * Functions + divideAndConquer, + product', + sum' +) where + +import Control.Arrow((***)) +import qualified Data.Monoid +import qualified Data.Ratio + +#if MIN_VERSION_parallel(3,0,0) +import qualified Control.Parallel.Strategies +#endif + +{- | + * The ratio of the original list-length at which to bisect. + + * CAVEAT: the value can overflow. +-} +type BisectionRatio = Data.Ratio.Ratio Int + +-- | The list-length beneath which to terminate bisection. +type MinLength = Int + +{- | + * Reduces a list to a single scalar encapsulated in a 'Data.Monoid.Monoid', + using a /divide-and-conquer/ strategy, + bisecting the list and recursively evaluating each part; . + + * By choosing a 'bisectionRatio' other than @(1 % 2)@, the bisection can be made asymmetrical. + The specified ratio represents the length of the left-hand portion, over the original list-length; + eg. @(1 % 3)@ results in the first part, half the length of the second. + + * This process of recursive bisection, is terminated beneath the specified minimum list-length, + after which the /monoid/'s binary operator is directly /folded/ over the list. + + * One can view this as a , + in which the list is exploded into a binary tree-structure + (each leaf of which contains a list of up to 'minLength' integers, and each node of which contains an associative binary operator), + and then collapsed to a scalar, by application of the operators. +-} +divideAndConquer :: Data.Monoid.Monoid monoid + => BisectionRatio -- ^ The ratio of the original list-length at which to bisect. + -> MinLength -- ^ For efficiency, the list will not be bisected, when it's length has been reduced to this value. + -> [monoid] -- ^ The list on which to operate. + -> monoid -- ^ The resulting scalar. +divideAndConquer bisectionRatio minLength l + | any ($ apportion minLength) [ + (< 1), --The left-hand list may be null. + (> pred minLength) --The right-hand list may be null. + ] = error $ "Factory.Math.DivideAndConquer.divideAndConquer:\tbisectionRatio='" ++ show bisectionRatio ++ "' is incompatible with minLength=" ++ show minLength ++ "." + | otherwise = slave (length l) l + where + apportion :: Int -> Int + apportion list = (list * Data.Ratio.numerator bisectionRatio) `div` Data.Ratio.denominator bisectionRatio + + slave len list + | len <= minLength = Data.Monoid.mconcat list --Fold the monoid's binary operator over the list. + | otherwise = uncurry Data.Monoid.mappend . +#if MIN_VERSION_parallel(3,0,0) + Control.Parallel.Strategies.withStrategy ( + Control.Parallel.Strategies.parTuple2 Control.Parallel.Strategies.rseq Control.Parallel.Strategies.rseq + ) . +#endif + (slave cut *** slave (len - cut)) $ splitAt cut list where --Apply the monoid's binary operator to the two operands resulting from bisection. + cut = apportion len + +{- | + * Multiplies the specified list of numbers. + + * Since the result can be large, 'divideAndConquer' is used in an attempt to form operands of a similar order of magnitude, + which creates scope for the use of more efficient multiplication-algorithms. +-} +product' :: Num n + => BisectionRatio -- ^ The ratio of the original list-length at which to bisect. + -> MinLength -- ^ For efficiency, the list will not be bisected, when it's length has been reduced to this value. + -> [n] -- ^ The numbers whose product is required. + -> n -- ^ The resulting product. +product' bisectionRatio minLength = Data.Monoid.getProduct . divideAndConquer bisectionRatio minLength . map Data.Monoid.Product + +{- | + * Sums the specified list of numbers. + + * Since the result can be large, 'divideAndConquer' is used in an attempt to form operands of a similar order of magnitude, + which creates scope for the use of more efficient multiplication-algorithms. + /Multiplication/ is required for the /addition/ of 'Data.Ratio.Rational' numbers by cross-multiplication; + this function is unlikely to be useful for other numbers. +-} +sum' :: Num n + => BisectionRatio -- ^ The ratio of the original list-length at which to bisect. + -> MinLength -- ^ For efficiency, the list will not be bisected, when it's length has been reduced to this value. + -> [n] -- ^ The numbers whose sum is required. + -> n -- ^ The resulting sum. +sum' bisectionRatio minLength = Data.Monoid.getSum . divideAndConquer bisectionRatio minLength . map Data.Monoid.Sum + diff --git a/src/Factory/Math/Factorial.hs b/src/Factory/Math/Factorial.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Factorial.hs @@ -0,0 +1,37 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] + + * Whilst this particular function is the subject of many introductory examples to Haskell, + the simple algorithms appropriate for that forum, leave a large margin for performance-improvement. + This module provides the interface for alternative algorithms. + + * . +-} + +module Factory.Math.Factorial( +-- * Type-classes + Algorithm(..) +) where + +-- | Defines the methods expected of a /factorial/-algorithm. +class Algorithm algorithm where + factorial :: Integral i => algorithm -> i -> i + diff --git a/src/Factory/Math/Fibonacci.hs b/src/Factory/Math/Fibonacci.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Fibonacci.hs @@ -0,0 +1,42 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] . +-} + +module Factory.Math.Fibonacci( +-- * Constants + fibonacci, + primeIndexedFibonacci +) where + +import qualified Data.Numbers.Primes + +-- | A constant ordered list of the /Fibonacci/-numbers. +fibonacci :: Integral i => [i] +fibonacci = 0 : scanl (+) 1 fibonacci + +{- | + * The subset of 'fibonacci', /indexed/ by a /prime/-number. + + * . +-} +primeIndexedFibonacci :: Integral i => [i] +primeIndexedFibonacci = map (fibonacci !!) Data.Numbers.Primes.primes + diff --git a/src/Factory/Math/Implementations/Factorial.hs b/src/Factory/Math/Implementations/Factorial.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Factorial.hs @@ -0,0 +1,138 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] + + * Provides implementations of the class 'Math.Factorial.Algorithm'. + + * Provides additional functions related to /factorials/, but which depends on a specific implementation, + and which therefore can't be accessed throught the class-interface. + + * . + + * . + + * . +-} + +module Factory.Math.Implementations.Factorial( +-- * Types +-- ** Data-types + Algorithm(..), +-- * Functions + primeFactors, +-- primeMultiplicity, + risingFactorial, + fallingFactorial, +-- ** Operators + (!/!) +) where + +import qualified Data.Numbers.Primes +import qualified Factory.Data.Bounds as Data.Bounds +import qualified Factory.Data.PrimeFactors as Data.PrimeFactors +import qualified Factory.Math.Factorial as Math.Factorial +import qualified ToolShed.Defaultable as Defaultable + +infixl 7 !/! --Same as (/). + +-- | The algorithms by which /factorial/ has been implemented. +data Algorithm = + Bisection -- ^ The integers from which the /factorial/ is composed, are multiplied using @Data.Bounds.product'@. + | PrimeFactorisation -- ^ The /prime factors/ of the /factorial/ are extracted, then raised to the appropriate power, before multiplication. + deriving (Eq, Read, Show) + +instance Defaultable.Defaultable Algorithm where + defaultValue = Bisection + +instance Math.Factorial.Algorithm Algorithm where + factorial algorithm n + | n < 2 = 1 + | otherwise = case algorithm of + Bisection -> risingFactorial 2 $ n - 1 + PrimeFactorisation -> Data.PrimeFactors.product' (recip 5) {-empirical-} 10 {-empirical-} $ primeFactors n + +{- | + * Returns the /prime factors/, of the /factorial/ of the specifed integer. + + * Precisely all the primes less than or equal to the specified integer /n/, are included in /n!/; + only the multiplicity of each of these known prime components need be determined. + + * . + + * CAVEAT: currently a hotspot. +-} +primeFactors :: Integral base + => base -- ^ The number, whose /factorial/ is to be factorised. + -> Data.PrimeFactors.Factors base base -- ^ The /base/ and /exponent/ of each /prime factor/ in the /factorial/, ordered by increasing /base/ (and decreasing /exponent/). +primeFactors n = takeWhile ((> 0) . snd) $ map (\prime -> (prime, primeMultiplicity prime n)) Data.Numbers.Primes.primes + +{- | + * The number of times a specific /prime/, can be factored from the /factorial/ of the specified integer. + + * General purpose /prime-factorisation/ has /exponential time-complexity/, + so use /Legendre's Theorem/, which relates only to the /prime factors/ of /factorials/. + + * . +-} +primeMultiplicity :: Integral i + => i -- ^ A prime number. + -> i -- ^ The integer, the factorial of which the prime is a factor. + -> i -- ^ The number of times the prime occurs in the factorial. +primeMultiplicity prime = sum . takeWhile (> 0) . tail . iterate (`div` prime) + +-- | Returns the /rising factorial/; +risingFactorial :: Integral i + => i -- ^ The lower bound of the integer-range, whose product is returned. + -> i -- ^ The number of integers in the range above. + -> i -- ^ The result. +risingFactorial _ 0 = 1 +risingFactorial 0 _ = 0 +risingFactorial x n = Data.Bounds.product' (recip 2) 64 $ Data.Bounds.normalise (x, (x + n) - 1) + +-- | Returns the /falling factorial/; +fallingFactorial :: Integral i + => i -- ^ The upper bound of the integer-range, whose product is returned. + -> i -- ^ The number of integers in the range beneath. + -> i -- ^ The result. +fallingFactorial _ 0 = 1 +fallingFactorial 0 _ = 0 +fallingFactorial x n = Data.Bounds.product' (recip 2) 64 $ Data.Bounds.normalise (x, (x - n) + 1) + +{- | + * Returns the ratio of two factorials. + + * It is more efficient than evaluating both factorials, and then dividing. + + * For more complex combinations of factorials, such as in the /Binomial coefficient/, + extract the /prime factors/ using 'primeFactors' + then manipulate them using the module "Data.PrimeFactors", + and evaluate it using by /Data.PrimeFactors.product'/. +-} +(!/!) :: (Integral i, Fractional f) + => i -- ^ The /numerator/. + -> i -- ^ The /denominator/. + -> f -- ^ The resulting fraction. +numerator !/! denominator + | numerator <= 1 = recip . fromIntegral $ Math.Factorial.factorial (Defaultable.defaultValue :: Algorithm) denominator + | denominator <= 1 = fromIntegral $ Math.Factorial.factorial (Defaultable.defaultValue :: Algorithm) numerator + | numerator == denominator = 1 + | numerator < denominator = recip $ denominator !/! numerator --Recurse. + | otherwise = fromIntegral $ Data.Bounds.product' (recip 2) 64 (succ denominator, numerator) + diff --git a/src/Factory/Math/Implementations/Pi/AGM/Algorithm.hs b/src/Factory/Math/Implementations/Pi/AGM/Algorithm.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Pi/AGM/Algorithm.hs @@ -0,0 +1,42 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines the set of /Arithmetic-geometric Mean/-type /Pi/-algorithms which have been implemented; currently just one. +-} + +module Factory.Math.Implementations.Pi.AGM.Algorithm( +-- * Types +-- ** Data-types + Algorithm(..) +) where + +import qualified Factory.Math.Implementations.Pi.AGM.BrentSalamin as Math.Implementations.Pi.AGM.BrentSalamin +import qualified Factory.Math.Pi as Math.Pi +import qualified Factory.Math.SquareRoot as Math.SquareRoot +import qualified ToolShed.Defaultable as Defaultable + +-- | Defines the available algorithms. +data Algorithm squareRootAlgorithm = BrentSalamin squareRootAlgorithm deriving (Eq, Read, Show) + +instance Defaultable.Defaultable squareRootAlgorithm => Defaultable.Defaultable (Algorithm squareRootAlgorithm) where + defaultValue = BrentSalamin Defaultable.defaultValue + +instance Math.SquareRoot.Algorithm squareRootAlgorithm => Math.Pi.Algorithm (Algorithm squareRootAlgorithm) where + openR (BrentSalamin squareRootAlgorithm) = Math.Implementations.Pi.AGM.BrentSalamin.openR squareRootAlgorithm + diff --git a/src/Factory/Math/Implementations/Pi/AGM/BrentSalamin.hs b/src/Factory/Math/Implementations/Pi/AGM/BrentSalamin.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Pi/AGM/BrentSalamin.hs @@ -0,0 +1,65 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] + + * Implements the /Brent-Salamin/ (AKA /Gauss-Legendre/) algorithm; + , + , + . + + * The precision of the result approximately doubles for each iteration. + + [@CAVEAT@] Assumptions on the convergence-rate result in rounding-errors, when only a small number of digits are requested. +-} + +module Factory.Math.Implementations.Pi.AGM.BrentSalamin( +-- * Functions + openR +) where + +import Control.Arrow((&&&)) +import qualified Data.Ratio +import qualified Factory.Math.ArithmeticGeometricMean as Math.ArithmeticGeometricMean +import qualified Factory.Math.Power as Math.Power +import qualified Factory.Math.Precision as Math.Precision +import qualified Factory.Math.SquareRoot as Math.SquareRoot + +{- | + * Returns /Pi/, accurate to the specified number of decimal digits. + + * This algorithm is based on the /arithmetic-geometric/ mean of @1@ and @(1 / sqrt 2)@, + but there are many confusingly similar formulations. + The algorithm I've used here, where @a@ is the /arithmetic mean/ and @g@ is the /geometric mean/, is equivalent to other common formulations: + +> pi = (a[N-1] + g[N-1])^2 / (1 - sum [2^n * (a[n] - g[n])^2]) where n = [0 .. N-1] +> => 4*a[N]^2 / (1 - sum [2^n * (a[n]^2 - 2*a[n]*g[n] + g[n]^2)]) +> => 4*a[N]^2 / (1 - sum [2^n * (a[n]^2 + 2*a[n]*g[n] + g[n]^2 - 4*a[n]*g[n])]) +> => 4*a[N]^2 / (1 - sum [2^n * ((a[n] + g[n])^2 - 4*a[n]*g[n])]) +> => 4*a[N]^2 / (1 - sum [2^(n-1) * 4 * (a[n-1]^2 - g[n-1]^2)]) where n = [1 .. N] +> => 4*a[N]^2 / (1 - sum [2^(n+1) * (a[n-1]^2 - g[n-1]^2)]) + +-} +openR :: Math.SquareRoot.Algorithm squareRootAlgorithm => squareRootAlgorithm -> Math.Precision.DecimalDigits -> Data.Ratio.Rational +openR squareRootAlgorithm decimalDigits = uncurry (/) . ( + Math.Power.square . uncurry (+) . last &&& negate . pred . sum . zipWith (*) (iterate (* 2) 1) . map (Math.Power.square . Math.ArithmeticGeometricMean.spread) + ) . take ( + Math.Precision.getIterationsRequired Math.Precision.quadraticConvergence 1 decimalDigits + ) $ Math.ArithmeticGeometricMean.convergeToAGM squareRootAlgorithm decimalDigits (1, Math.SquareRoot.squareRoot squareRootAlgorithm decimalDigits (recip 2 :: Data.Ratio.Rational)) + diff --git a/src/Factory/Math/Implementations/Pi/BBP/Algorithm.hs b/src/Factory/Math/Implementations/Pi/BBP/Algorithm.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Pi/BBP/Algorithm.hs @@ -0,0 +1,47 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines the set of /Bailey-Borwein-Plouffe/-type formulae which have been implemented. +-} + +module Factory.Math.Implementations.Pi.BBP.Algorithm( +-- * Types +-- ** Data-types + Algorithm(..) +) where + +import qualified Factory.Math.Implementations.Pi.BBP.Base65536 as Math.Implementations.Pi.BBP.Base65536 +import qualified Factory.Math.Implementations.Pi.BBP.Bellard as Math.Implementations.Pi.BBP.Bellard +import qualified Factory.Math.Implementations.Pi.BBP.Implementation as Math.Implementations.Pi.BBP.Implementation +import qualified Factory.Math.Pi as Math.Pi +import qualified ToolShed.Defaultable as Defaultable + +-- | Defines those /BBP/-type series which have been implemented. +data Algorithm = + Base65536 -- ^ A /base/-@2^16@ version of the formula. + | Bellard -- ^ A /nega-base/ @2^10@ version of the formula. + deriving (Eq, Read, Show) + +instance Defaultable.Defaultable Algorithm where + defaultValue = Base65536 + +instance Math.Pi.Algorithm Algorithm where + openR Base65536 = Math.Implementations.Pi.BBP.Implementation.openR Math.Implementations.Pi.BBP.Base65536.series + openR Bellard = Math.Implementations.Pi.BBP.Implementation.openR Math.Implementations.Pi.BBP.Bellard.series + diff --git a/src/Factory/Math/Implementations/Pi/BBP/Base65536.hs b/src/Factory/Math/Implementations/Pi/BBP/Base65536.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Pi/BBP/Base65536.hs @@ -0,0 +1,38 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines a specific base-@2^16@ /BBP/-formula; + +-} + +module Factory.Math.Implementations.Pi.BBP.Base65536( +-- * Constants + series +) where + +import qualified Factory.Math.Implementations.Pi.BBP.Series as Math.Implementations.Pi.BBP.Series + +-- | Defines the parameters of this specific series. +series :: Math.Implementations.Pi.BBP.Series.Series +series = Math.Implementations.Pi.BBP.Series.MkSeries { + Math.Implementations.Pi.BBP.Series.numerators = zipWith ($) (concat $ repeat [id, id, id, negate]) $ map (2 ^) [15 :: Integer, 14, 14, 12, 11, 10, 10, 8, 7, 6, 6, 4, 3, 2, 2, 0], + Math.Implementations.Pi.BBP.Series.getDenominators = \i -> map (32 * fromIntegral i +) [2, 3, 4, 7, 10, 11, 12, 15, 18, 19, 20, 23, 26, 27, 28, 31], + Math.Implementations.Pi.BBP.Series.seriesScalingFactor = recip $ 2 ^ (13 :: Int), + Math.Implementations.Pi.BBP.Series.base = 2 ^ (16 :: Int) +} diff --git a/src/Factory/Math/Implementations/Pi/BBP/Bellard.hs b/src/Factory/Math/Implementations/Pi/BBP/Bellard.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Pi/BBP/Bellard.hs @@ -0,0 +1,41 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines /Bellard/'s nega-base-@2^10@ /BBP/-formula; +-} + +module Factory.Math.Implementations.Pi.BBP.Bellard( +-- * Constants + series +) where + +import Control.Arrow((&&&)) +import qualified Factory.Math.Implementations.Pi.BBP.Series as Math.Implementations.Pi.BBP.Series + +-- | Defines the parameters of this specific series. +series :: Math.Implementations.Pi.BBP.Series.Series +series = Math.Implementations.Pi.BBP.Series.MkSeries { + Math.Implementations.Pi.BBP.Series.numerators = zipWith ($) [negate, negate, id, negate, negate, negate, id] $ map (2 ^) [5 :: Integer, 0, 8, 6, 2, 2, 0], + Math.Implementations.Pi.BBP.Series.getDenominators = \i -> let + f, t :: Integer + (f, t) = (4 *) &&& (10 *) $ fromIntegral i + in [f + 1, f + 3, t + 1, t + 3, t + 5, t + 7, t + 9], + Math.Implementations.Pi.BBP.Series.seriesScalingFactor = recip $ 2 ^ (6 :: Int), + Math.Implementations.Pi.BBP.Series.base = negate $ 2 ^ (10 :: Int) +} diff --git a/src/Factory/Math/Implementations/Pi/BBP/Implementation.hs b/src/Factory/Math/Implementations/Pi/BBP/Implementation.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Pi/BBP/Implementation.hs @@ -0,0 +1,58 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] + + * Implements a /Bailey-Borwein-Plouffe/ formula; + + * Surprisingly, because of the huge size of the 'Data.Ratio.Rational' quantities, + it is a /single/ call to @Factory.Math.Summation.sum'@, rather than the calculation of the many terms in the series, which is the performance-bottleneck. +-} + +module Factory.Math.Implementations.Pi.BBP.Implementation( +-- * Functions + openR +) where + +import Data.Ratio((%)) +import qualified Data.Ratio +import qualified Factory.Math.Implementations.Pi.BBP.Series as Math.Implementations.Pi.BBP.Series +import qualified Factory.Math.Precision as Math.Precision +import qualified Factory.Math.Summation as Math.Summation + +-- | Returns /Pi/, accurate to the specified number of decimal digits. +openR :: + Math.Implementations.Pi.BBP.Series.Series -- ^ This /Pi/-algorithm is parameterised by the type of other algorithms to use. + -> Math.Precision.DecimalDigits -- ^ The number of decimal digits required. + -> Data.Ratio.Rational +openR Math.Implementations.Pi.BBP.Series.MkSeries { + Math.Implementations.Pi.BBP.Series.numerators = numerators, + Math.Implementations.Pi.BBP.Series.getDenominators = getDenominators, + Math.Implementations.Pi.BBP.Series.seriesScalingFactor = seriesScalingFactor, + Math.Implementations.Pi.BBP.Series.base = base +} decimalDigits = (seriesScalingFactor *) . Math.Summation.sum' 8 . take ( + Math.Precision.getTermsRequired ( + recip . fromIntegral $ abs {-potentially negative-} base --The convergence-rate. + ) decimalDigits + ) . zipWith (*) ( + iterate (/ fromIntegral base) 1 --Generate the scaling-ratio, between successive terms. + ) $ map ( + sum . zipWith (%) numerators . getDenominators + ) [0 ..] + diff --git a/src/Factory/Math/Implementations/Pi/BBP/Series.hs b/src/Factory/Math/Implementations/Pi/BBP/Series.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Pi/BBP/Series.hs @@ -0,0 +1,38 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines a /Bailey-Borwein-Plouffe/ formula; +-} + +module Factory.Math.Implementations.Pi.BBP.Series( +-- * Types +-- ** Data-types + Series(..) +) where + +import qualified Data.Ratio + +-- | Defines a series corresponding to a specific /BBP/-formula. +data Series = MkSeries { + numerators :: [Integer], -- ^ The constant numerators from which each term in the series is composed. + getDenominators :: Int -> [Integer], -- ^ Generates the term-dependent denominators from which each term in the series is composed. + seriesScalingFactor :: Data.Ratio.Rational, -- ^ The ratio by which the sum to infinity of the series, must be scaled to result in /Pi/. + base :: Integer -- ^ The geometric ratio, by which successive terms are scaled. +} + diff --git a/src/Factory/Math/Implementations/Pi/Borwein/Algorithm.hs b/src/Factory/Math/Implementations/Pi/Borwein/Algorithm.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Pi/Borwein/Algorithm.hs @@ -0,0 +1,56 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines the set of /Borwein/-type algorithms (currently only one) which have been implemented; . +-} + +module Factory.Math.Implementations.Pi.Borwein.Algorithm( +-- * Types +-- ** Data-types + Algorithm(..) +) where + +import qualified Factory.Math.Factorial as Math.Factorial +import qualified Factory.Math.Implementations.Pi.Borwein.Borwein1993 as Math.Implementations.Pi.Borwein.Borwein1993 +import qualified Factory.Math.Implementations.Pi.Borwein.Implementation as Math.Implementations.Pi.Borwein.Implementation +import qualified Factory.Math.Pi as Math.Pi +import qualified Factory.Math.SquareRoot as Math.SquareRoot +import qualified ToolShed.Defaultable as Defaultable + +{- | + * Define those /Borwein/-series which have been implemented. + + * Though currently there's only one, provision has been made for the addition of more. +-} +data Algorithm squareRootAlgorithm factorialAlgorithm = + Borwein1993 squareRootAlgorithm factorialAlgorithm -- ^ . + deriving (Eq, Read, Show) + +instance ( + Defaultable.Defaultable squareRootAlgorithm, + Defaultable.Defaultable factorialAlgorithm + ) => Defaultable.Defaultable (Algorithm squareRootAlgorithm factorialAlgorithm) where + defaultValue = Borwein1993 Defaultable.defaultValue Defaultable.defaultValue + +instance ( + Math.SquareRoot.Algorithm squareRootAlgorithm, + Math.Factorial.Algorithm factorialAlgorithm + ) => Math.Pi.Algorithm (Algorithm squareRootAlgorithm factorialAlgorithm) where + openR (Borwein1993 squareRootAlgorithm factorialAlgorithm) = Math.Implementations.Pi.Borwein.Implementation.openR Math.Implementations.Pi.Borwein.Borwein1993.series squareRootAlgorithm factorialAlgorithm + diff --git a/src/Factory/Math/Implementations/Pi/Borwein/Borwein1993.hs b/src/Factory/Math/Implementations/Pi/Borwein/Borwein1993.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Pi/Borwein/Borwein1993.hs @@ -0,0 +1,74 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines the /Borwein/ series for /Pi/; +-} + +module Factory.Math.Implementations.Pi.Borwein.Borwein1993( +-- * Constants + series +) where + +--import Control.Arrow((***)) +import qualified Data.Ratio +import Data.Ratio((%)) +--import Factory.Data.PrimeFactors((>*<), (>/<), (>^)) +--import qualified Factory.Data.PrimeFactors as Data.PrimeFactors +import qualified Factory.Math.Factorial as Math.Factorial +import qualified Factory.Math.Implementations.Factorial as Math.Implementations.Factorial +import qualified Factory.Math.Implementations.Pi.Borwein.Series as Math.Implementations.Pi.Borwein.Series +import qualified Factory.Math.Power as Math.Power +import qualified Factory.Math.Precision as Math.Precision +import qualified Factory.Math.SquareRoot as Math.SquareRoot + +-- | Defines the parameters of the /Borwein/ series. +series :: (Math.SquareRoot.Algorithm squareRootAlgorithm, Math.Factorial.Algorithm factorialAlgorithm) => Math.Implementations.Pi.Borwein.Series.Series squareRootAlgorithm factorialAlgorithm +series = Math.Implementations.Pi.Borwein.Series.MkSeries { + Math.Implementations.Pi.Borwein.Series.terms = \squareRootAlgorithm factorialAlgorithm decimalDigits -> let + simplify, squareRoot :: Data.Ratio.Rational -> Data.Ratio.Rational + simplify = Math.Precision.simplify (decimalDigits - 1 {-ignore single integral digit-}) --This makes a gigantic difference to performance. + squareRoot = simplify . Math.SquareRoot.squareRoot squareRootAlgorithm decimalDigits + + sqrt5, a, b, c3 :: Data.Ratio.Rational + sqrt5 = squareRoot 5 + + a = 63365028312971999585426220 + sqrt5 * (28337702140800842046825600 + 384 * squareRoot (10891728551171178200467436212395209160385656017 + 4870929086578810225077338534541688721351255040 * sqrt5)) + b = 7849910453496627210289749000 + 3510586678260932028965606400 * sqrt5 + 2515968 * squareRoot (3110 * (6260208323789001636993322654444020882161 + 2799650273060444296577206890718825190235 * sqrt5)) + c3 = simplify . Math.Power.cube $ negate 214772995063512240 - sqrt5 * (96049403338648032 + 1296 * squareRoot (10985234579463550323713318473 + 4912746253692362754607395912 * sqrt5)) + in ( + squareRoot $ negate c3, --The factor into which the series must be divided, to yield Pi. + zipWith ( +{- + \n power -> let + product' = Data.PrimeFactors.product' (recip 2) 10 + in uncurry (/) . ( + (* (a + b * fromIntegral n)) . fromIntegral . product' *** (* power) . fromIntegral . product' + ) $ Math.Implementations.Factorial.primeFactors (6 * n) >/< ( + Math.Implementations.Factorial.primeFactors (3 * n) >*< Math.Implementations.Factorial.primeFactors n >^ 3 + ) +-} + \n power -> ( + Math.Implementations.Factorial.risingFactorial (3 * n + 1) (3 * n) % Math.Power.cube (Math.Factorial.factorial factorialAlgorithm n) + ) * ( + (a + b * fromIntegral n) / power + ) + ) [0 :: Integer ..] $ iterate (* c3) 1 + ), + Math.Implementations.Pi.Borwein.Series.convergenceRate = 10 ** negate 50 --Empirical. +} diff --git a/src/Factory/Math/Implementations/Pi/Borwein/Implementation.hs b/src/Factory/Math/Implementations/Pi/Borwein/Implementation.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Pi/Borwein/Implementation.hs @@ -0,0 +1,57 @@ +{-# LANGUAGE CPP #-} +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines /Borwein/ series for /Pi/; +-} + +module Factory.Math.Implementations.Pi.Borwein.Implementation( +-- * Functions + openR +) where + +import qualified Control.Arrow +import qualified Data.Ratio +import qualified Factory.Math.Implementations.Pi.Borwein.Series as Math.Implementations.Pi.Borwein.Series +import qualified Factory.Math.Precision as Math.Precision + +#if MIN_VERSION_parallel(3,0,0) +import qualified Control.Parallel.Strategies +#endif + +-- | Returns /Pi/, accurate to the specified number of decimal digits. +openR :: + Math.Implementations.Pi.Borwein.Series.Series squareRootAlgorithm factorialAlgorithm -- ^ This /Pi/-algorithm is parameterised by the type of other algorithms to use. + -> squareRootAlgorithm -- ^ The specific /square-root/ algorithm to apply to the above series. + -> factorialAlgorithm -- ^ The specific /factorial/-algorithm to apply to the above series. + -> Math.Precision.DecimalDigits -- ^ The number of decimal digits required. + -> Data.Ratio.Rational +openR Math.Implementations.Pi.Borwein.Series.MkSeries { + Math.Implementations.Pi.Borwein.Series.terms = terms, + Math.Implementations.Pi.Borwein.Series.convergenceRate = convergenceRate +} squareRootAlgorithm factorialAlgorithm decimalDigits = uncurry (/) +#if MIN_VERSION_parallel(3,0,0) + . Control.Parallel.Strategies.withStrategy (Control.Parallel.Strategies.parTuple2 Control.Parallel.Strategies.rdeepseq Control.Parallel.Strategies.rdeepseq) +#endif + . Control.Arrow.second ( + sum . take ( + Math.Precision.getTermsRequired convergenceRate decimalDigits + ) + ) $ terms squareRootAlgorithm factorialAlgorithm decimalDigits + diff --git a/src/Factory/Math/Implementations/Pi/Borwein/Series.hs b/src/Factory/Math/Implementations/Pi/Borwein/Series.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Pi/Borwein/Series.hs @@ -0,0 +1,44 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines a -type series for /Pi/. +-} + +module Factory.Math.Implementations.Pi.Borwein.Series( +-- * Types +-- ** Data-types + Series(..) +) where + +import qualified Data.Ratio +import qualified Factory.Math.Precision as Math.Precision + +-- | Defines a series corresponding to a specific /Borwein/-formula. +data Series squareRootAlgorithm factorialAlgorithm = MkSeries { + terms :: + squareRootAlgorithm + -> factorialAlgorithm + -> Math.Precision.DecimalDigits + -> ( + Data.Ratio.Rational, --The factor into which the sum to infinity of the sequence, must be divided to result in /Pi/ + [Data.Ratio.Rational] --The sequence of terms, the sum to infinity of which defines the series. + ), + convergenceRate :: Math.Precision.ConvergenceRate -- ^ The expected number of digits of /Pi/, per term in the series. +} + diff --git a/src/Factory/Math/Implementations/Pi/Ramanujan/Algorithm.hs b/src/Factory/Math/Implementations/Pi/Ramanujan/Algorithm.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Pi/Ramanujan/Algorithm.hs @@ -0,0 +1,55 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines the set of /Ramanujan/-type algorithms which have been implemented; . +-} + +module Factory.Math.Implementations.Pi.Ramanujan.Algorithm( +-- * Types +-- ** Data-types + Algorithm(..) +) where + +import qualified Factory.Math.Factorial as Math.Factorial +import qualified Factory.Math.Implementations.Pi.Ramanujan.Chudnovsky as Math.Implementations.Pi.Ramanujan.Chudnovsky +import qualified Factory.Math.Implementations.Pi.Ramanujan.Classic as Math.Implementations.Pi.Ramanujan.Classic +import qualified Factory.Math.Implementations.Pi.Ramanujan.Implementation as Math.Implementations.Pi.Ramanujan.Implementation +import qualified Factory.Math.Pi as Math.Pi +import qualified Factory.Math.SquareRoot as Math.SquareRoot +import qualified ToolShed.Defaultable as Defaultable + +-- | Define those /Ramanujan/-series which have been implemented. +data Algorithm squareRootAlgorithm factorialAlgorithm = + Classic squareRootAlgorithm factorialAlgorithm -- ^ The original version. + | Chudnovsky squareRootAlgorithm factorialAlgorithm -- ^ A variant found by the /Chudnovsky brothers/. + deriving (Eq, Read, Show) + +instance ( + Defaultable.Defaultable squareRootAlgorithm, + Defaultable.Defaultable factorialAlgorithm + ) => Defaultable.Defaultable (Algorithm squareRootAlgorithm factorialAlgorithm) where + defaultValue = Chudnovsky Defaultable.defaultValue Defaultable.defaultValue + +instance ( + Math.SquareRoot.Algorithm squareRootAlgorithm, + Math.Factorial.Algorithm factorialAlgorithm + ) => Math.Pi.Algorithm (Algorithm squareRootAlgorithm factorialAlgorithm) where + openR (Classic squareRootAlgorithm factorialAlgorithm) = Math.Implementations.Pi.Ramanujan.Implementation.openR Math.Implementations.Pi.Ramanujan.Classic.series squareRootAlgorithm factorialAlgorithm + openR (Chudnovsky squareRootAlgorithm factorialAlgorithm) = Math.Implementations.Pi.Ramanujan.Implementation.openR Math.Implementations.Pi.Ramanujan.Chudnovsky.series squareRootAlgorithm factorialAlgorithm + diff --git a/src/Factory/Math/Implementations/Pi/Ramanujan/Chudnovsky.hs b/src/Factory/Math/Implementations/Pi/Ramanujan/Chudnovsky.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Pi/Ramanujan/Chudnovsky.hs @@ -0,0 +1,63 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines the /Chudnovsky/ series for /Pi/; . +-} + +module Factory.Math.Implementations.Pi.Ramanujan.Chudnovsky( +-- * Constants + series +) where + +--import Control.Arrow((***)) +import Data.Ratio((%)) +--import Factory.Data.PrimeFactors((>/<), (>*<), (>^)) +--import qualified Factory.Data.PrimeFactors as Data.PrimeFactors +import qualified Factory.Math.Factorial as Math.Factorial +import qualified Factory.Math.Implementations.Factorial as Math.Implementations.Factorial +import qualified Factory.Math.Implementations.Pi.Ramanujan.Series as Math.Implementations.Pi.Ramanujan.Series +import qualified Factory.Math.Power as Math.Power +import qualified Factory.Math.SquareRoot as Math.SquareRoot + +-- | Defines the parameters of the /Chudnovsky/ series. +series :: ( + Math.SquareRoot.Algorithm squareRootAlgorithm, + Math.Factorial.Algorithm factorialAlgorithm + ) => Math.Implementations.Pi.Ramanujan.Series.Series squareRootAlgorithm factorialAlgorithm +series = Math.Implementations.Pi.Ramanujan.Series.MkSeries { + Math.Implementations.Pi.Ramanujan.Series.terms = \factorialAlgorithm -> zipWith ( +{- + \n power -> let + product' = Data.PrimeFactors.product' (recip 2) 10 + in uncurry (%) . ( + (* (13591409 + 545140134 * n)) . product' *** (* power) . product' + ) $ Math.Implementations.Factorial.primeFactors (6 * n) >/< ( + Math.Implementations.Factorial.primeFactors (3 * n) >*< Math.Implementations.Factorial.primeFactors n >^ 3 + ) +-} + \n power -> ( + Math.Implementations.Factorial.risingFactorial (3 * n + 1) (3 * n) % Math.Power.cube (Math.Factorial.factorial factorialAlgorithm n) + ) * ( + (13591409 + 545140134 * n) % power + ) -- CAVEAT: the order in which these terms are evaluated radically affects performance. + ) [0 ..] $ iterate (* (Math.Power.cube $ negate 640320 :: Integer)) 1, + Math.Implementations.Pi.Ramanujan.Series.getSeriesScalingFactor = \squareRootAlgorithm decimalDigits -> 426880 * Math.SquareRoot.squareRoot squareRootAlgorithm decimalDigits (10005 :: Integer), + Math.Implementations.Pi.Ramanujan.Series.convergenceRate = 10 ** negate 14.0 --Empirical. +} + diff --git a/src/Factory/Math/Implementations/Pi/Ramanujan/Classic.hs b/src/Factory/Math/Implementations/Pi/Ramanujan/Classic.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Pi/Ramanujan/Classic.hs @@ -0,0 +1,60 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines the /Ramanujan/ series for /Pi/; . +-} + +module Factory.Math.Implementations.Pi.Ramanujan.Classic( +-- * Constants + series +) where + +--import Control.Arrow((***)) +import Data.Ratio((%)) +--import Factory.Data.PrimeFactors((>/<), (>^)) +--import qualified Factory.Data.PrimeFactors as Data.PrimeFactors +import qualified Factory.Math.Factorial as Math.Factorial +import qualified Factory.Math.Implementations.Factorial as Math.Implementations.Factorial +import qualified Factory.Math.Implementations.Pi.Ramanujan.Series as Math.Implementations.Pi.Ramanujan.Series +import qualified Factory.Math.Power as Math.Power +import qualified Factory.Math.SquareRoot as Math.SquareRoot + +-- | Defines the parameters of the /Ramanujan/ series. +series :: (Math.SquareRoot.Algorithm squareRootAlgorithm, Math.Factorial.Algorithm factorialAlgorithm) => Math.Implementations.Pi.Ramanujan.Series.Series squareRootAlgorithm factorialAlgorithm +series = Math.Implementations.Pi.Ramanujan.Series.MkSeries { + Math.Implementations.Pi.Ramanujan.Series.terms = \factorialAlgorithm -> let + toFourthPower = (^ (4 :: Int)) + in zipWith ( +{- + \n power -> let + product' = Data.PrimeFactors.product' (recip 2) 10 + in uncurry (%) . ( + (* (1103 + 26390 * n)) . product' *** (* power) . product' + ) $ Math.Implementations.Factorial.primeFactors (4 * n) >/< Math.Implementations.Factorial.primeFactors n >^ 4 +-} + \n power -> ( + Math.Implementations.Factorial.risingFactorial (n + 1) (3 * n) % Math.Power.cube (Math.Factorial.factorial factorialAlgorithm n) + ) * ( + (1103 + 26390 * n) % power + ) -- CAVEAT: the order in which these terms are evaluated radically affects performance. + ) [0 ..] $ iterate (* toFourthPower 396) 1, + Math.Implementations.Pi.Ramanujan.Series.getSeriesScalingFactor = \squareRootAlgorithm decimalDigits -> 9801 / Math.SquareRoot.squareRoot squareRootAlgorithm decimalDigits (8 :: Integer), + Math.Implementations.Pi.Ramanujan.Series.convergenceRate = 10 ** negate 7.9 --Empirical. +} + diff --git a/src/Factory/Math/Implementations/Pi/Ramanujan/Implementation.hs b/src/Factory/Math/Implementations/Pi/Ramanujan/Implementation.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Pi/Ramanujan/Implementation.hs @@ -0,0 +1,59 @@ +{-# LANGUAGE CPP #-} +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Implements a /Ramanujan/-type series for /Pi/; . +-} + +module Factory.Math.Implementations.Pi.Ramanujan.Implementation( +-- * Functions + openR +) where + +import qualified Data.Ratio +import qualified Factory.Math.Implementations.Pi.Ramanujan.Series as Math.Implementations.Pi.Ramanujan.Series +import qualified Factory.Math.Precision as Math.Precision +import qualified Factory.Math.Summation as Math.Summation + +#if MIN_VERSION_parallel(3,0,0) +import qualified Control.Parallel.Strategies +#endif + +-- | Returns /Pi/, accurate to the specified number of decimal digits. +openR :: + Math.Implementations.Pi.Ramanujan.Series.Series squareRootAlgorithm factorialAlgorithm -- ^ This /Pi/-algorithm is parameterised by the type of other algorithms to use. + -> squareRootAlgorithm -- ^ The specific /square-root/ algorithm to apply to the above series. + -> factorialAlgorithm -- ^ The specific /factorial/-algorithm to apply to the above series. + -> Math.Precision.DecimalDigits -- ^ The number of decimal digits required. + -> Data.Ratio.Rational +openR Math.Implementations.Pi.Ramanujan.Series.MkSeries { + Math.Implementations.Pi.Ramanujan.Series.terms = terms, + Math.Implementations.Pi.Ramanujan.Series.getSeriesScalingFactor = getSeriesScalingFactor, + Math.Implementations.Pi.Ramanujan.Series.convergenceRate = convergenceRate +} squareRootAlgorithm factorialAlgorithm decimalDigits = uncurry (/) +#if MIN_VERSION_parallel(3,0,0) + $ Control.Parallel.Strategies.withStrategy (Control.Parallel.Strategies.parTuple2 Control.Parallel.Strategies.rdeepseq Control.Parallel.Strategies.rdeepseq) +#endif + ( + getSeriesScalingFactor squareRootAlgorithm decimalDigits, + Math.Summation.sumR 64 . take ( + Math.Precision.getTermsRequired convergenceRate decimalDigits + ) $ terms factorialAlgorithm + ) + diff --git a/src/Factory/Math/Implementations/Pi/Ramanujan/Series.hs b/src/Factory/Math/Implementations/Pi/Ramanujan/Series.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Pi/Ramanujan/Series.hs @@ -0,0 +1,38 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines a -type series for /Pi/. +-} + +module Factory.Math.Implementations.Pi.Ramanujan.Series( +-- * Types +-- ** Data-types + Series(..) +) where + +import qualified Data.Ratio +import qualified Factory.Math.Precision as Math.Precision + +-- | Defines a series corresponding to a specific /Ramanujan/-formula. +data Series squareRootAlgorithm factorialAlgorithm = MkSeries { + terms :: factorialAlgorithm -> [Data.Ratio.Rational], -- ^ The sequence of terms, the sum to infinity of which defines the series. + getSeriesScalingFactor :: squareRootAlgorithm -> Math.Precision.DecimalDigits -> Data.Ratio.Rational, -- ^ The ratio by which the sum to infinity of the sequence, must be scaled to result in /Pi/. + convergenceRate :: Math.Precision.ConvergenceRate -- ^ The expected number of digits of /Pi/, per term in the series. +} + diff --git a/src/Factory/Math/Implementations/Pi/Spigot/Algorithm.hs b/src/Factory/Math/Implementations/Pi/Spigot/Algorithm.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Pi/Spigot/Algorithm.hs @@ -0,0 +1,50 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines the set of /Spigot/-algorithms which have been implemented. +-} + +module Factory.Math.Implementations.Pi.Spigot.Algorithm( +-- * Types +-- ** Data-types + Algorithm(..) +) where + +import Data.Ratio((%)) +import qualified Factory.Math.Implementations.Pi.Spigot.Gosper as Math.Implementations.Pi.Spigot.Gosper +import qualified Factory.Math.Implementations.Pi.Spigot.RabinowitzWagon as Math.Implementations.Pi.Spigot.RabinowitzWagon +import qualified Factory.Math.Implementations.Pi.Spigot.Spigot as Math.Implementations.Pi.Spigot.Spigot +import qualified Factory.Math.Pi as Math.Pi +import qualified ToolShed.Defaultable as Defaultable + +-- | Define those /Spigot/-algorithms which have been implemented. +data Algorithm = + Gosper -- ^ A /continued fraction/ discovered by /Gosper/. + | RabinowitzWagon -- ^ A /continued fraction/ discovered by /Rabinowitz/ and /Wagon/. + deriving (Eq, Read, Show) + +instance Defaultable.Defaultable Algorithm where + defaultValue = Gosper + +instance Math.Pi.Algorithm Algorithm where + openI Gosper = Math.Implementations.Pi.Spigot.Spigot.openI Math.Implementations.Pi.Spigot.Gosper.series + openI RabinowitzWagon = Math.Implementations.Pi.Spigot.Spigot.openI Math.Implementations.Pi.Spigot.RabinowitzWagon.series + + openR algorithm decimalDigits = Math.Pi.openI algorithm decimalDigits % (10 ^ (decimalDigits - 1)) + diff --git a/src/Factory/Math/Implementations/Pi/Spigot/Gosper.hs b/src/Factory/Math/Implementations/Pi/Spigot/Gosper.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Pi/Spigot/Gosper.hs @@ -0,0 +1,39 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines the /Gosper/ series; +-} + +module Factory.Math.Implementations.Pi.Spigot.Gosper( +-- * Constants + series +) where + +import qualified Factory.Math.Implementations.Pi.Spigot.Series as Math.Implementations.Pi.Spigot.Series +import qualified Factory.Math.Precision as Math.Precision + +-- | Defines a series which converges to /Pi/. +series :: Integral i => Math.Implementations.Pi.Spigot.Series.Series i +series = Math.Implementations.Pi.Spigot.Series.MkSeries { + Math.Implementations.Pi.Spigot.Series.baseNumerators = map (\i -> i * (2 * i - 1)) [1 ..], + Math.Implementations.Pi.Spigot.Series.baseDenominators = map ((* 3) . (\i -> (i + 1) * (i + 2))) [3, 6 ..], + Math.Implementations.Pi.Spigot.Series.coefficients = [3, 8 ..], --5n - 2 + Math.Implementations.Pi.Spigot.Series.nTerms = Math.Precision.getTermsRequired $ 1 / 13 {-empirical convergence-rate-} +} + diff --git a/src/Factory/Math/Implementations/Pi/Spigot/RabinowitzWagon.hs b/src/Factory/Math/Implementations/Pi/Spigot/RabinowitzWagon.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Pi/Spigot/RabinowitzWagon.hs @@ -0,0 +1,40 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines the /Rabinowitz-Wagon/ series; + + . +-} + +module Factory.Math.Implementations.Pi.Spigot.RabinowitzWagon( +-- * Constants + series +) where + +import qualified Factory.Math.Implementations.Pi.Spigot.Series as Math.Implementations.Pi.Spigot.Series +import qualified Factory.Math.Precision as Math.Precision + +-- | Defines a series which converges to /Pi/. +series :: Integral i => Math.Implementations.Pi.Spigot.Series.Series i +series = Math.Implementations.Pi.Spigot.Series.MkSeries { + Math.Implementations.Pi.Spigot.Series.baseNumerators = [1 ..], + Math.Implementations.Pi.Spigot.Series.baseDenominators = [3, 5 ..], + Math.Implementations.Pi.Spigot.Series.coefficients = repeat 2, + Math.Implementations.Pi.Spigot.Series.nTerms = Math.Precision.getTermsRequired $ 10 ** negate (3 / 10) {-convergence-rate-} +} diff --git a/src/Factory/Math/Implementations/Pi/Spigot/Series.hs b/src/Factory/Math/Implementations/Pi/Spigot/Series.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Pi/Spigot/Series.hs @@ -0,0 +1,53 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines the parameters of a series used in a /Spigot/-table to generate /Pi/. +-} + +module Factory.Math.Implementations.Pi.Spigot.Series( +-- * Types +-- ** Data-types + Series(..), +-- * Functions + bases +) where + +import Data.Ratio((%)) +import qualified Data.Ratio +import qualified Factory.Math.Precision as Math.Precision + +{- | + * Defines a series composed from a sum of terms, each one of which is the product of a coefficient and a base. + + * The coefficents and bases of the series are described in /Horner form/; @Pi = c1 + (b1 * (c2 + b2 * (c3 + b3 * (...))))@. +-} +data Series i = MkSeries { + coefficients :: [i], + baseNumerators :: [i], + baseDenominators :: [i], + nTerms :: Math.Precision.DecimalDigits -> Int -- ^ The width of the spigot-table, required to accurately generate the requested number of digits. +} + +-- | Combines 'baseNumerators' and 'baseDenominators', and as a side-effect, expresses the ratio in lowest terms. +bases :: Integral i => Series i -> [Data.Ratio.Ratio i] +bases MkSeries { + baseNumerators = n, + baseDenominators = d +} = zipWith (%) n d + diff --git a/src/Factory/Math/Implementations/Pi/Spigot/Spigot.hs b/src/Factory/Math/Implementations/Pi/Spigot/Spigot.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Pi/Spigot/Spigot.hs @@ -0,0 +1,153 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] + + * Implements a /spigot/-algorithm; . + + * Uses the traditional algorithm, rather than the /unbounded/ algorithm described by . +-} + +module Factory.Math.Implementations.Pi.Spigot.Spigot( +-- * Types +-- ** Type-synonyms +-- Base, +-- Coefficients, +-- I, +-- Pi, +-- PreDigits, +-- QuotRem, +-- * Constants + decimal, +-- * Functions +-- carryAndDivide, +-- mkRow, +-- processColumns, + openI, +-- ** Accessors +-- getQuotient, +-- getRemainder +) where + +import Data.Ratio((%)) +import qualified Control.Arrow +import qualified Data.Char +import qualified Data.Ratio +import qualified Factory.Math.Implementations.Pi.Spigot.Series as Math.Implementations.Pi.Spigot.Series +import qualified Factory.Math.Precision as Math.Precision + +{- | + * The type in which all arithmetic is performed. + + * A small dynamic range, 32 bits or more, is typically adequate. +-} +type I = Int + +-- | The constant base in which we want the resulting value of /Pi/ to be expressed. +decimal :: I +decimal = 10 + +-- | Coerce the polymorphic type 'Data.Ratio.Ratio' to suit the base used in our series. +type Base = Data.Ratio.Ratio I + +-- | Coerce the polymorphic type returned by 'quotRem' to our specific requirements. +type QuotRem = (I, I) + +-- Accessors. +getQuotient, getRemainder :: QuotRem -> I +getQuotient = fst +getRemainder = snd + +type PreDigits = [I] +type Pi = [I] +type Coefficients = [I] + +{- | + * For a digit on one row of the spigot-table, add any numerator carried from the similar calculation one column to the right. + + * Divide the result of this summation, by the denominator of the base, to get the quotient and remainder. + + * Determine the quantity to carry to the similar calculation one column to the left, by multiplying the quotient by the numerator of the base. +-} +carryAndDivide :: (Base, I) -> QuotRem -> QuotRem +carryAndDivide (base, lhs) rhs + | n < d = (0, n) --In some degenerate cases, the result of the subsequent calculation can be more simply determined. + | otherwise = Control.Arrow.first (* Data.Ratio.numerator base) $ n `quotRem` d + where + d, n :: I + d = Data.Ratio.denominator base + n = lhs + getQuotient rhs --Carry numerator from the column to the right and add it to the current digit. + +{- | + * Fold 'carryAndDivide', from right to left, over the columns of a row in the spigot-table, continuously checking for overflow. + + * Release any previously withheld result-digits, after any adjustment for overflow in the current result-digit. + + * Withhold the current result-digit until the risk of overflow in subsequent result-digits has been assessed. + + * Call 'mkRow'. +-} +processColumns + :: Math.Implementations.Pi.Spigot.Series.Series I + -> PreDigits + -> [(Base, I)] -- ^ Data-row. + -> Pi +processColumns series preDigits l + | overflowMargin > 1 = preDigits ++ nextRow [digit] --There's neither overflow, nor risk of impact from subsequent overflow. + | overflowMargin == 1 = nextRow $ preDigits ++ [digit] --There's no overflow, but risk of impact from subsequent overflow. + | otherwise = map ((`mod` decimal) . (+ 1)) preDigits ++ nextRow [0] --Overflow => propagate the excess to previously withheld preDigits. + where + results :: [QuotRem] + results = init $ scanr carryAndDivide (0, undefined) l + + digit :: I + digit = getQuotient $ head results + + overflowMargin :: I + overflowMargin = decimal - digit + + nextRow :: [I] -> [I] + nextRow preDigits' = mkRow series preDigits' $ map getRemainder results + +{- | + * Multiply the remainders from the previous row. + + * Zip them with the constant bases, with an addition one stuck on the front to perform the conversion to decimal, to create a new row of the spigot-table. + + * Call 'processColumns'. +-} +mkRow :: Math.Implementations.Pi.Spigot.Series.Series I -> PreDigits -> Coefficients -> Pi +mkRow series preDigits = processColumns series preDigits . zip (1 % decimal : Math.Implementations.Pi.Spigot.Series.bases series) . map (* decimal) + +{- | + * Initialises a /spigot/-table with the row of 'Math.Implementations.Pi.Spigot.Series.coefficients'. + + * Ensures that the row has suffient terms to accurately generate the required number of digits. + + * Extracts only those digits which are guaranteed to be accurate. + + * CAVEAT: the result is returned as an 'Integer', i.e. without any decimal point. +-} +openI :: Math.Implementations.Pi.Spigot.Series.Series I -> Math.Precision.DecimalDigits -> Integer +openI series decimalDigits = read . map ( + Data.Char.intToDigit . fromIntegral + ) . take decimalDigits . mkRow series [] . take ( + Math.Implementations.Pi.Spigot.Series.nTerms series decimalDigits + ) $ Math.Implementations.Pi.Spigot.Series.coefficients series + diff --git a/src/Factory/Math/Implementations/Primality.hs b/src/Factory/Math/Implementations/Primality.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/Primality.hs @@ -0,0 +1,220 @@ +{-# LANGUAGE CPP #-} +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] + + * Determines whether an integer is prime. + + * . + + * + + * CAVEAT: it doesn't determine the prime-factors of composite numbers, just that they exist. +-} + +module Factory.Math.Implementations.Primality( +-- * Types +-- ** Data-types + Algorithm(..) +-- * Functions +-- ** Predicates +-- isPrimeByAKS, +-- isPrimeByMillerRabin, +-- witnessesCompositeness +) where + +import Control.Arrow((&&&)) +import qualified Control.DeepSeq +import qualified Data.Numbers.Primes +import qualified Factory.Data.MonicPolynomial as Data.MonicPolynomial +import qualified Factory.Data.Polynomial as Data.Polynomial +import qualified Factory.Data.QuotientRing as Data.QuotientRing +import qualified Factory.Math.MultiplicativeOrder as Math.MultiplicativeOrder +import qualified Factory.Math.Power as Math.Power +import qualified Factory.Math.Primality as Math.Primality +import qualified Factory.Math.PrimeFactorisation as Math.PrimeFactorisation +import qualified ToolShed.Defaultable as Defaultable + +#if MIN_VERSION_parallel(3,0,0) +import qualified Control.Parallel.Strategies +#endif + +-- | The algorithms by which /primality/-testing has been implemented. +data Algorithm factorisationAlgorithm = + AKS factorisationAlgorithm -- ^ . + | MillerRabin -- ^ . + deriving (Eq, Read, Show) + +instance Defaultable.Defaultable (Algorithm factorisationAlgorithm) where + defaultValue = MillerRabin + +instance Math.PrimeFactorisation.Algorithm factorisationAlgorithm => Math.Primality.Algorithm (Algorithm factorisationAlgorithm) where + isPrime _ 2 = True --The only even prime. + isPrime algorithm candidate + | candidate < 2 || ( + any ( + (== 0) . (candidate `rem`) --The candidate has a small prime-factor, and is therefore composite. + ) . filter ( + (candidate >=) . (* 2) --The candidate must be at least double the small prime, for it to be a potential factor. + ) . take 5 {-arbitrarily-} $ Data.Numbers.Primes.primes --Excludes even numbers, provided at least the 1st prime is tested. + ) = False + | otherwise = ( + case algorithm of + AKS factorisationAlgorithm -> isPrimeByAKS factorisationAlgorithm + MillerRabin -> isPrimeByMillerRabin + ) candidate + +{- | + * An implementation of the /Agrawal-Kayal-Saxena/ primality-test; , + using the /Lenstra/ and /Pomerance/ algorithm. + + * CAVEAT: this deterministic algorithm has a theoretical time-complexity of @O(log^6)@, + and therefore can't compete with the performance of probabilistic ones. + + * The /formal polynomials/ used in this algorithm, are conceptually different from /polynomial functions/; + the /indeterminate/ and its powers, are merely used to name a sequence of pigeon-holes in which /coefficients/ are stored, + and is never substituted for a specific value. + This mind-shift, allows one to introduce concepts like /modular/ arithmetic on polynomials, + which merely represent an operation on their coefficients and the pigeon-hole in which they're placed. + + [@Manindra Agrawal, Neeraj Kayal and Nitin Saxena@] . + + [@H. W. Lenstra, Jr. and Carl Pomerance@] . + + [@Salembier and Southerington@] , + + [@R. Crandall and J. Papadopoulos@] , + + [@Andreas Klappenecker@] , + + [@Vibhor Bhatt and G. K. Patra@] , +-} +isPrimeByAKS :: (Math.PrimeFactorisation.Algorithm factorisationAlgorithm, Control.DeepSeq.NFData i, Integral i) => factorisationAlgorithm -> i -> Bool +isPrimeByAKS factorisationAlgorithm n = and [ + not $ Math.Power.isPerfectPower n, --Step 1. + Math.Primality.areCoprime n `all` filter (/= n) [2 .. r], --Step 3. +#if MIN_VERSION_parallel(3,0,0) + and $ Control.Parallel.Strategies.parMap Control.Parallel.Strategies.rdeepseq --Benefits from '+RTS -H100M', which reduces garbage-collections. +#else + all +#endif + ( + \a -> let +-- lhs, rhs :: Data.Polynomial.Polynomial i i + lhs = Data.Polynomial.raiseModulo (Data.Polynomial.mkLinear 1 a) n {-power-} n {-modulus-} + rhs = Data.Polynomial.mod' (Data.Polynomial.mkPolynomial [(1, n), (a, 0)]) n + in Data.QuotientRing.areCongruentModulo ( + Data.MonicPolynomial.mkMonicPolynomial lhs + ) ( + Data.MonicPolynomial.mkMonicPolynomial rhs + ) ( + Data.MonicPolynomial.mkMonicPolynomial modulus + ) -- Because all these polynomials are /monic/, one can establish /congruence/ using /integer/-division. + ) [ + 1 .. floor . (* lg) . sqrt $ fromIntegral r + ] --Step 4; (x + a)^n ~ x^n + a mod (x^r - 1, n). + ] where + lg :: Double + lg = logBase 2 $ fromIntegral n + +-- r :: i + r = fst . head . dropWhile ( + (<= floor (Math.Power.square lg)) . snd + ) . map ( + id &&& Math.MultiplicativeOrder.multiplicativeOrder factorisationAlgorithm n + ) $ Math.Primality.areCoprime n `filter` [2 ..] --Step 2. + +-- modulus :: Data.Polynomial.Polynomial i i + modulus = Data.Polynomial.mkPolynomial [(1, r), (negate 1, 0)] + +{- | + * Uses the specified 'base' in an attempt to prove the /compositeness/ of an integer. + + * This is the opposite of the /Miller Test/; . + + * If the result is 'True', then the candidate is /composite/; regrettably the converse isn't true. + Amongst the set of possible bases, over three-quarters are /witnesses/ to the compositeness of a /composite/ candidate, + the remainder belong to the subset of /liars/. + In consequence, many false results must be accumulated for different bases, to convincingly identify a prime. +-} +witnessesCompositeness :: Integral i + => i -- ^ Candidate integer. + -> i + -> Int + -> i -- ^ Base. + -> Bool +witnessesCompositeness candidate oddRemainder nPowersOfTwo base = all ( + $ ((`mod` candidate) . Math.Power.square) `iterate` Math.Power.raiseModulo base oddRemainder candidate --Repeatedly modulo-square. + ) [ + (/= 1) . head, --Check whether the zeroeth modulo-power is incongruent to one. + all (/= pred candidate) . take nPowersOfTwo --Check whether any modulo-power is incongruent to -1. + ] + +{- | + * Repeatedly calls 'witnessesCompositeness', to progressively increase the probability of detecting a /composite/ number, + until ultimately the candidate integer is proven to be prime. + + * Should all bases be tested, then the test is deterministic, but at an efficiency /lower/ than performing prime-factorisation. + + * The test becomes deterministic, for any candidate integer, when the number of tests reaches the limit defined by /Eric Bach/. + + * A testing of smaller set of bases, is sufficient for candidates smaller than various thresholds; . + + * . + + * + + * . + + * , . +-} +isPrimeByMillerRabin :: Integral i => i -> Bool +isPrimeByMillerRabin primeCandidate = not $ witnessesCompositeness primeCandidate ( + fst $ last binaryFactors --Odd-remainder. + ) ( + length binaryFactors --The number of times that 'two' can be factored-out from 'predecessor'. + ) `any` testBases where + predecessor = primeCandidate - 1 + binaryFactors = takeWhile ((== 0) . snd) . tail {-drop the original-} $ iterate ((`quotRem` 2) . fst) (predecessor, 0) --Factor-out powers of two. + testBases + | null fewestPrimeBases = let + millersTestSet = floor . (* 2 {-Eric Bach-}) . Math.Power.square . toRational {-avoid premature rounding-} $ log (fromIntegral primeCandidate :: Double {-overflows at 10^851-}) + in [2 .. predecessor `min` millersTestSet] + | otherwise = head fewestPrimeBases `take` Data.Numbers.Primes.primes + where + fewestPrimeBases = map fst $ dropWhile ((primeCandidate >=) . snd) [ + (0, 9), --All odd integers less this, are prime, and require no further verification. + (1, 2047), + (2, 1373653), + (3, 25326001), + (4, 3215031751), + (5, 2152302898747), --Jaeschke ... + (6, 3474749660383), + (8, 341550071728321), + (11, 3825123056546413051), --Zhang ... + (12, 318665857834031151167461), + (13, 3317044064679887385961981), + (14, 6003094289670105800312596501), + (15, 59276361075595573263446330101), + (17, 564132928021909221014087501701), + (19, 1543267864443420616877677640751301), + (20, 10 ^ (36 :: Int)) --At least. + ] + diff --git a/src/Factory/Math/Implementations/PrimeFactorisation.hs b/src/Factory/Math/Implementations/PrimeFactorisation.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/PrimeFactorisation.hs @@ -0,0 +1,150 @@ +{-# LANGUAGE CPP #-} +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] + + * Implements several different prime-factorisation algorithms. + + * . +-} + +module Factory.Math.Implementations.PrimeFactorisation( +-- * Types +-- ** Data-types + Algorithm( +-- DixonsMethod, + FermatsMethod, + TrialDivision + ) +-- * Functions +-- factoriseByDixonsMethod +-- factoriseByFermatsMethod +-- factoriseByTrialDivision +) where + +import Control.Arrow((&&&), (***)) +import qualified Control.Arrow +import qualified Control.DeepSeq +import qualified Data.Maybe +import qualified Data.Numbers.Primes +import qualified Factory.Data.Exponential as Data.Exponential +import Factory.Data.Exponential((<^)) +import qualified Factory.Data.PrimeFactors as Data.PrimeFactors +import qualified Factory.Math.Power as Math.Power +import qualified Factory.Math.PrimeFactorisation as Math.PrimeFactorisation +import qualified ToolShed.Defaultable as Defaultable + +#if MIN_VERSION_parallel(3,0,0) +import qualified Control.Parallel.Strategies +#endif + +-- | The algorithms by which prime-factorisation has been implemented. +data Algorithm + = DixonsMethod -- ^ . + | FermatsMethod -- ^ . + | TrialDivision -- ^ . + deriving (Eq, Read, Show) + +instance Defaultable.Defaultable Algorithm where + defaultValue = TrialDivision + +instance Math.PrimeFactorisation.Algorithm Algorithm where + primeFactors algorithm = case algorithm of + DixonsMethod -> factoriseByDixonsMethod + FermatsMethod -> Data.PrimeFactors.reduce . factoriseByFermatsMethod + TrialDivision -> factoriseByTrialDivision + +{- | + * . +-} +factoriseByDixonsMethod :: Integral base => base -> Data.PrimeFactors.Factors base exponent +factoriseByDixonsMethod = undefined + +{- | + * . + + * . + + * . + + * @i = f1 * f2@ Assume a non-trivial factorisation, ie. one in which both factors exceed one. + => @i = (larger + smaller) * (larger - smaller)@ Represent the co-factors as a sum and difference. + => @i = larger^2 - smaller^2@ Which has an integral solution if @i@ is neither /even/ nor a /perfect square/. + => @sqrt (larger^2 - i) = smaller@ Search for /larger/, which results in an integral value for /smaller/. + + * Given that the smaller factor /f2/, can't be less than 3 (/i/ isn't /even/), then the larger /f1/, can't be greater than @(i `div` 3)@. + So: @(f2 >= 3) && (f1 <= i `div` 3)@ Two equations which can be used to solve for /larger/. + => @(larger - smaller >= 3) && (larger + smaller <= i `div` 3)@ Add these to eliminate /smaller/. + => @larger <= (i + 9) `div` 6@ The upper bound of the search-space. + + * This algorithm works best when there's a factor close to the /square-root/. +-} +factoriseByFermatsMethod :: (Control.DeepSeq.NFData base, Integral base, Control.DeepSeq.NFData exponent, Num exponent) => base -> Data.PrimeFactors.Factors base exponent +factoriseByFermatsMethod i + | i <= 3 = [Data.Exponential.rightIdentity i] + | even i = Data.Exponential.rightIdentity 2 : factoriseByFermatsMethod (i `div` 2) {-recurse-} + | Data.Maybe.isJust maybeSquareNumber = (<^ 2) `map` factoriseByFermatsMethod (Data.Maybe.fromJust maybeSquareNumber) {-recurse-} + | null factors = [Data.Exponential.rightIdentity i] --Prime. + | otherwise = uncurry (++) . +#if MIN_VERSION_parallel(3,0,0) + Control.Parallel.Strategies.withStrategy ( + Control.Parallel.Strategies.parTuple2 Control.Parallel.Strategies.rdeepseq Control.Parallel.Strategies.rdeepseq --CAVEAT: unproductive on the size of integers tested so far. + ) . +#endif + ( + factoriseByFermatsMethod *** factoriseByFermatsMethod --Divide and conquer. + ) $ head factors + where +-- maybeSquareNumber :: Integral i => Maybe i + maybeSquareNumber = Math.Power.maybeSquareNumber i + +-- factors :: Integral i => [i] + factors = map ( + ( + uncurry (+) &&& uncurry (-) --Construct the co-factors as the sum and difference of /larger/ and /smaller/. + ) . Control.Arrow.second Data.Maybe.fromJust + ) . filter ( + Data.Maybe.isJust . snd --Search for a perfect square. + ) . map ( + Control.Arrow.second $ Math.Power.maybeSquareNumber {-hotspot-} . (+ negate i) --Associate the corresponding value of /smaller/. + ) . takeWhile ( + (<= (i + 9) `div` 6) . fst --Terminate the search at the maximum value of /larger/. + ) . Math.Power.squaresFrom {-hotspot-} . ceiling $ sqrt (fromIntegral i :: Double) --Start the search at the minimum value of /larger/. + +{- | + * Decomposes the specified integer, into a product of /prime/-factors, + using , AKA . + + * This works best when the factors are small. +-} +factoriseByTrialDivision :: (Integral base, Num exponent) => base -> Data.PrimeFactors.Factors base exponent +factoriseByTrialDivision = slave Data.Numbers.Primes.primes where + slave primes i + | null primeCandidates = [Data.Exponential.rightIdentity i] + | otherwise = Data.Exponential.rightIdentity lowestPrimeFactor `Data.PrimeFactors.insert'` slave primeCandidates (i `quot` lowestPrimeFactor) + where + primeCandidates = dropWhile ( + (/= 0) . (i `rem`) + ) $ takeWhile ( + <= Math.PrimeFactorisation.maxBoundPrimeFactor i + ) primes + + lowestPrimeFactor = head primeCandidates + diff --git a/src/Factory/Math/Implementations/SquareRoot.hs b/src/Factory/Math/Implementations/SquareRoot.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Implementations/SquareRoot.hs @@ -0,0 +1,192 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Implements 'Math.SquareRoot.Algorithm' by a variety of methods. + + [@CAVEAT@] + + Caller may benefit from application of 'Math.Precision.simplify' before operating on the result; + which though of the required accuracy, may not be the most concise rational number satisfying that criterion. +-} +module Factory.Math.Implementations.SquareRoot( +-- * Types +-- ** Type-synonyms +-- ProblemSpecification, + Terms, +-- ** Data-types + Algorithm(..) +-- * Functions +-- squareRootByContinuedFraction, +-- squareRootByIteration, +-- squareRootByTaylorSeries, +-- taylorSeriesCoefficients +) where + +import Control.Arrow((***)) +import Factory.Data.PrimeFactors((>/<), (>^)) +import qualified Factory.Data.PrimeFactors as Data.PrimeFactors +import qualified Factory.Math.Implementations.Factorial as Math.Implementations.Factorial +import qualified Factory.Math.Power as Math.Power +import qualified Factory.Math.Precision as Math.Precision +import qualified Factory.Math.SquareRoot as Math.SquareRoot +import qualified Factory.Math.Summation as Math.Summation +import qualified ToolShed.Defaultable as Defaultable + +-- | The number of terms in a series. +type Terms = Int + +-- | The algorithms by which the /square-root/ has been implemented. +data Algorithm + = BakhshaliApproximation -- ^ + | ContinuedFraction -- ^ . + | HalleysMethod -- ^ . + | NewtonRaphsonIteration -- ^ . + | TaylorSeries Terms -- ^ . + deriving (Eq, Read, Show) + +instance Defaultable.Defaultable Algorithm where + defaultValue = NewtonRaphsonIteration + +-- | Returns an improved estimate for the /square-root/ of the specified value, to the required precision, using the supplied initial estimate.. +type ProblemSpecification operand + = Math.SquareRoot.Estimate + -> Math.Precision.DecimalDigits -- ^ The required precision. + -> operand -- ^ The value for which to find the /square-root/. + -> Math.SquareRoot.Result + +instance Math.SquareRoot.Algorithm Algorithm where + squareRootFrom _ _ _ 0 = 0 + squareRootFrom _ _ _ 1 = 1 + squareRootFrom algorithm estimate@(x, decimalDigits) requiredDecimalDigits y + | decimalDigits >= requiredDecimalDigits = x + | requiredDecimalDigits <= 0 = error $ "Factory.Math.Implementations.SquareRoot.squareRootFrom:\tinvalid number of required decimal digits; " ++ show requiredDecimalDigits + | y < 0 = error $ "Factory.Math.Implementations.SquareRoot.squareRootFrom:\tthere's no real square-root of " ++ show y + | otherwise = ( + case algorithm of + ContinuedFraction -> squareRootByContinuedFraction + _ -> squareRootByIteration algorithm + ) estimate requiredDecimalDigits y + +instance Math.SquareRoot.Iterator Algorithm where + step BakhshaliApproximation y x + | dy == 0 = x --The estimate was precise. + | otherwise = x' - dx' --Correct the estimate. + where + dy, dydx, dx, x', dydx', dx' :: Math.SquareRoot.Result + dy = Math.SquareRoot.getDiscrepancy y x + dydx = 2 * x + dx = dy / dydx + x' = x + dx --Identical to Newton-Raphson iteration. + dydx' = 2 * x' + dx' = Math.Power.square dx / dydx' + +{- + * /Halley's/ method; + +> X(n+1) = Xn - f(Xn) / [f'(Xn) - f''(Xn) * f(Xn) / 2 * f'(Xn)] +> => Xn - (Xn^2 - Y) / [2Xn - 2 * (Xn^2 - Y) / 2 * 2Xn] where Y = X^2, f(X) = X^2 - Y, f'(X) = 2X, f''(X) = 2 +> => Xn - 1 / [2Xn / (Xn^2 - Y) - 1 / 2Xn] +-} + step HalleysMethod y x + | dy == 0 = x --The estimate was precise. + | otherwise = x - dx --Correct the estimate. + where + dy, dydx, dx :: Math.SquareRoot.Result + dy = negate $ Math.SquareRoot.getDiscrepancy y x --Use the estimate to determine the error in 'y'. + dydx = 2 * x --The gradient, at the estimated value 'x'. + dx = recip $ dydx / dy - recip dydx + +-- step NewtonRaphsonIteration y x = (x + realToFrac y / x) / 2 --This is identical to the /Babylonian Method/. +-- step NewtonRaphsonIteration y x = x / 2 + realToFrac y / (2 * x) --Faster. + step NewtonRaphsonIteration y x = x / 2 + (realToFrac y / 2) / x --Faster still. + + step (TaylorSeries terms) y x = squareRootByTaylorSeries terms y x + + step algorithm _ _ = error $ "Factory.Math.Implementations.SquareRoot.step:\tinappropriate algorithm; " ++ show algorithm + + convergenceOrder BakhshaliApproximation = Math.Precision.quarticConvergence + convergenceOrder ContinuedFraction = Math.Precision.linearConvergence + convergenceOrder HalleysMethod = Math.Precision.cubicConvergence + convergenceOrder NewtonRaphsonIteration = Math.Precision.quadraticConvergence + convergenceOrder (TaylorSeries terms) = terms --The order of convergence, per iteration, equals the number of terms in the series on each iteration. + +{- | + * Uses /continued-fractions/, to iterate towards the principal /square-root/ of the specified positive integer; + , + , + . + + + * The convergence of the /continued-fraction/ is merely /1st order/ (linear). +-} +squareRootByContinuedFraction :: Real operand => ProblemSpecification operand +squareRootByContinuedFraction (initialEstimate, initialDecimalDigits) requiredDecimalDigits y = initialEstimate + (convergents initialEstimate !! Math.Precision.getTermsRequired (10 ^^ negate initialDecimalDigits) requiredDecimalDigits) where + convergents :: Math.SquareRoot.Result -> [Math.SquareRoot.Result] + convergents x = iterate ((Math.SquareRoot.getDiscrepancy y x /) . ((2 * x) +)) 0 + +{- | + * The constant coefficients of the /Taylor-series/ for a /square-root/; . + + * @ ((-1)^n * factorial(2*n)) / ((1 - 2*n) * 4^n * factorial(n^2)) @. +-} +taylorSeriesCoefficients :: Fractional f => [f] +taylorSeriesCoefficients = zipWith ( + \powers n -> let + doubleN = 2 * n + product' = Data.PrimeFactors.product' (recip 2) {-arbitrary-} 10 {-arbitrary-} + in uncurry (/) . ( + fromIntegral . product' *** fromIntegral . (* ((1 - doubleN) * powers)) . product' + ) $ Math.Implementations.Factorial.primeFactors doubleN >/< Math.Implementations.Factorial.primeFactors n >^ 2 + ) ( + iterate (* negate 4) 1 -- (-4)^n + ) [0 :: Integer ..] -- n + +{- | + * Returns the /Taylor-series/ for the /square-root/ of the specified value, to any requested number of terms. + + * . + + * The convergence of the series is merely /linear/, + in that each term increases the precision, by a constant number of decimal places, equal to the those in the original estimate. + + * By feeding-back the improved estimate, to form a new series, the order of convergence, on each successive iteration, + becomes proportional to the number of terms; + +> Terms Convergence +> ===== =========== +> 2 terms /quadratic/ +> 3 terms /cubic/ +-} +squareRootByTaylorSeries :: Real operand + => Terms -- ^ The number of terms of the infinite series, to evaluate. + -> operand -- ^ The value for which the /square-root/ is required. + -> Math.SquareRoot.Result -- ^ An initial estimate. + -> Math.SquareRoot.Result +squareRootByTaylorSeries _ _ 0 = error "Factory.Math.Implementations.SquareRoot.squareRootByTaylorSeries:\talgorithm can't cope with estimated value of zero." +squareRootByTaylorSeries terms y x + | terms < 2 = error $ "Factory.Math.Implementations.SquareRoot.squareRootByTaylorSeries:\tinvalid number of terms; " ++ show terms + | otherwise = Math.Summation.sumR' . take terms . zipWith (*) taylorSeriesCoefficients $ iterate (* relativeError) x + where + relativeError :: Math.SquareRoot.Result + relativeError = (realToFrac y / Math.Power.square x) - 1 --Pedantically, this is the error in y, which is twice the magnitude of the error in x. + +-- | Iterates from the estimated value, towards the /square-root/, a sufficient number of times to achieve the required accuracy. +squareRootByIteration :: Real operand => Algorithm -> ProblemSpecification operand +squareRootByIteration algorithm (initialEstimate, initialDecimalDigits) requiredDecimalDigits y = iterate (Math.SquareRoot.step algorithm y) initialEstimate !! Math.Precision.getIterationsRequired (Math.SquareRoot.convergenceOrder algorithm) initialDecimalDigits requiredDecimalDigits + diff --git a/src/Factory/Math/MultiplicativeOrder.hs b/src/Factory/Math/MultiplicativeOrder.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/MultiplicativeOrder.hs @@ -0,0 +1,66 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Exports the /Multiplicative Order/ of an integer, in a specific /modular/ arithmetic. + +-} + +module Factory.Math.MultiplicativeOrder( +-- * Functions + multiplicativeOrder +) where + +import qualified Control.DeepSeq +import qualified Factory.Data.Exponential as Data.Exponential +import qualified Factory.Math.Power as Math.Power +import qualified Factory.Math.Primality as Math.Primality +import qualified Factory.Math.PrimeFactorisation as Math.PrimeFactorisation + +{- | + * The smallest positive integral power to which the specified integral base must be raised, + to be congruent with one, in the specified /modular/ arithmetic. + + * Based on . + + * . + + * . +-} +multiplicativeOrder :: (Math.PrimeFactorisation.Algorithm primeFactorisationAlgorithm, Control.DeepSeq.NFData i, Integral i) + => primeFactorisationAlgorithm + -> i -- ^ Base. + -> i -- ^ Modulus. + -> i -- ^ Result. +multiplicativeOrder primeFactorisationAlgorithm base modulus + | modulus < 2 = error $ "Factory.Math.MultiplicativeOrder.multiplicativeOrder:\tinvalid modulus; " ++ show modulus + | not $ Math.Primality.areCoprime base modulus = error $ "Factory.Math.MultiplicativeOrder.multiplicativeOrder:\targuments aren't coprime; " ++ show (base, modulus) + | otherwise = foldr (lcm . multiplicativeOrder') 1 $ Math.PrimeFactorisation.primeFactors primeFactorisationAlgorithm modulus --Combine the /multiplicative order/ of the constituent /prime-factors/. + where +-- multiplicativeOrder' :: (Control.DeepSeq.NFData i, Integral i) => Data.Exponential.Exponential i -> i + multiplicativeOrder' e = product . map ( + \e' -> let + d :: Int + d = length . takeWhile (/= 1) . iterate ( + \y -> Math.Power.raiseModulo y (Data.Exponential.getBase e') pk + ) $ Math.Power.raiseModulo base (totient `div` Data.Exponential.evaluate e') pk + in Data.Exponential.getBase e' ^ d + ) $ Math.PrimeFactorisation.primeFactors primeFactorisationAlgorithm totient where + pk = Data.Exponential.evaluate e + totient = Math.PrimeFactorisation.primePowerTotient e + diff --git a/src/Factory/Math/Pi.hs b/src/Factory/Math/Pi.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Pi.hs @@ -0,0 +1,101 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines the classes of /Pi/-algorithm which have been implemented. +-} + +module Factory.Math.Pi( +-- * Type-classes + Algorithm(..), +-- * Types +-- ** Data-types + Category(..) +) where + +import qualified Data.Ratio +import qualified Factory.Math.Precision as Math.Precision +import qualified ToolShed.Defaultable as Defaultable + +{- | + * Defines the methods expected of a /Pi/-algorithm. + + * Most of the implementations naturally return a 'Rational', but the spigot-algorithms naturally produce a @[Int]@; + though representing /Pi/ as a big integer with the decimal point removed is clearly incorrect. + + * Since representing /Pi/ as either a 'Rational' or promoted to an 'Integer', is inconvenient, an alternative decimal 'String'-representation is provided. +-} +class Algorithm algorithm where + openR :: algorithm -> Math.Precision.DecimalDigits -> Data.Ratio.Rational -- ^ Returns the value of /Pi/ as a 'Rational'. + + openI :: algorithm -> Math.Precision.DecimalDigits -> Integer -- ^ Returns the value of /Pi/, promoted by the required precision to form an integer. + openI _ 1 = 3 + openI algorithm decimalDigits + | decimalDigits <= 0 = error $ "Factory.Math.Pi.openI:\tinsufficient decimalDigits=" ++ show decimalDigits + | otherwise = round . Math.Precision.promote (openR algorithm decimalDigits) $ decimalDigits - 1 + + openS :: algorithm -> Math.Precision.DecimalDigits -> String -- ^ Returns the value of /Pi/ as a decimal 'String'. + openS _ 1 = "3" + openS algorithm decimalDigits + | decimalDigits <= 0 = "" + | decimalDigits <= 16 = take (decimalDigits + 1) $ show (pi :: Double) + | otherwise = "3." ++ tail (show $ openI algorithm decimalDigits) --Insert a decimal point. + +-- | Categorises the various algorithms. +data Category agm bbp borwein ramanujan spigot + = AGM agm -- ^ Algorithms based on the /Arithmetic-geometric Mean/. + | BBP bbp -- ^ . + | Borwein borwein -- ^ . + | Ramanujan ramanujan -- ^ . + | Spigot spigot -- ^ Algorithms from which the digits of /Pi/ slowly drip, one by one. + deriving (Eq, Read, Show) + +instance ( + Defaultable.Defaultable agm, + Defaultable.Defaultable bbp, + Defaultable.Defaultable borwein, + Defaultable.Defaultable ramanujan, + Defaultable.Defaultable spigot + ) => Defaultable.Defaultable (Category agm bbp borwein ramanujan spigot) where + defaultValue = BBP Defaultable.defaultValue + +instance ( + Algorithm agm, + Algorithm bbp, + Algorithm borwein, + Algorithm ramanujan, + Algorithm spigot + ) => Algorithm (Category agm bbp borwein ramanujan spigot) where + openR algorithm decimalDigits + | decimalDigits <= 0 = error $ "Factory.Math.Pi.openR:\tinsufficient decimalDigits=" ++ show decimalDigits + | decimalDigits <= 16 = Math.Precision.simplify (decimalDigits - 1) (pi :: Double) + | otherwise = ( + case algorithm of + AGM agm -> openR agm + BBP bbp -> openR bbp + Borwein borwein -> openR borwein + Ramanujan ramanujan -> openR ramanujan + Spigot spigot -> openR spigot + ) decimalDigits + + openI _ 1 = 3 + openI (Spigot spigot) decimalDigits = openI spigot decimalDigits + openI algorithm decimalDigits + | decimalDigits <= 0 = error $ "Factory.Math.Pi.openI:\tinsufficient decimalDigits=" ++ show decimalDigits + | otherwise = round . Math.Precision.promote (openR algorithm decimalDigits) $ decimalDigits - 1 + diff --git a/src/Factory/Math/Power.hs b/src/Factory/Math/Power.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Power.hs @@ -0,0 +1,137 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Exports functions involving integral powers. +-} + +module Factory.Math.Power( +-- * Functions + square, + squaresFrom, + maybeSquareNumber, + cube, + cubeRoot, + raiseModulo, +-- ** Predicates + isPerfectPower +) where + +import qualified Data.Set + +-- | Mainly for convenience. +{-# INLINE square #-} +square :: Num n => n -> n +square = (^ (2 :: Int)) + +-- | Just for convenience. +cube :: Num n => n -> n +cube = (^ (3 :: Int)) + +{- | + * Iteratively generate sequential /squares/, from the specified initial value, + based on the fact that @(x + 1)^2 = x^2 + 2 * x + 1@. + + * The initial value doesn't need to be either positive or integral. +-} +squaresFrom :: Num n => n -> [(n, n)] +squaresFrom from = iterate (\(x, y) -> (x + 1, y + 2 * x + 1)) (from, square from) + +-- | Just for convenience. +cubeRoot :: Double -> Double +cubeRoot = (** recip 3) + +{- | + * Raise an arbitrary number to the specified positive integral power, using /modular/ arithmetic. + + * Implements exponentiation as a sequence of either /squares/ or multiplications by the base; + . + + * . +-} +raiseModulo :: (Integral i, Integral power) + => i -- ^ Base. + -> power + -> i -- ^ Modulus. + -> i -- ^ Result. +raiseModulo _ _ 0 = error "Factory.Math.Power.raiseModulo:\tzero modulus." +raiseModulo _ _ 1 = 0 +raiseModulo _ 0 modulus = 1 `mod` modulus +raiseModulo base power modulus + | base < 0 = (`mod` modulus) . (if odd power then negate else id) $ raiseModulo (negate base) power modulus --Recurse. + | power < 0 = error $ "Factory.Math.Power.raiseModulo:\tnegative power; " ++ show power + | first `elem` [0, 1] = first + | otherwise = slave power + where + first = base `mod` modulus + + slave 1 = first + slave e = (`mod` modulus) . (if r == 0 {-even-} then id else (* base)) . square $ slave q {-recurse-} where + (q, r) = e `quotRem` 2 + +{- | + * Returns @(Just . sqrt)@ if the specified integer is a /square number/ (AKA /perfect square/). + + * . + + * . + + * @(square . sqrt)@ is expensive, so the modulus of the operand is tested first, in an attempt to prove it isn't a /perfect square/. + The set of tests, and the valid moduli within each test, are ordered to maximize the rate of failure-detection. +-} +maybeSquareNumber :: Integral i => i -> Maybe i +maybeSquareNumber i +-- | i < 0 = Nothing --This function is performance-sensitive, but this test is neither strictly nor frequently required. + | all (\(modulus, valid) -> mod i modulus `elem` valid) [ +-- --Distribution of moduli amongst perfect squares Cumulative failure-detection. + (16, [0,1,4,9]), --All moduli are equally likely. 75% + (9, [0,1,4,7]), --Zero occurs 33%, the others only 22%. 88% + (17, [1,2,4,8,9,13,15,16,0]), --Zero only occurs 5.8%, the others 11.8%. 94% +-- These additional tests, aren't always cost-effective. + (13, [1,3,4,9,10,12,0]), --Zero only occurs 7.7%, the others 15.4%. 97% + (7, [1,2,4,0]), --Zero only occurs 14.3%, the others 28.6%. 98% + (5, [1,4,0]) --Zero only occurs 20%, the others 40%. 99% + +-- ] && fromIntegral iSqrt == sqrt' = Just iSqrt --CAVEAT: erroneously True for 187598574531033120, whereas 187598574531033121 is square. + ] && square iSqrt == i = Just iSqrt + | otherwise = Nothing + where + sqrt' :: Double + sqrt' = sqrt $ fromIntegral i + + iSqrt = round sqrt' + +{- | + * An integer @(> 1)@ which can be expressed as an integral power @(> 1)@ of a smaller /natural/ number. + + * CAVEAT: /zero/ and /one/ are normally excluded from this set. + + * . + + * . +-} +isPerfectPower :: Integral i => i -> Bool +isPerfectPower i + | i < square 2 = False + | otherwise = i `Data.Set.member` foldr ( + \n set -> if n `Data.Set.member` set + then set +-- else Data.Set.union set . Data.Set.fromList . takeWhile (<= i) . iterate (* n) $ square n --TODO: test relative speed. + else foldr Data.Set.insert set . takeWhile (<= i) . iterate (* n) $ square n + ) Data.Set.empty [2 .. round $ sqrt (fromIntegral i :: Double)] + diff --git a/src/Factory/Math/Precision.hs b/src/Factory/Math/Precision.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Precision.hs @@ -0,0 +1,118 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines the unit with which precision is measured, and operations on it. +-} +module Factory.Math.Precision( +-- * Types +-- ** Type-synonyms + ConvergenceOrder, + ConvergenceRate, + DecimalDigits, +-- * Constants + linearConvergence, + quadraticConvergence, + cubicConvergence, + quarticConvergence, +-- * Functions + getIterationsRequired, + getTermsRequired, + promote, + simplify +) where + +import qualified Data.Ratio + +-- | The /order of convergence/; . +type ConvergenceOrder = Int + +-- | The /rate of convergence/; . +type ConvergenceRate = Double + +-- | A number of decimal digits. +type DecimalDigits = Int + +-- | /Linear/ convergence-rate; which may be qualified by the /rate of convergence/. +linearConvergence :: ConvergenceOrder +linearConvergence = 1 + +-- | /Quadratic/ convergence-rate. +quadraticConvergence :: ConvergenceOrder +quadraticConvergence = 2 + +-- | /Cubic/ convergence-rate. +cubicConvergence :: ConvergenceOrder +cubicConvergence = 3 + +-- | /Quartic/ convergence-rate. +quarticConvergence :: ConvergenceOrder +quarticConvergence = 4 + +-- | The predicted number of iterations, required to achieve a specific accuracy, at a given /order of convergence/. +getIterationsRequired :: Integral i + => ConvergenceOrder + -> DecimalDigits -- ^ The precision of the initial estimate. + -> DecimalDigits -- ^ The required precision. + -> i +getIterationsRequired convergenceOrder initialDecimalDigits requiredDecimalDigits + | initialDecimalDigits <= 0 = error $ "Factory.Math.Precision.getIterationsRequired:\tinsufficient 'initialDecimalDigits'; " ++ show initialDecimalDigits + | precisionRatio <= 1 = 0 + | otherwise = ceiling $ fromIntegral convergenceOrder `logBase` precisionRatio + where + precisionRatio :: Double + precisionRatio = fromIntegral requiredDecimalDigits / fromIntegral initialDecimalDigits + +{- | + * The predicted number of terms which must be extracted from a series, + if it is to converge to the required accuracy, + at the specified linear /convergence-rate/. + + * The /convergence-rate/ of a series, is the error in the series after summation of @(n+1)th@ terms, + divided by the error after only @n@ terms, as the latter tends to infinity. + As such, for a /convergent/ series (in which the error get smaller with successive terms), it's value lies in the range @0 .. 1@. + + * . +-} +getTermsRequired :: Integral i + => ConvergenceRate + -> DecimalDigits -- ^ The additional number of correct decimal digits. + -> i +getTermsRequired _ 0 = 0 +getTermsRequired convergenceRate requiredDecimalDigits + | convergenceRate <= 0 || convergenceRate >= 1 = error $ "Factory.Math.Precision.getTermsRequired:\t (0 < convergence-rate < 1); " ++ show convergenceRate + | requiredDecimalDigits < 0 = error $ "Factory.Math.Precision.getTermsRequired:\t'requiredDecimalDigits' must be positive; " ++ show requiredDecimalDigits + | otherwise = ceiling $ fromIntegral requiredDecimalDigits / negate (logBase 10 convergenceRate) + +-- | Promotes the specified number, by a number of 'DecimalDigits'. +promote :: Num n => n -> DecimalDigits -> n +promote x = (* x) . (10 ^) + +{- | + * Reduces a 'Data.Ratio.Rational' to the minimal form required for the specified number of /fractional/ decimal places; + irrespective of the number of integral decimal places. + + * A 'Data.Ratio.Rational' approximation to an irrational number, may be very long, and provide an unknown excess precision. + Whilst this doesn't sound harmful, it costs in performance and memory-requirement, and being unpredictable isn't actually useful. +-} +simplify :: RealFrac operand + => DecimalDigits -- ^ The number of places after the decimal point, which are required. + -> operand + -> Data.Ratio.Rational +simplify decimalDigits operand = Data.Ratio.approxRational operand . recip $ 4 * 10 ^ (decimalDigits + 1) --Tolerate any error less than half the least significant digit required. + diff --git a/src/Factory/Math/Primality.hs b/src/Factory/Math/Primality.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Primality.hs @@ -0,0 +1,92 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] + + * Exports a common interface for primality-implementations. + + * Provides utilities for these implementations. +-} + +module Factory.Math.Primality( +-- * Type-classes + Algorithm(..), +-- * Functions + carmichaelNumbers, +-- ** Predicates + areCoprime, + isFermatWitness, + isCarmichaelNumber +) where + +import qualified Control.DeepSeq +import qualified Factory.Math.Power as Math.Power + +-- | Defines the methods expected of a primality-algorithm. +class Algorithm algorithm where + isPrime :: (Control.DeepSeq.NFData i, Integral i) => algorithm -> i -> Bool + +{- | + 'True' if the two specified integers are /relatively prime/, + i.e. if they share no common positive factors except one. + + * @1@ and @-1@ are the only numbers which are /coprime/ to themself. + + * . + + * . +-} +areCoprime :: Integral i => i -> i -> Bool +areCoprime i = (== 1) . gcd i + +{- | + * Tests /Fermat's Little Theorem/ for all applicable values, as a probabilistic primality-test. + + * . + + * . + + * . + + * CAVEAT: this primality-test fails for the /Carmichael numbers/. + + * TODO: confirm that all values must be tested. +-} +isFermatWitness :: Integral i => i -> Bool +isFermatWitness i = not . all isFermatPseudoPrime $ filter (areCoprime i) [2 .. i - 1] where + isFermatPseudoPrime base = Math.Power.raiseModulo base (i - 1) i == 1 --CAVEAT: a /Fermat Pseudo-prime/ must also be a /composite/ number. + +{- | + * A /Carmichael number/ is an odd /composite/ number which satisfies /Fermat's little theorem/. + + * . + + * . +-} +isCarmichaelNumber :: (Algorithm algorithm, Control.DeepSeq.NFData i, Integral i) => algorithm -> i -> Bool +isCarmichaelNumber algorithm i = not $ or [ + i <= 2, + even i, + isFermatWitness i, + isPrime algorithm i + ] + +-- | An ordered list of the /Carmichael/ numbers; . +carmichaelNumbers :: (Algorithm algorithm, Control.DeepSeq.NFData i, Integral i) => algorithm -> [i] +carmichaelNumbers algorithm = isCarmichaelNumber algorithm `filter` [3, 5 ..] diff --git a/src/Factory/Math/PrimeFactorisation.hs b/src/Factory/Math/PrimeFactorisation.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/PrimeFactorisation.hs @@ -0,0 +1,135 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] + + * . + + * Exports a common interface to permit decomposition of positive integers, + into the unique combination of /prime/-factors known to exist according to the /Fundamental Theorem of Arithmetic/; . + + * Leveraging this abstract capability, it derives the /smoothness/, /power-smoothness/, and /omega/-numbers. + + * Filters the list of /regular-numbers/ from the list of /smoothness/. + + * CAVEAT: to avoid wasting time, it may be advantageous to check /Factory.Math.Primality.isPrime/ first. +-} + +module Factory.Math.PrimeFactorisation( +-- * Type-classes + Algorithm(..), +-- * Functions + maxBoundPrimeFactor, + smoothness, + powerSmoothness, + regularNumbers, + primePowerTotient, + eulersTotient, + omega +) where + +import qualified Control.DeepSeq +import qualified Data.List +import qualified Factory.Data.Exponential as Data.Exponential +import qualified Factory.Data.PrimeFactors as Data.PrimeFactors + +-- | Defines the methods expected of a /factorisation/-algorithm. +class Algorithm algorithm where + primeFactors :: (Control.DeepSeq.NFData base, Integral base) + => algorithm + -> base -- ^ The operand + -> Data.PrimeFactors.Factors base Int {-arbitrarily-} + +{- | + * The upper limit for a prime to be considered as a candidate factor of the specified number. + + * One might naively think that this limit is @(x `div` 2)@ for an even number, + but though a prime-factor /greater/ than the /square-root/ of the number can exist, + its smaller /cofactor/ decomposes to a prime which must be less than the /square-root/. + + * NB: rather using @(primeFactor <= sqrt numerator)@ to filter the candidate prime-factors of a given numerator, + one can alternatively use @(numerator >= primeFactor ^ 2)@ to filter what can potentially be factored by a given prime-factor. +-} +maxBoundPrimeFactor :: Integral i => i -> i +maxBoundPrimeFactor = floor . (sqrt :: Double -> Double) . fromIntegral + +{- | + * A constant, zero-indexed, conceptually infinite, list, of the /smooth/ness of all positive integers. + + * . + + * . +-} +smoothness :: (Algorithm algorithm, Control.DeepSeq.NFData base, Integral base) => algorithm -> [base] +smoothness algorithm = 0 : map (Data.Exponential.getBase . last . primeFactors algorithm) [1 ..] + +{- | + * A constant, zero-indexed, conceptually infinite, list of the /power-smooth/ness of all positive integers. + + * . +-} +powerSmoothness :: (Algorithm algorithm, Control.DeepSeq.NFData base, Integral base) => algorithm -> [base] +powerSmoothness algorithm = 0 : map (maximum . map Data.Exponential.evaluate . primeFactors algorithm) [1 ..] + +{- | + * Filters 'smoothness', to derive the constant list of /Hamming-numbers/. + + * . +-} +regularNumbers :: (Algorithm algorithm, Control.DeepSeq.NFData base, Integral base) => algorithm -> [base] +regularNumbers algorithm = map fst . filter ((<= (5 :: Integer)) . snd) . zip [1 ..] . tail $ smoothness algorithm + +{- | + * /Euler's Totient/ for a /power/ of a /prime/-number. + + * By /Olofsson/; @(phi(n^k) = n^(k - 1) * phi(n))@ + and since @(phi(prime) = prime - 1)@ + + * CAVEAT: checks neither the primality nor the bounds of the specified value; therefore for internal use only. +-} +primePowerTotient :: (Integral base, Integral exponent) => Data.Exponential.Exponential base exponent -> base +primePowerTotient (base, exponent') = pred base * base ^ pred exponent' + +{- | + * The number of /coprimes/ less than or equal to the specified positive integer. + + * . + + * . + + * AKA /EulerPhi/. +-} +eulersTotient :: (Algorithm algorithm, Control.DeepSeq.NFData i, Integral i) => algorithm -> i -> i +eulersTotient _ 1 = 1 +eulersTotient algorithm i + | i <= 0 = error $ "Factory.Math.PrimeFactorisation.eulersTotient:\tundefined for; " ++ show i + | otherwise = product . map primePowerTotient $ primeFactors algorithm i + +{- | + * A constant, zero-indexed, conceptually infinite, list of the /small omega/ numbers, the number of /distinct/ prime factors; cf. /big omega/. + + * . + + * + + * . +-} +omega :: (Algorithm algorithm, Control.DeepSeq.NFData i, Integral i) => algorithm -> [i] +omega algorithm = map (Data.List.genericLength . primeFactors algorithm) [0 :: Integer ..] + diff --git a/src/Factory/Math/Radix.hs b/src/Factory/Math/Radix.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Radix.hs @@ -0,0 +1,118 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Facilitates representation of 'Integral' values in alternative 'Integral' bases. +-} + +module Factory.Math.Radix( +-- * Constants +-- decodes, +-- digits, +-- encodes, +-- * Functions + digitSum, + digitalRoot, + fromBase, + toBase +) where + +import Data.Array((!)) +import qualified Data.Array +import qualified Data.Char +import qualified Data.List +import qualified Data.Maybe + +-- | Characters used to represent the digits of numbers in @(-36 <= base <= 36)@. +digits :: String +digits = ['0' .. '9'] ++ ['a' .. 'z'] + +-- | Constant random-access lookup for 'digits'. +encodes :: (Data.Array.Ix index, Integral index) => Data.Array.Array index Char +encodes = Data.Array.listArray (0, fromIntegral $ length digits - 1) digits + +-- | Constant reverse-lookup for 'digits'. +decodes :: Integral i => [(Char, i)] +decodes = zip digits [0 ..] + +{- | + * Convert the specified integral decimal quantity, to an alternative base, and represent the result as a 'String'. + + * Both negative decimals and negative bases are permissible. + + * The conversion to 'Char' can only succeed where printable and intelligible characters exist to represent all digits in the chosen base; + which in practice means @(-36 <= base <= 36)@. +-} +toBase :: (Integral base, Integral decimal) => base -> decimal -> String +toBase 10 decimal = show decimal --Base unchanged. +toBase _ 0 = "0" --Zero has the same representation in any base. +toBase base decimal + | abs base < 2 = error $ "Factory.Math.Radix.toBase:\tan arbitrary integer can't be represented in base " ++ show base + | abs base > Data.List.genericLength digits = error $ "Factory.Math.Radix.toBase:\tunable to clearly represent the complete set of digits in base " ++ show base + | base > 0 && decimal < 0 = '-' : map (toDigit . fromIntegral) (fromDecimal (negate decimal) []) + | otherwise = (toDigit . fromIntegral) `map` fromDecimal decimal [] + where + fromDecimal 0 = id + fromDecimal n + | remainder < 0 = fromDecimal (quotient + 1) . ((remainder - fromIntegral base) :) --This can only occur when base is negative; cf. 'divMod'. + | otherwise = fromDecimal quotient . (remainder :) + where + (quotient, remainder) = n `quotRem` fromIntegral base + + toDigit :: Int -> Char + toDigit n + | n >&< encodes = encodes ! n + | otherwise = error $ "Factory.Math.Radix.toBase.toDigit:\tno suitable character-representation for integer " ++ show n + where + (>&<) :: Int -> Data.Array.Array Int Char -> Bool + index >&< array = ($ index) `all` [(>= lower), (<= upper)] where + (lower, upper) = Data.Array.bounds array + +{- | + * Convert the 'String'-representation of a number in the specified base, to a decimal integer. + + * Both negative numbers and negative bases are permissible. +-} +fromBase :: (Integral base, Integral decimal, Read decimal) => base -> String -> decimal +fromBase 10 s = read s --Base unchanged. +fromBase _ "0" = 0 --Zero has the same representation in any base. +fromBase base s + | abs base < 2 = error $ "Factory.Math.Radix.fromBase:\tan arbitrary integer can't be represented in base " ++ show base + | abs base > Data.List.genericLength digits = error $ "Factory.Math.Radix.fromBase:\tunable to clearly represent the complete set of digits in base " ++ show base + | base > 0 && head s == '-' = negate . fromBase base $ tail s --Recurse. + | otherwise = Data.List.foldl' (\l -> ((l * fromIntegral base) +) . fromIntegral . fromDigit) 0 s where + fromDigit :: Char -> Int + fromDigit c = case c `lookup` decodes of + Just i + | i >= abs (fromIntegral base) -> error $ "Factory.Math.Radix.fromBase.fromDigit:\tillegal char " ++ show c ++ ", for base " ++ show base + | otherwise -> i + _ -> error $ "Factory.Math.Radix.fromBase.fromDigit:\tunrecognised char " ++ show c + +{- | + * . + + * . +-} +digitSum :: (Integral base, Integral decimal) => base -> decimal -> decimal +digitSum 10 = fromIntegral . foldr ((+) . Data.Char.digitToInt) 0 . show +digitSum base = sum . Data.Maybe.mapMaybe (`lookup` decodes) . toBase base + +-- | . +digitalRoot :: Integral decimal => decimal -> decimal +digitalRoot = head . dropWhile (> 9) . iterate (digitSum (10 :: Int)) + diff --git a/src/Factory/Math/SquareRoot.hs b/src/Factory/Math/SquareRoot.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/SquareRoot.hs @@ -0,0 +1,121 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] + + * Exports a common interface for /square-root/ implementations. + + * Provides utilities for these implementations. +-} + +module Factory.Math.SquareRoot( +-- * Type-classes + Algorithm(..), + Iterator(..), +-- * Types +-- ** Type-synonyms + Result, + Estimate, +-- * Functions + getAccuracy, + getDiscrepancy, + getEstimate, +-- rSqrt, +-- ** Predicates + isPrecise +) where + +import qualified Data.Ratio +import qualified Factory.Math.Power as Math.Power +import qualified Factory.Math.Precision as Math.Precision + +-- | The result-type; actually, only the concrete return-type of 'Math.Precision.simplify', stops it being a polymorphic instance of 'Fractional'. +type Result = Data.Ratio.Rational + +-- | Contains an estimate for the /square-root/ of a value, and its accuracy. +type Estimate = (Result, Math.Precision.DecimalDigits) + +-- | Defines the methods expected of a /square-root/ algorithm. +class Algorithm algorithm where + squareRootFrom :: Real operand + => algorithm + -> Estimate -- ^ An initial estimate from which to start. + -> Math.Precision.DecimalDigits -- ^ The required precision. + -> operand -- ^ The value for which to find the /square-root/. + -> Result -- ^ Returns an improved estimate of the /square-root/, found using the specified algorithm, accurate to at least the required number of decimal digits. + + squareRoot :: Real operand + => algorithm + -> Math.Precision.DecimalDigits -- ^ The required precision. + -> operand -- ^ The value for which to find the /square-root/. + -> Result -- ^ Returns an estimate of the /square-root/, found using the specified algorithm, accurate to at least the required number of decimal digits. + squareRoot algorithm decimalDigits operand = squareRootFrom algorithm (getEstimate operand) decimalDigits operand --Default implementation + +-- | The interface required to iterate, from an estimate of the required value, to the next approximation. +class Iterator algorithm where + step :: Real operand + => algorithm + -> operand -- ^ The value for which the /square-root/ is required; @y@. + -> Result -- ^ The current estimate; @x(n)@. + -> Result -- ^ An improved estimate; @x(n+1)@. + + convergenceOrder :: algorithm -> Math.Precision.ConvergenceOrder -- ^ The ultimate ratio of successive terms as the iteration converges. + +-- | Generalise 'sqrt' to operate on any 'Real' operand. +rSqrt :: Real operand => operand -> Double +rSqrt = sqrt . realToFrac + +-- | Uses 'Double'-precision floating-point arithmetic, to obtain an initial estimate for the /square-root/, and its accuracy. +getEstimate :: Real operand => operand -> Estimate +getEstimate y + | y < 0 = error $ "Factory.Math.SquareRoot.getEstimate:\tthere's no real square-root of " ++ show y + | otherwise = (Math.Precision.simplify decimalDigits {-doubles performance by roughly halving number's length-} . toRational $ rSqrt y, decimalDigits) + where + decimalDigits :: Math.Precision.DecimalDigits + decimalDigits = 16 -- . + +{- | + * The signed difference between the square of an estimate for the /square-root/ of a value, and that value. + + * Positive when the estimate is too low. + + * CAVEAT: the magnitude is twice the error in the /square-root/. +-} +getDiscrepancy :: Real operand => operand -> Result -> Result +getDiscrepancy y x = realToFrac y - Math.Power.square x + +-- | 'True' if the specified estimate for the /square-root/, is precise. +isPrecise :: Real operand => operand -> Result -> Bool +isPrecise y x = getDiscrepancy y x == 0 + +{- | + * For a given value and an estimate of its /square-root/, + returns the number of decimals digits to which the /square-root/ is accurate; including the integral digits. + + * CAVEAT: the result returned for an exact match has been bodged. +-} +getAccuracy :: Real operand => operand -> Result -> Math.Precision.DecimalDigits +getAccuracy y x + | absoluteError == 0 = maxBound --Bodge. +-- | otherwise = length . takeWhile (< 1) $ iterate (* 10) relativeError --CAVEAT: too slow. + | otherwise = length $ show (round $ realToFrac y / absoluteError :: Integer) + where + absoluteError :: Result + absoluteError = abs (getDiscrepancy y x) / 2 --NB: the magnitude of the error in 'y', is twice the error in its square-root, 'x'. + diff --git a/src/Factory/Math/Statistics.hs b/src/Factory/Math/Statistics.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Statistics.hs @@ -0,0 +1,79 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Miscellaneous statistical functions. +-} + +module Factory.Math.Statistics( +-- * Functions + mean, + nCr, + nPr +) where + +import Control.Arrow((***)) +import Control.Parallel(par, pseq) +import qualified Data.List +--import qualified Factory.Data.PrimeFactors as Data.PrimeFactors +--import Factory.Data.PrimeFactors((>/<), (>*<)) +import qualified Factory.Math.Factorial as Math.Factorial +import qualified Factory.Math.Implementations.Factorial as Math.Implementations.Factorial + +-- | Determines the of the supplied numbers. +mean :: (Real r, Fractional f) => [r] -> f +mean [] = error "Factory.Math.Statistics.mean:\tundefined result for specified null-list" +mean l = uncurry (/) . (realToFrac *** fromIntegral) $ foldr (\s -> (+ s) *** succ) (0, 0 :: Int) l + +-- | The number of unordered combinations of /r/ objects taken from /n/; . +nCr :: (Math.Factorial.Algorithm factorialAlgorithm, Integral i) + => factorialAlgorithm + -> i -- ^ The total number of items from which to select. + -> i -- ^ The number of iterms in a sample. + -> i -- ^ The number of combinations. +nCr _ 0 _ = 1 +nCr _ _ 0 = 1 +nCr factorialAlgorithm n r + | n < 0 = error $ "Factory.Math.Statistics.nCr:\tinvalid n; " ++ show n + | r < 0 = error $ "Factory.Math.Statistics.nCr:\tinvalid r; " ++ show r + | n < r = 0 +{- + | otherwise = uncurry div $ product' *** product' $ Math.Implementations.Factorial.primeFactors n >/< ( + Math.Implementations.Factorial.primeFactors r >*< Math.Implementations.Factorial.primeFactors (n - r) + ) where + product' = Data.PrimeFactors.product' (recip 2) 10 +-} + | otherwise = numerator `par` (denominator `pseq` numerator `div` denominator) + where + [smaller, bigger] = Data.List.sort [r, n - r] + numerator = Math.Implementations.Factorial.risingFactorial (bigger + 1) (n - bigger) + denominator = Math.Factorial.factorial factorialAlgorithm smaller + +-- | The number of permutations of /r/ objects taken from /n/; . +nPr :: Integral i + => i -- ^ The total number of items from which to select. + -> i -- ^ The number of items in a sample. + -> i -- ^ The number of permutations. +nPr 0 _ = 1 +nPr _ 0 = 1 +nPr n r + | n < 0 = error $ "Factory.Math.Statistics.nPr:\tinvalid n; " ++ show n + | r < 0 = error $ "Factory.Math.Statistics.nPr:\tinvalid r; " ++ show r + | n < r = 0 + | otherwise = Math.Implementations.Factorial.fallingFactorial n r + diff --git a/src/Factory/Math/Summation.hs b/src/Factory/Math/Summation.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Math/Summation.hs @@ -0,0 +1,113 @@ +{-# LANGUAGE CPP #-} +{- + Copyright (C) 2010 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Provides an alternative algorithm for the summation of /rational/ numbers. +-} + +module Factory.Math.Summation( +-- * Functions + sum', + sumR', + sumR +) where + +import qualified Data.List +import qualified Data.Ratio +import Data.Ratio((%)) +import qualified Control.DeepSeq + +#if MIN_VERSION_parallel(3,0,0) +import qualified Control.Parallel.Strategies +import qualified ToolShed.ListPlus as ListPlus +#endif + +{- | + * Sums a list of numbers of arbitrary type. + + * Sparks the summation of @(list-length / chunk-size)@ chunks from the list, each of the specified size (thought the last chunk may be smaller), + then recursively sums the list of results from each spark. + + * CAVEAT: unless the numbers are large, 'Data.Ratio.Rational' (requiring /cross-multiplication/), or the list long, + 'sum' is too light-weight for sparking to be productive, + therefore it is more likely to be the parallelised deep /evaluation/ of list-elements which saves time. +-} +sum' :: (Num n, Control.DeepSeq.NFData n) +#if MIN_VERSION_toolshed(11,1,1) + => ListPlus.ChunkLength +#else + => Int -- ^ The Chunk-length. +#endif + -> [n] + -> n +#if MIN_VERSION_parallel(3,0,0) +sum' chunkLength + | chunkLength <= 1 = error $ "Factory.Math.Summation.sum':\tinvalid chunk-size; " ++ show chunkLength + | otherwise = slave + where + slave :: (Num n, Control.DeepSeq.NFData n) => [n] -> n + slave [] = 0 + slave [x] = x + slave l = slave {-recurse-} . Control.Parallel.Strategies.parMap Control.Parallel.Strategies.rdeepseq sum $ ListPlus.chunk chunkLength l +#else +sum' _ = sum +#endif + +{- + * Sums a list of /rational/ type numbers. + + * CAVEAT: though faster than 'Data.List.sum', this algorithm has poor space-complexity, making it unsuitable for unrestricted use. +-} +{-# INLINE sumR' #-} --This makes a staggering difference. +sumR' :: Integral i => [Data.Ratio.Ratio i] -> Data.Ratio.Ratio i +sumR' l = foldr (\ratio -> ((Data.Ratio.numerator ratio * (commonDenominator `div` Data.Ratio.denominator ratio)) +)) 0 l % commonDenominator where +-- commonDenominator = foldr (lcm . Data.Ratio.denominator) 1 l + commonDenominator = Data.List.foldl' (\multiple -> lcm multiple . Data.Ratio.denominator) 1 l --Slightly faster. + +{- | + * Sums a list of /rational/ numbers. + + * Sparks the summation of @(list-length / chunk-length)@ chunks from the list, each of the specified size (thought the last chunk may be smaller), + then recursively sums the list of results from each spark. + + * CAVEAT: memory-use is proportional to chunk-size. +-} +{-# INLINE sumR #-} --This makes a staggering difference to calls from other modules. +sumR :: (Integral i, Control.DeepSeq.NFData i) +#if MIN_VERSION_toolshed(11,1,1) + => ListPlus.ChunkLength +#else + => Int -- ^ The Chunk-length. +#endif + -> [Data.Ratio.Ratio i] + -> Data.Ratio.Ratio i +sumR chunkLength + | chunkLength <= 1 = error $ "Factory.Math.Summation.sumR:\tinvalid chunk-size; " ++ show chunkLength + | otherwise = slave + where + slave :: (Integral i, Control.DeepSeq.NFData i) => [Data.Ratio.Ratio i] -> Data.Ratio.Ratio i + slave l + | length l <= chunkLength = sumR' l + | otherwise = slave {-recurse-} . +#if MIN_VERSION_parallel(3,0,0) + Control.Parallel.Strategies.parMap Control.Parallel.Strategies.rdeepseq +#else + map +#endif + sumR' $ ListPlus.chunk chunkLength l diff --git a/src/Factory/Test/CommandOptions.hs b/src/Factory/Test/CommandOptions.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/CommandOptions.hs @@ -0,0 +1,48 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines the available set of command-line options; of which there's currently only one. +-} + +module Factory.Test.CommandOptions( +-- * Types +-- ** Data-types + CommandOptions(..), +-- * Functions +-- ** Mutators + setVerbose +) where + +import ToolShed.Defaultable as Defaultable + +-- | Declare a record used to contain command-line options. +data CommandOptions = MkCommandOptions { + verbose :: Bool -- ^ Whether additional informative output should be generated, where applicable. +} + +instance Defaultable CommandOptions where + defaultValue = MkCommandOptions { verbose = False } + +-- | Mutator. +setVerbose :: CommandOptions -> CommandOptions +setVerbose commandOptions = commandOptions { + verbose = True +} + + diff --git a/src/Factory/Test/Performance/Factorial.hs b/src/Factory/Test/Performance/Factorial.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/Performance/Factorial.hs @@ -0,0 +1,68 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Times functions exported from module "Math.Factorial". +-} + +module Factory.Test.Performance.Factorial( +-- * Functions + factorialPerformance, + factorialPerformanceControl, + factorialPerformanceGraph, + factorialPerformanceGraphControl +) where + +import qualified Control.DeepSeq +import qualified Data.List +import qualified Factory.Math.Factorial as Math.Factorial +import qualified ToolShed.TimePure as TimePure + +-- | Measures the CPU-time required by 'Math.Factorial.factorial'. +factorialPerformance :: (Math.Factorial.Algorithm algorithm, Control.DeepSeq.NFData i, Integral i) => algorithm -> i -> IO (Double, i) +factorialPerformance algorithm = TimePure.getCPUSeconds . Math.Factorial.factorial algorithm + +-- | Measures the CPU-time required by a naive implementation. +factorialPerformanceControl :: (Control.DeepSeq.NFData i, Integral i) => i -> IO (Double, i) +--factorialPerformanceControl i = TimePure.getCPUSeconds $ product [1 .. i] --CAVEAT: too lazy. +factorialPerformanceControl i = TimePure.getCPUSeconds $ Data.List.foldl' (*) 1 [2 .. i] + +{- | + * Measure the CPU-time required by 'Math.Factorial.factorial', against an exponentially increasing operand. + + * CAVEAT: nothing is returned, since the result is printed ... and it never terminates. +-} +factorialPerformanceGraph :: Math.Factorial.Algorithm algorithm => Bool -> algorithm -> IO () +factorialPerformanceGraph verbose algorithm = mapM_ ( + \operand -> factorialPerformance algorithm operand >>= putStrLn . shows operand . showChar '\t' . ( + if verbose + then (`shows` "") + else (`shows` "") . fst + ) + ) $ iterate (* 2) (1 :: Integer) + +-- | Graphs the CPU-time required by a naive implementation, against an exponentially increasing operand. +factorialPerformanceGraphControl :: Bool -> IO () +factorialPerformanceGraphControl verbose = mapM_ ( + \operand -> factorialPerformanceControl operand >>= putStrLn . shows operand . showChar '\t' . ( + if verbose + then (`shows` "") + else (`shows` "") . fst + ) + ) $ iterate (* 2) (1 :: Integer) + diff --git a/src/Factory/Test/Performance/Pi.hs b/src/Factory/Test/Performance/Pi.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/Performance/Pi.hs @@ -0,0 +1,81 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Times functions exported from module "Math.Pi". +-} + +module Factory.Test.Performance.Pi( +-- * Types +-- ** Type-synonyms + Category, +-- * Functions + piPerformance, + piPerformanceGraph +) where + +import qualified Factory.Math.Factorial as Math.Factorial +import qualified Factory.Math.Implementations.Pi.AGM.Algorithm as Math.Implementations.Pi.AGM.Algorithm +import qualified Factory.Math.Implementations.Pi.BBP.Algorithm as Math.Implementations.Pi.BBP.Algorithm +import qualified Factory.Math.Implementations.Pi.Borwein.Algorithm as Math.Implementations.Pi.Borwein.Algorithm +import qualified Factory.Math.Implementations.Pi.Ramanujan.Algorithm as Math.Implementations.Pi.Ramanujan.Algorithm +import qualified Factory.Math.Implementations.Pi.Spigot.Algorithm as Math.Implementations.Pi.Spigot.Algorithm +import qualified Factory.Math.Pi as Math.Pi +import qualified Factory.Math.Precision as Math.Precision +import qualified Factory.Math.SquareRoot as Math.SquareRoot +import qualified ToolShed.TimePure as TimePure + +-- | The type of a /Pi/-algorithm, including where required, the algorithm for /square-root/s and /factorial/s. +type Category squareRootAlgorithm factorialAlgorithm = Math.Pi.Category ( + Math.Implementations.Pi.AGM.Algorithm.Algorithm squareRootAlgorithm + ) Math.Implementations.Pi.BBP.Algorithm.Algorithm ( + Math.Implementations.Pi.Borwein.Algorithm.Algorithm squareRootAlgorithm factorialAlgorithm + ) ( + Math.Implementations.Pi.Ramanujan.Algorithm.Algorithm squareRootAlgorithm factorialAlgorithm + ) Math.Implementations.Pi.Spigot.Algorithm.Algorithm + +-- | Measures the CPU-time required to find Pi to the required precision. +piPerformance :: ( + Math.SquareRoot.Algorithm squareRootAlgorithm, + Math.Factorial.Algorithm factorialAlgorithm + ) => Category squareRootAlgorithm factorialAlgorithm -> Math.Precision.DecimalDigits -> IO (Double, String) +piPerformance category = TimePure.getCPUSeconds . Math.Pi.openS category + +{- | + * Measures the CPU-time required to determine /Pi/ to an exponentially increasing precision-requirement. + + * CAVEAT: nothing is returned, since the result is printed ... and it never terminates. +-} +piPerformanceGraph :: ( + Math.SquareRoot.Algorithm squareRootAlgorithm, + Show squareRootAlgorithm, + Math.Factorial.Algorithm factorialAlgorithm, + Show factorialAlgorithm + ) => RealFrac i + => Category squareRootAlgorithm factorialAlgorithm -- ^ The algorithm. + -> i -- ^ The factor by which the precision is increased on each iteration. + -> Math.Precision.DecimalDigits -- ^ The maximum precision required. + -> Bool -- ^ Whether to return the digits of /Pi/. + -> IO () +piPerformanceGraph category factor maxDecimalDigits verbose = mapM_ ( + \decimalDigits -> piPerformance category decimalDigits >>= putStrLn . shows decimalDigits . showChar '\t' . ( + if verbose + then (`shows` "") + else (`shows` "") . fst + ) + ) . takeWhile (<= maxDecimalDigits) . map round $ iterate (* factor) 1 diff --git a/src/Factory/Test/Performance/Primality.hs b/src/Factory/Test/Performance/Primality.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/Performance/Primality.hs @@ -0,0 +1,53 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Times functions exported from module "Math.Primality". +-} + +module Factory.Test.Performance.Primality( +-- * Functions + carmichaelNumbersPerformance, + isPrimePerformance, + isPrimePerformanceGraph +) where + +import qualified Factory.Math.Fibonacci as Math.Fibonacci +import qualified Factory.Math.Primality as Math.Primality +import qualified ToolShed.TimePure as TimePure + +-- | Measures the CPU-time required to find the specified number of /Carmichael/ numbers, which is returned together with the requested list. +carmichaelNumbersPerformance :: Math.Primality.Algorithm primalityAlgorithm => primalityAlgorithm -> Int -> IO (Double, [Integer]) +carmichaelNumbersPerformance primalityAlgorithm i + | i < 0 = error $ "Factory.Test.Performance.Primality.carmichaelNumbersPerformance:\tnegative number; " ++ show i + | otherwise = TimePure.getCPUSeconds . take i $ Math.Primality.carmichaelNumbers primalityAlgorithm + +-- | Measures the CPU-time required to determine whether the specified integer is prime, which is returned together with the Boolean result. +isPrimePerformance :: Math.Primality.Algorithm primalityAlgorithm => primalityAlgorithm -> Integer -> IO (Double, Bool) +isPrimePerformance primalityAlgorithm = TimePure.getCPUSeconds . Math.Primality.isPrime primalityAlgorithm + +{- | + * Measures the CPU-time required to determine whether /prime-indexed Fibonacci-numbers/ are actually /prime/. + + * CAVEAT: nothing is returned, since the result is printed ... and it never terminates. +-} +isPrimePerformanceGraph :: Math.Primality.Algorithm primalityAlgorithm => primalityAlgorithm -> IO () +isPrimePerformanceGraph primalityAlgorithm = mapM_ ( + \operand -> isPrimePerformance primalityAlgorithm operand >>= putStrLn . shows operand . showChar '\t' . (`shows` "") + ) Math.Fibonacci.primeIndexedFibonacci + diff --git a/src/Factory/Test/Performance/PrimeFactorisation.hs b/src/Factory/Test/Performance/PrimeFactorisation.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/Performance/PrimeFactorisation.hs @@ -0,0 +1,50 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Times functions exported by module "Math.PrimeFactorisation". +-} + +module Factory.Test.Performance.PrimeFactorisation( +-- * Functions + primeFactorsPerformance, + primeFactorsPerformanceGraph +) where + +import qualified Factory.Data.PrimeFactors as Data.PrimeFactors +import qualified Factory.Math.Fibonacci as Math.Fibonacci +import qualified Factory.Math.PrimeFactorisation as Math.PrimeFactorisation +import qualified ToolShed.TimePure as TimePure + +-- | Measures the CPU-time required to prime-factorise the specified integer, which is returned together with the resulting list of factors. +primeFactorsPerformance :: Math.PrimeFactorisation.Algorithm algorithm => algorithm -> Integer -> IO (Double, Data.PrimeFactors.Factors Integer Int) +primeFactorsPerformance algorithm = TimePure.getCPUSeconds . Math.PrimeFactorisation.primeFactors algorithm + +{- | + * Measure the CPU-time required by 'Math.PrimeFactorisation.primeFactors', + arbitrarily against the /Fibonacci/-numbers (which seemed to fit the requirements). + + * CAVEAT: nothing is returned, since the result is printed ... and it never terminates. +-} +primeFactorsPerformanceGraph :: Math.PrimeFactorisation.Algorithm algorithm => algorithm -> Int -> IO () +primeFactorsPerformanceGraph algorithm tests + | tests < 0 = error $ "Factory.Test.Performance.PrimeFactorisation.primeFactorsPerformanceGraph:\tnegative number; " ++ show tests + | otherwise = mapM_ ( + \operand -> primeFactorsPerformance algorithm operand >>= putStrLn . shows operand . showChar '\t' . (`shows` "") + ) . take tests . dropWhile (< 2) $ Math.Fibonacci.fibonacci + diff --git a/src/Factory/Test/Performance/SquareRoot.hs b/src/Factory/Test/Performance/SquareRoot.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/Performance/SquareRoot.hs @@ -0,0 +1,54 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Times functions exported from module "Math.SquareRoot". +-} + +module Factory.Test.Performance.SquareRoot( +-- * Functions + squareRootPerformance, + squareRootPerformanceGraph +) where + +import qualified Control.Arrow +import qualified Factory.Math.Precision as Math.Precision +import qualified Factory.Math.SquareRoot as Math.SquareRoot +import qualified ToolShed.TimePure as TimePure + +-- | Measures the CPU-time required by 'Math.SquareRoot.squareRootFrom', which is returned together with the approximate rational result. +squareRootPerformance :: (Math.SquareRoot.Algorithm algorithm, Real operand) => algorithm -> operand -> Math.Precision.DecimalDigits -> IO (Double, Math.SquareRoot.Result) +squareRootPerformance algorithm operand requiredDecimalDigits = TimePure.getCPUSeconds $ Math.SquareRoot.squareRoot algorithm requiredDecimalDigits operand + +{- | + * Measures the CPU-time required by 'Math.SquareRoot.squareRootFrom', and the resulting accuracy, + using the specified algorithm, to an exponentially increasing precision-requirement. + + * CAVEAT: nothing is returned, since the result is printed ... and it never terminates. +-} +squareRootPerformanceGraph :: ( + Math.SquareRoot.Algorithm algorithm, + Math.SquareRoot.Iterator algorithm, + Show algorithm, + Real operand + ) => algorithm -> operand -> IO () +squareRootPerformanceGraph algorithm operand = mapM_ ( + \requiredDecimalDigits -> putStrLn . ( + \(cpuSeconds, actualDecimalDigits) -> shows algorithm . showChar '\t' . shows requiredDecimalDigits . showChar '\t' . shows actualDecimalDigits . showChar '\t' $ shows cpuSeconds "" + ) . Control.Arrow.second (Math.SquareRoot.getAccuracy operand) =<< squareRootPerformance algorithm operand requiredDecimalDigits + ) $ iterate (* max 2 (Math.SquareRoot.convergenceOrder algorithm)) 16 diff --git a/src/Factory/Test/Performance/Statistics.hs b/src/Factory/Test/Performance/Statistics.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/Performance/Statistics.hs @@ -0,0 +1,40 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Times functions exported from module "Math.Statistics". +-} + +module Factory.Test.Performance.Statistics( +-- * Functions + nCrPerformance +) where + +import qualified Control.DeepSeq +import qualified Factory.Math.Factorial as Math.Factorial +import qualified Factory.Math.Statistics as Math.Statistics +import qualified ToolShed.TimePure as TimePure + +-- | Measures the CPU-time required by 'Math.Statistics.nCr'. +nCrPerformance :: (Math.Factorial.Algorithm factorialAlgorithm, Control.DeepSeq.NFData i, Integral i) + => factorialAlgorithm + -> i -- ^ The total number from which to select. + -> i -- ^ The number of items in a sample. + -> IO (Double, i) +nCrPerformance factorialAlgorithm n r = TimePure.getCPUSeconds $ Math.Statistics.nCr factorialAlgorithm n r + diff --git a/src/Factory/Test/QuickCheck/ArithmeticGeometricMean.hs b/src/Factory/Test/QuickCheck/ArithmeticGeometricMean.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/QuickCheck/ArithmeticGeometricMean.hs @@ -0,0 +1,65 @@ +{-# LANGUAGE CPP #-} +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines /QuickCheck/-properties for "Math.ArithmeticGeometricMean". +-} + +module Factory.Test.QuickCheck.ArithmeticGeometricMean( +-- * Types +-- ** Type-synonyms +-- Testable, +-- * Functions + quickChecks +) where + +import qualified Factory.Math.ArithmeticGeometricMean as Math.ArithmeticGeometricMean +import qualified Factory.Math.Implementations.SquareRoot as Math.Implementations.SquareRoot +import qualified Factory.Math.Precision as Math.Precision +import Factory.Test.QuickCheck.SquareRoot() +import qualified Test.QuickCheck +import Test.QuickCheck((==>)) + +#if MIN_VERSION_base(4,3,0) +import Data.Tuple(swap) +#else +-- | Swap the components of a pair. +swap :: (a, b) -> (b, a) +swap (a, b) = (b, a) +#endif + +type Testable = Math.Implementations.SquareRoot.Algorithm -> Math.Precision.DecimalDigits -> Math.ArithmeticGeometricMean.AGM -> Int -> Test.QuickCheck.Property + +-- | Defines invariant properties. +quickChecks :: IO () +quickChecks = Test.QuickCheck.quickCheck `mapM_` [prop_symmetrical, prop_bounds] where + prop_symmetrical, prop_bounds :: Testable + prop_symmetrical squareRootAlgorithm decimalDigits agm index = Math.ArithmeticGeometricMean.isValid agm ==> Test.QuickCheck.label "prop_symmetrical" . and . tail . take index' $ zipWith (==) ( + Math.ArithmeticGeometricMean.convergeToAGM squareRootAlgorithm decimalDigits' agm + ) ( + Math.ArithmeticGeometricMean.convergeToAGM squareRootAlgorithm decimalDigits' $ swap agm + ) where + decimalDigits' = 1 + (decimalDigits `mod` 64) + index' = 1 + (index `mod` 8) + + prop_bounds squareRootAlgorithm decimalDigits agm index = all ($ agm) [Math.ArithmeticGeometricMean.isValid, uncurry (/=)] ==> Test.QuickCheck.label "prop_bounds" . all (uncurry (>=)) . tail . take index' $ Math.ArithmeticGeometricMean.convergeToAGM squareRootAlgorithm decimalDigits' agm + where + decimalDigits' = 33 {-test is sensitive to rounding-errors-} + (decimalDigits `mod` 96) + index' = 1 + (index `mod` 5) + diff --git a/src/Factory/Test/QuickCheck/Bounds.hs b/src/Factory/Test/QuickCheck/Bounds.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/QuickCheck/Bounds.hs @@ -0,0 +1,43 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines /QuickCheck/-properties for "Data.Bounds". +-} + +module Factory.Test.QuickCheck.Bounds( +-- * Functions + quickChecks +) where + +import qualified Data.Ratio +import qualified Factory.Data.Bounds as Data.Bounds +import qualified Test.QuickCheck + +-- | Defines invariant properties. +quickChecks :: IO () +quickChecks = Test.QuickCheck.quickCheckWith Test.QuickCheck.stdArgs {Test.QuickCheck.maxSuccess = 1000} prop_product where + prop_product :: Data.Ratio.Ratio Integer -> Integer -> Data.Bounds.Bounds Integer -> Test.QuickCheck.Property + prop_product ratio minLength bounds = Test.QuickCheck.label "prop_product" $ Data.Bounds.product' ratio' minLength' bounds' == product (Data.Bounds.toList bounds') where + bounds' = Data.Bounds.normalise bounds + minLength' = 1 + minLength `mod` 1000 + ratio' = if r > 1 + then recip r + else r + where + r = abs ratio diff --git a/src/Factory/Test/QuickCheck/Factorial.hs b/src/Factory/Test/QuickCheck/Factorial.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/QuickCheck/Factorial.hs @@ -0,0 +1,72 @@ +{-# LANGUAGE CPP #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines /QuickCheck/-properties for "Math.Implementations.Factorial". +-} + +module Factory.Test.QuickCheck.Factorial( +-- * Types +-- ** Type-synonyms +-- Testable, +-- * Functions + quickChecks +) where + +import Data.Ratio((%)) +import qualified Factory.Math.Factorial as Math.Factorial +import qualified Factory.Math.Implementations.Factorial as Math.Implementations.Factorial +import Factory.Math.Implementations.Factorial((!/!)) +import qualified Test.QuickCheck +import Test.QuickCheck((==>)) + +instance Test.QuickCheck.Arbitrary Math.Implementations.Factorial.Algorithm where + arbitrary = Test.QuickCheck.elements [Math.Implementations.Factorial.Bisection, Math.Implementations.Factorial.PrimeFactorisation] +#if !(MIN_VERSION_QuickCheck(2,1,0)) + coarbitrary = undefined --CAVEAT: stops warnings from ghc. +#endif + +type Testable = Integer -> Integer -> Test.QuickCheck.Property + +-- | Defines invariant properties. +quickChecks :: IO () +quickChecks = Test.QuickCheck.quickCheck `mapM_` [prop_equivalence, prop_symmetry, prop_x0, prop_0n] >> Test.QuickCheck.quickCheck prop_ratio >> Test.QuickCheck.quickCheck prop_consistency where + prop_equivalence, prop_symmetry, prop_x0, prop_0n :: Testable + prop_equivalence x n = Test.QuickCheck.label "prop_equivalence" $ Math.Implementations.Factorial.risingFactorial x n == sign * Math.Implementations.Factorial.fallingFactorial (negate x) n && Math.Implementations.Factorial.fallingFactorial x n == sign * Math.Implementations.Factorial.risingFactorial (negate x) n where + sign :: Integer + sign + | odd n = negate 1 + | otherwise = 1 + + prop_symmetry x n = Test.QuickCheck.label "prop_symmetry" $ Math.Implementations.Factorial.risingFactorial x n == Math.Implementations.Factorial.fallingFactorial (x + n - 1) n + + prop_x0 x _ = Test.QuickCheck.label "prop_x0" $ all (== 1) $ map ($ 0) [Math.Implementations.Factorial.risingFactorial x, Math.Implementations.Factorial.fallingFactorial x] + + prop_0n _ n = Test.QuickCheck.label "prop_0n" $ all (== if n == 0 then 1 else 0) $ map ($ n) [Math.Implementations.Factorial.risingFactorial 0, Math.Implementations.Factorial.fallingFactorial 0] + + prop_ratio :: Math.Implementations.Factorial.Algorithm -> Integer -> Integer -> Test.QuickCheck.Property + prop_ratio algorithm i j = Test.QuickCheck.label "prop_ratio" $ n !/! d == Math.Factorial.factorial algorithm n % Math.Factorial.factorial algorithm d where + n = (i `mod` 100000) - 1 + d = (j `mod` 100000) - 1 + + prop_consistency :: Math.Implementations.Factorial.Algorithm -> Math.Implementations.Factorial.Algorithm -> Integer -> Test.QuickCheck.Property + prop_consistency l r i = l /= r ==> Test.QuickCheck.label "prop_consistency" $ Math.Factorial.factorial l n == Math.Factorial.factorial r n where + n = (i `mod` 100000) - 1 + diff --git a/src/Factory/Test/QuickCheck/MonicPolynomial.hs b/src/Factory/Test/QuickCheck/MonicPolynomial.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/QuickCheck/MonicPolynomial.hs @@ -0,0 +1,74 @@ +{-# LANGUAGE CPP #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Implements 'Test.QuickCheck.Arbitrary' and defines /QuickCheck/-properties for "Data.MonicPolynomial". +-} + +module Factory.Test.QuickCheck.MonicPolynomial( +-- * Types +-- ** Type-synonyms +-- P +-- * Functions + quickChecks +) where + +import Factory.Data.Ring((=*=), (=+=), (=^)) +import Factory.Test.QuickCheck.Polynomial() +import qualified Factory.Data.MonicPolynomial as Data.MonicPolynomial +import qualified Factory.Data.Polynomial as Data.Polynomial +import qualified Factory.Data.QuotientRing as Data.QuotientRing +import qualified Factory.Data.Ring as Data.Ring +import qualified Test.QuickCheck + +instance ( + Test.QuickCheck.Arbitrary c, + Integral c, + Test.QuickCheck.Arbitrary e, + Integral e + ) => Test.QuickCheck.Arbitrary (Data.MonicPolynomial.MonicPolynomial c e) where + arbitrary = do + polynomial <- Test.QuickCheck.arbitrary + + return . Data.MonicPolynomial.mkMonicPolynomial $ ((1, Data.Polynomial.getDegree polynomial + 1) :) `Data.Polynomial.lift` polynomial +#if !(MIN_VERSION_QuickCheck(2,1,0)) + coarbitrary = undefined --CAVEAT: stops warnings from ghc. +#endif + +type P = Data.MonicPolynomial.MonicPolynomial Integer Integer + +-- | Defines invariant properties. +quickChecks :: IO () +quickChecks = Test.QuickCheck.quickCheck `mapM_` [prop_quotRem, prop_quotientRingNormalised] >> Test.QuickCheck.quickCheck prop_perfectPower >> Test.QuickCheck.quickCheck prop_isDivisibleBy where + prop_quotRem, prop_quotientRingNormalised :: P -> P -> Test.QuickCheck.Property + prop_quotRem numerator denominator = Test.QuickCheck.label "prop_quotRem" $ numerator == denominator =*= quotient =+= remainder where + (quotient, remainder) = numerator `Data.QuotientRing.quotRem'` denominator + + prop_quotientRingNormalised numerator denominator = Test.QuickCheck.label "prop_quotientRingNormalised" $ all (Data.Polynomial.isNormalised . Data.MonicPolynomial.getPolynomial) [numerator `Data.QuotientRing.quot'` denominator, numerator `Data.QuotientRing.rem'` denominator] + + prop_perfectPower :: P -> Int -> Test.QuickCheck.Property + prop_perfectPower polynomial power = Test.QuickCheck.label "prop_perfectPower" $ iterate (`Data.QuotientRing.quot'` polynomial) (polynomial =^ power') !! pred power' == polynomial where + power' :: Int + power' = 1 + power `mod` 100 + + prop_isDivisibleBy :: [P] -> Test.QuickCheck.Property + prop_isDivisibleBy monicPolynomials = Test.QuickCheck.label "prop_isDivisibleBy" $ all (Data.QuotientRing.isDivisibleBy (Data.Ring.product' (recip 2) {-TODO-} 10 monicPolynomials)) monicPolynomials + + diff --git a/src/Factory/Test/QuickCheck/Pi.hs b/src/Factory/Test/QuickCheck/Pi.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/QuickCheck/Pi.hs @@ -0,0 +1,132 @@ +{-# LANGUAGE CPP #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines /QuickCheck/-properties for "Math.Pi". +-} + +module Factory.Test.QuickCheck.Pi( +-- * Types +-- ** Type-synonyms +-- Testable, +-- * Functions + quickChecks +) where + +import Control.Applicative((<$>), (<*>)) +import Factory.Test.QuickCheck.Factorial() +import Factory.Test.QuickCheck.SquareRoot() +import qualified Factory.Math.Factorial as Math.Factorial +import qualified Factory.Math.Implementations.Factorial as Math.Implementations.Factorial +import qualified Factory.Math.Implementations.Pi.AGM.Algorithm as Math.Implementations.Pi.AGM.Algorithm +import qualified Factory.Math.Implementations.Pi.BBP.Algorithm as Math.Implementations.Pi.BBP.Algorithm +import qualified Factory.Math.Implementations.Pi.Borwein.Algorithm as Math.Implementations.Pi.Borwein.Algorithm +import qualified Factory.Math.Implementations.Pi.Ramanujan.Algorithm as Math.Implementations.Pi.Ramanujan.Algorithm +import qualified Factory.Math.Implementations.Pi.Spigot.Algorithm as Math.Implementations.Pi.Spigot.Algorithm +import qualified Factory.Math.Implementations.SquareRoot as Math.Implementations.SquareRoot +import qualified Factory.Math.Pi as Math.Pi +import qualified Factory.Math.Precision as Math.Precision +import qualified Factory.Math.SquareRoot as Math.SquareRoot +import qualified Test.QuickCheck +import Test.QuickCheck((==>)) + +instance ( + Test.QuickCheck.Arbitrary squareRootAlgorithm, + Math.SquareRoot.Algorithm squareRootAlgorithm + ) => Test.QuickCheck.Arbitrary (Math.Implementations.Pi.AGM.Algorithm.Algorithm squareRootAlgorithm) where + arbitrary = Math.Implementations.Pi.AGM.Algorithm.BrentSalamin <$> Test.QuickCheck.arbitrary +#if !(MIN_VERSION_QuickCheck(2,1,0)) + coarbitrary = undefined --CAVEAT: stops warnings from ghc. +#endif + +instance Test.QuickCheck.Arbitrary Math.Implementations.Pi.BBP.Algorithm.Algorithm where + arbitrary = Test.QuickCheck.elements [Math.Implementations.Pi.BBP.Algorithm.Bellard, Math.Implementations.Pi.BBP.Algorithm.Base65536] +#if !(MIN_VERSION_QuickCheck(2,1,0)) + coarbitrary = undefined --CAVEAT: stops warnings from ghc. +#endif + +instance ( + Test.QuickCheck.Arbitrary squareRootAlgorithm, + Math.SquareRoot.Algorithm squareRootAlgorithm, + Test.QuickCheck.Arbitrary factorialAlgorithm, + Math.Factorial.Algorithm factorialAlgorithm + ) => Test.QuickCheck.Arbitrary (Math.Implementations.Pi.Borwein.Algorithm.Algorithm squareRootAlgorithm factorialAlgorithm) where + arbitrary = Test.QuickCheck.oneof [ + Math.Implementations.Pi.Borwein.Algorithm.Borwein1993 <$> Test.QuickCheck.arbitrary <*> Test.QuickCheck.arbitrary + ] +#if !(MIN_VERSION_QuickCheck(2,1,0)) + coarbitrary = undefined --CAVEAT: stops warnings from ghc. +#endif + +instance ( + Test.QuickCheck.Arbitrary squareRootAlgorithm, + Math.SquareRoot.Algorithm squareRootAlgorithm, + Test.QuickCheck.Arbitrary factorialAlgorithm, + Math.Factorial.Algorithm factorialAlgorithm + ) => Test.QuickCheck.Arbitrary (Math.Implementations.Pi.Ramanujan.Algorithm.Algorithm squareRootAlgorithm factorialAlgorithm) where + arbitrary = Test.QuickCheck.oneof [ + Math.Implementations.Pi.Ramanujan.Algorithm.Classic <$> Test.QuickCheck.arbitrary <*> Test.QuickCheck.arbitrary, + Math.Implementations.Pi.Ramanujan.Algorithm.Chudnovsky <$> Test.QuickCheck.arbitrary <*> Test.QuickCheck.arbitrary + ] +#if !(MIN_VERSION_QuickCheck(2,1,0)) + coarbitrary = undefined --CAVEAT: stops warnings from ghc. +#endif + +instance Test.QuickCheck.Arbitrary Math.Implementations.Pi.Spigot.Algorithm.Algorithm where + arbitrary = Test.QuickCheck.elements [Math.Implementations.Pi.Spigot.Algorithm.RabinowitzWagon, Math.Implementations.Pi.Spigot.Algorithm.Gosper] +#if !(MIN_VERSION_QuickCheck(2,1,0)) + coarbitrary = undefined --CAVEAT: stops warnings from ghc. +#endif + +instance ( + Test.QuickCheck.Arbitrary agm, + Test.QuickCheck.Arbitrary bbp, + Test.QuickCheck.Arbitrary borwein, + Test.QuickCheck.Arbitrary ramanujan, + Test.QuickCheck.Arbitrary spigot + ) => Test.QuickCheck.Arbitrary (Math.Pi.Category agm bbp borwein ramanujan spigot) where + arbitrary = Test.QuickCheck.oneof [ + Math.Pi.AGM <$> Test.QuickCheck.arbitrary, + Math.Pi.BBP <$> Test.QuickCheck.arbitrary, + Math.Pi.Borwein <$> Test.QuickCheck.arbitrary, + Math.Pi.Ramanujan <$> Test.QuickCheck.arbitrary, + Math.Pi.Spigot <$> Test.QuickCheck.arbitrary + ] +#if !(MIN_VERSION_QuickCheck(2,1,0)) + coarbitrary = undefined --CAVEAT: stops warnings from ghc. +#endif + +type Category = Math.Pi.Category ( + Math.Implementations.Pi.AGM.Algorithm.Algorithm Math.Implementations.SquareRoot.Algorithm + ) Math.Implementations.Pi.BBP.Algorithm.Algorithm ( + Math.Implementations.Pi.Borwein.Algorithm.Algorithm Math.Implementations.SquareRoot.Algorithm Math.Implementations.Factorial.Algorithm + ) ( + Math.Implementations.Pi.Ramanujan.Algorithm.Algorithm Math.Implementations.SquareRoot.Algorithm Math.Implementations.Factorial.Algorithm + ) Math.Implementations.Pi.Spigot.Algorithm.Algorithm + +type Testable = Category -> Category -> Math.Precision.DecimalDigits -> Test.QuickCheck.Property + +-- | Defines invariant properties. +quickChecks :: IO () +quickChecks = Test.QuickCheck.quickCheck prop_consistency where + prop_consistency :: Testable + prop_consistency l r decimalDigits = l /= r ==> Test.QuickCheck.label "prop_consistency" $ Math.Pi.openI l decimalDigits' - Math.Pi.openI r decimalDigits' <= 1 {-rounding error-} where + decimalDigits' = 1 + (decimalDigits `mod` 250) + diff --git a/src/Factory/Test/QuickCheck/Polynomial.hs b/src/Factory/Test/QuickCheck/Polynomial.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/QuickCheck/Polynomial.hs @@ -0,0 +1,122 @@ +{-# LANGUAGE CPP #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Implements 'Test.QuickCheck.Arbitrary' and defines /QuickCheck/-properties for "Data.Polynomial". +-} + +module Factory.Test.QuickCheck.Polynomial( +-- * Functions + quickChecks +) where + +import Control.Applicative((<$>)) +import Control.Arrow((***)) +import Factory.Data.Ring((=*=), (=+=), (=-=), (=^)) +import qualified Data.Numbers.Primes +import qualified Data.Ratio +import qualified Factory.Data.Polynomial as Data.Polynomial +import qualified Factory.Data.QuotientRing as Data.QuotientRing +import qualified Factory.Data.Ring as Data.Ring +import qualified Test.QuickCheck +import Test.QuickCheck((==>)) + +instance ( + Test.QuickCheck.Arbitrary c, + Integral c, + Test.QuickCheck.Arbitrary e, + Integral e + ) => Test.QuickCheck.Arbitrary (Data.Polynomial.Polynomial c e) where + arbitrary = Data.Polynomial.mkPolynomial . map ((+ negate 4) . (`mod` 8) *** (`mod` 8)) <$> Test.QuickCheck.arbitrary +#if !(MIN_VERSION_QuickCheck(2,1,0)) + coarbitrary = undefined --CAVEAT: stops warnings from ghc. +#endif + +-- | Defines invariant properties. +quickChecks :: IO () +quickChecks + = Test.QuickCheck.quickCheck prop_congruence + >> Test.QuickCheck.quickCheck `mapM_` [prop_quotRem, prop_degree, prop_ringNormalised, prop_quotientRingNormalised] + >> Test.QuickCheck.quickCheck `mapM_` [prop_power, prop_perfectPower, prop_normalised] + >> Test.QuickCheck.quickCheck prop_raiseModuloNormalised + >> Test.QuickCheck.quickCheck `mapM_` [prop_integralDomain, prop_isDivisibleBy] + where + prop_congruence :: Int -> Test.QuickCheck.Property + prop_congruence i = Test.QuickCheck.label "prop_congruence" $ Data.Polynomial.areCongruentModulo (Data.Polynomial.mkLinear 1 (negate 1) =^ prime) (Data.Polynomial.mkPolynomial [(1, prime), (negate 1, 0)]) prime where + prime :: Integer + prime = Data.Numbers.Primes.primes !! mod i 100 + + prop_quotRem, prop_degree, prop_ringNormalised, prop_quotientRingNormalised :: Data.Polynomial.Polynomial Integer Integer -> Data.Polynomial.Polynomial Integer Integer -> Test.QuickCheck.Property + prop_quotRem numerator denominator = denominator' /= Data.Polynomial.zero ==> Test.QuickCheck.label "prop_quotRem" $ numerator' == denominator' =*= quotient =+= remainder where + numerator', denominator' :: Data.Polynomial.Polynomial Data.Ratio.Rational Integer + numerator' = Data.Polynomial.realCoefficientsToFrac numerator + denominator' = Data.Polynomial.realCoefficientsToFrac denominator + + (quotient, remainder) = numerator' `Data.QuotientRing.quotRem'` denominator' + + prop_degree numerator denominator = denominator' /= Data.Polynomial.zero ==> Test.QuickCheck.label "prop_degree" $ remainder == Data.Polynomial.zero || Data.Polynomial.getDegree remainder < Data.Polynomial.getDegree denominator' where + numerator', denominator' :: Data.Polynomial.Polynomial Data.Ratio.Rational Integer + numerator' = Data.Polynomial.realCoefficientsToFrac numerator + denominator' = Data.Polynomial.realCoefficientsToFrac denominator + + remainder = numerator' `Data.QuotientRing.rem'` denominator' + + prop_ringNormalised l r = Test.QuickCheck.label "prop_ringNormalised" $ all Data.Polynomial.isNormalised [l =*= r, l =+= r, l =-= r] + + prop_quotientRingNormalised numerator denominator = denominator' /= Data.Polynomial.zero ==> Test.QuickCheck.label "prop_quotientRingNormalised" $ all Data.Polynomial.isNormalised [numerator' `Data.QuotientRing.quot'` denominator', numerator' `Data.QuotientRing.rem'` denominator'] where + numerator', denominator' :: Data.Polynomial.Polynomial Data.Ratio.Rational Integer + numerator' = Data.Polynomial.realCoefficientsToFrac numerator + denominator' = Data.Polynomial.realCoefficientsToFrac denominator + + prop_power, prop_perfectPower, prop_normalised :: Data.Polynomial.Polynomial Integer Integer -> Int -> Test.QuickCheck.Property + prop_power polynomial power = Test.QuickCheck.label "prop_power" $ polynomial =^ power' == iterate (=*= polynomial) polynomial !! pred power' where + power' :: Int + power' = 1 + power `mod` 100 + + prop_perfectPower polynomial power = polynomial' /= Data.Polynomial.zero ==> Test.QuickCheck.label "prop_perfectPower" $ iterate (`Data.QuotientRing.quot'` polynomial') (polynomial' =^ power') !! pred power' == polynomial' where + polynomial' :: Data.Polynomial.Polynomial Data.Ratio.Rational Integer + polynomial' = Data.Polynomial.realCoefficientsToFrac polynomial + + power' :: Int + power' = 1 + power `mod` 100 + + prop_normalised polynomial i = Test.QuickCheck.label "prop_normalised" $ all Data.Polynomial.isNormalised [ + polynomial =^ power', + polynomial `Data.Polynomial.mod'` modulus' + ] where + power' :: Int + power' = 1 + i `mod` 100 + + modulus' :: Integer + modulus' = 1 + fromIntegral i `mod` 100 + + prop_raiseModuloNormalised :: Data.Polynomial.Polynomial Integer Integer -> Integer -> Integer -> Test.QuickCheck.Property + prop_raiseModuloNormalised polynomial power modulus = Test.QuickCheck.label "prop_raiseModuloNormalised" . Data.Polynomial.isNormalised $ Data.Polynomial.raiseModulo polynomial power' modulus' where + power', modulus' :: Integer + power' = 1 + power `mod` 100 + modulus' = 1 + modulus `mod` 100 + + prop_integralDomain, prop_isDivisibleBy :: [Data.Polynomial.Polynomial Integer Integer] -> Test.QuickCheck.Property + prop_integralDomain polynomials = Data.Polynomial.zero `notElem` polynomials ==> Test.QuickCheck.label "prop_integralDomain" $ Data.Ring.product' (recip 2) {-TODO-} 10 polynomials /= Data.Polynomial.zero + + prop_isDivisibleBy polynomials = Test.QuickCheck.label "prop_isDivisibleBy" . all (Data.QuotientRing.isDivisibleBy (Data.Ring.product' (recip 2) {-TODO-} 10 polynomials')) $ filter (/= Data.Polynomial.zero) polynomials' where + polynomials' :: [Data.Polynomial.Polynomial Data.Ratio.Rational Integer] + polynomials' = map Data.Polynomial.realCoefficientsToFrac polynomials + diff --git a/src/Factory/Test/QuickCheck/Power.hs b/src/Factory/Test/QuickCheck/Power.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/QuickCheck/Power.hs @@ -0,0 +1,53 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines /QuickCheck/-properties for "Math.Power". +-} + +module Factory.Test.QuickCheck.Power( +-- * Functions + quickChecks +) where + +import qualified Data.List +import qualified Factory.Math.Power as Math.Power +import qualified Test.QuickCheck +import Test.QuickCheck((==>)) + +-- | Defines invariant properties. +quickChecks :: IO () +quickChecks = + Test.QuickCheck.quickCheck prop_maybeSquareNumber + >> Test.QuickCheck.quickCheckWith Test.QuickCheck.stdArgs {Test.QuickCheck.maxSuccess = 10000} prop_notSquare + >> Test.QuickCheck.quickCheck prop_squaresFrom + >> Test.QuickCheck.quickCheck prop_raiseModulo + where + prop_maybeSquareNumber, prop_notSquare :: Integer -> Test.QuickCheck.Property + prop_maybeSquareNumber i = Test.QuickCheck.label "prop_maybeSquareNumber" $ Math.Power.maybeSquareNumber (Math.Power.square i) == Just (abs i) + + prop_notSquare i = abs i > 0 ==> Test.QuickCheck.label "prop_notSquare" $ Math.Power.maybeSquareNumber (i ^ (10 {-promote rounding-error using big number-} :: Int) + 1) == Nothing + + prop_squaresFrom :: Integer -> Integer -> Test.QuickCheck.Property + prop_squaresFrom from l = Test.QuickCheck.label "prop_squaresFrom" . (\(x, y) -> y == Math.Power.square x) . Data.List.genericIndex (Math.Power.squaresFrom from) $ abs l + + prop_raiseModulo :: Integer -> Integer -> Integer -> Test.QuickCheck.Property + prop_raiseModulo b e m = m /= 0 ==> Test.QuickCheck.label "prop_raiseModulo" $ Math.Power.raiseModulo b e' m == (b ^ e') `mod` m where + e' :: Integer + e' = abs e + diff --git a/src/Factory/Test/QuickCheck/Primality.hs b/src/Factory/Test/QuickCheck/Primality.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/QuickCheck/Primality.hs @@ -0,0 +1,77 @@ +{-# LANGUAGE CPP #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Implements 'Test.QuickCheck.Arbitrary' and defines /QuickCheck/-properties for "Math.Implementations.Primality". +-} + +module Factory.Test.QuickCheck.Primality( +-- * Functions + quickChecks +) where + +import Control.Applicative((<$>)) +import Factory.Test.QuickCheck.PrimeFactorisation() +import qualified Data.List +import qualified Data.Numbers.Primes +import qualified Factory.Math.Implementations.Primality as Math.Implementations.Primality +import qualified Factory.Math.Implementations.PrimeFactorisation as Math.Implementations.PrimeFactorisation +import qualified Factory.Math.Primality as Math.Primality +import qualified Test.QuickCheck +import Test.QuickCheck((==>)) + +instance Test.QuickCheck.Arbitrary factorisationAlgorithm => Test.QuickCheck.Arbitrary (Math.Implementations.Primality.Algorithm factorisationAlgorithm) where + arbitrary = Test.QuickCheck.oneof [ + Math.Implementations.Primality.AKS <$> Test.QuickCheck.arbitrary, + return Math.Implementations.Primality.MillerRabin + ] +#if !(MIN_VERSION_QuickCheck(2,1,0)) + coarbitrary = undefined --CAVEAT: stops warnings from ghc. +#endif + +-- | Defines invariant properties. +quickChecks :: IO () +quickChecks + = Test.QuickCheck.quickCheck prop_prime + >> Test.QuickCheck.quickCheck prop_composite + >> Test.QuickCheck.quickCheck prop_consistency + where + prop_prime :: Math.Implementations.Primality.Algorithm Math.Implementations.PrimeFactorisation.Algorithm -> Integer -> Test.QuickCheck.Property + prop_prime primalityAlgorithm i = Test.QuickCheck.label "prop_prime" $ Math.Primality.isPrime primalityAlgorithm prime where + normalise n + | primalityAlgorithm == Math.Implementations.Primality.MillerRabin = n `mod` 1000000 --Limited by the efficiency of 'Data.Numbers.Primes.primes'. + | otherwise = n `mod` 59 + + prime :: Integer + prime = Data.List.genericIndex Data.Numbers.Primes.primes $ normalise i + + prop_composite :: Math.Implementations.Primality.Algorithm Math.Implementations.PrimeFactorisation.Algorithm -> [Integer] -> Test.QuickCheck.Property + prop_composite primalityAlgorithm l = length l > 1 ==> Test.QuickCheck.label "prop_composite" . not $ Math.Primality.isPrime primalityAlgorithm composite where + normalise n + | primalityAlgorithm == Math.Implementations.Primality.MillerRabin = n `mod` 1000000 + | otherwise = n `mod` 10 + + composite :: Integer + composite = product . map (Data.List.genericIndex Data.Numbers.Primes.primes . normalise) $ take 8 l + + prop_consistency :: Math.Implementations.Primality.Algorithm Math.Implementations.PrimeFactorisation.Algorithm -> Math.Implementations.Primality.Algorithm Math.Implementations.PrimeFactorisation.Algorithm -> Integer -> Test.QuickCheck.Property + prop_consistency l r i = l /= r ==> Test.QuickCheck.label "prop_consistency" $ Math.Primality.isPrime l i' == Math.Primality.isPrime r i' where + i' = i `mod` 512 + diff --git a/src/Factory/Test/QuickCheck/PrimeFactorisation.hs b/src/Factory/Test/QuickCheck/PrimeFactorisation.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/QuickCheck/PrimeFactorisation.hs @@ -0,0 +1,98 @@ +{-# LANGUAGE CPP #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Implements 'Test.QuickCheck.Arbitrary' and defines /QuickCheck/-properties for "Math.PrimeFactorisation". +-} + +module Factory.Test.QuickCheck.PrimeFactorisation( +-- * Functions + quickChecks +) where + +import qualified Data.List +import qualified Data.Numbers.Primes +import qualified Factory.Data.PrimeFactors as Data.PrimeFactors +import qualified Factory.Data.Exponential as Data.Exponential +import qualified Factory.Math.Implementations.PrimeFactorisation as Math.Implementations.PrimeFactorisation +import qualified Factory.Math.MultiplicativeOrder as Math.MultiplicativeOrder +import qualified Factory.Math.PrimeFactorisation as Math.PrimeFactorisation +import qualified Test.QuickCheck +import Test.QuickCheck((==>)) + +instance Test.QuickCheck.Arbitrary Math.Implementations.PrimeFactorisation.Algorithm where + arbitrary = Test.QuickCheck.oneof [ + Test.QuickCheck.elements [ + Math.Implementations.PrimeFactorisation.TrialDivision, + Math.Implementations.PrimeFactorisation.FermatsMethod + ] + ] +#if !(MIN_VERSION_QuickCheck(2,1,0)) + coarbitrary = undefined --CAVEAT: stops warnings from ghc. +#endif + +-- | Defines invariant properties. +quickChecks :: IO () +quickChecks = + Test.QuickCheck.quickCheck prop_consistency + >> Test.QuickCheck.quickCheck `mapM_` [prop_primeFactors, prop_smoothness, prop_eulersTotientP, prop_eulersTotientInequality] + >> Test.QuickCheck.quickCheck `mapM_` [prop_eulersTotient, prop_lagrange, prop_multiplicativeOrder, prop_perfectPower] + where + prop_consistency :: Integer -> Test.QuickCheck.Property + prop_consistency i = Test.QuickCheck.label "prop_consistency" $ (Math.PrimeFactorisation.primeFactors Math.Implementations.PrimeFactorisation.TrialDivision i' :: Data.PrimeFactors.Factors Integer Int) == Math.PrimeFactorisation.primeFactors Math.Implementations.PrimeFactorisation.FermatsMethod i' where + i' :: Integer + i' = 1 + (i `mod` 1000000) + + prop_primeFactors, prop_smoothness, prop_eulersTotientP, prop_eulersTotientInequality :: Math.Implementations.PrimeFactorisation.Algorithm -> Integer -> Test.QuickCheck.Property + prop_primeFactors algorithm i = Test.QuickCheck.label "prop_primeFactors" $ Data.PrimeFactors.product' (recip 2) {-TODO-} 10 (Math.PrimeFactorisation.primeFactors algorithm i') == i' where + i' :: Integer + i' = 1 + (i `mod` 1000000) + + prop_smoothness algorithm i = Test.QuickCheck.label "prop_smoothness" $ (Math.PrimeFactorisation.smoothness algorithm !! (2 ^ i')) <= (2 :: Integer) where + i' :: Integer + i' = i `mod` 20 + + prop_eulersTotientP algorithm i = Test.QuickCheck.label "prop_eulersTotient" $ Math.PrimeFactorisation.eulersTotient algorithm prime == prime - 1 where + prime :: Integer + prime = Data.List.genericIndex Data.Numbers.Primes.primes (i `mod` 10000) + + prop_eulersTotientInequality algorithm i = i `notElem` [2, 6] ==> Test.QuickCheck.label "prop_eulersTotientInequality" $ Math.PrimeFactorisation.eulersTotient algorithm i' >= floor (sqrt $ fromIntegral i' :: Double) where + i' = 1 + (i `mod` 100000) + + prop_eulersTotient, prop_lagrange, prop_multiplicativeOrder, prop_perfectPower :: Math.Implementations.PrimeFactorisation.Algorithm -> Integer -> Integer -> Test.QuickCheck.Property + prop_eulersTotient algorithm i power = Test.QuickCheck.label "prop_eulersTotient" $ Math.PrimeFactorisation.eulersTotient algorithm (base ^ power') == (base ^ (power' - 1)) * (base - 1) where + base :: Integer + base = Data.List.genericIndex Data.Numbers.Primes.primes (i `mod` 8) + + power' = 1 + (power `mod` 5) + + prop_lagrange algorithm base modulus = gcd base modulus' == 1 ==> Test.QuickCheck.label "prop_lagrange" $ (Math.PrimeFactorisation.eulersTotient algorithm modulus' `rem` Math.MultiplicativeOrder.multiplicativeOrder algorithm base modulus') == 0 where + modulus' :: Integer + modulus' = 2 + abs modulus + + prop_multiplicativeOrder algorithm base modulus = gcd base modulus' == 1 ==> Test.QuickCheck.label "prop_multiplicativeOrder" $ ( + base ^ Math.MultiplicativeOrder.multiplicativeOrder algorithm base modulus' + ) `mod` modulus' == 1 where + modulus' :: Integer + modulus' = 2 + abs modulus + + prop_perfectPower algorithm b e = Test.QuickCheck.label "prop_perfectPower" $ foldr1 gcd ( + map Data.Exponential.getExponent . Math.PrimeFactorisation.primeFactors algorithm $ (2 + b `mod` 10 :: Integer) ^ (2 + e `mod` 5) + ) > 1 diff --git a/src/Factory/Test/QuickCheck/QuickChecks.hs b/src/Factory/Test/QuickCheck/QuickChecks.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/QuickCheck/QuickChecks.hs @@ -0,0 +1,58 @@ +{-# LANGUAGE CPP #-} +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Calls the /quickChecks/-functions for modules supporting this feature. +-} + +module Factory.Test.QuickCheck.QuickChecks( +-- * Functions + run +) where + +import qualified Factory.Test.QuickCheck.ArithmeticGeometricMean +import qualified Factory.Test.QuickCheck.Bounds +import qualified Factory.Test.QuickCheck.Factorial +import qualified Factory.Test.QuickCheck.MonicPolynomial +import qualified Factory.Test.QuickCheck.Pi +import qualified Factory.Test.QuickCheck.Polynomial +import qualified Factory.Test.QuickCheck.Power +import qualified Factory.Test.QuickCheck.Primality +import qualified Factory.Test.QuickCheck.PrimeFactorisation +import qualified Factory.Test.QuickCheck.Radix +import qualified Factory.Test.QuickCheck.SquareRoot +import qualified Factory.Test.QuickCheck.Statistics +import qualified Factory.Test.QuickCheck.Summation + +-- | Run the /quickChecks/-functions for modules supporting this feature. +run :: IO () +run = putStrLn "ArithmeticGeometricMean" >> Factory.Test.QuickCheck.ArithmeticGeometricMean.quickChecks + >> putStrLn "Bounds" >> Factory.Test.QuickCheck.Bounds.quickChecks + >> putStrLn "Factorial" >> Factory.Test.QuickCheck.Factorial.quickChecks + >> putStrLn "MonicPolynomial" >> Factory.Test.QuickCheck.MonicPolynomial.quickChecks + >> putStrLn "Pi" >> Factory.Test.QuickCheck.Pi.quickChecks + >> putStrLn "Polynomial" >> Factory.Test.QuickCheck.Polynomial.quickChecks + >> putStrLn "Power" >> Factory.Test.QuickCheck.Power.quickChecks + >> putStrLn "Primality" >> Factory.Test.QuickCheck.Primality.quickChecks + >> putStrLn "PrimeFactorisation" >> Factory.Test.QuickCheck.PrimeFactorisation.quickChecks + >> putStrLn "Radix" >> Factory.Test.QuickCheck.Radix.quickChecks + >> putStrLn "SquareRoot" >> Factory.Test.QuickCheck.SquareRoot.quickChecks + >> putStrLn "Statistics" >> Factory.Test.QuickCheck.Statistics.quickChecks + >> putStrLn "Summation" >> Factory.Test.QuickCheck.Summation.quickChecks + diff --git a/src/Factory/Test/QuickCheck/Radix.hs b/src/Factory/Test/QuickCheck/Radix.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/QuickCheck/Radix.hs @@ -0,0 +1,46 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines /QuickCheck/-properties for "Math.Radix". +-} + +module Factory.Test.QuickCheck.Radix( +-- * Types +-- ** Type-synonyms +-- Testable, +-- * Functions + quickChecks +) where + +import qualified Factory.Math.Radix as Math.Radix +import qualified Test.QuickCheck +import Test.QuickCheck((==>)) + +type Testable = (Int, Integer) -> Test.QuickCheck.Property + +-- | Defines invariant properties. +quickChecks :: IO () +quickChecks = Test.QuickCheck.quickCheck `mapM_` [prop_reversable, prop_digitalRoot] where + prop_reversable, prop_digitalRoot :: Testable + prop_reversable (b, n) = abs base > 1 ==> Test.QuickCheck.label "prop_reversable" $ Math.Radix.fromBase base (Math.Radix.toBase base n) == n where + base = (b `mod` 73) - 36 + + prop_digitalRoot (_, n) = Test.QuickCheck.label "prop_digitalRoot" $ Math.Radix.digitalRoot n' == 9 where + n' = 9 * (1 + abs n) + diff --git a/src/Factory/Test/QuickCheck/SquareRoot.hs b/src/Factory/Test/QuickCheck/SquareRoot.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/QuickCheck/SquareRoot.hs @@ -0,0 +1,91 @@ +{-# LANGUAGE CPP #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Implements 'Test.QuickCheck.Arbitrary' and defines /QuickCheck/-properties for "Math.Implementations.SquareRoot". +-} + +module Factory.Test.QuickCheck.SquareRoot( +-- * Types +-- ** Type-synonyms +-- Testable, +-- * Functions + quickChecks +) where + +import Control.Applicative((<$>)) +import Data.Ratio((%)) +import qualified Data.Ratio +import qualified Factory.Math.Implementations.SquareRoot as Math.Implementations.SquareRoot +import qualified Factory.Math.Power as Math.Power +import qualified Factory.Math.Precision as Math.Precision +import qualified Factory.Math.SquareRoot as Math.SquareRoot +import qualified Test.QuickCheck + +instance Test.QuickCheck.Arbitrary (Math.Implementations.SquareRoot.Algorithm) where + arbitrary = Test.QuickCheck.oneof [ + Test.QuickCheck.elements [ + Math.Implementations.SquareRoot.BakhshaliApproximation, + Math.Implementations.SquareRoot.ContinuedFraction, + Math.Implementations.SquareRoot.HalleysMethod, + Math.Implementations.SquareRoot.NewtonRaphsonIteration + ], + Math.Implementations.SquareRoot.TaylorSeries <$> Test.QuickCheck.elements [2 .. 32] + ] +#if !(MIN_VERSION_QuickCheck(2,1,0)) + coarbitrary = undefined --CAVEAT: stops warnings from ghc. +#endif + +type Testable = (Math.Implementations.SquareRoot.Algorithm, Math.Precision.DecimalDigits, Data.Ratio.Rational) -> Test.QuickCheck.Property + +-- | Defines invariant properties. +quickChecks :: IO () +quickChecks = Test.QuickCheck.quickCheck `mapM_` [prop_accuracy, prop_factorable, prop_perfectSquare] where + prop_accuracy, prop_factorable, prop_perfectSquare :: Testable + prop_accuracy (algorithm, decimalDigits, operand) = Test.QuickCheck.label "prop_accuracy" . (>= requiredDecimalDigits) . Math.SquareRoot.getAccuracy operand' $ Math.SquareRoot.squareRoot algorithm requiredDecimalDigits operand' where + requiredDecimalDigits :: Math.Precision.DecimalDigits + requiredDecimalDigits = 1 + (decimalDigits `mod` 1024) + + operand' :: Data.Ratio.Rational + operand' = abs operand + + prop_factorable (algorithm, decimalDigits, operand) = Test.QuickCheck.label "prop_factorable" . (<= 5) . ( + * 10 ^ requiredDecimalDigits --Promote the relative error. + ) . abs $ 1 - ( + Math.SquareRoot.squareRoot algorithm requiredDecimalDigits ( + toRational $ Data.Ratio.numerator operand' + ) / Math.SquareRoot.squareRoot algorithm requiredDecimalDigits ( + toRational $ Data.Ratio.denominator operand' + ) + ) / Math.SquareRoot.squareRoot algorithm requiredDecimalDigits operand' where + requiredDecimalDigits :: Math.Precision.DecimalDigits + requiredDecimalDigits = 1 + (decimalDigits `mod` 1024) + + operand' :: Data.Ratio.Rational + operand' = 1 + abs operand + + prop_perfectSquare (algorithm, decimalDigits, operand) = Test.QuickCheck.label "prop_perfectSquare" . Math.SquareRoot.isPrecise perfectSquare $ Math.SquareRoot.squareRoot algorithm requiredDecimalDigits perfectSquare where + requiredDecimalDigits :: Math.Precision.DecimalDigits + requiredDecimalDigits = 1 + (decimalDigits `mod` 32768) + + operand', perfectSquare :: Data.Ratio.Rational + operand' = (abs (Data.Ratio.numerator operand) `min` (2 ^ (32 :: Int))) % (abs (Data.Ratio.denominator operand) `min` (2 ^ (32 :: Int))) --Avoid floating-point rounding-errors in 'Math.SquareRoot.rSqrt'. + perfectSquare = Math.Power.square operand' + diff --git a/src/Factory/Test/QuickCheck/Statistics.hs b/src/Factory/Test/QuickCheck/Statistics.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/QuickCheck/Statistics.hs @@ -0,0 +1,68 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines /QuickCheck/-properties for "Math.Statistics". +-} + +module Factory.Test.QuickCheck.Statistics( +-- * Functions + quickChecks +) where + +import qualified Data.List +import qualified Data.Numbers.Primes +import qualified Factory.Math.Implementations.Factorial as Math.Implementations.Factorial +import qualified Factory.Math.Statistics as Math.Statistics +import Factory.Test.QuickCheck.Factorial() +import qualified Test.QuickCheck +import Test.QuickCheck((==>)) + +-- | Defines invariant properties. +quickChecks :: IO () +quickChecks = + Test.QuickCheck.quickCheck `mapM_` [prop_nC0, prop_nC1, prop_sum] + >> Test.QuickCheck.quickCheck `mapM_` [prop_symmetry, prop_prime] + >> Test.QuickCheck.quickCheck `mapM_` [prop_nP0, prop_nP1] + >> Test.QuickCheck.quickCheck prop_balance where + prop_nC0, prop_nC1, prop_sum :: Math.Implementations.Factorial.Algorithm -> Integer -> Test.QuickCheck.Property + prop_nC0 algorithm n = Test.QuickCheck.label "prop_nC0" $ Math.Statistics.nCr algorithm (abs n) 0 == 1 + + prop_nC1 algorithm i = Test.QuickCheck.label "prop_nC1" $ Math.Statistics.nCr algorithm n 1 == n where + n = 1 + abs i + + prop_sum algorithm i = Test.QuickCheck.label "prop_sum" $ sum (Math.Statistics.nCr algorithm n `map` [0 .. n]) == 2 ^ n where + n = 1 + abs i + + prop_symmetry, prop_prime :: Math.Implementations.Factorial.Algorithm -> (Integer, Integer) -> Test.QuickCheck.Property + prop_symmetry algorithm (i, j) = Test.QuickCheck.label "prop_symmetry" $ Math.Statistics.nCr algorithm n r == Math.Statistics.nCr algorithm n (n - r) where + [r, n] = Data.List.sort $ map abs [i, j] + + prop_prime algorithm (i, j) = r `notElem` [0, n] ==> Test.QuickCheck.label "prop_prime" $ (Math.Statistics.nCr algorithm n r `mod` n) == 0 where + n = Data.Numbers.Primes.primes !! fromIntegral (i `mod` 500000) + r = j `mod` n --Ensure r is smaller than n. + + prop_nP0, prop_nP1 :: Integer -> Test.QuickCheck.Property + prop_nP0 n = Test.QuickCheck.label "prop_nP0" $ Math.Statistics.nPr (abs n) 0 == 1 + + prop_nP1 i = Test.QuickCheck.label "prop_nP1" $ Math.Statistics.nPr n 1 == n where + n = 1 + abs i + + prop_balance :: [Integer] -> Test.QuickCheck.Property + prop_balance l = not (null l) ==> Test.QuickCheck.label "prop_balance" . (< 1e-11 {-rounding errors-}) . abs . sum $ map (\i -> fromIntegral i - (Math.Statistics.mean l :: Double)) l + diff --git a/src/Factory/Test/QuickCheck/Summation.hs b/src/Factory/Test/QuickCheck/Summation.hs new file mode 100644 --- /dev/null +++ b/src/Factory/Test/QuickCheck/Summation.hs @@ -0,0 +1,43 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] Defines /QuickCheck/-properties for "Math.Summation". +-} + +module Factory.Test.QuickCheck.Summation( +-- * Functions + quickChecks +) where + +import qualified Data.Ratio +import qualified Factory.Math.Summation as Math.Summation +import qualified Test.QuickCheck +import Test.QuickCheck((==>)) + +-- | Defines invariant properties. +quickChecks :: IO () +quickChecks = Test.QuickCheck.quickCheck `mapM_` [prop_sum, prop_sumR] where + prop_sum, prop_sumR :: Int -> [Data.Ratio.Rational] -> Test.QuickCheck.Property + prop_sum chunkSize l = not (null l) ==> Test.QuickCheck.label "prop_sum" $ Math.Summation.sum' chunkSize' l == sum l where + chunkSize' = 2 + (chunkSize `mod` length l) + + prop_sumR chunkSize l = not (null l) ==> Test.QuickCheck.label "prop_sumR" $ Math.Summation.sumR chunkSize' l == sum l where + chunkSize' = 2 + (chunkSize `mod` length l) + + diff --git a/src/Main.hs b/src/Main.hs new file mode 100644 --- /dev/null +++ b/src/Main.hs @@ -0,0 +1,165 @@ +{- + Copyright (C) 2011 Dr. Alistair Ward + + 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 . +-} +{- | + [@AUTHOR@] Dr. Alistair Ward + + [@DESCRIPTION@] + + * Contains the entry-point to the program. + + * Facilitates testing. +-} + +module Main( +-- * Type-classes +-- CommandLineAction, +-- * Functions + main +) where + +import qualified Data.List +import qualified Data.Ratio +import qualified Distribution.Package +import qualified Distribution.Text +import qualified Distribution.Version +import qualified Factory.Math.Implementations.Factorial as Math.Implementations.Factorial +import qualified Factory.Math.Implementations.Primality as Math.Implementations.Primality +import qualified Factory.Math.Implementations.PrimeFactorisation as Math.Implementations.PrimeFactorisation +import qualified Factory.Math.Implementations.SquareRoot as Math.Implementations.SquareRoot +import qualified Factory.Test.CommandOptions as Test.CommandOptions +import qualified Factory.Test.Performance.Factorial as Test.Performance.Factorial +import qualified Factory.Test.Performance.Pi as Test.Performance.Pi +import qualified Factory.Test.Performance.Primality as Test.Performance.Primality +import qualified Factory.Test.Performance.PrimeFactorisation as Test.Performance.PrimeFactorisation +import qualified Factory.Test.Performance.SquareRoot as Test.Performance.SquareRoot +import qualified Factory.Test.Performance.Statistics as Test.Performance.Statistics +import qualified Factory.Test.QuickCheck.QuickChecks as Test.QuickCheck.QuickChecks +import qualified System +import qualified System.Console.GetOpt as G +import qualified System.IO +import qualified System.IO.Error +import qualified ToolShed.Defaultable as Defaultable + +-- Local convenience definitions. +type PrimalityAlgorithm = Math.Implementations.Primality.Algorithm Math.Implementations.PrimeFactorisation.Algorithm +type PiCategory = Test.Performance.Pi.Category Math.Implementations.SquareRoot.Algorithm Math.Implementations.Factorial.Algorithm + +-- | Used to thread user-defined command-line options, though the list of functions which implement them. +type CommandLineAction = Test.CommandOptions.CommandOptions -> IO Test.CommandOptions.CommandOptions --Supplied as the type-argument to 'G.OptDescr'. + +-- | Parses the command-line arguments, to determine 'Test.CommandOptions.CommandOptions'. +main :: IO () +main = do + progName <- System.getProgName + args <- System.getArgs + + let + usage :: String + usage = "Usage:\t" ++ G.usageInfo progName optDescrList + +--Define the command-line options, and the 'CommandLineAction's used to handle them. + optDescrList :: [G.OptDescr CommandLineAction] + optDescrList = [ +-- String [String] (G.ArgDescr CommandLineAction) String + G.Option "" ["carmichaelNumbersPerformance"] (carmichaelNumbersPerformance `G.ReqArg` "(Math.Implementations.Primality.Algorithm, Int)") "Test the performance of 'Math.Primality.carmichaelNumbers'.", + G.Option "" ["factorialPerformance"] (factorialPerformance `G.ReqArg` "(Math.Implementations.Factorial.Algorithm, Integer)") "Test the performance of 'Math.Factorial.factorial'.", + G.Option "" ["factorialPerformanceGraph"] (factorialPerformanceGraph `G.ReqArg` "Math.Implementations.Factorial.Algorithm") "Test the performance of 'Math.Factorial.factorial', with an exponentially increasing operand.", + G.Option "" ["factorialPerformanceGraphControl"] (G.NoArg factorialPerformanceGraphControl) "Test the performance of a naive factorial-implementation, with an exponentially increasing operand.", + G.Option "" ["isPrimePerformance"] (isPrimePerformance `G.ReqArg` "(Math.Implementations.Primality.Algorithm, Integer)") "Test the performance of 'Math.Primality.isPrime'.", + G.Option "" ["isPrimePerformanceGraph"] (isPrimePerformanceGraph `G.ReqArg` "Math.Implementations.Primality.Algorithm") "Test the performance of 'Math.Primality.isPrime', against the prime-indexed Fibonacci-numbers.", + G.Option "" ["nCrPerformance"] (nCrPerformance `G.ReqArg` "(Math.Implementations.Factorial.Algorithm, Integer, Integer)") "Test the performance of 'Math.Factorial.factorial'.", + G.Option "" ["piPerformance"] (piPerformance `G.ReqArg` "(Math.Pi.Category, Math.Precision.DecimalDigits)") "Test the performance of 'Math.Pi.openI'.", + G.Option "" ["piPerformanceGraph"] (piPerformanceGraph `G.ReqArg` "(Math.Pi.Category, Double, Math.Precision.DecimalDigits)") "Test the performance of 'Math.Pi.openI', with an exponential precision-requirement (of the specified exponent), up to the specified limit.", + G.Option "" ["primeFactorsPerformance"] (primeFactorsPerformance `G.ReqArg` "(Math.Implementations.PrimeFactorisation.Algorithm, Integer)") "Test the performance of 'Math.PrimeFactorisation.primeFactors'.", + G.Option "" ["primeFactorsPerformanceGraph"] (primeFactorsPerformanceGraph `G.ReqArg` "(Math.Implementations.PrimeFactorisation.Algorithm, Int)") "Test the performance of 'Math.PrimeFactorisation.primeFactors', on the odd integers from the Fibonacci sequence.", + G.Option "" ["squareRootPerformance"] (squareRootPerformance `G.ReqArg` "(Math.Implementations.SquareRoot.Algorithm, Data.Ratio.Rational, DecimalDigits)") "Test 'Math.SquareRoot.squareRoot'.", + G.Option "" ["squareRootPerformanceGraph"] (squareRootPerformanceGraph `G.ReqArg` "(Math.Implementations.SquareRoot.Algorithm, Data.Ratio.Rational)") "Test the performance of 'Math.SquareRoot.squareRoot', with an exponentially increasing precision-requirement.", + G.Option "" ["verbose"] (G.NoArg $ return {-to IO-monad-} . Test.CommandOptions.setVerbose) ("Provide additional information where available; default '" ++ show (Test.CommandOptions.verbose Defaultable.defaultValue) ++ "'."), + G.Option "" ["version"] (G.NoArg $ const printVersion) "Print version-information & then exit.", + G.Option "q" ["runQuickChecks"] (G.NoArg $ const runQuickChecks) "Run Quick-checks using arbitrary data & then exit.", + G.Option "?" ["help"] (G.NoArg $ const printUsage) "Display this help-text & then exit." + ] where + printVersion, printUsage, runQuickChecks :: IO Test.CommandOptions.CommandOptions + printVersion = System.IO.hPutStrLn System.IO.stderr (Distribution.Text.display packageIdentifier ++ "\n\nCopyright (C) 2011 Dr. Alistair Ward.\nThis program comes with ABSOLUTELY NO WARRANTY.\nThis is free software, and you are welcome to redistribute it under certain conditions.\n\nWritten by Dr. Alistair Ward.") >> System.exitWith System.ExitSuccess where + packageIdentifier :: Distribution.Package.PackageIdentifier + packageIdentifier = Distribution.Package.PackageIdentifier { + Distribution.Package.pkgName = Distribution.Package.PackageName "factory", + Distribution.Package.pkgVersion = Distribution.Version.Version [0, 0, 0, 2] [] + } + + printUsage = System.IO.hPutStrLn System.IO.stderr usage >> System.exitWith System.ExitSuccess + runQuickChecks = Test.QuickCheck.QuickChecks.run >> System.exitWith System.ExitSuccess + + factorialPerformanceGraphControl :: Test.CommandOptions.CommandOptions -> IO Test.CommandOptions.CommandOptions + factorialPerformanceGraphControl commandOptions = Test.Performance.Factorial.factorialPerformanceGraphControl (Test.CommandOptions.verbose commandOptions) >> System.exitWith (System.ExitFailure 1) + + carmichaelNumbersPerformance, factorialPerformance, factorialPerformanceGraph, isPrimePerformance, isPrimePerformanceGraph, piPerformance, piPerformanceGraph, primeFactorsPerformance, primeFactorsPerformanceGraph, squareRootPerformance, squareRootPerformanceGraph :: String -> CommandLineAction + + carmichaelNumbersPerformance arg _ = Test.Performance.Primality.carmichaelNumbersPerformance algorithm i >>= print >> System.exitWith System.ExitSuccess where + algorithm :: PrimalityAlgorithm + (algorithm, i) = read arg + + factorialPerformance arg _ = Test.Performance.Factorial.factorialPerformance algorithm i >>= print >> System.exitWith System.ExitSuccess where + algorithm :: Math.Implementations.Factorial.Algorithm + i :: Integer + (algorithm, i) = read arg + + factorialPerformanceGraph arg commandOptions = Test.Performance.Factorial.factorialPerformanceGraph (Test.CommandOptions.verbose commandOptions) (read arg :: Math.Implementations.Factorial.Algorithm) >> System.exitWith (System.ExitFailure 1) + + isPrimePerformance arg _ = Test.Performance.Primality.isPrimePerformance algorithm i >>= print >> System.exitWith System.ExitSuccess where + algorithm :: PrimalityAlgorithm + (algorithm, i) = read arg + + isPrimePerformanceGraph arg _ = Test.Performance.Primality.isPrimePerformanceGraph (read arg :: Math.Implementations.Primality.Algorithm Math.Implementations.PrimeFactorisation.Algorithm) >> System.exitWith (System.ExitFailure 1) + + nCrPerformance arg _ = Test.Performance.Statistics.nCrPerformance algorithm n r >>= print >> System.exitWith System.ExitSuccess where + algorithm :: Math.Implementations.Factorial.Algorithm + n, r :: Integer + (algorithm, n, r) = read arg + + piPerformance arg _ = Test.Performance.Pi.piPerformance category decimalDigits >>= print >> System.exitWith System.ExitSuccess where + category :: PiCategory + (category, decimalDigits) = read arg + + piPerformanceGraph arg commandOptions = Test.Performance.Pi.piPerformanceGraph category factor maxDecimalDigits (Test.CommandOptions.verbose commandOptions) >> System.exitWith (System.ExitFailure 1) where + category :: PiCategory + factor :: Double + (category, factor, maxDecimalDigits) = read arg + + primeFactorsPerformance arg _ = Test.Performance.PrimeFactorisation.primeFactorsPerformance algorithm i >>= print >> System.exitWith System.ExitSuccess where + algorithm :: Math.Implementations.PrimeFactorisation.Algorithm + (algorithm, i) = read arg + + primeFactorsPerformanceGraph arg _ = Test.Performance.PrimeFactorisation.primeFactorsPerformanceGraph algorithm i >> System.exitWith (System.ExitFailure 1) where + algorithm :: Math.Implementations.PrimeFactorisation.Algorithm + (algorithm, i) = read arg + + squareRootPerformance arg _ = Test.Performance.SquareRoot.squareRootPerformance algorithm operand decimalDigits >>= print >> System.exitWith System.ExitSuccess where + algorithm :: Math.Implementations.SquareRoot.Algorithm + operand :: Data.Ratio.Rational + (algorithm, operand, decimalDigits) = read arg + + squareRootPerformanceGraph arg _ = Test.Performance.SquareRoot.squareRootPerformanceGraph algorithm operand >> System.exitWith (System.ExitFailure 1) where + algorithm :: Math.Implementations.SquareRoot.Algorithm + operand :: Data.Ratio.Rational + (algorithm, operand) = read arg + +-- G.getOpt :: G.ArgOrder CommandLineAction -> [G.OptDescr Action] -> [String] -> ([Action], [String], [String]) + case G.getOpt G.RequireOrder optDescrList args of + (commandLineActions, _, []) -> Data.List.foldl' (>>=) (return {-to IO-monad-} Defaultable.defaultValue) commandLineActions >> System.exitWith System.ExitSuccess + (_, _, errors) -> System.IO.Error.ioError . System.IO.Error.userError $ concat errors ++ usage --Throw. +