packages feed

factory (empty) → 0.0.0.2

raw patch · 75 files changed

+6962/−0 lines, 75 filesdep +Cabaldep +QuickCheckdep +arraysetup-changed

Dependencies added: Cabal, QuickCheck, array, base, containers, deepseq, haskell98, parallel, primes, toolshed

Files

+ LICENSE view
@@ -0,0 +1,675 @@+                    GNU GENERAL PUBLIC LICENSE+                       Version 3, 29 June 2007++ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>+ Everyone is permitted to copy and distribute verbatim copies+ of this license document, but changing it is not allowed.++                            Preamble++  The GNU General Public License is a free, copyleft license for+software and other kinds of works.++  The licenses for most software and other practical works are designed+to take away your freedom to share and change the works.  By contrast,+the GNU General Public License is intended to guarantee your freedom to+share and change all versions of a program--to make sure it remains free+software for all its users.  We, the Free Software Foundation, use the+GNU General Public License for most of our software; it applies also to+any other work released this way by its authors.  You can apply it to+your programs, too.++  When we speak of free software, we are referring to freedom, not+price.  Our General Public Licenses are designed to make sure that you+have the freedom to distribute copies of free software (and charge for+them if you wish), that you receive source code or can get it if you+want it, that you can change the software or use pieces of it in new+free programs, and that you know you can do these things.++  To protect your rights, we need to prevent others from denying you+these rights or asking you to surrender the rights.  Therefore, you have+certain responsibilities if you distribute copies of the software, or if+you modify it: responsibilities to respect the freedom of others.++  For example, if you distribute copies of such a program, whether+gratis or for a fee, you must pass on to the recipients the same+freedoms that you received.  You must make sure that they, too, receive+or can get the source code.  And you must show them these terms so they+know their rights.++  Developers that use the GNU GPL protect your rights with two steps:+(1) assert copyright on the software, and (2) offer you this License+giving you legal permission to copy, distribute and/or modify it.++  For the developers' and authors' protection, the GPL clearly explains+that there is no warranty for this free software.  For both users' and+authors' sake, the GPL requires that modified versions be marked as+changed, so that their problems will not be attributed erroneously to+authors of previous versions.++  Some devices are designed to deny users access to install or run+modified versions of the software inside them, although the manufacturer+can do so.  This is fundamentally incompatible with the aim of+protecting users' freedom to change the software.  The systematic+pattern of such abuse occurs in the area of products for individuals to+use, which is precisely where it is most unacceptable.  Therefore, we+have designed this version of the GPL to prohibit the practice for those+products.  If such problems arise substantially in other domains, we+stand ready to extend this provision to those domains in future versions+of the GPL, as needed to protect the freedom of users.++  Finally, every program is threatened constantly by software patents.+States should not allow patents to restrict development and use of+software on general-purpose computers, but in those that do, we wish to+avoid the special danger that patents applied to a free program could+make it effectively proprietary.  To prevent this, the GPL assures that+patents cannot be used to render the program non-free.++  The precise terms and conditions for copying, distribution and+modification follow.++                       TERMS AND CONDITIONS++  0. Definitions.++  "This License" refers to version 3 of the GNU General Public License.++  "Copyright" also means copyright-like laws that apply to other kinds of+works, such as semiconductor masks.++  "The Program" refers to any copyrightable work licensed under this+License.  Each licensee is addressed as "you".  "Licensees" and+"recipients" may be individuals or organizations.++  To "modify" a work means to copy from or adapt all or part of the work+in a fashion requiring copyright permission, other than the making of an+exact copy.  The resulting work is called a "modified version" of the+earlier work or a work "based on" the earlier work.++  A "covered work" means either the unmodified Program or a work based+on the Program.++  To "propagate" a work means to do anything with it that, without+permission, would make you directly or secondarily liable for+infringement under applicable copyright law, except executing it on a+computer or modifying a private copy.  Propagation includes copying,+distribution (with or without modification), making available to the+public, and in some countries other activities as well.++  To "convey" a work means any kind of propagation that enables other+parties to make or receive copies.  Mere interaction with a user through+a computer network, with no transfer of a copy, is not conveying.++  An interactive user interface displays "Appropriate Legal Notices"+to the extent that it includes a convenient and prominently visible+feature that (1) displays an appropriate copyright notice, and (2)+tells the user that there is no warranty for the work (except to the+extent that warranties are provided), that licensees may convey the+work under this License, and how to view a copy of this License.  If+the interface presents a list of user commands or options, such as a+menu, a prominent item in the list meets this criterion.++  1. Source Code.++  The "source code" for a work means the preferred form of the work+for making modifications to it.  "Object code" means any non-source+form of a work.++  A "Standard Interface" means an interface that either is an official+standard defined by a recognized standards body, or, in the case of+interfaces specified for a particular programming language, one that+is widely used among developers working in that language.++  The "System Libraries" of an executable work include anything, other+than the work as a whole, that (a) is included in the normal form of+packaging a Major Component, but which is not part of that Major+Component, and (b) serves only to enable use of the work with that+Major Component, or to implement a Standard Interface for which an+implementation is available to the public in source code form.  A+"Major Component", in this context, means a major essential component+(kernel, window system, and so on) of the specific operating system+(if any) on which the executable work runs, or a compiler used to+produce the work, or an object code interpreter used to run it.++  The "Corresponding Source" for a work in object code form means all+the source code needed to generate, install, and (for an executable+work) run the object code and to modify the work, including scripts to+control those activities.  However, it does not include the work's+System Libraries, or general-purpose tools or generally available free+programs which are used unmodified in performing those activities but+which are not part of the work.  For example, Corresponding Source+includes interface definition files associated with source files for+the work, and the source code for shared libraries and dynamically+linked subprograms that the work is specifically designed to require,+such as by intimate data communication or control flow between those+subprograms and other parts of the work.++  The Corresponding Source need not include anything that users+can regenerate automatically from other parts of the Corresponding+Source.++  The Corresponding Source for a work in source code form is that+same work.++  2. Basic Permissions.++  All rights granted under this License are granted for the term of+copyright on the Program, and are irrevocable provided the stated+conditions are met.  This License explicitly affirms your unlimited+permission to run the unmodified Program.  The output from running a+covered work is covered by this License only if the output, given its+content, constitutes a covered work.  This License acknowledges your+rights of fair use or other equivalent, as provided by copyright law.++  You may make, run and propagate covered works that you do not+convey, without conditions so long as your license otherwise remains+in force.  You may convey covered works to others for the sole purpose+of having them make modifications exclusively for you, or provide you+with facilities for running those works, provided that you comply with+the terms of this License in conveying all material for which you do+not control copyright.  Those thus making or running the covered works+for you must do so exclusively on your behalf, under your direction+and control, on terms that prohibit them from making any copies of+your copyrighted material outside their relationship with you.++  Conveying under any other circumstances is permitted solely under+the conditions stated below.  Sublicensing is not allowed; section 10+makes it unnecessary.++  3. Protecting Users' Legal Rights From Anti-Circumvention Law.++  No covered work shall be deemed part of an effective technological+measure under any applicable law fulfilling obligations under article+11 of the WIPO copyright treaty adopted on 20 December 1996, or+similar laws prohibiting or restricting circumvention of such+measures.++  When you convey a covered work, you waive any legal power to forbid+circumvention of technological measures to the extent such circumvention+is effected by exercising rights under this License with respect to+the covered work, and you disclaim any intention to limit operation or+modification of the work as a means of enforcing, against the work's+users, your or third parties' legal rights to forbid circumvention of+technological measures.++  4. Conveying Verbatim Copies.++  You may convey verbatim copies of the Program's source code as you+receive it, in any medium, provided that you conspicuously and+appropriately publish on each copy an appropriate copyright notice;+keep intact all notices stating that this License and any+non-permissive terms added in accord with section 7 apply to the code;+keep intact all notices of the absence of any warranty; and give all+recipients a copy of this License along with the Program.++  You may charge any price or no price for each copy that you convey,+and you may offer support or warranty protection for a fee.++  5. Conveying Modified Source Versions.++  You may convey a work based on the Program, or the modifications to+produce it from the Program, in the form of source code under the+terms of section 4, provided that you also meet all of these conditions:++    a) The work must carry prominent notices stating that you modified+    it, and giving a relevant date.++    b) The work must carry prominent notices stating that it is+    released under this License and any conditions added under section+    7.  This requirement modifies the requirement in section 4 to+    "keep intact all notices".++    c) You must license the entire work, as a whole, under this+    License to anyone who comes into possession of a copy.  This+    License will therefore apply, along with any applicable section 7+    additional terms, to the whole of the work, and all its parts,+    regardless of how they are packaged.  This License gives no+    permission to license the work in any other way, but it does not+    invalidate such permission if you have separately received it.++    d) If the work has interactive user interfaces, each must display+    Appropriate Legal Notices; however, if the Program has interactive+    interfaces that do not display Appropriate Legal Notices, your+    work need not make them do so.++  A compilation of a covered work with other separate and independent+works, which are not by their nature extensions of the covered work,+and which are not combined with it such as to form a larger program,+in or on a volume of a storage or distribution medium, is called an+"aggregate" if the compilation and its resulting copyright are not+used to limit the access or legal rights of the compilation's users+beyond what the individual works permit.  Inclusion of a covered work+in an aggregate does not cause this License to apply to the other+parts of the aggregate.++  6. Conveying Non-Source Forms.++  You may convey a covered work in object code form under the terms+of sections 4 and 5, provided that you also convey the+machine-readable Corresponding Source under the terms of this License,+in one of these ways:++    a) Convey the object code in, or embodied in, a physical product+    (including a physical distribution medium), accompanied by the+    Corresponding Source fixed on a durable physical medium+    customarily used for software interchange.++    b) Convey the object code in, or embodied in, a physical product+    (including a physical distribution medium), accompanied by a+    written offer, valid for at least three years and valid for as+    long as you offer spare parts or customer support for that product+    model, to give anyone who possesses the object code either (1) a+    copy of the Corresponding Source for all the software in the+    product that is covered by this License, on a durable physical+    medium customarily used for software interchange, for a price no+    more than your reasonable cost of physically performing this+    conveying of source, or (2) access to copy the+    Corresponding Source from a network server at no charge.++    c) Convey individual copies of the object code with a copy of the+    written offer to provide the Corresponding Source.  This+    alternative is allowed only occasionally and noncommercially, and+    only if you received the object code with such an offer, in accord+    with subsection 6b.++    d) Convey the object code by offering access from a designated+    place (gratis or for a charge), and offer equivalent access to the+    Corresponding Source in the same way through the same place at no+    further charge.  You need not require recipients to copy the+    Corresponding Source along with the object code.  If the place to+    copy the object code is a network server, the Corresponding Source+    may be on a different server (operated by you or a third party)+    that supports equivalent copying facilities, provided you maintain+    clear directions next to the object code saying where to find the+    Corresponding Source.  Regardless of what server hosts the+    Corresponding Source, you remain obligated to ensure that it is+    available for as long as needed to satisfy these requirements.++    e) Convey the object code using peer-to-peer transmission, provided+    you inform other peers where the object code and Corresponding+    Source of the work are being offered to the general public at no+    charge under subsection 6d.++  A separable portion of the object code, whose source code is excluded+from the Corresponding Source as a System Library, need not be+included in conveying the object code work.++  A "User Product" is either (1) a "consumer product", which means any+tangible personal property which is normally used for personal, family,+or household purposes, or (2) anything designed or sold for incorporation+into a dwelling.  In determining whether a product is a consumer product,+doubtful cases shall be resolved in favor of coverage.  For a particular+product received by a particular user, "normally used" refers to a+typical or common use of that class of product, regardless of the status+of the particular user or of the way in which the particular user+actually uses, or expects or is expected to use, the product.  A product+is a consumer product regardless of whether the product has substantial+commercial, industrial or non-consumer uses, unless such uses represent+the only significant mode of use of the product.++  "Installation Information" for a User Product means any methods,+procedures, authorization keys, or other information required to install+and execute modified versions of a covered work in that User Product from+a modified version of its Corresponding Source.  The information must+suffice to ensure that the continued functioning of the modified object+code is in no case prevented or interfered with solely because+modification has been made.++  If you convey an object code work under this section in, or with, or+specifically for use in, a User Product, and the conveying occurs as+part of a transaction in which the right of possession and use of the+User Product is transferred to the recipient in perpetuity or for a+fixed term (regardless of how the transaction is characterized), the+Corresponding Source conveyed under this section must be accompanied+by the Installation Information.  But this requirement does not apply+if neither you nor any third party retains the ability to install+modified object code on the User Product (for example, the work has+been installed in ROM).++  The requirement to provide Installation Information does not include a+requirement to continue to provide support service, warranty, or updates+for a work that has been modified or installed by the recipient, or for+the User Product in which it has been modified or installed.  Access to a+network may be denied when the modification itself materially and+adversely affects the operation of the network or violates the rules and+protocols for communication across the network.++  Corresponding Source conveyed, and Installation Information provided,+in accord with this section must be in a format that is publicly+documented (and with an implementation available to the public in+source code form), and must require no special password or key for+unpacking, reading or copying.++  7. Additional Terms.++  "Additional permissions" are terms that supplement the terms of this+License by making exceptions from one or more of its conditions.+Additional permissions that are applicable to the entire Program shall+be treated as though they were included in this License, to the extent+that they are valid under applicable law.  If additional permissions+apply only to part of the Program, that part may be used separately+under those permissions, but the entire Program remains governed by+this License without regard to the additional permissions.++  When you convey a copy of a covered work, you may at your option+remove any additional permissions from that copy, or from any part of+it.  (Additional permissions may be written to require their own+removal in certain cases when you modify the work.)  You may place+additional permissions on material, added by you to a covered work,+for which you have or can give appropriate copyright permission.++  Notwithstanding any other provision of this License, for material you+add to a covered work, you may (if authorized by the copyright holders of+that material) supplement the terms of this License with terms:++    a) Disclaiming warranty or limiting liability differently from the+    terms of sections 15 and 16 of this License; or++    b) Requiring preservation of specified reasonable legal notices or+    author attributions in that material or in the Appropriate Legal+    Notices displayed by works containing it; or++    c) Prohibiting misrepresentation of the origin of that material, or+    requiring that modified versions of such material be marked in+    reasonable ways as different from the original version; or++    d) Limiting the use for publicity purposes of names of licensors or+    authors of the material; or++    e) Declining to grant rights under trademark law for use of some+    trade names, trademarks, or service marks; or++    f) Requiring indemnification of licensors and authors of that+    material by anyone who conveys the material (or modified versions of+    it) with contractual assumptions of liability to the recipient, for+    any liability that these contractual assumptions directly impose on+    those licensors and authors.++  All other non-permissive additional terms are considered "further+restrictions" within the meaning of section 10.  If the Program as you+received it, or any part of it, contains a notice stating that it is+governed by this License along with a term that is a further+restriction, you may remove that term.  If a license document contains+a further restriction but permits relicensing or conveying under this+License, you may add to a covered work material governed by the terms+of that license document, provided that the further restriction does+not survive such relicensing or conveying.++  If you add terms to a covered work in accord with this section, you+must place, in the relevant source files, a statement of the+additional terms that apply to those files, or a notice indicating+where to find the applicable terms.++  Additional terms, permissive or non-permissive, may be stated in the+form of a separately written license, or stated as exceptions;+the above requirements apply either way.++  8. Termination.++  You may not propagate or modify a covered work except as expressly+provided under this License.  Any attempt otherwise to propagate or+modify it is void, and will automatically terminate your rights under+this License (including any patent licenses granted under the third+paragraph of section 11).++  However, if you cease all violation of this License, then your+license from a particular copyright holder is reinstated (a)+provisionally, unless and until the copyright holder explicitly and+finally terminates your license, and (b) permanently, if the copyright+holder fails to notify you of the violation by some reasonable means+prior to 60 days after the cessation.++  Moreover, your license from a particular copyright holder is+reinstated permanently if the copyright holder notifies you of the+violation by some reasonable means, this is the first time you have+received notice of violation of this License (for any work) from that+copyright holder, and you cure the violation prior to 30 days after+your receipt of the notice.++  Termination of your rights under this section does not terminate the+licenses of parties who have received copies or rights from you under+this License.  If your rights have been terminated and not permanently+reinstated, you do not qualify to receive new licenses for the same+material under section 10.++  9. Acceptance Not Required for Having Copies.++  You are not required to accept this License in order to receive or+run a copy of the Program.  Ancillary propagation of a covered work+occurring solely as a consequence of using peer-to-peer transmission+to receive a copy likewise does not require acceptance.  However,+nothing other than this License grants you permission to propagate or+modify any covered work.  These actions infringe copyright if you do+not accept this License.  Therefore, by modifying or propagating a+covered work, you indicate your acceptance of this License to do so.++  10. Automatic Licensing of Downstream Recipients.++  Each time you convey a covered work, the recipient automatically+receives a license from the original licensors, to run, modify and+propagate that work, subject to this License.  You are not responsible+for enforcing compliance by third parties with this License.++  An "entity transaction" is a transaction transferring control of an+organization, or substantially all assets of one, or subdividing an+organization, or merging organizations.  If propagation of a covered+work results from an entity transaction, each party to that+transaction who receives a copy of the work also receives whatever+licenses to the work the party's predecessor in interest had or could+give under the previous paragraph, plus a right to possession of the+Corresponding Source of the work from the predecessor in interest, if+the predecessor has it or can get it with reasonable efforts.++  You may not impose any further restrictions on the exercise of the+rights granted or affirmed under this License.  For example, you may+not impose a license fee, royalty, or other charge for exercise of+rights granted under this License, and you may not initiate litigation+(including a cross-claim or counterclaim in a lawsuit) alleging that+any patent claim is infringed by making, using, selling, offering for+sale, or importing the Program or any portion of it.++  11. Patents.++  A "contributor" is a copyright holder who authorizes use under this+License of the Program or a work on which the Program is based.  The+work thus licensed is called the contributor's "contributor version".++  A contributor's "essential patent claims" are all patent claims+owned or controlled by the contributor, whether already acquired or+hereafter acquired, that would be infringed by some manner, permitted+by this License, of making, using, or selling its contributor version,+but do not include claims that would be infringed only as a+consequence of further modification of the contributor version.  For+purposes of this definition, "control" includes the right to grant+patent sublicenses in a manner consistent with the requirements of+this License.++  Each contributor grants you a non-exclusive, worldwide, royalty-free+patent license under the contributor's essential patent claims, to+make, use, sell, offer for sale, import and otherwise run, modify and+propagate the contents of its contributor version.++  In the following three paragraphs, a "patent license" is any express+agreement or commitment, however denominated, not to enforce a patent+(such as an express permission to practice a patent or covenant not to+sue for patent infringement).  To "grant" such a patent license to a+party means to make such an agreement or commitment not to enforce a+patent against the party.++  If you convey a covered work, knowingly relying on a patent license,+and the Corresponding Source of the work is not available for anyone+to copy, free of charge and under the terms of this License, through a+publicly available network server or other readily accessible means,+then you must either (1) cause the Corresponding Source to be so+available, or (2) arrange to deprive yourself of the benefit of the+patent license for this particular work, or (3) arrange, in a manner+consistent with the requirements of this License, to extend the patent+license to downstream recipients.  "Knowingly relying" means you have+actual knowledge that, but for the patent license, your conveying the+covered work in a country, or your recipient's use of the covered work+in a country, would infringe one or more identifiable patents in that+country that you have reason to believe are valid.++  If, pursuant to or in connection with a single transaction or+arrangement, you convey, or propagate by procuring conveyance of, a+covered work, and grant a patent license to some of the parties+receiving the covered work authorizing them to use, propagate, modify+or convey a specific copy of the covered work, then the patent license+you grant is automatically extended to all recipients of the covered+work and works based on it.++  A patent license is "discriminatory" if it does not include within+the scope of its coverage, prohibits the exercise of, or is+conditioned on the non-exercise of one or more of the rights that are+specifically granted under this License.  You may not convey a covered+work if you are a party to an arrangement with a third party that is+in the business of distributing software, under which you make payment+to the third party based on the extent of your activity of conveying+the work, and under which the third party grants, to any of the+parties who would receive the covered work from you, a discriminatory+patent license (a) in connection with copies of the covered work+conveyed by you (or copies made from those copies), or (b) primarily+for and in connection with specific products or compilations that+contain the covered work, unless you entered into that arrangement,+or that patent license was granted, prior to 28 March 2007.++  Nothing in this License shall be construed as excluding or limiting+any implied license or other defenses to infringement that may+otherwise be available to you under applicable patent law.++  12. No Surrender of Others' Freedom.++  If conditions are imposed on you (whether by court order, agreement or+otherwise) that contradict the conditions of this License, they do not+excuse you from the conditions of this License.  If you cannot convey a+covered work so as to satisfy simultaneously your obligations under this+License and any other pertinent obligations, then as a consequence you may+not convey it at all.  For example, if you agree to terms that obligate you+to collect a royalty for further conveying from those to whom you convey+the Program, the only way you could satisfy both those terms and this+License would be to refrain entirely from conveying the Program.++  13. Use with the GNU Affero General Public License.++  Notwithstanding any other provision of this License, you have+permission to link or combine any covered work with a work licensed+under version 3 of the GNU Affero General Public License into a single+combined work, and to convey the resulting work.  The terms of this+License will continue to apply to the part which is the covered work,+but the special requirements of the GNU Affero General Public License,+section 13, concerning interaction through a network will apply to the+combination as such.++  14. Revised Versions of this License.++  The Free Software Foundation may publish revised and/or new versions of+the GNU General Public License from time to time.  Such new versions will+be similar in spirit to the present version, but may differ in detail to+address new problems or concerns.++  Each version is given a distinguishing version number.  If the+Program specifies that a certain numbered version of the GNU General+Public License "or any later version" applies to it, you have the+option of following the terms and conditions either of that numbered+version or of any later version published by the Free Software+Foundation.  If the Program does not specify a version number of the+GNU General Public License, you may choose any version ever published+by the Free Software Foundation.++  If the Program specifies that a proxy can decide which future+versions of the GNU General Public License can be used, that proxy's+public statement of acceptance of a version permanently authorizes you+to choose that version for the Program.++  Later license versions may give you additional or different+permissions.  However, no additional obligations are imposed on any+author or copyright holder as a result of your choosing to follow a+later version.++  15. Disclaimer of Warranty.++  THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY+APPLICABLE LAW.  EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR+PURPOSE.  THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM+IS WITH YOU.  SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.++  16. Limitation of Liability.++  IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF+SUCH DAMAGES.++  17. Interpretation of Sections 15 and 16.++  If the disclaimer of warranty and limitation of liability provided+above cannot be given local legal effect according to their terms,+reviewing courts shall apply local law that most closely approximates+an absolute waiver of all civil liability in connection with the+Program, unless a warranty or assumption of liability accompanies a+copy of the Program in return for a fee.++                     END OF TERMS AND CONDITIONS++            How to Apply These Terms to Your New Programs++  If you develop a new program, and you want it to be of the greatest+possible use to the public, the best way to achieve this is to make it+free software which everyone can redistribute and change under these terms.++  To do so, attach the following notices to the program.  It is safest+to attach them to the start of each source file to most effectively+state the exclusion of warranty; and each file should have at least+the "copyright" line and a pointer to where the full notice is found.++    <one line to give the program's name and a brief idea of what it does.>+    Copyright (C) <year>  <name of author>++    This program is free software: you can redistribute it and/or modify+    it under the terms of the GNU General Public License as published by+    the Free Software Foundation, either version 3 of the License, or+    (at your option) any later version.++    This program is distributed in the hope that it will be useful,+    but WITHOUT ANY WARRANTY; without even the implied warranty of+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the+    GNU General Public License for more details.++    You should have received a copy of the GNU General Public License+    along with this program.  If not, see <http://www.gnu.org/licenses/>.++Also add information on how to contact you by electronic and paper mail.++  If the program does terminal interaction, make it output a short+notice like this when it starts in an interactive mode:++    <program>  Copyright (C) <year>  <name of author>+    This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.+    This is free software, and you are welcome to redistribute it+    under certain conditions; type `show c' for details.++The hypothetical commands `show w' and `show c' should show the appropriate+parts of the General Public License.  Of course, your program's commands+might be different; for a GUI interface, you would use an "about box".++  You should also get your employer (if you work as a programmer) or school,+if any, to sign a "copyright disclaimer" for the program, if necessary.+For more information on this, and how to apply and follow the GNU GPL, see+<http://www.gnu.org/licenses/>.++  The GNU General Public License does not permit incorporating your program+into proprietary programs.  If your program is a subroutine library, you+may consider it more useful to permit linking proprietary applications with+the library.  If this is what you want to do, use the GNU Lesser General+Public License instead of this License.  But first, please read+<http://www.gnu.org/philosophy/why-not-lgpl.html>.+
+ Setup.hs view
@@ -0,0 +1,5 @@+#!/usr/bin/env runhaskell++import qualified	Distribution.Simple++main	= Distribution.Simple.defaultMain
+ changelog view
@@ -0,0 +1,16 @@+2011-03-01 Dr. Alistair Ward <factory at functionalley dot eu>++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
@@ -0,0 +1,11 @@+Author:+	Dr. Alistair Ward <factory at functionalley dot eu>.++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'.
+ factory.cabal view
@@ -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 <at> functionalley <dot> eu+Bug-reports:		factory <at> functionalley <dot> 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+
+ makefile view
@@ -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 <http://www.gnu.org/licenses/>.+ +.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/:.*//'+
+ src/Factory/Data/Bounds.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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; <http://en.wikipedia.org/wiki/Divide_and_conquer_algorithm>.++	* 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 <http://en.wikipedia.org/wiki/Hylomorphism_%28computer_science%29>,+	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+
+ src/Factory/Data/Exponential.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]++	* Describes a simple numeric type, designed to contain an /exponential/ number.++	* <http://en.wikipedia.org/wiki/Exponentiation>.+-}++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'; <http://en.wikipedia.org/wiki/Identity_element>.+-}+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+
+ src/Factory/Data/MonicPolynomial.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]++	* Describes a /monic polynomial; <http://en.wikipedia.org/wiki/Monic_polynomial#Classifications>;+	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)+
+ src/Factory/Data/Monomial.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]++	* Describes a <http://en.wikipedia.org/wiki/Monomial> 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+
+ src/Factory/Data/Polynomial.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]++	* Describes a <http://en.wikipedia.org/wiki/Univariate> polynomial and operations on it.++	* <http://en.wikipedia.org/wiki/Polynomial>.++	* <http://mathworld.wolfram.com/Polynomial.html>.+-}++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 (<http://en.wikipedia.org/wiki/Laurent_polynomial>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 <http://en.wikipedia.org/wiki/Monomial_order>);+	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/; <http://en.wikipedia.org/wiki/Polynomial_ring>.+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/.+	<http://en.wikipedia.org/wiki/Polynomial_long_division>.+	<http://demonstrations.wolfram.com/PolynomialLongDivision/>.+-}+	_ `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.++	* <http://en.wikipedia.org/wiki/Monic_polynomial#Classifications>.+-}+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.++	* <http://planetmath.org/encyclopedia/PolynomialCongruence.html>.+-}+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/.++	* <http://en.wikipedia.org/wiki/Degree_of_a_polynomial>.++	* <http://mathworld.wolfram.com/PolynomialDegree.html>.+-}+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.++	* <http://en.wikipedia.org/wiki/Scalar_multiplication>.+-}+(*=) :: (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`)+
+ src/Factory/Data/PrimeFactors.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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+
+ src/Factory/Data/QuotientRing.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]++	* Describes a /Quotient Ring/; <http://en.wikipedia.org/wiki/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'.++	* <http://www.usna.edu/Users/math/wdj/book/node74.html>.+-}+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+
+ src/Factory/Data/Ring.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]++	* Describes a /ring/ and operations on its members.++	* <http://en.wikipedia.org/wiki/Ring_%28mathematics%29>.++	* <http://www.numericana.com/answer/rings.htm>.+-}++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/; <http://en.wikipedia.org/wiki/Commutativity>.+	(=*=)			:: r -> r -> r	-- ^ Multiplication of two members.+	additiveInverse		:: r -> r	-- ^ The operand required to yield /zero/ under addition; <http://en.wikipedia.org/wiki/Additive_inverse>.+	multiplicativeIdentity	:: r		-- ^ The /identity/-member under multiplication; <http://mathworld.wolfram.com/MultiplicativeIdentity.html>.+	additiveIdentity	:: r		-- ^ The /identity/-member under addition (AKA /zero/); <http://en.wikipedia.org/wiki/Additive_identity>.++	(=-=) :: 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;+	<http://en.wikipedia.org/wiki/Exponentiation_by_squaring>.+-}+(=^) :: (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+
+ src/Factory/Math/ArithmeticGeometricMean.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]	Determines the /Arithmetic-geometric mean/; <http://en.wikipedia.org/wiki/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/; <http://en.wikipedia.org/wiki/Arithmetic_mean>.+type ArithmeticMean	= Data.Ratio.Rational++-- | The type of the /geometric mean/; <http://en.wikipedia.org/wiki/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]+
+ src/Factory/Math/DivideAndConquer.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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; <http://en.wikipedia.org/wiki/Divide_and_conquer_algorithm>.++	* 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 <http://en.wikipedia.org/wiki/Hylomorphism_%28computer_science%29>,+	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+
+ src/Factory/Math/Factorial.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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.++	* <http://mathworld.wolfram.com/Factorial.html>.+-}++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+
+ src/Factory/Math/Fibonacci.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]	<http://en.wikipedia.org/wiki/Fibonacci_number>.+-}++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.++	* <http://primes.utm.edu/glossary/page.php?sort=FibonacciPrime>.+-}+primeIndexedFibonacci :: Integral i => [i]+primeIndexedFibonacci	= map (fibonacci !!) Data.Numbers.Primes.primes+
+ src/Factory/Math/Implementations/Factorial.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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.++	* <http://en.wikipedia.org/wiki/Factorial>.++	* <http://mathworld.wolfram.com/Factorial.html>.++	* <http://www.luschny.de/math/factorial/FastFactorialFunctions.htm>.+-}++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.++	* <http://en.wikipedia.org/wiki/Factorial#Number_theory>.++	* 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/.++	* <http://www.proofwiki.org/wiki/Multiplicity_of_Prime_Factor_in_Factorial>.+-}+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/; <http://mathworld.wolfram.com/RisingFactorial.html>+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/; <http://mathworld.wolfram.com/FallingFactorial.html>+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)+
+ src/Factory/Math/Implementations/Pi/AGM/Algorithm.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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+
+ src/Factory/Math/Implementations/Pi/AGM/BrentSalamin.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]++	* Implements the /Brent-Salamin/ (AKA /Gauss-Legendre/) algorithm;+		<http://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_algorithm>,+		<http://mathworld.wolfram.com/Brent-SalaminFormula.html>,+		<http://www.pi314.net/eng/salamin.php>.++	* 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))+
+ src/Factory/Math/Implementations/Pi/BBP/Algorithm.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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+
+ src/Factory/Math/Implementations/Pi/BBP/Base65536.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]	Defines a specific base-@2^16@ /BBP/-formula; <http://mathworld.wolfram.com/PiFormulas.html>++-}++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)+}
+ src/Factory/Math/Implementations/Pi/BBP/Bellard.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]	Defines /Bellard/'s nega-base-@2^10@ /BBP/-formula; <http://en.wikipedia.org/wiki/Bellard%27s_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)+}
+ src/Factory/Math/Implementations/Pi/BBP/Implementation.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]++	* Implements a /Bailey-Borwein-Plouffe/ formula; <http://mathworld.wolfram.com/PiFormulas.html>++	* 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 ..]+
+ src/Factory/Math/Implementations/Pi/BBP/Series.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]	Defines a /Bailey-Borwein-Plouffe/ formula; <http://mathworld.wolfram.com/PiFormulas.html>+-}++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.+}+
+ src/Factory/Math/Implementations/Pi/Borwein/Algorithm.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]	Defines the set of /Borwein/-type algorithms (currently only one) which have been implemented; <http://www.pi314.net/eng/borwein.php>.+-}++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	-- ^ <http://en.wikipedia.org/wiki/Borwein%27s_algorithm>.+	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+
+ src/Factory/Math/Implementations/Pi/Borwein/Borwein1993.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]	Defines the /Borwein/ series for /Pi/; <http://en.wikipedia.org/wiki/Borwein%27s_algorithm#Jonathan_Borwein_and_Peter_Borwein.27s_Version_.281993.29>+-}++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.+}
+ src/Factory/Math/Implementations/Pi/Borwein/Implementation.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]	Defines /Borwein/ series for /Pi/; <http://en.wikipedia.org/wiki/Borwein%27s_algorithm>+-}++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 +
+ src/Factory/Math/Implementations/Pi/Borwein/Series.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]	Defines a <http://en.wikipedia.org/wiki/Srinivasa_Borwein>-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.+}+
+ src/Factory/Math/Implementations/Pi/Ramanujan/Algorithm.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]	Defines the set of /Ramanujan/-type algorithms which have been implemented; <http://en.wikipedia.org/wiki/Pi>.+-}++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+
+ src/Factory/Math/Implementations/Pi/Ramanujan/Chudnovsky.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]	Defines the /Chudnovsky/ series for /Pi/; <http://en.wikipedia.org/wiki/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.+}+
+ src/Factory/Math/Implementations/Pi/Ramanujan/Classic.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]	Defines the /Ramanujan/ series for /Pi/; <http://en.wikipedia.org/wiki/Srinivasa_Ramanujan>.+-}++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.+}+
+ src/Factory/Math/Implementations/Pi/Ramanujan/Implementation.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]	Implements a /Ramanujan/-type series for /Pi/; <http://en.wikipedia.org/wiki/Srinivasa_Ramanujan>.+-}++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+	)+
+ src/Factory/Math/Implementations/Pi/Ramanujan/Series.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]	Defines a <http://en.wikipedia.org/wiki/Srinivasa_Ramanujan>-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.+}+
+ src/Factory/Math/Implementations/Pi/Spigot/Algorithm.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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))+
+ src/Factory/Math/Implementations/Pi/Spigot/Gosper.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]	Defines the /Gosper/ series; <http://www.pi314.net/eng/goutte.php>+-}++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-}+}+
+ src/Factory/Math/Implementations/Pi/Spigot/RabinowitzWagon.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]	Defines the /Rabinowitz-Wagon/ series;+	<http://web.comlab.ox.ac.uk/oucl/work/jeremy.gibbons/publications/spigot.pdf>+	<http://www.mathpropress.com/stan/bibliography/spigot.pdf>.+-}++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-}+}
+ src/Factory/Math/Implementations/Pi/Spigot/Series.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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+
+ src/Factory/Math/Implementations/Pi/Spigot/Spigot.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]++	* Implements a /spigot/-algorithm; <http://en.wikipedia.org/wiki/Spigot_algorithm>.++	* Uses the traditional algorithm, rather than the /unbounded/ algorithm described by <http://www.comlab.ox.ac.uk/jeremy.gibbons/publications/spigot.pdf>.+-}++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+
+ src/Factory/Math/Implementations/Primality.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]++	* Determines whether an integer is prime.++	* <http://en.wikipedia.org/wiki/Primality_test>.++	* <http://primes.utm.edu/index.html>++	* 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	-- ^ <http://en.wikipedia.org/wiki/AKS_primality_test>.+	| MillerRabin			-- ^ <http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test>.+	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; <http://en.wikipedia.org/wiki/AKS_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@]	<http://www.cse.iitk.ac.in/users/manindra/algebra/primality_v6.pdf>.++	[@H. W. Lenstra, Jr. and Carl Pomerance@]		<http://www.math.dartmouth.edu/~carlp/PDF/complexity12.pdf>.++	[@Salembier and Southerington@]				<http://ece.gmu.edu/courses/ECE746/project/F06_Project_resources/Salembier_Southerington_AKS.pdf>,++	[@R. Crandall and J. Papadopoulos@]			<http://images.apple.com/acg/pdf/aks3.pdf>,++	[@Andreas Klappenecker@]				<http://faculty.cs.tamu.edu/klappi/629/aks.ps>,++	[@Vibhor Bhatt and G. K. Patra@]			<http://www.cmmacs.ernet.in/cmmacs/Publications/resch_rep/rrcm0307.pdf>,+-}+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/; <http://mathworld.wolfram.com/MillersPrimalityTest.html>.++	* 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; <http://primes.utm.edu/prove/prove2_3.html>.++	* <http://en.wikipedia.org/wiki/Miller-Rabin_primality_test>.++	* <http://mathworld.wolfram.com/Rabin-MillerStrongPseudoprimeTest.html>++	* <http://mathworld.wolfram.com/StrongPseudoprime.html>.++	* <http://oeis.org/A014233>, <http://oeis.org/A006945>.+-}+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.+			 ]+
+ src/Factory/Math/Implementations/PrimeFactorisation.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]++	* Implements several different prime-factorisation algorithms.++	* <http://www.tug.org/texinfohtml/coreutils.html#factor-invocation>.+-}++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	-- ^ <http://en.wikipedia.org/wiki/Dixon%27s_factorization_method>.+	| FermatsMethod	-- ^ <http://en.wikipedia.org/wiki/Fermat%27s_factorization_method>.+	| TrialDivision	-- ^ <http://en.wikipedia.org/wiki/Trial_division>.+	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++{- |+	* <http://en.wikipedia.org/wiki/Dixon%27s_factorization_method>.+-}+factoriseByDixonsMethod :: Integral base => base -> Data.PrimeFactors.Factors base exponent+factoriseByDixonsMethod	= undefined++{- |+	* <http://en.wikipedia.org/wiki/Fermat%27s_factorization_method>.++	* <http://mathworld.wolfram.com/FermatsFactorizationMethod.html>.++	* <http://en.wikipedia.org/wiki/Congruence_of_squares>.++	*	@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 <http://mathworld.wolfram.com/DirectSearchFactorization.html>, AKA <http://en.wikipedia.org/wiki/Trial_division>.++	* 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+
+ src/Factory/Math/Implementations/SquareRoot.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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	-- ^ <http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Bakhshali_approximation>+	| ContinuedFraction		-- ^ <http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Continued_fraction_expansion>.+	| HalleysMethod			-- ^ <http://en.wikipedia.org/wiki/Halley%27s_method>.+	| NewtonRaphsonIteration	-- ^ <http://en.wikipedia.org/wiki/Newton%27s_method>.+	| TaylorSeries Terms		-- ^ <http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Taylor_series>.+	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; <http://mathworld.wolfram.com/HalleysMethod.html>++>		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;+	<http://en.wikipedia.org/wiki/Solving_quadratic_equations_with_continued_fractions>,+	<http://en.wikipedia.org/wiki/Generalized_continued_fraction#Roots_of_positive_numbers>,+	<http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Continued_fraction_expansion>.+	<http://www.myreckonings.com/Dead_Reckoning/Online/Materials/General%20Method%20for%20Extracting%20Roots.pdf>++	* The convergence <http://en.wikipedia.org/wiki/Rate_of_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/; <http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Taylor_series>.++	* @ ((-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.++	* <http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Taylor_series>.++	* 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+
+ src/Factory/Math/MultiplicativeOrder.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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 <http://rosettacode.org/wiki/Multiplicative_order#Haskell>.++	* <http://en.wikipedia.org/wiki/Multiplicative_order>.++	* <http://mathworld.wolfram.com/MultiplicativeOrder.html>.+-}+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+
+ src/Factory/Math/Pi.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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		-- ^ <http://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula>.+	| Borwein borwein	-- ^ <http://en.wikipedia.org/wiki/Borwein%27s_algorithm>.+	| Ramanujan ramanujan	-- ^ <http://www.pi314.net/eng/ramanujan.php>.+	| 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+
+ src/Factory/Math/Power.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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;+	<http://en.wikipedia.org/wiki/Exponentiation_by_squaring>.++	* <http://en.wikipedia.org/wiki/Modular_exponentiation>.+-}+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/).++	* <http://en.wikipedia.org/wiki/Square_number>.++	* <http://mathworld.wolfram.com/SquareNumber.html>.++	* @(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.++	* <http://en.wikipedia.org/wiki/Perfect_power>.++	* <http://mathworld.wolfram.com/PerfectPower.html>.+-}+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)]+
+ src/Factory/Math/Precision.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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/; <http://en.wikipedia.org/wiki/Rate_of_convergence>.+type ConvergenceOrder	= Int++-- | The /rate of convergence/; <http://en.wikipedia.org/wiki/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@.++	* <http://en.wikipedia.org/wiki/Rate_of_convergence>.+-}+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.+
+ src/Factory/Math/Primality.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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.++	* <http://en.wikipedia.org/wiki/Coprime>.++	* <http://mathworld.wolfram.com/RelativelyPrime.html>.+-}+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.++	* <http://en.wikipedia.org/wiki/Fermat%27s_little_theorem>.++	* <http://en.wikipedia.org/wiki/Fermat_primality_test>.++	* <http://en.wikipedia.org/wiki/Fermat_pseudoprime>.++	* 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/.++	* <http://en.wikipedia.org/wiki/Carmichael_number>.++	* <http://mathworld.wolfram.com/CarmichaelNumber.html>.+-}+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; <http://en.wikipedia.org/wiki/Carmichael_number>.+carmichaelNumbers :: (Algorithm algorithm, Control.DeepSeq.NFData i, Integral i) => algorithm -> [i]+carmichaelNumbers algorithm	= isCarmichaelNumber algorithm `filter` [3, 5 ..]
+ src/Factory/Math/PrimeFactorisation.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@AUTHOR@]	Dr. Alistair Ward++ [@DESCRIPTION@]++	* <http://en.wikipedia.org/wiki/Integer_factorization>.++	* 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/; <http://en.wikipedia.org/wiki/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.++	* <http://en.wikipedia.org/wiki/Smooth_number>.++	* <http://mathworld.wolfram.com/SmoothNumber.html>.+-}+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.++	* <http://en.wikipedia.org/wiki/Smooth_number#Powersmooth_numbers>.+-}+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/.++	* <http://en.wikipedia.org/wiki/Regular_number>.+-}+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.++	* <http://en.wikipedia.org/wiki/Euler%27s_totient_function>.++	* <http://mathworld.wolfram.com/TotientFunction.html>.++	* 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/.++	* <http://oeis.org/wiki/Omega%28n%29,_number_of_distinct_primes_dividing_n>.++	* <http://mathworld.wolfram.com/DistinctPrimeFactors.html>++	* <http://planetmath.org/encyclopedia/NumberOfDistinctPrimeFactorsFunction.html>.+-}+omega :: (Algorithm algorithm, Control.DeepSeq.NFData i, Integral i) => algorithm -> [i]+omega algorithm	= map (Data.List.genericLength . primeFactors algorithm) [0 :: Integer ..]+
+ src/Factory/Math/Radix.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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++{- |+	* <http://mathworld.wolfram.com/DigitSum.html>.++	* <http://en.wikipedia.org/wiki/Digit_sum>.+-}+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++-- | <http://en.wikipedia.org/wiki/Digital_root>.+digitalRoot :: Integral decimal => decimal -> decimal+digitalRoot	= head . dropWhile (> 9) . iterate (digitSum (10 :: Int))+
+ src/Factory/Math/SquareRoot.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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	-- <http://en.wikipedia.org/wiki/IEEE_floating_point>.++{- |+	* 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'.+
+ src/Factory/Math/Statistics.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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 <http://en.wikipedia.org/wiki/Mean> 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/; <http://en.wikipedia.org/wiki/Combination>.+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/; <http://en.wikipedia.org/wiki/Permutations>.+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+
+ src/Factory/Math/Summation.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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
+ src/Factory/Test/CommandOptions.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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+}++
+ src/Factory/Test/Performance/Factorial.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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)+
+ src/Factory/Test/Performance/Pi.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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
+ src/Factory/Test/Performance/Primality.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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+
+ src/Factory/Test/Performance/PrimeFactorisation.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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+
+ src/Factory/Test/Performance/SquareRoot.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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
+ src/Factory/Test/Performance/Statistics.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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+
+ src/Factory/Test/QuickCheck/ArithmeticGeometricMean.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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)+
+ src/Factory/Test/QuickCheck/Bounds.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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
+ src/Factory/Test/QuickCheck/Factorial.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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+
+ src/Factory/Test/QuickCheck/MonicPolynomial.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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++
+ src/Factory/Test/QuickCheck/Pi.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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)+
+ src/Factory/Test/QuickCheck/Polynomial.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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+
+ src/Factory/Test/QuickCheck/Power.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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+
+ src/Factory/Test/QuickCheck/Primality.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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+
+ src/Factory/Test/QuickCheck/PrimeFactorisation.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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
+ src/Factory/Test/QuickCheck/QuickChecks.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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+
+ src/Factory/Test/QuickCheck/Radix.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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)+
+ src/Factory/Test/QuickCheck/SquareRoot.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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'+
+ src/Factory/Test/QuickCheck/Statistics.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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+
+ src/Factory/Test/QuickCheck/Summation.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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)++
+ src/Main.hs view
@@ -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 <http://www.gnu.org/licenses/>.+-}+{- |+ [@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.+