diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,674 @@ + GNU GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The GNU General Public License is a free, copyleft license for +software and other kinds of works. + + The licenses for most software and other practical works are designed +to take away your freedom to share and change the works. By contrast, +the GNU General Public License is intended to guarantee your freedom to +share and change all versions of a program--to make sure it remains free +software for all its users. We, the Free Software Foundation, use the +GNU General Public License for most of our software; it applies also to +any other work released this way by its authors. You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +them if you wish), that you receive source code or can get it if you +want it, that you can change the software or use pieces of it in new +free programs, and that you know you can do these things. + + To protect your rights, we need to prevent others from denying you +these rights or asking you to surrender the rights. Therefore, you have +certain responsibilities if you distribute copies of the software, or if +you modify it: responsibilities to respect the freedom of others. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must pass on to the recipients the same +freedoms that you received. You must make sure that they, too, receive +or can get the source code. And you must show them these terms so they +know their rights. + + Developers that use the GNU GPL protect your rights with two steps: +(1) assert copyright on the software, and (2) offer you this License +giving you legal permission to copy, distribute and/or modify it. + + For the developers' and authors' protection, the GPL clearly explains +that there is no warranty for this free software. For both users' and +authors' sake, the GPL requires that modified versions be marked as +changed, so that their problems will not be attributed erroneously to +authors of previous versions. + + Some devices are designed to deny users access to install or run +modified versions of the software inside them, although the manufacturer +can do so. This is fundamentally incompatible with the aim of +protecting users' freedom to change the software. The systematic +pattern of such abuse occurs in the area of products for individuals to +use, which is precisely where it is most unacceptable. Therefore, we +have designed this version of the GPL to prohibit the practice for those +products. If such problems arise substantially in other domains, we +stand ready to extend this provision to those domains in future versions +of the GPL, as needed to protect the freedom of users. + + Finally, every program is threatened constantly by software patents. +States should not allow patents to restrict development and use of +software on general-purpose computers, but in those that do, we wish to +avoid the special danger that patents applied to a free program could +make it effectively proprietary. To prevent this, the GPL assures that +patents cannot be used to render the program non-free. + + The precise terms and conditions for copying, distribution and +modification follow. + + TERMS AND CONDITIONS + + 0. Definitions. + + "This License" refers to version 3 of the GNU General Public License. + + "Copyright" also means copyright-like laws that apply to other kinds of +works, such as semiconductor masks. + + "The Program" refers to any copyrightable work licensed under this +License. Each licensee is addressed as "you". "Licensees" and +"recipients" may be individuals or organizations. + + To "modify" a work means to copy from or adapt all or part of the work +in a fashion requiring copyright permission, other than the making of an +exact copy. The resulting work is called a "modified version" of the +earlier work or a work "based on" the earlier work. + + A "covered work" means either the unmodified Program or a work based +on the Program. + + To "propagate" a work means to do anything with it that, without +permission, would make you directly or secondarily liable for +infringement under applicable copyright law, except executing it on a +computer or modifying a private copy. Propagation includes copying, +distribution (with or without modification), making available to the +public, and in some countries other activities as well. + + To "convey" a work means any kind of propagation that enables other +parties to make or receive copies. Mere interaction with a user through +a computer network, with no transfer of a copy, is not conveying. + + An interactive user interface displays "Appropriate Legal Notices" +to the extent that it includes a convenient and prominently visible +feature that (1) displays an appropriate copyright notice, and (2) +tells the user that there is no warranty for the work (except to the +extent that warranties are provided), that licensees may convey the +work under this License, and how to view a copy of this License. If +the interface presents a list of user commands or options, such as a +menu, a prominent item in the list meets this criterion. + + 1. Source Code. + + The "source code" for a work means the preferred form of the work +for making modifications to it. "Object code" means any non-source +form of a work. + + A "Standard Interface" means an interface that either is an official +standard defined by a recognized standards body, or, in the case of +interfaces specified for a particular programming language, one that +is widely used among developers working in that language. + + The "System Libraries" of an executable work include anything, other +than the work as a whole, that (a) is included in the normal form of +packaging a Major Component, but which is not part of that Major +Component, and (b) serves only to enable use of the work with that +Major Component, or to implement a Standard Interface for which an +implementation is available to the public in source code form. A +"Major Component", in this context, means a major essential component +(kernel, window system, and so on) of the specific operating system +(if any) on which the executable work runs, or a compiler used to +produce the work, or an object code interpreter used to run it. + + The "Corresponding Source" for a work in object code form means all +the source code needed to generate, install, and (for an executable +work) run the object code and to modify the work, including scripts to +control those activities. However, it does not include the work's +System Libraries, or general-purpose tools or generally available free +programs which are used unmodified in performing those activities but +which are not part of the work. For example, Corresponding Source +includes interface definition files associated with source files for +the work, and the source code for shared libraries and dynamically +linked subprograms that the work is specifically designed to require, +such as by intimate data communication or control flow between those +subprograms and other parts of the work. + + The Corresponding Source need not include anything that users +can regenerate automatically from other parts of the Corresponding +Source. + + The Corresponding Source for a work in source code form is that +same work. + + 2. Basic Permissions. + + All rights granted under this License are granted for the term of +copyright on the Program, and are irrevocable provided the stated +conditions are met. This License explicitly affirms your unlimited +permission to run the unmodified Program. The output from running a +covered work is covered by this License only if the output, given its +content, constitutes a covered work. This License acknowledges your +rights of fair use or other equivalent, as provided by copyright law. + + You may make, run and propagate covered works that you do not +convey, without conditions so long as your license otherwise remains +in force. You may convey covered works to others for the sole purpose +of having them make modifications exclusively for you, or provide you +with facilities for running those works, provided that you comply with +the terms of this License in conveying all material for which you do +not control copyright. Those thus making or running the covered works +for you must do so exclusively on your behalf, under your direction +and control, on terms that prohibit them from making any copies of +your copyrighted material outside their relationship with you. + + Conveying under any other circumstances is permitted solely under +the conditions stated below. Sublicensing is not allowed; section 10 +makes it unnecessary. + + 3. Protecting Users' Legal Rights From Anti-Circumvention Law. + + No covered work shall be deemed part of an effective technological +measure under any applicable law fulfilling obligations under article +11 of the WIPO copyright treaty adopted on 20 December 1996, or +similar laws prohibiting or restricting circumvention of such +measures. + + When you convey a covered work, you waive any legal power to forbid +circumvention of technological measures to the extent such circumvention +is effected by exercising rights under this License with respect to +the covered work, and you disclaim any intention to limit operation or +modification of the work as a means of enforcing, against the work's +users, your or third parties' legal rights to forbid circumvention of +technological measures. + + 4. Conveying Verbatim Copies. + + You may convey verbatim copies of the Program's source code as you +receive it, in any medium, provided that you conspicuously and +appropriately publish on each copy an appropriate copyright notice; +keep intact all notices stating that this License and any +non-permissive terms added in accord with section 7 apply to the code; +keep intact all notices of the absence of any warranty; and give all +recipients a copy of this License along with the Program. + + You may charge any price or no price for each copy that you convey, +and you may offer support or warranty protection for a fee. + + 5. Conveying Modified Source Versions. + + You may convey a work based on the Program, or the modifications to +produce it from the Program, in the form of source code under the +terms of section 4, provided that you also meet all of these conditions: + + a) The work must carry prominent notices stating that you modified + it, and giving a relevant date. + + b) The work must carry prominent notices stating that it is + released under this License and any conditions added under section + 7. This requirement modifies the requirement in section 4 to + "keep intact all notices". + + c) You must license the entire work, as a whole, under this + License to anyone who comes into possession of a copy. This + License will therefore apply, along with any applicable section 7 + additional terms, to the whole of the work, and all its parts, + regardless of how they are packaged. This License gives no + permission to license the work in any other way, but it does not + invalidate such permission if you have separately received it. + + d) If the work has interactive user interfaces, each must display + Appropriate Legal Notices; however, if the Program has interactive + interfaces that do not display Appropriate Legal Notices, your + work need not make them do so. + + A compilation of a covered work with other separate and independent +works, which are not by their nature extensions of the covered work, +and which are not combined with it such as to form a larger program, +in or on a volume of a storage or distribution medium, is called an +"aggregate" if the compilation and its resulting copyright are not +used to limit the access or legal rights of the compilation's users +beyond what the individual works permit. Inclusion of a covered work +in an aggregate does not cause this License to apply to the other +parts of the aggregate. + + 6. Conveying Non-Source Forms. + + You may convey a covered work in object code form under the terms +of sections 4 and 5, provided that you also convey the +machine-readable Corresponding Source under the terms of this License, +in one of these ways: + + a) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by the + Corresponding Source fixed on a durable physical medium + customarily used for software interchange. + + b) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by a + written offer, valid for at least three years and valid for as + long as you offer spare parts or customer support for that product + model, to give anyone who possesses the object code either (1) a + copy of the Corresponding Source for all the software in the + product that is covered by this License, on a durable physical + medium customarily used for software interchange, for a price no + more than your reasonable cost of physically performing this + conveying of source, or (2) access to copy the + Corresponding Source from a network server at no charge. + + c) Convey individual copies of the object code with a copy of the + written offer to provide the Corresponding Source. This + alternative is allowed only occasionally and noncommercially, and + only if you received the object code with such an offer, in accord + with subsection 6b. + + d) Convey the object code by offering access from a designated + place (gratis or for a charge), and offer equivalent access to the + Corresponding Source in the same way through the same place at no + further charge. You need not require recipients to copy the + Corresponding Source along with the object code. If the place to + copy the object code is a network server, the Corresponding Source + may be on a different server (operated by you or a third party) + that supports equivalent copying facilities, provided you maintain + clear directions next to the object code saying where to find the + Corresponding Source. Regardless of what server hosts the + Corresponding Source, you remain obligated to ensure that it is + available for as long as needed to satisfy these requirements. + + e) Convey the object code using peer-to-peer transmission, provided + you inform other peers where the object code and Corresponding + Source of the work are being offered to the general public at no + charge under subsection 6d. + + A separable portion of the object code, whose source code is excluded +from the Corresponding Source as a System Library, need not be +included in conveying the object code work. + + A "User Product" is either (1) a "consumer product", which means any +tangible personal property which is normally used for personal, family, +or household purposes, or (2) anything designed or sold for incorporation +into a dwelling. In determining whether a product is a consumer product, +doubtful cases shall be resolved in favor of coverage. For a particular +product received by a particular user, "normally used" refers to a +typical or common use of that class of product, regardless of the status +of the particular user or of the way in which the particular user +actually uses, or expects or is expected to use, the product. A product +is a consumer product regardless of whether the product has substantial +commercial, industrial or non-consumer uses, unless such uses represent +the only significant mode of use of the product. + + "Installation Information" for a User Product means any methods, +procedures, authorization keys, or other information required to install +and execute modified versions of a covered work in that User Product from +a modified version of its Corresponding Source. The information must +suffice to ensure that the continued functioning of the modified object +code is in no case prevented or interfered with solely because +modification has been made. + + If you convey an object code work under this section in, or with, or +specifically for use in, a User Product, and the conveying occurs as +part of a transaction in which the right of possession and use of the +User Product is transferred to the recipient in perpetuity or for a +fixed term (regardless of how the transaction is characterized), the +Corresponding Source conveyed under this section must be accompanied +by the Installation Information. But this requirement does not apply +if neither you nor any third party retains the ability to install +modified object code on the User Product (for example, the work has +been installed in ROM). + + The requirement to provide Installation Information does not include a +requirement to continue to provide support service, warranty, or updates +for a work that has been modified or installed by the recipient, or for +the User Product in which it has been modified or installed. Access to a +network may be denied when the modification itself materially and +adversely affects the operation of the network or violates the rules and +protocols for communication across the network. + + Corresponding Source conveyed, and Installation Information provided, +in accord with this section must be in a format that is publicly +documented (and with an implementation available to the public in +source code form), and must require no special password or key for +unpacking, reading or copying. + + 7. Additional Terms. + + "Additional permissions" are terms that supplement the terms of this +License by making exceptions from one or more of its conditions. +Additional permissions that are applicable to the entire Program shall +be treated as though they were included in this License, to the extent +that they are valid under applicable law. If additional permissions +apply only to part of the Program, that part may be used separately +under those permissions, but the entire Program remains governed by +this License without regard to the additional permissions. + + When you convey a copy of a covered work, you may at your option +remove any additional permissions from that copy, or from any part of +it. (Additional permissions may be written to require their own +removal in certain cases when you modify the work.) You may place +additional permissions on material, added by you to a covered work, +for which you have or can give appropriate copyright permission. + + Notwithstanding any other provision of this License, for material you +add to a covered work, you may (if authorized by the copyright holders of +that material) supplement the terms of this License with terms: + + a) Disclaiming warranty or limiting liability differently from the + terms of sections 15 and 16 of this License; or + + b) Requiring preservation of specified reasonable legal notices or + author attributions in that material or in the Appropriate Legal + Notices displayed by works containing it; or + + c) Prohibiting misrepresentation of the origin of that material, or + requiring that modified versions of such material be marked in + reasonable ways as different from the original version; or + + d) Limiting the use for publicity purposes of names of licensors or + authors of the material; or + + e) Declining to grant rights under trademark law for use of some + trade names, trademarks, or service marks; or + + f) Requiring indemnification of licensors and authors of that + material by anyone who conveys the material (or modified versions of + it) with contractual assumptions of liability to the recipient, for + any liability that these contractual assumptions directly impose on + those licensors and authors. + + All other non-permissive additional terms are considered "further +restrictions" within the meaning of section 10. If the Program as you +received it, or any part of it, contains a notice stating that it is +governed by this License along with a term that is a further +restriction, you may remove that term. If a license document contains +a further restriction but permits relicensing or conveying under this +License, you may add to a covered work material governed by the terms +of that license document, provided that the further restriction does +not survive such relicensing or conveying. + + If you add terms to a covered work in accord with this section, you +must place, in the relevant source files, a statement of the +additional terms that apply to those files, or a notice indicating +where to find the applicable terms. + + Additional terms, permissive or non-permissive, may be stated in the +form of a separately written license, or stated as exceptions; +the above requirements apply either way. + + 8. Termination. + + You may not propagate or modify a covered work except as expressly +provided under this License. Any attempt otherwise to propagate or +modify it is void, and will automatically terminate your rights under +this License (including any patent licenses granted under the third +paragraph of section 11). + + However, if you cease all violation of this License, then your +license from a particular copyright holder is reinstated (a) +provisionally, unless and until the copyright holder explicitly and +finally terminates your license, and (b) permanently, if the copyright +holder fails to notify you of the violation by some reasonable means +prior to 60 days after the cessation. + + Moreover, your license from a particular copyright holder is +reinstated permanently if the copyright holder notifies you of the +violation by some reasonable means, this is the first time you have +received notice of violation of this License (for any work) from that +copyright holder, and you cure the violation prior to 30 days after +your receipt of the notice. + + Termination of your rights under this section does not terminate the +licenses of parties who have received copies or rights from you under +this License. If your rights have been terminated and not permanently +reinstated, you do not qualify to receive new licenses for the same +material under section 10. + + 9. Acceptance Not Required for Having Copies. + + You are not required to accept this License in order to receive or +run a copy of the Program. Ancillary propagation of a covered work +occurring solely as a consequence of using peer-to-peer transmission +to receive a copy likewise does not require acceptance. However, +nothing other than this License grants you permission to propagate or +modify any covered work. These actions infringe copyright if you do +not accept this License. Therefore, by modifying or propagating a +covered work, you indicate your acceptance of this License to do so. + + 10. Automatic Licensing of Downstream Recipients. + + Each time you convey a covered work, the recipient automatically +receives a license from the original licensors, to run, modify and +propagate that work, subject to this License. You are not responsible +for enforcing compliance by third parties with this License. + + An "entity transaction" is a transaction transferring control of an +organization, or substantially all assets of one, or subdividing an +organization, or merging organizations. If propagation of a covered +work results from an entity transaction, each party to that +transaction who receives a copy of the work also receives whatever +licenses to the work the party's predecessor in interest had or could +give under the previous paragraph, plus a right to possession of the +Corresponding Source of the work from the predecessor in interest, if +the predecessor has it or can get it with reasonable efforts. + + You may not impose any further restrictions on the exercise of the +rights granted or affirmed under this License. For example, you may +not impose a license fee, royalty, or other charge for exercise of +rights granted under this License, and you may not initiate litigation +(including a cross-claim or counterclaim in a lawsuit) alleging that +any patent claim is infringed by making, using, selling, offering for +sale, or importing the Program or any portion of it. + + 11. Patents. + + A "contributor" is a copyright holder who authorizes use under this +License of the Program or a work on which the Program is based. The +work thus licensed is called the contributor's "contributor version". + + A contributor's "essential patent claims" are all patent claims +owned or controlled by the contributor, whether already acquired or +hereafter acquired, that would be infringed by some manner, permitted +by this License, of making, using, or selling its contributor version, +but do not include claims that would be infringed only as a +consequence of further modification of the contributor version. For +purposes of this definition, "control" includes the right to grant +patent sublicenses in a manner consistent with the requirements of +this License. + + Each contributor grants you a non-exclusive, worldwide, royalty-free +patent license under the contributor's essential patent claims, to +make, use, sell, offer for sale, import and otherwise run, modify and +propagate the contents of its contributor version. + + In the following three paragraphs, a "patent license" is any express +agreement or commitment, however denominated, not to enforce a patent +(such as an express permission to practice a patent or covenant not to +sue for patent infringement). To "grant" such a patent license to a +party means to make such an agreement or commitment not to enforce a +patent against the party. + + If you convey a covered work, knowingly relying on a patent license, +and the Corresponding Source of the work is not available for anyone +to copy, free of charge and under the terms of this License, through a +publicly available network server or other readily accessible means, +then you must either (1) cause the Corresponding Source to be so +available, or (2) arrange to deprive yourself of the benefit of the +patent license for this particular work, or (3) arrange, in a manner +consistent with the requirements of this License, to extend the patent +license to downstream recipients. "Knowingly relying" means you have +actual knowledge that, but for the patent license, your conveying the +covered work in a country, or your recipient's use of the covered work +in a country, would infringe one or more identifiable patents in that +country that you have reason to believe are valid. + + If, pursuant to or in connection with a single transaction or +arrangement, you convey, or propagate by procuring conveyance of, a +covered work, and grant a patent license to some of the parties +receiving the covered work authorizing them to use, propagate, modify +or convey a specific copy of the covered work, then the patent license +you grant is automatically extended to all recipients of the covered +work and works based on it. + + A patent license is "discriminatory" if it does not include within +the scope of its coverage, prohibits the exercise of, or is +conditioned on the non-exercise of one or more of the rights that are +specifically granted under this License. You may not convey a covered +work if you are a party to an arrangement with a third party that is +in the business of distributing software, under which you make payment +to the third party based on the extent of your activity of conveying +the work, and under which the third party grants, to any of the +parties who would receive the covered work from you, a discriminatory +patent license (a) in connection with copies of the covered work +conveyed by you (or copies made from those copies), or (b) primarily +for and in connection with specific products or compilations that +contain the covered work, unless you entered into that arrangement, +or that patent license was granted, prior to 28 March 2007. + + Nothing in this License shall be construed as excluding or limiting +any implied license or other defenses to infringement that may +otherwise be available to you under applicable patent law. + + 12. No Surrender of Others' Freedom. + + If conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot convey a +covered work so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you may +not convey it at all. For example, if you agree to terms that obligate you +to collect a royalty for further conveying from those to whom you convey +the Program, the only way you could satisfy both those terms and this +License would be to refrain entirely from conveying the Program. + + 13. Use with the GNU Affero General Public License. + + Notwithstanding any other provision of this License, you have +permission to link or combine any covered work with a work licensed +under version 3 of the GNU Affero General Public License into a single +combined work, and to convey the resulting work. The terms of this +License will continue to apply to the part which is the covered work, +but the special requirements of the GNU Affero General Public License, +section 13, concerning interaction through a network will apply to the +combination as such. + + 14. Revised Versions of this License. + + The Free Software Foundation may publish revised and/or new versions of +the GNU General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + + Each version is given a distinguishing version number. If the +Program specifies that a certain numbered version of the GNU General +Public License "or any later version" applies to it, you have the +option of following the terms and conditions either of that numbered +version or of any later version published by the Free Software +Foundation. If the Program does not specify a version number of the +GNU General Public License, you may choose any version ever published +by the Free Software Foundation. + + If the Program specifies that a proxy can decide which future +versions of the GNU General Public License can be used, that proxy's +public statement of acceptance of a version permanently authorizes you +to choose that version for the Program. + + Later license versions may give you additional or different +permissions. However, no additional obligations are imposed on any +author or copyright holder as a result of your choosing to follow a +later version. + + 15. Disclaimer of Warranty. + + THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY +APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT +HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY +OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, +THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM +IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF +ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. Limitation of Liability. + + IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS +THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY +GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE +USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF +DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD +PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), +EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF +SUCH DAMAGES. + + 17. Interpretation of Sections 15 and 16. + + If the disclaimer of warranty and limitation of liability provided +above cannot be given local legal effect according to their terms, +reviewing courts shall apply local law that most closely approximates +an absolute waiver of all civil liability in connection with the +Program, unless a warranty or assumption of liability accompanies a +copy of the Program in return for a fee. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +state the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + +Also add information on how to contact you by electronic and paper mail. + + If the program does terminal interaction, make it output a short +notice like this when it starts in an interactive mode: + + Copyright (C) + This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, your program's commands +might be different; for a GUI interface, you would use an "about box". + + You should also get your employer (if you work as a programmer) or school, +if any, to sign a "copyright disclaimer" for the program, if necessary. +For more information on this, and how to apply and follow the GNU GPL, see +. + + The GNU General Public License does not permit incorporating your program +into proprietary programs. If your program is a subroutine library, you +may consider it more useful to permit linking proprietary applications with +the library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. But first, please read +. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/cpp_bits/interface.cpp b/cpp_bits/interface.cpp new file mode 100644 --- /dev/null +++ b/cpp_bits/interface.cpp @@ -0,0 +1,754 @@ + +#include + +// ----------------------------------------------------------------------------- +// types + +typedef void Var; // Variable +typedef void CF; // CanonicalForm +typedef void Fac; // CFFactor = Factor +typedef void FacList; // CFFList = List +typedef void Iter; // ListIterator + +// ----------------------------------------------------------------------------- +// constants + +extern "C" int level_base() { return LEVELBASE; } +extern "C" int level_trans() { return LEVELTRANS; } +extern "C" int level_quot() { return LEVELQUOT; } +extern "C" int level_expr() { return LEVELEXPR; } + +// ----------------------------------------------------------------------------- +// version info + +extern "C" +void get_factory_version(char *str, int n) +{ + strncpy(str,FACTORYVERSION,n); + str[n-1]=0; +} + +extern "C" +void get_package_version(char *str, int n) +{ + strncpy(str,PACKAGE_VERSION,n); + str[n-1]=0; +} + +extern "C" int have_FLINT() { +#ifdef HAVE_FLINT + return 1; +#else + return 0; +#endif +} + +extern "C" int have_NTL() { +#ifdef HAVE_NTL + return 1; +#else + return 0; +#endif +} + +extern "C" int have_GMP() { +#ifdef HAVE_GMP + return 1; +#else + return 0; +#endif +} + +// ----------------------------------------------------------------------------- +// config + +extern "C" +void set_default_switches() +{ +#ifdef HAVE_NTL + // printf("we have NTL\n"); + On(SW_USE_CHINREM_GCD); + On(SW_USE_NTL_SORT); +#endif +#ifdef HAVE_FLINT + // printf("we have FLINT\n"); + #if defined(SW_USE_FL_GCD_P) + On(SW_USE_FL_GCD_P); // apparently 4.0.3 (on Debian 9) does not yet have these... + #endif + #if defined(SW_USE_FL_GCD_0) + On(SW_USE_FL_GCD_0); + #endif +#endif + On(SW_USE_EZGCD); + On(SW_USE_EZGCD_P); + On(SW_USE_QGCD); +} + +// ----------------------------------------------------------------------------- +// Variable + +extern "C" +void free_var (Var *ptr) +{ + Variable *varp = (Variable*) ptr; + delete varp; +} + +extern "C" +Var *new_var_level(int level) +{ + Variable *varp = new Variable(level); + return (void*)varp; +} + +extern "C" +Var *new_var_name(char name) +{ + Variable *varp = new Variable(name); + return (void*)varp; +} + +extern "C" +Var *new_var_level_name(int level, char name) +{ + Variable *varp = new Variable(level,name); + return (void*)varp; +} + +extern "C" +Var *root_of(CF *ptr) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + Variable *varp = new Variable( rootOf( *cfp ) ); + return (void*)varp; +} + +// ------------------------------------- + +extern "C" +int get_var_level(Var *ptr) +{ + Variable *varp = (Variable*) ptr; + return varp->level(); +} + +extern "C" +char get_var_name(Var *ptr) +{ + Variable *varp = (Variable*) ptr; + return varp->name(); +} + +extern "C" +int has_mipo(Var *ptr) +{ + Variable *varp = (Variable*) ptr; + return hasMipo( *varp ); +} + +extern "C" +CF *get_mipo(Var *ptr1 , Var *ptr2) +{ + Variable *varp1 = (Variable*) ptr1; + Variable *varp2 = (Variable*) ptr2; + CanonicalForm *cfp = new CanonicalForm( getMipo( *varp1 , *varp2 ) ); + return cfp; +} + +extern "C" +void set_mipo(Var *ptr1 , CF *ptr2) +{ + Variable *varp = (Variable *) ptr1; + CanonicalForm *cfp = (CanonicalForm*) ptr2; + setMipo( *varp, *cfp ); +} + +extern "C" +void set_reduce(Var *ptr, int flag) +{ + Variable *varp = (Variable*) ptr; + setReduce( *varp , flag); +} + +// ----------------------------------------------------------------------------- +// CanonicalForm + +extern "C" +void free_cf (CF *ptr) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + delete cfp; +} + +extern "C" +CF *copy_cf(CF *ptr) +{ + CanonicalForm *cfp1 = (CanonicalForm*) ptr; + CanonicalForm *cfp2 = new CanonicalForm(*cfp1); + return (void*)cfp2; +} + +// ----------------------------------------------------------------------------- +// Factor + +extern "C" +void free_fac (CFFactor *ptr) +{ + CFFactor *facp = (CFFactor*) ptr; + delete facp; +} + +extern "C" +Fac *copy_fac(Fac *ptr) +{ + CFFactor *facp1 = (CFFactor*) ptr; + CFFactor *facp2 = new CFFactor(*facp1); + return (void*)facp2; +} + +extern "C" +int get_fac_expo(Fac *ptr) +{ + CFFactor *facp = (CFFactor*) ptr; + return facp->exp(); +} + +extern "C" +CF *get_factor(Fac *ptr) +{ + CFFactor *facp = (CFFactor*) ptr; + CanonicalForm *cfp = new CanonicalForm( facp->factor() ); + return (void*)cfp; +} + +// ----------------------------------------------------------------------------- +// Iterator + +extern "C" +void free_iter(Iter *ptr) +{ + ListIterator *iterp = (ListIterator *) ptr; + delete iterp; +} + +extern "C" +Iter *next_iter(Iter *ptr) +{ + ListIterator *iterp1 = (ListIterator *) ptr; + (*iterp1) ++; + return iterp1; +} + +// ----------------------------------------------------------------------------- +// List + +extern "C" +void free_faclist (FacList *ptr) +{ + CFFList *listp = (CFFList*) ptr; + delete listp; +} + +extern "C" +int get_list_length(FacList *ptr) +{ + CFFList *listp = (CFFList*) ptr; + return listp->length(); +} + +extern "C" +void flatten_faclist (FacList *ptr, Fac **tgt_arr) +{ + CFFList *listp = (CFFList*) ptr; + ListIterator iter = ListIterator(*listp); + int n = listp->length(); + for(int i=0 ; iisZero()); +} + +extern "C" +int is_one(CF *ptr) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + return (cfp->isOne()); +} + +extern "C" +int is_imm(CF *ptr) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + return (cfp->isImm()); +} + +extern "C" +int is_univariate(CF *ptr) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + return (cfp->isUnivariate()); +} + +// ----------------------------------------------------------------------------- +// base domain predicates + +extern "C" +int in_ZZ(CF *ptr) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + return (cfp->inZ()); +} + +extern "C" +int in_QQ(CF *ptr) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + return (cfp->inQ()); +} + +extern "C" +int in_GF(CF *ptr) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + return (cfp->inGF()); +} + +extern "C" +int in_FF(CF *ptr) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + return (cfp->inFF()); +} + +// ----------------------------------------------------------------------------- +// domain predicates + +extern "C" +int in_BaseDomain(CF *ptr) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + return (cfp->inBaseDomain()); +} + +extern "C" +int in_CoeffDomain(CF *ptr) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + return (cfp->inCoeffDomain()); +} + +extern "C" +int in_PolyDomain(CF *ptr) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + return (cfp->inPolyDomain()); +} + +extern "C" +int in_Extension(CF *ptr) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + return (cfp->inExtension()); +} + +extern "C" +int in_QuotDomain(CF *ptr) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + return (cfp->inQuotDomain()); +} + +// ----------------------------------------------------------------------------- +// imm_value, degree, level + +extern "C" +long int smallint_value(CF *ptr) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + return (cfp->intval()); // returns a long!!! +} + +extern "C" +int degree_of(CF *ptr) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + return (cfp->degree()); +} + +extern "C" +int level_of(CF *ptr) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + return (cfp->level()); +} + +extern "C" +Var *mvar_of(CF *ptr) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + Variable *varp = new Variable( cfp->mvar() ); + return varp; +} + +// ----------------------------------------------------------------------------- + +extern "C" +CF *index_poly(CF *ptr, int idx) +{ + CanonicalForm *cfp1 = (CanonicalForm*) ptr; + CanonicalForm *cfp2 = new CanonicalForm( (*cfp1)[idx] ); + return cfp2; +} + +extern "C" +CF *map_into(CF *ptr) +{ + CanonicalForm *cfp1 = (CanonicalForm*) ptr; + CanonicalForm *cfp2 = new CanonicalForm( cfp1->mapinto() ); + return cfp2; +} + +extern "C" +CF *substitute(CF *ptr1, Var *ptrv, CF *ptr2) +{ + CanonicalForm *cfp1 = (CanonicalForm*) ptr1; + Variable *varp = (Variable* ) ptrv; + CanonicalForm *cfp2 = (CanonicalForm*) ptr2; + CanonicalForm *cfp3 = new CanonicalForm( (*cfp1)( *cfp2 , *varp ) ); + return cfp3; +} + +// ----------------------------------------------------------------------------- +// numerator and denominator + +extern "C" +CF *numer(CF *ptr) +{ + CanonicalForm *cfp1 = (CanonicalForm*) ptr; + CanonicalForm *cfp2 = new CanonicalForm( cfp1->num() ); + return cfp2; +} + +extern "C" +CF *denom(CF *ptr) +{ + CanonicalForm *cfp1 = (CanonicalForm*) ptr; + CanonicalForm *cfp2 = new CanonicalForm( cfp1->den() ); + return cfp2; +} + +// ----------------------------------------------------------------------------- +// binary operations + +extern "C" +int is_equal(CF *ptr1, CF* ptr2) +{ + CanonicalForm *cfp1 = (CanonicalForm*) ptr1; + CanonicalForm *cfp2 = (CanonicalForm*) ptr2; + return ( (*cfp1) == (*cfp2) ); +} + +extern "C" +CF *plus_cf(CF *ptr1, CF* ptr2) +{ + CanonicalForm *cfp1 = (CanonicalForm*) ptr1; + CanonicalForm *cfp2 = (CanonicalForm*) ptr2; + CanonicalForm *cfp3 = new CanonicalForm( (*cfp1) + (*cfp2) ); + return cfp3; +} + +extern "C" +CF *minus_cf(CF *ptr1, CF* ptr2) +{ + CanonicalForm *cfp1 = (CanonicalForm*) ptr1; + CanonicalForm *cfp2 = (CanonicalForm*) ptr2; + CanonicalForm *cfp3 = new CanonicalForm( (*cfp1) - (*cfp2) ); + return cfp3; +} + +extern "C" +CF *times_cf(CF *ptr1, CF* ptr2) +{ + CanonicalForm *cfp1 = (CanonicalForm*) ptr1; + CanonicalForm *cfp2 = (CanonicalForm*) ptr2; + CanonicalForm *cfp3 = new CanonicalForm( (*cfp1) * (*cfp2) ); + return cfp3; +} + +extern "C" +CF *pow_cf(CF *ptr1, int n) +{ + CanonicalForm *cfp1 = (CanonicalForm*) ptr1; + CanonicalForm *cfp3 = new CanonicalForm( power( *cfp1 , n ) ); + return cfp3; +} + +extern "C" +CF *div_cf(CF *ptr1, CF* ptr2) +{ + CanonicalForm *cfp1 = (CanonicalForm*) ptr1; + CanonicalForm *cfp2 = (CanonicalForm*) ptr2; + CanonicalForm *cfp3 = new CanonicalForm( div( *cfp1 , *cfp2 ) ); + return cfp3; +} + +extern "C" +CF *mod_cf(CF *ptr1, CF* ptr2) +{ + CanonicalForm *cfp1 = (CanonicalForm*) ptr1; + CanonicalForm *cfp2 = (CanonicalForm*) ptr2; + CanonicalForm *cfp3 = new CanonicalForm( mod( *cfp1 , *cfp2 ) ); + return cfp3; +} + +extern "C" +CF *gcd_poly_cf(CF *ptr1, CF* ptr2) +{ + CanonicalForm *cfp1 = (CanonicalForm*) ptr1; + CanonicalForm *cfp2 = (CanonicalForm*) ptr2; + CanonicalForm *cfp3 = new CanonicalForm( gcd_poly( *cfp1 , *cfp2 ) ); + return cfp3; +} + +extern "C" +CF *reduce_cf(CF *ptr1, CF* ptr2) +{ + CanonicalForm *cfp1 = (CanonicalForm*) ptr1; + CanonicalForm *cfp2 = (CanonicalForm*) ptr2; + CanonicalForm *cfp3 = new CanonicalForm( reduce( *cfp1 , *cfp2 ) ); + return cfp3; +} + +// ----------------------------------------------------------------------------- +// GMP conversion + +extern "C" +int get_gf_value(CF *ptr) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + return gf_value( *cfp ); +} + +extern "C" +void get_gmp_numerator(CF *ptr, mpz_ptr tgt) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + gmp_numerator( *cfp , tgt ); +} + +extern "C" +void get_gmp_denominator(CF *ptr, mpz_ptr tgt) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + gmp_denominator( *cfp , tgt ); +} + +extern "C" +CF *make_ZZ_from_gmp(mpz_srcptr src) +{ + mpz_t copy; + mpz_init_set(copy,src); + // make_cf DOES NOT COPY, just copies the header (mpz_t structure)... + // ...so we have to copy it here because our src is ephemeral. + CanonicalForm *cfp = new CanonicalForm( make_cf(copy) ); + return cfp; +} + +extern "C" +CF *make_QQ_from_gmp(mpz_srcptr src1, mpz_srcptr src2, int normalize_flag ) +{ + mpz_t copy1; + mpz_t copy2; + mpz_init_set(copy1,src1); + mpz_init_set(copy2,src2); + CanonicalForm *cfp = new CanonicalForm( make_cf( copy1 , copy2 , normalize_flag ) ); + return cfp; +} + +extern "C" +CF *make_GF(int z) +{ + CanonicalForm *cfp = new CanonicalForm( make_cf_from_gf(z) ); + return cfp; +} + + +/* +void gmp_numerator ( const CanonicalForm & f, mpz_ptr result ); + +void gmp_denominator ( const CanonicalForm & f, mpz_ptr result ); + +int gf_value (const CanonicalForm & f ); + +CanonicalForm make_cf ( const mpz_ptr n ); + +CanonicalForm make_cf ( const mpz_ptr n, const mpz_ptr d, bool normalize ); + +CanonicalForm make_cf_from_gf ( const int z ); +*/ + + +// ----------------------------------------------------------------------------- + +extern "C" +int get_characteristic() +{ + return getCharacteristic(); +} + +// prime fields and QQ +extern "C" +void set_characteristic1(int c) +{ + setCharacteristic(c); +} + +/* +// prime power fields (meaning FF?) +// apparently this does not exist anymore?!?! +extern "C" +void set_characteristic2(int c, int n) +{ + setCharacteristic(c,n); +} +*/ + +// Galois fields +extern "C" +void set_characteristic3(int c, int n, char var) +{ + setCharacteristic(c,n,var); +} + +/* +// already defined above +extern "C" +int get_gf_value(CF *ptr) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + return gf_value( cfp ); +} +*/ + +extern "C" +int is_FF_in_GF(CF *ptr) +{ + CanonicalForm *cfp = (CanonicalForm*) ptr; + return cfp -> isFFinGF(); +} + +extern "C" +int get_GF_degree() +{ + return getGFDegree(); +} + +extern "C" +CF *get_GF_generator() +{ + return new CanonicalForm ( getGFGenerator() ); +} + + +// ----------------------------------------------------------------------------- diff --git a/singular-factory.cabal b/singular-factory.cabal new file mode 100644 --- /dev/null +++ b/singular-factory.cabal @@ -0,0 +1,82 @@ +Name: singular-factory +Version: 0.1 +Synopsis: Multivariate polynomial factorization via bindings to Singular-factory +Description: Bindings to Singular-factory. Factory is the subset of the Singular + computer algebra system which does factorization of multivariate + polynomials (over integers, rationals, and finite fields) +License: GPL +License-file: LICENSE +Author: Balazs Komuves +Copyright: (c) 2019 Balazs Komuves +Maintainer: bkomuves (plus) hackage (at) gmail (dot) com +homepage: https://github.com/bkomuves/singular-factory +Stability: Experimental +Category: Math +Tested-With: GHC == 8.6.5 +Cabal-Version: 2.0 +Build-Type: Simple + +Source-repository head + type: git + location: https://github.com/bkomuves/singular-factory + +flag pkg-cfg-prefix + description: whether the pkg-config factory package name has a prefix + manual: False + +-------------------------------------------------------------------------------- + +Library + Build-Depends: base >= 4 && < 5, containers >= 0.5, text >= 1.0, + filepath >= 1.0, directory >= 1.3, process >= 1.2, + hgmp >= 0.1.1 + + Exposed-Modules: Math.Singular.Factory + Math.Singular.Factory.Domains + Math.Singular.Factory.Variables + Math.Singular.Factory.Polynomial + Math.Singular.Factory.Expr + Math.Singular.Factory.Parser + Math.Singular.Factory.Counting + Math.Singular.Factory.GFTables + Math.Singular.Factory.Internal.CanonicalForm + Math.Singular.Factory.Internal.Factory + Math.Singular.Factory.Internal.DList + + Default-Language: Haskell2010 + + Default-extensions: BangPatterns, TypeSynonymInstances, FlexibleInstances, FlexibleContexts, + KindSignatures, DataKinds, EmptyDataDecls, ScopedTypeVariables + + Other-extensions: ForeignFunctionInterface, CPP, PatternSynonyms + + Hs-Source-Dirs: src + + cc-options: -x c++ + c-sources: cpp_bits/interface.cpp + extra-libraries: stdc++ + + if flag(pkg-cfg-prefix) + pkgconfig-depends: singular-factory + else + pkgconfig-depends: factory + +-------------------------------------------------------------------------------- + +test-suite singular-factory-tests + + type: exitcode-stdio-1.0 + + hs-source-dirs: test + main-is: TestSuite.hs + + build-depends: base >= 4 && < 5, random, + singular-factory >= 0.1, + QuickCheck >= 2, + tasty, tasty-quickcheck, tasty-hunit + + Default-Language: Haskell2010 + Default-Extensions: CPP, BangPatterns + +-------------------------------------------------------------------------------- + diff --git a/src/Math/Singular/Factory.hs b/src/Math/Singular/Factory.hs new file mode 100644 --- /dev/null +++ b/src/Math/Singular/Factory.hs @@ -0,0 +1,88 @@ +{- | + +Haskell bindings to the Singular-Factory multivariate polynomial factorization engine. + +Singular is a computer algebra system developed at the University of Kaiserslautern, +and Factory is the polynomial factorization engine of Singular. + +It handles multivariate polynomials over the integers, rationals, and finite fields +(and also algebraic and trancendental extensions of these, but those are not yet +supported). + +Links: + +* Singular: https://www.singular.uni-kl.de/ + +* Factory: https://www.singular.uni-kl.de/dox/html/factory_page.html + +-} + +module Math.Singular.Factory + ( -- * Initialization + initialize + -- * Configuration + , factoryVersion + , factoryConfig , FactoryConfig(..) + , printVersion + , printConfig + -- * Re-exported modules + , module Math.Singular.Factory.Domains + , module Math.Singular.Factory.Variables + , module Math.Singular.Factory.Polynomial + ) + where + +-------------------------------------------------------------------------------- + +import Math.Singular.Factory.Domains +import Math.Singular.Factory.Variables +import Math.Singular.Factory.Polynomial + +import Math.Singular.Factory.GFTables +import Math.Singular.Factory.Internal.Factory +import Math.Singular.Factory.Internal.CanonicalForm + +-------------------------------------------------------------------------------- + +-- | You are supposed to call 'initialize' before doing anything +-- (though factorization over characteristic zeros should work without it) +initialize = do + initGFTables + setDefaultSwitches + +-------------------------------------------------------------------------------- + +data FactoryConfig = FactoryConfig + { withFLINT :: Bool + , withNTL :: Bool + , withGMP :: Bool + } + deriving (Eq,Show) + +factoryConfig :: FactoryConfig +factoryConfig = FactoryConfig + { withFLINT = haveFLINT + , withNTL = haveNTL + , withGMP = haveGMP + } + +-------------------------------------------------------------------------------- + +-- | Prints factory\'s version +printVersion = do + putStrLn $ "factory version = " ++ factoryVersion + +-- | Prints some more detailed configuration info +printConfig = do + putStrLn $ "factory version = " ++ factoryVersion + putStrLn $ "package version = " ++ packageVersion + mb_gftables <- getGFTablesDir + putStrLn $ "gftables folder = " ++ show mb_gftables + putStrLn $ "have FLINT = " ++ show haveFLINT + putStrLn $ "have NTL = " ++ show haveNTL + putStrLn $ "have GMP = " ++ show haveGMP + +-------------------------------------------------------------------------------- + + + diff --git a/src/Math/Singular/Factory/Counting.hs b/src/Math/Singular/Factory/Counting.hs new file mode 100644 --- /dev/null +++ b/src/Math/Singular/Factory/Counting.hs @@ -0,0 +1,75 @@ + +-- | Brute-force counting of solutions of polynomial equations over (small) finite field + +{-# LANGUAGE BangPatterns, TypeApplications, ScopedTypeVariables, DataKinds #-} +module Math.Singular.Factory.Counting where + +-------------------------------------------------------------------------------- + +import Data.Proxy + +import Math.Singular.Factory.Domains +import Math.Singular.Factory.Polynomial +import Math.Singular.Factory.GFTables + +-------------------------------------------------------------------------------- +-- * Hypersurfaces + +-- | Count points of a hypersurface in an affine space +-- over a finite field. +-- +-- The @Int@ input is the number of variables (that is, the dimension). +-- +countAffineHypersurface :: FiniteDomain domain => Int -> Polynomial vars domain -> Int +countAffineHypersurface = go where + go !1 !p = countTrues [ polyIsZero (substitute1 1 (konst a) p) | a <- enumerateDomain ] + go !k !p = sum [ go (k-1) (substitute1 k (konst a) p) | a <- enumerateDomain ] + +-- | Count points of a hypersurface in an projective space +-- over a finite field. +-- +-- The @Int@ input is the number of /variables/ (that is, the dimension plus one)! +-- +-- NOTE: We assume that the input is a homogeneous polynomial, but this is not checked!! +countProjectiveHypersurface :: FiniteDomain domain => Int -> Polynomial vars domain -> Int +countProjectiveHypersurface = go False where + go !isPos !1 !p = countTrues [ polyIsZero $ substitute1 1 (konst a) p | a <- projEnumerateDomain isPos , let b = isPos || not (isZero a) , b ] + go !isPos !k !p = sum [ go b (k-1) $ substitute1 k (konst a) p | a <- projEnumerateDomain isPos , let b = isPos || not (isZero a) ] + +-------------------------------------------------------------------------------- +-- * Systems of equations (varieties) + +-- | Count solutions of a system of polynomial equations in an affine space +-- over a finite field. +-- +-- The @Int@ input is the number of variables (that is, the dimension). +countAffineSolutions :: FiniteDomain domain => Int -> [Polynomial vars domain] -> Int +countAffineSolutions = go where + go !1 !ps = countTrues [ and $ map polyIsZero $ map (substitute1 1 (konst a)) ps | a <- enumerateDomain ] + go !k !ps = sum [ go (k-1) $ map (substitute1 k (konst a)) ps | a <- enumerateDomain ] + +-- | Count solutions of a system of polynomial equations in a projective space +-- over a finite field. +-- +-- The @Int@ input is the number of /variables/ (that is, the dimension plus one)! +-- +-- NOTE: We assume that the input is a list of homogeneous polynomial, but this is not checked!! +countProjectiveSolutions :: FiniteDomain domain => Int -> [Polynomial vars domain] -> Int +countProjectiveSolutions = go False where + go !isPos !1 !ps = countTrues [ and $ map polyIsZero $ map (substitute1 1 (konst a)) ps | a <- projEnumerateDomain isPos , let b = isPos || not (isZero a) , b ] + go !isPos !k !ps = sum [ go b (k-1) $ map (substitute1 k (konst a)) ps | a <- projEnumerateDomain isPos , let b = isPos || not (isZero a) ] + +-------------------------------------------------------------------------------- +-- * Misc helpers + +-- | if the input is False, return [0,1], otherwise all elements of the domain +projEnumerateDomain :: FiniteDomain domain => Bool -> [domain] +projEnumerateDomain True = enumerateDomain +projEnumerateDomain False = [0,1] + +countTrues :: [Bool] -> Int +countTrues = foldr f 0 where + f False !c = c + f True !c = c+1 + +-------------------------------------------------------------------------------- diff --git a/src/Math/Singular/Factory/Domains.hs b/src/Math/Singular/Factory/Domains.hs new file mode 100644 --- /dev/null +++ b/src/Math/Singular/Factory/Domains.hs @@ -0,0 +1,324 @@ + +-- | Base domains. +-- +-- These are the base rings and fields Factory can work with, namely: +-- +-- * the ring integers +-- * the field of rationals +-- * finite fields (prime fields and Galois fields) +-- +-- Another representation for finite fields are explicit algebraic extensions +-- of prime fields. This has less limitations (does not rely on precomputed +-- tables), but it is not implemented yet. +-- +-- Note1: non-prime order Galois fields are supported only for small orders! +-- (this is a limitation by singular-factory). +-- Also for them to work, we need to be able to figure out the location of the +-- \"gftables\" directory first. +-- +-- Note2: as Factory has the base domain as a global state (...), this whole +-- library is not at all thread safe! +-- + +{-# LANGUAGE + BangPatterns, PatternSynonyms, KindSignatures, DataKinds, + FlexibleInstances, TypeSynonymInstances + #-} +module Math.Singular.Factory.Domains where + +-------------------------------------------------------------------------------- + +import Data.Ratio + +import GHC.TypeLits +import Data.Proxy +import Data.IORef + +import Control.Monad +import System.IO.Unsafe as Unsafe + +import Math.Singular.Factory.Internal.CanonicalForm +import Math.Singular.Factory.Internal.Factory hiding ( FF , GF ) + +-------------------------------------------------------------------------------- +-- * The global characteristics + +type Prime = Int + +data FactoryChar + = CharZero -- ^ QQ + | CharFp !Prime -- ^ prime field + | CharGF !Prime !Int -- ^ Galois field + deriving (Eq,Show) + +-- | Unfortunately, Factory maintains a global state... +theFactoryChar :: IORef FactoryChar +theFactoryChar = Unsafe.unsafePerformIO $ newIORef CharZero + +setFactoryChar :: FactoryChar -> IO () +setFactoryChar new = do + old <- readIORef theFactoryChar + when (new /= old) $ do + writeIORef theFactoryChar new + case new of + CharZero -> setCharacteristic1 0 + CharFp p -> setCharacteristic1 p + CharGF p n -> setCharacteristic3 p n '@' -- we use '@' for the Galois field variable for now + +mapIntoCF :: FactoryChar -> CF -> CF +mapIntoCF char cf = Unsafe.unsafePerformIO $ do + setFactoryChar char + mapIntoIO cf + +-------------------------------------------------------------------------------- +-- * Prime fields + +-- | Haskell prime fields (this is provided only for completeness) +newtype Fp (p :: Nat) + = Fp Int + deriving (Eq) + +fpPrime :: KnownNat p => Fp p -> Int +fpPrime = fromIntegral . natVal . proxyP where + proxyP :: Fp p -> Proxy p + proxyP _ = Proxy + +modp :: (KnownNat p, Integral a) => a -> Fp p +modp k = fp where + fp = Fp $ mod (fromIntegral k) (fpPrime fp) + +-- proxyFp :: Fp p -> Proxy p +-- proxyFp _ = Proxy + +instance Show (Fp p) where + show (Fp k) = show k + +instance KnownNat p => Num (Fp p) where + fromInteger = modp + negate fp@(Fp k) = if k==0 then Fp 0 else Fp (fpPrime fp - k) + (Fp a) + (Fp b) = modp (a+b) + (Fp a) - (Fp b) = modp (a-b) + (Fp a) * (Fp b) = modp (a*b) + abs = id + signum = const (Fp 1) + +fpToFF :: KnownNat p => Fp p -> FF p +fpToFF (Fp k) = mkFF k + +-------------------------------------------------------------------------------- +-- * Finite fields + +-- | Factory prime fields +newtype FF (p :: Nat) + = FF { unFF :: CF } + +mkFF :: (KnownNat p, Integral a) => a -> FF p +mkFF !k = ff where + ff = Unsafe.unsafePerformIO $ do + setBaseDomain (mkProxy ff) + cf <- (mapIntoIO =<< makeIntegerCF (fromIntegral k)) + return (FF cf) + +instance Eq (FF p) where + (FF a) == (FF b) = nativeEqCF a b + +instance Show (FF p) where + show (FF cf) = show (valueFF cf) -- "[" ++ show (valueFF cf) ++ "]" + +ffPrime :: KnownNat p => FF p -> Int +ffPrime = fromIntegral . natVal . proxyP where + proxyP :: FF p -> Proxy p + proxyP _ = Proxy + +instance KnownNat p => Num (FF p) where + fromInteger = mkFF + negate (FF a) = FF (negate a) + (FF a) + (FF b) = FF (a + b) + (FF a) - (FF b) = FF (a - b) + (FF a) * (FF b) = FF (a * b) + abs = id + signum = const 1 + +instance (KnownNat p) => Fractional (FF p) where + fromRational q = mkFF (numerator q) / mkFF (denominator q) + (FF a) / (FF b) = FF (divCF a b) + +-------------------------------------------------------------------------------- +-- * Galois fields + +-- | Galois fields @GF(p^n)@. +-- +-- The (nonzero) elements are represented as powers of the canonical generator. +-- +-- The symbol is the name of the canonical generator (used for pretty-printing). +-- +-- Note: because of how Factory is implemented, it is /required/ that @n >= 2@... +-- (use 'FF' for prime fields). +-- +-- Also, the sizes are really limited (because they rely on tables, and only +-- small tables are included): @p < 256@ and @p^n < 65536@ +-- +newtype GF (p :: Nat) (n :: Nat) (x :: Symbol) + = GF { unGF :: CF } + +-- | Create elements of the prime subfield. For the rest, you can use the +-- powers of the generator. +mkGF :: (KnownNat p, KnownNat n, KnownSymbol x, Integral a) => a -> GF p n x +mkGF !k = gf where + gf = Unsafe.unsafePerformIO $ do + setBaseDomain (mkProxy gf) + cf <- (mapIntoIO =<< makeIntegerCF (fromIntegral k)) + return (GF cf) + +-- | The canonical generator of the (multiplicative group of the) Galois field +genGF :: (KnownNat p, KnownNat n, KnownSymbol x) => GF p n x +genGF = gf where + gf = Unsafe.unsafePerformIO $ do + setBaseDomain (mkProxy gf) + cf <- getGFGenerator + return (GF cf) + +-- | A power of the canonical generator +genPowGF :: (KnownNat p, KnownNat n, KnownSymbol x) => Int -> GF p n x +genPowGF e = powGF genGF e + +powGF :: (KnownNat p, KnownNat n, KnownSymbol x) => GF p n x -> Int -> GF p n x +powGF (GF cf) e = GF (powCF cf e) + +instance Eq (GF p n x) where + (GF a) == (GF b) = nativeEqCF a b + +instance (KnownNat p, KnownNat n, KnownSymbol x) => Show (GF p n x) where + show gf@(GF cf) = showGFValue1 (gfSymbol gf) (valueGF cf) + +instance (KnownNat p, KnownNat n, KnownSymbol x) => Num (GF p n x) where + fromInteger = mkGF + negate (GF a) = GF (negate a) + (GF a) + (GF b) = GF (a + b) + (GF a) - (GF b) = GF (a - b) + (GF a) * (GF b) = GF (a * b) + abs = id + signum = const 1 + +instance (KnownNat p, KnownNat n, KnownSymbol x) => Fractional (GF p n x) where + fromRational q = mkGF (numerator q) / mkGF (denominator q) + (GF a) / (GF b) = GF (divCF a b) + +gfPrime :: KnownNat p => GF p n x -> Int +gfPrime = fromIntegral . natVal . proxyP where + proxyP :: GF p n x -> Proxy p + proxyP _ = Proxy + +gfExponent :: KnownNat n => GF p n x -> Int +gfExponent = fromIntegral . natVal . proxyE where + proxyE :: GF p n x -> Proxy n + proxyE _ = Proxy + +gfSymbol :: KnownSymbol x => GF p n x -> String +gfSymbol = symbolVal . proxyX where + proxyX :: GF p n x -> Proxy x + proxyX _ = Proxy + + +-------------------------------------------------------------------------------- +-- * Base domains + +class (Eq a, Show a, Num a) => BaseDomain a where + characteristic :: Proxy a -> Int + charExponent :: Proxy a -> Int + baseDomainName :: Proxy a -> String + factoryChar :: Proxy a -> FactoryChar + baseToCF :: a -> CF + unsafeCfToBase :: CF -> a + isZero :: a -> Bool + isOne :: a -> Bool + +setBaseDomain :: BaseDomain a => Proxy a -> IO () +setBaseDomain = setFactoryChar . factoryChar + +instance BaseDomain Integer where + characteristic _ = 0 + charExponent _ = 1 + baseDomainName _ = "ZZ" + factoryChar _ = CharZero + baseToCF x = Unsafe.unsafePerformIO (makeIntegerCF x) + unsafeCfToBase = valueZZ + isZero n = (n == 0) + isOne n = (n == 1) + +instance BaseDomain Rational where + characteristic _ = 0 + charExponent _ = 1 + baseDomainName _ = "QQ" + factoryChar _ = CharZero + baseToCF x = Unsafe.unsafePerformIO (makeRationalCF x) + unsafeCfToBase = valueQQ + isZero q = (q == 0) + isOne q = (q == 1) + +instance KnownNat p => BaseDomain (Fp p) where + characteristic pxy = (fpPrime $ proxyUndef pxy) + charExponent _ = 1 + baseDomainName pxy = "F_" ++ show (characteristic pxy) + factoryChar pxy = CharFp (characteristic pxy) + baseToCF x = baseToCF (fpToFF x) + unsafeCfToBase cf = Fp (valueFF cf) + isZero (Fp k) = (k == 0) + isOne (Fp k) = (k == 1) + +instance KnownNat p => BaseDomain (FF p) where + characteristic pxy = (ffPrime $ proxyUndef pxy) + charExponent pxy = 1 + baseDomainName pxy = "FF(" ++ show (characteristic pxy) ++ ")" where + factoryChar pxy = CharFp (characteristic pxy) + baseToCF (FF cf) = cf + unsafeCfToBase = FF + isZero (FF cf) = isZeroCF cf + isOne (FF cf) = isOneCF cf + +instance (KnownNat p, KnownNat n, KnownSymbol x) => BaseDomain (GF p n x) where + characteristic pxy = (gfPrime $ proxyUndef pxy) + charExponent pxy = (gfExponent $ proxyUndef pxy) + baseDomainName pxy = "GF(" ++ show (characteristic pxy) ++ ")" where + factoryChar pxy = CharGF (characteristic pxy) (charExponent pxy) + baseToCF (GF cf) = cf + unsafeCfToBase = GF + isZero (GF cf) = isZeroCF cf + isOne (GF cf) = isOneCF cf + +-------------------------------------------------------------------------------- +-- * Finite domains + +class BaseDomain domain => FiniteDomain domain where + domainSize :: Proxy domain -> Int + enumerateDomain :: [domain] + +instance KnownNat p => FiniteDomain (Fp p) where + domainSize pxy = characteristic pxy + enumerateDomain = list where + list = [ Fp i | i<-[0..p-1] ] + p = characteristic $ mkProxy (head list) + +instance KnownNat p => FiniteDomain (FF p) where + domainSize pxy = characteristic pxy + enumerateDomain = list where + list = [ mkFF i | i<-[0..p-1] ] + p = characteristic $ mkProxy (head list) + +instance (KnownNat p, KnownNat n, KnownSymbol x) => FiniteDomain (GF p n x) where + domainSize pxy = characteristic pxy ^ charExponent pxy + enumerateDomain = list where + list = 0 : [ genPowGF i | i<-[0..n-2] ] where + pxy = mkProxy (head list) + n = domainSize pxy + +-------------------------------------------------------------------------------- +-- * Proxy + +mkProxy :: a -> Proxy a +mkProxy _ = Proxy + +proxyUndef :: Proxy a -> a +proxyUndef _ = undefined + +-------------------------------------------------------------------------------- diff --git a/src/Math/Singular/Factory/Expr.hs b/src/Math/Singular/Factory/Expr.hs new file mode 100644 --- /dev/null +++ b/src/Math/Singular/Factory/Expr.hs @@ -0,0 +1,111 @@ + +-- | Polynomial expressions (used for parsing) + +{-# LANGUAGE BangPatterns, DeriveFunctor, DeriveFoldable, DeriveTraversable #-} +module Math.Singular.Factory.Expr where + +-------------------------------------------------------------------------------- + +import Data.Char +import Data.List + +import Control.Applicative +import Control.Monad + +import Data.Text.Lazy ( Text ) +import qualified Data.Text.Lazy as T +import qualified Data.Text.Lazy.Read as T + +import Math.Singular.Factory.Internal.DList as DList + +-------------------------------------------------------------------------------- +-- * Types + +data Sign + = Plus + | Minus + deriving (Eq,Ord,Show) + +negateIfMinus :: Num a => Sign -> a -> a +negateIfMinus Plus = id +negateIfMinus Minus = negate + +-------------------------------------------------------------------------------- + +-- | Monomials +newtype Monom var + = Monom [(var,Int)] + deriving (Eq,Ord,Show,Functor,Foldable,Traversable) + +-- | A monomial multiplied by a constant +data Term coeff var + = Term !coeff !(Monom var) + deriving (Eq,Ord,Show,Functor,Foldable,Traversable) + +-- | Polynomials as linear combination of monomials +newtype GenPoly coeff var + = GenPoly [Term coeff var] + deriving (Eq,Ord,Show,Functor,Foldable,Traversable) + +-- | Polynomial expressions +data Expr v + = VarE !v + | KstE !Integer + | NegE (Expr v) + | LinE [(Sign,Expr v)] + | MulE [Expr v] + | PowE (Expr v) !Int + deriving (Eq,Ord,Show,Functor,Foldable,Traversable) + +-------------------------------------------------------------------------------- +-- * Evaluation to @Num@ + +evalSign :: Num c => Sign -> c -> c +evalSign Plus = id +evalSign Minus = negate + +evalMonom :: Num c => (var -> c) -> Monom var -> c +evalMonom f (Monom list) = product (map g list) where + g (var,expo) = (f var)^expo + +evalTerm :: Num c => (coeff -> c) -> (var -> c) -> Term coeff var -> c +evalTerm f g (Term coeff monom) = f coeff * evalMonom g monom + +evalGenPoly :: Num c => (coeff -> c) -> (var -> c) -> GenPoly coeff var -> c +evalGenPoly f g (GenPoly terms) = sum (map (evalTerm f g) terms) + +evalExpr :: Num c => (var -> c) -> Expr var -> c +evalExpr evalVar = go where + go expr = case expr of + VarE v -> evalVar v + KstE k -> fromInteger k + NegE e -> negate (go e) + LinE xs -> sum [ evalSign pm (go x) | (pm,x) <- xs ] + MulE xs -> product (map go xs) + PowE e k -> (go e)^k + +-------------------------------------------------------------------------------- + +{- + +prettExpr :: Expr -> String +prettExpr expr = DList.toList (prettPrecExpr 0 expr) + +prettyPrecExpr :: Int -> Expr -> DList Char +prettyPrecExpr = go where + + chr c = DList.singleton c + str s = DList.append s + + go !prec !expr = case expr of + VarE !Var + KstE !Integer + NegE Expr + LinE [(Sign,Expr)] + MulE [Expr] + PowE Expr !Int + +-} + +-------------------------------------------------------------------------------- + diff --git a/src/Math/Singular/Factory/GFTables.hs b/src/Math/Singular/Factory/GFTables.hs new file mode 100644 --- /dev/null +++ b/src/Math/Singular/Factory/GFTables.hs @@ -0,0 +1,145 @@ + +-- | Guessing the location of factory\'s gftables. +-- +-- For the Galois field domains to work (for non-primer order), you have to +-- first initialize singular-factory by setting this directory. +-- +-- This can be done either by `initGFTables` or `initGFTables'`. +-- + +{-# LANGUAGE CPP, BangPatterns, ForeignFunctionInterface #-} +module Math.Singular.Factory.GFTables where + +------------------------------------------------------------------------------------------ + +import Control.Monad + +import Foreign.Ptr +import Foreign.C +import Foreign.C.String + +import System.FilePath +import System.Directory +import System.Process + +import Data.IORef +import System.IO.Unsafe as Unsafe + +-------------------------------------------------------------------------------- +-- * A global variable + +{-# NOINLINE theGFTablesDir #-} +theGFTablesDir :: IORef (Maybe FilePath) +theGFTablesDir = Unsafe.unsafePerformIO $ newIORef Nothing + +getGFTablesDir :: IO (Maybe FilePath) +getGFTablesDir = readIORef theGFTablesDir + +-------------------------------------------------------------------------------- +-- * Initialization + +-- | We try to guess the location. +initGFTables :: IO () +initGFTables = initGFTables' Nothing + +-- | Set the location of the small finite field table files. +-- +-- If you know where they are located, please set it. +-- If you don't know, we try to guess it, but I have no idea how +-- to figure this out in general (pkg-config does not seem to have this information...) +-- +initGFTables' :: Maybe FilePath -> IO () +initGFTables' mbdir = case mbdir of + Just fpath -> setGFTablesDir fpath + Nothing -> guessGFTablesDir >>= \d -> case d of + Just fpath -> do + -- putStrLn $ "gftables dir = " ++ (fpath "gftables") + setGFTablesDir fpath + Nothing -> do + writeIORef theGFTablesDir Nothing + putStrLn "WARNING: cannot find factory's gftables" + +setGFTablesDir :: FilePath -> IO () +setGFTablesDir fpath0 = do + fpath1 <- canonicalizePath fpath0 + withCString (fpath1 ++ "/") $ \ptr -> set_gftable_dir ptr + writeIORef theGFTablesDir (Just fpath1) + +-- | Apparently we need to manually find the directory containing the GF tables... +-- +-- On my debian install it is at @\/usr\/share\/singular\/factory\/gftables\/@, but how +-- to figure that out??? +foreign import ccall "set_gftable_dir" set_gftable_dir :: Ptr CChar -> IO () + +------------------------------------------------------------------------------------------ +-- * Guessing + +guessGFTablesDir :: IO (Maybe FilePath) +guessGFTablesDir = do + +#if defined(linux_HOST_OS) + guessLinux +#elif defined(darwin_HOST_OS) + guessHomebrew +#elif defined(mingw32_HOST_OS) || defined(mingw64_HOST_OS) + guessCygwin +#else + return Nothing +#endif + +------------------------------------------------------------------------------------------ + +infixr 5 >>> + +(>>>) :: IO (Maybe a) -> IO (Maybe a) -> IO (Maybe a) +(>>>) action1 action2 = do + mb <- action1 + case mb of + Just x -> return (Just x) + Nothing -> action2 + +testDir :: FilePath -> IO (Maybe FilePath) +testDir dir = doesFileExist (dir "gftables/361") >>= \b -> if b + then return (Just dir) + else return Nothing + +------------------------------------------------------------------------------------------ + +-- macOS w/ Homebrew +-- on my machine it looks like this: /usr/local/Cellar/singular/4.1.2p1_2/share/factory +guessHomebrew :: IO (Maybe FilePath) +guessHomebrew = do + let brew_prefix = "/usr/local" -- TODO: maybe use "brew config" (but it's rather slow) + let cellar = brew_prefix "Cellar" + let sing_root = cellar "singular" + entries <- map (sing_root ) <$> listDirectory sing_root + subdirs <- filterM doesDirectoryExist entries + -- print entries + -- print subdirs + foldl (>>>) (return Nothing) [ testDir (sing_root d "share/factory") | d <- subdirs] + +-- generic Linux +-- I guess typically it looks like this: +guessLinux :: IO (Maybe FilePath) +guessLinux = + testDir "/usr/share/singular/factory" >>> + testDir "/usr/share/singular" >>> + testDir "/usr/share/factory" >>> + testDir "/usr/local/share/singular/factory" >>> + testDir "/usr/local/share/singular" >>> + testDir "/usr/local/share/factory" + +-- cygwin +guessCygwin ::IO (Maybe FilePath) +guessCygwin = do + cygwin_root <- readCreateProcess (shell "cygpath -w /") "" + let test1 dir = testDir (cygwin_root dir) + id <$> + test1 "/usr/share/singular/factory" >>> + test1 "/usr/share/singular" >>> + test1 "/usr/share/factory" >>> + test1 "/usr/local/share/singular/factory" >>> + test1 "/usr/local/share/singular" >>> + test1 "/usr/local/share/factory" + +------------------------------------------------------------------------------------------ diff --git a/src/Math/Singular/Factory/Internal/CanonicalForm.hs b/src/Math/Singular/Factory/Internal/CanonicalForm.hs new file mode 100644 --- /dev/null +++ b/src/Math/Singular/Factory/Internal/CanonicalForm.hs @@ -0,0 +1,511 @@ + +-- | Low-level bindings to singular-factory + +{-# LANGUAGE BangPatterns, EmptyDataDecls, ForeignFunctionInterface #-} +module Math.Singular.Factory.Internal.CanonicalForm where + +-------------------------------------------------------------------------------- + +import Data.Word +import Data.Ratio + +import Foreign.C +import Foreign.Ptr +import Foreign.ForeignPtr +import Foreign.Marshal.Array + +import Numeric.GMP.Types +import qualified Numeric.GMP.Utils as GMP + +-------------------------------------------------------------------------------- +-- * types + +data Variable +data CanonicalForm +data Factor -- CFFactor +data FactorList -- List + +type Var = ForeignPtr Variable +type CF = ForeignPtr CanonicalForm +type Fac = ForeignPtr Factor +type FacList = ForeignPtr FactorList + +-------------------------------------------------------------------------------- +-- * constants + +foreign import ccall "level_trans" c_level_trans :: CInt +foreign import ccall "level_base" c_level_base :: CInt +foreign import ccall "level_quot" c_level_quot :: CInt +foreign import ccall "level_expr" c_level_expr :: CInt + +-------------------------------------------------------------------------------- +-- * versions + +foreign import ccall "get_factory_version" c_get_factory_version :: Ptr CChar -> Int -> IO () +foreign import ccall "get_package_version" c_get_package_version :: Ptr CChar -> Int -> IO () + +getFactoryVersion :: IO String +getFactoryVersion = do + let n = 16 + allocaArray n $ \ptr -> do + c_get_factory_version ptr n + peekCString ptr + +getPackageVersion :: IO String +getPackageVersion = do + let n = 16 + allocaArray n $ \ptr -> do + c_get_package_version ptr n + peekCString ptr + +foreign import ccall "have_FLINT" c_have_FLINT :: IO CInt +foreign import ccall "have_NTL" c_have_NTL :: IO CInt +foreign import ccall "have_GMP" c_have_GMP :: IO CInt + +-------------------------------------------------------------------------------- +-- * config + +foreign import ccall "set_default_switches" c_set_default_switches :: IO () + +setDefaultSwitches :: IO () +setDefaultSwitches = do + c_set_default_switches + +-------------------------------------------------------------------------------- +-- * memory management + +foreign import ccall "&free_var" + varFinalizerPtr :: FunPtr (Ptr Variable -> IO ()) + +foreign import ccall "&free_cf" + cfFinalizerPtr :: FunPtr (Ptr CanonicalForm -> IO ()) + +foreign import ccall "&free_fac" + facFinalizerPtr :: FunPtr (Ptr Factor -> IO ()) + +foreign import ccall "&free_faclist" + faclistFinalizerPtr :: FunPtr (Ptr FactorList -> IO ()) + +makeVar :: Ptr Variable -> IO Var +makeVar = newForeignPtr varFinalizerPtr + +makeCF :: Ptr CanonicalForm -> IO CF +makeCF = newForeignPtr cfFinalizerPtr + +makeFac :: Ptr Factor -> IO Fac +makeFac = newForeignPtr facFinalizerPtr + +makeFacList :: Ptr FactorList -> IO FacList +makeFacList = newForeignPtr faclistFinalizerPtr + +-------------------------------------------------------------------------------- +-- * variables + +foreign import ccall "new_var_level" c_new_var_level :: CInt -> IO (Ptr Variable) +foreign import ccall "new_var_name" c_new_var_name :: CChar -> IO (Ptr Variable) +foreign import ccall "new_var_level_name" c_new_var_level_name :: CInt -> CChar -> IO (Ptr Variable) + +newVarL :: Int -> IO Var +newVarL level = makeVar =<< c_new_var_level (fromIntegral level) + +newVarN :: Char -> IO Var +newVarN name = makeVar =<< c_new_var_name (castCharToCChar name) + +newVarLN :: Int -> Char -> IO Var +newVarLN level name = makeVar =<< c_new_var_level_name (fromIntegral level) (castCharToCChar name) + +newTransVar :: IO Var +newTransVar = makeVar =<< c_new_var_level (c_level_trans) + +foreign import ccall "root_of" c_root_of :: Ptr CanonicalForm -> IO (Ptr Variable) + +newRootOf :: CF -> IO Var +newRootOf cf = withForeignPtr cf $ \ptr -> makeVar =<< c_root_of ptr + +foreign import ccall "get_var_level" c_get_var_level :: Ptr Variable -> IO CInt +foreign import ccall "get_var_name" c_get_var_name :: Ptr Variable -> IO CChar + +getVarLevel :: Var -> IO Int +getVarLevel var = withForeignPtr var $ \ptr -> fromIntegral <$> (c_get_var_level ptr) + +getVarName :: Var -> IO Char +getVarName var = withForeignPtr var $ \ptr -> castCCharToChar <$> (c_get_var_name ptr) + +foreign import ccall "has_mipo" c_has_mipo :: Ptr Variable -> IO CInt +foreign import ccall "get_mipo" c_get_mipo :: Ptr Variable -> Ptr Variable -> IO (Ptr CanonicalForm) +foreign import ccall "set_mipo" c_set_mipo :: Ptr Variable -> Ptr CanonicalForm -> IO () +foreign import ccall "set_reduce" c_set_reduce :: Ptr Variable -> CInt -> IO () + +hasMinimalPoly :: Var -> IO Bool +hasMinimalPoly var = withForeignPtr var $ \ptr -> liftBool (c_has_mipo ptr) + +getMinimalPoly :: Var -> Var -> IO CF +getMinimalPoly var1 var2 = + withForeignPtr var1 $ \ptr1 -> + withForeignPtr var2 $ \ptr2 -> + makeCF =<< c_get_mipo ptr1 ptr2 + +setReduceFlag :: Var -> Bool -> IO () +setReduceFlag var flag = withForeignPtr var $ \ptr -> c_set_reduce ptr (if flag then 1 else 0) + +-------------------------------------------------------------------------------- +-- * factors + +foreign import ccall "get_factor" c_get_factor :: Ptr Factor -> IO (Ptr CanonicalForm) + +getFactorCF :: Fac -> IO CF +getFactorCF fac = withForeignPtr fac $ \ptr -> makeCF =<< c_get_factor ptr + +foreign import ccall "get_fac_expo" c_get_fac_expo :: Ptr Factor -> IO CInt + +getFactorExpo :: Fac -> IO Int +getFactorExpo fac = withForeignPtr fac $ \ptr -> fromIntegral <$> c_get_fac_expo ptr + +unpackFactor :: Fac -> IO (CF,Int) +unpackFactor fac = do + cf <- getFactorCF fac + expo <- getFactorExpo fac + return (cf,expo) + +-------------------------------------------------------------------------------- +-- * lists + +foreign import ccall "get_list_length" c_get_list_length :: Ptr FactorList -> IO CInt + +getListLength :: FacList -> IO Int +getListLength faclist = withForeignPtr faclist $ \ptr -> fromIntegral <$> c_get_list_length ptr + +-------------------------------------------------------------------------------- +-- * factorization + +foreign import ccall "flatten_faclist" c_flatten_faclist :: Ptr FactorList -> Ptr (Ptr Factor) -> IO () + +flattenFactorList :: FacList -> IO [Fac] +flattenFactorList faclist = withForeignPtr faclist $ \listptr -> do + n <- fromIntegral <$> c_get_list_length listptr + allocaArray n $ \arrptr -> do + c_flatten_faclist listptr arrptr + facptrs <- peekArray n arrptr + mapM makeFac facptrs + +unpackFactorList :: FacList -> IO [(CF,Int)] +unpackFactorList faclist = mapM unpackFactor =<< flattenFactorList faclist + +foreign import ccall "hs_factorize" c_hs_factorize :: Ptr CanonicalForm -> IO (Ptr FactorList) + +factorizeIO' :: CF -> IO FacList +factorizeIO' cf = withForeignPtr cf $ \ptr -> makeFacList =<< (c_hs_factorize ptr) + +factorizeIO :: CF -> IO [(CF,Int)] +factorizeIO cf = unpackFactorList =<< factorizeIO' cf + +{- +foreign import ccall "rat_factorize" c_rat_factorize :: Ptr CanonicalForm -> Ptr Variable -> CInt -> IO (Ptr FactorList) +foreign import ccall "fq_factorize" c_fq_factorize :: Ptr CanonicalForm -> Ptr Variable -> CInt -> IO (Ptr FactorList) +foreign import ccall "fp_factorize" c_fp_factorize :: Ptr CanonicalForm -> CInt -> IO (Ptr FactorList) +foreign import ccall "gf_factorize" c_gf_factorize :: Ptr CanonicalForm -> CInt -> IO (Ptr FactorList) + +fpFactorizeIO' :: CF -> Bool -> IO FacList +fpFactorizeIO' cf subst_check_flag = + withForeignPtr cf $ \ptr -> + makeFacList =<< (c_fp_factorize ptr (bool2cint subst_check_flag)) + +gfFactorizeIO' :: CF -> Bool -> IO FacList +gfFactorizeIO' cf subst_check_flag = + withForeignPtr cf $ \ptr -> + makeFacList =<< (c_gf_factorize ptr (bool2cint subst_check_flag)) + +ratFactorizeIO' :: CF -> Var -> Bool -> IO FacList +ratFactorizeIO' cf var subst_check_flag = + withForeignPtr cf $ \ptr -> withForeignPtr var $ \varptr -> + makeFacList =<< (c_rat_factorize ptr varptr (bool2cint subst_check_flag)) + +fqFactorizeIO' :: CF -> Var -> Bool -> IO FacList +fqFactorizeIO' cf var subst_check_flag = + withForeignPtr cf $ \ptr -> withForeignPtr var $ \varptr -> + makeFacList =<< (c_fq_factorize ptr varptr (bool2cint subst_check_flag)) + +fpFactorizeIO :: CF -> Bool -> IO [(CF,Int)] +fpFactorizeIO cf flag = unpackFactorList =<< fpFactorizeIO' cf flag + +gfFactorizeIO :: CF -> Bool -> IO [(CF,Int)] +gfFactorizeIO cf flag = unpackFactorList =<< gfFactorizeIO' cf flag + +ratFactorizeIO :: CF -> Var -> Bool -> IO [(CF,Int)] +ratFactorizeIO cf var flag = unpackFactorList =<< ratFactorizeIO' cf var flag + +fqFactorizeIO :: CF -> Var -> Bool -> IO [(CF,Int)] +fqFactorizeIO cf var flag = unpackFactorList =<< fqFactorizeIO' cf var flag +-} + +-------------------------------------------------------------------------------- +-- * basic CFs + +foreign import ccall "empty_cf" c_empty_cf :: IO (Ptr CanonicalForm) + +newEmptyCF :: IO CF +newEmptyCF = makeCF =<< c_empty_cf + +foreign import ccall "const_cf" c_const_cf :: CInt -> IO (Ptr CanonicalForm) + +newSmallConstCF :: Int -> IO CF +newSmallConstCF k = makeCF =<< c_const_cf (fromIntegral k) + +foreign import ccall "var_cf" c_var_cf :: Ptr Variable -> IO (Ptr CanonicalForm) +foreign import ccall "var_pow_cf" c_var_pow_cf :: Ptr Variable -> Int -> IO (Ptr CanonicalForm) + +varIO :: Var -> IO CF +varIO var = withForeignPtr var $ \vptr -> makeCF =<< c_var_cf vptr + +varPowIO :: Var -> Int -> IO CF +varPowIO var expo = withForeignPtr var $ \vptr -> makeCF =<< c_var_pow_cf vptr expo + +-------------------------------------------------------------------------------- +-- * basic CF predicates + +bool2cint :: Bool -> CInt +bool2cint True = 1 +bool2cint False = 0 + +cint2bool :: CInt -> Bool +cint2bool k = (k/=0) + +liftBool :: IO CInt -> IO Bool +liftBool action = action >>= \k -> return (cint2bool k) + +---------------------------------------- + +foreign import ccall "is_zero" c_is_zero :: Ptr CanonicalForm -> IO CInt +foreign import ccall "is_one" c_is_one :: Ptr CanonicalForm -> IO CInt + +isZeroIO :: CF -> IO Bool +isZeroIO cf = withForeignPtr cf $ \ptr -> liftBool (c_is_zero ptr) + +isOneIO :: CF -> IO Bool +isOneIO cf = withForeignPtr cf $ \ptr -> liftBool (c_is_one ptr) + +foreign import ccall "is_imm" c_is_imm :: Ptr CanonicalForm -> IO CInt +foreign import ccall "is_univariate" c_is_univariate :: Ptr CanonicalForm -> IO CInt + +isImmediateIO :: CF -> IO Bool +isImmediateIO cf = withForeignPtr cf $ \ptr -> liftBool (c_is_imm ptr) + +isUnivariateIO :: CF -> IO Bool +isUnivariateIO cf = withForeignPtr cf $ \ptr -> liftBool (c_is_univariate ptr) + +foreign import ccall "in_ZZ" c_in_ZZ :: Ptr CanonicalForm -> IO CInt +foreign import ccall "in_QQ" c_in_QQ :: Ptr CanonicalForm -> IO CInt +foreign import ccall "in_GF" c_in_GF :: Ptr CanonicalForm -> IO CInt +foreign import ccall "in_FF" c_in_FF :: Ptr CanonicalForm -> IO CInt + +isInZZ_IO :: CF -> IO Bool +isInZZ_IO cf = withForeignPtr cf $ \ptr -> liftBool (c_in_ZZ ptr) + +isInQQ_IO :: CF -> IO Bool +isInQQ_IO cf = withForeignPtr cf $ \ptr -> liftBool (c_in_QQ ptr) + +isInGF_IO :: CF -> IO Bool +isInGF_IO cf = withForeignPtr cf $ \ptr -> liftBool (c_in_GF ptr) + +isInFF_IO :: CF -> IO Bool +isInFF_IO cf = withForeignPtr cf $ \ptr -> liftBool (c_in_FF ptr) + +foreign import ccall "in_BaseDomain" c_in_BaseDomain :: Ptr CanonicalForm -> IO CInt +foreign import ccall "in_CoeffDomain" c_in_CoeffDomain :: Ptr CanonicalForm -> IO CInt +foreign import ccall "in_PolyDomain" c_in_PolyDomain :: Ptr CanonicalForm -> IO CInt +foreign import ccall "in_Extension" c_in_Extension :: Ptr CanonicalForm -> IO CInt +foreign import ccall "in_QuotDomain" c_in_QuotDomain :: Ptr CanonicalForm -> IO CInt + +isInBaseDomainIO :: CF -> IO Bool +isInBaseDomainIO cf = withForeignPtr cf $ \ptr -> liftBool (c_in_BaseDomain ptr) + +isInCoeffDomainIO :: CF -> IO Bool +isInCoeffDomainIO cf = withForeignPtr cf $ \ptr -> liftBool (c_in_CoeffDomain ptr) + +isInPolyDomainIO :: CF -> IO Bool +isInPolyDomainIO cf = withForeignPtr cf $ \ptr -> liftBool (c_in_PolyDomain ptr) + +isInExtensionIO :: CF -> IO Bool +isInExtensionIO cf = withForeignPtr cf $ \ptr -> liftBool (c_in_Extension ptr) + +isInQuotDomainIO :: CF -> IO Bool +isInQuotDomainIO cf = withForeignPtr cf $ \ptr -> liftBool (c_in_QuotDomain ptr) + +-------------------------------------------------------------------------------- +-- * basic properties + +foreign import ccall "degree_of" c_degree_of :: Ptr CanonicalForm -> IO CInt +foreign import ccall "level_of" c_level_of :: Ptr CanonicalForm -> IO CInt +foreign import ccall "mvar_of" c_mvar_of :: Ptr CanonicalForm -> IO (Ptr Variable) + +getDegree :: CF -> IO Int +getDegree cf = withForeignPtr cf $ \ptr -> fromIntegral <$> (c_degree_of ptr) + +getLevel :: CF -> IO Int +getLevel cf = withForeignPtr cf $ \ptr -> fromIntegral <$> (c_level_of ptr) + +getMainVar :: CF -> IO Var +getMainVar cf = withForeignPtr cf $ \ptr -> makeVar =<< (c_mvar_of ptr) + +-------------------------------------------------------------------------------- +-- * small values + +foreign import ccall "smallint_value" c_smallint_value :: Ptr CanonicalForm -> IO CLong -- !!! + +getSmallIntValue :: CF -> IO Int +getSmallIntValue cf = withForeignPtr cf $ \ptr -> fromIntegral <$> (c_smallint_value ptr) + +foreign import ccall "numer" c_numer :: Ptr CanonicalForm -> IO (Ptr CanonicalForm) +foreign import ccall "denom" c_denom :: Ptr CanonicalForm -> IO (Ptr CanonicalForm) + +getNumer :: CF -> IO CF +getNumer cf = withForeignPtr cf $ \ptr -> makeCF =<< (c_numer ptr) + +getDenom :: CF -> IO CF +getDenom cf = withForeignPtr cf $ \ptr -> makeCF =<< (c_denom ptr) + +foreign import ccall "index_poly" c_index_poly :: Ptr CanonicalForm -> Int -> IO (Ptr CanonicalForm) +foreign import ccall "map_into" c_map_into :: Ptr CanonicalForm -> IO (Ptr CanonicalForm) +foreign import ccall "substitute" c_substitute :: Ptr CanonicalForm -> Ptr Variable -> Ptr CanonicalForm -> IO (Ptr CanonicalForm) + +-------------------------------------------------------------------------------- +-- * Polynomial operations + +-- | Get the given degree part +getCfAtIndex :: CF -> Int -> IO CF +getCfAtIndex cf idx = withForeignPtr cf $ \ptr -> makeCF =<< (c_index_poly ptr $ fromIntegral idx) + +-- | Map into the current base domain +mapIntoIO :: CF -> IO CF +mapIntoIO cf = withForeignPtr cf $ \ptr -> makeCF =<< (c_map_into ptr) + +-- | Map into the current base domain +substituteIO :: Var -> CF -> (CF -> IO CF) +substituteIO var what cf = + withForeignPtr var $ \pvar -> + withForeignPtr what $ \pwhat -> + withForeignPtr cf $ \pcf -> + makeCF =<< (c_substitute pcf pvar pwhat) + +-------------------------------------------------------------------------------- +-- * Binary operations + +foreign import ccall "is_equal" c_is_equal :: Ptr CanonicalForm -> Ptr CanonicalForm -> IO CInt + +isEqualIO :: CF -> CF -> IO Bool +isEqualIO x y = withForeignPtr x $ \p -> withForeignPtr y $ \q -> liftBool (c_is_equal p q) + +foreign import ccall "plus_cf" c_plus_cf :: Ptr CanonicalForm -> Ptr CanonicalForm -> IO (Ptr CanonicalForm) +foreign import ccall "minus_cf" c_minus_cf :: Ptr CanonicalForm -> Ptr CanonicalForm -> IO (Ptr CanonicalForm) +foreign import ccall "times_cf" c_times_cf :: Ptr CanonicalForm -> Ptr CanonicalForm -> IO (Ptr CanonicalForm) +foreign import ccall "pow_cf" c_pow_cf :: Ptr CanonicalForm -> CInt -> IO (Ptr CanonicalForm) +foreign import ccall "div_cf" c_div_cf :: Ptr CanonicalForm -> Ptr CanonicalForm -> IO (Ptr CanonicalForm) +foreign import ccall "mod_cf" c_mod_cf :: Ptr CanonicalForm -> Ptr CanonicalForm -> IO (Ptr CanonicalForm) + +foreign import ccall "gcd_poly_cf" c_gcd_poly_cf :: Ptr CanonicalForm -> Ptr CanonicalForm -> IO (Ptr CanonicalForm) +foreign import ccall "reduce_cf" c_reduce_cf :: Ptr CanonicalForm -> Ptr CanonicalForm -> IO (Ptr CanonicalForm) + +plusIO :: CF -> CF -> IO CF +plusIO x y = withForeignPtr x $ \p -> withForeignPtr y $ \q -> makeCF =<< (c_plus_cf p q) + +minusIO :: CF -> CF -> IO CF +minusIO x y = withForeignPtr x $ \p -> withForeignPtr y $ \q -> makeCF =<< (c_minus_cf p q) + +timesIO :: CF -> CF -> IO CF +timesIO x y = withForeignPtr x $ \p -> withForeignPtr y $ \q -> makeCF =<< (c_times_cf p q) + +powIO :: CF -> Int -> IO CF +powIO x n = withForeignPtr x $ \p -> makeCF =<< (c_pow_cf p $ fromIntegral n) + +divIO :: CF -> CF -> IO CF +divIO x y = withForeignPtr x $ \p -> withForeignPtr y $ \q -> makeCF =<< (c_div_cf p q) + +modIO :: CF -> CF -> IO CF +modIO x y = withForeignPtr x $ \p -> withForeignPtr y $ \q -> makeCF =<< (c_mod_cf p q) + +gcdPolyIO :: CF -> CF -> IO CF +gcdPolyIO x y = withForeignPtr x $ \p -> withForeignPtr y $ \q -> makeCF =<< (c_gcd_poly_cf p q) + +reduceIO :: CF -> CF -> IO CF +reduceIO x y = withForeignPtr x $ \p -> withForeignPtr y $ \q -> makeCF =<< (c_reduce_cf p q) + +-------------------------------------------------------------------------------- +-- * GMP compatibility layer + +foreign import ccall "get_gmp_numerator" c_get_gmp_numerator :: Ptr CanonicalForm -> Ptr MPZ -> IO () +foreign import ccall "get_gmp_denominator" c_get_gmp_denominator :: Ptr CanonicalForm -> Ptr MPZ -> IO () + +getGmpNumerator :: CF -> IO Integer +getGmpNumerator cf = + withForeignPtr cf $ \cfp -> + GMP.withOutInteger_ (\mpz_ptr -> c_get_gmp_numerator cfp mpz_ptr) + +getGmpDenominator :: CF -> IO Integer +getGmpDenominator cf = + withForeignPtr cf $ \cfp -> + GMP.withOutInteger_ (\mpz_ptr -> c_get_gmp_denominator cfp mpz_ptr) + +foreign import ccall "make_ZZ_from_gmp" c_make_ZZ_from_gmp :: Ptr MPZ -> IO (Ptr CanonicalForm) +foreign import ccall "make_QQ_from_gmp" c_make_QQ_from_gmp :: Ptr MPZ -> Ptr MPZ -> CInt -> IO (Ptr CanonicalForm) + +makeIntegerCF :: Integer -> IO CF +makeIntegerCF n = + GMP.withInInteger n $ \mpz_ptr -> + makeCF =<< (c_make_ZZ_from_gmp mpz_ptr) + +makeRationalCF :: Rational -> IO CF +makeRationalCF n = + GMP.withInInteger (numerator n) $ \mpz_ptr1 -> + GMP.withInInteger (denominator n) $ \mpz_ptr2 -> + makeCF =<< (c_make_QQ_from_gmp mpz_ptr1 mpz_ptr2 0) + +-------------------------------------------------------------------------------- +-- * Base domain characteristic + +foreign import ccall "get_characteristic" c_get_characteristic :: IO CInt +foreign import ccall "set_characteristic1" c_set_characteristic1 :: CInt -> IO () +-- foreign import ccall "set_characteristic2" c_set_characteristic2 :: CInt -> CInt -> IO () +foreign import ccall "set_characteristic3" c_set_characteristic3 :: CInt -> CInt -> CChar -> IO () + +getCharacteristic :: IO Int +getCharacteristic = fromIntegral <$> c_get_characteristic + +-- | prime fields and QQ +setCharacteristic1 :: Int -> IO () +setCharacteristic1 = c_set_characteristic1 . fromIntegral + +{- +-- | prime power fields (meaning FF?) +-- apparently they removed this ???? +setCharacteristic2 :: Int -> Int -> IO () +setCharacteristic2 p n = c_set_characteristic2 (fromIntegral p) (fromIntegral n) +-} + +-- | Galois fields +setCharacteristic3 :: Int -> Int -> Char -> IO () +setCharacteristic3 p n ch = c_set_characteristic3 (fromIntegral p) (fromIntegral n) (castCharToCChar ch) + +foreign import ccall "get_gf_value" c_get_gf_value :: Ptr CanonicalForm -> IO CInt +foreign import ccall "is_FF_in_GF" c_is_FF_in_GF :: Ptr CanonicalForm -> IO CInt +foreign import ccall "get_GF_degree" c_get_GF_degree :: IO CInt +foreign import ccall "get_GF_generator" c_get_GF_generator :: IO (Ptr CanonicalForm) + +-- | This returns the exponent of the canonical generator. +-- If the input is zero, it appears to return the order of the field, q, but don't rely on this... +getGFValue :: CF -> IO Int +getGFValue cf = withForeignPtr cf $ \ptr -> fromIntegral <$> (c_get_gf_value ptr) + +-- | True if element of the prime subfield +isFFinGF_IO :: CF -> IO Bool +isFFinGF_IO cf = withForeignPtr cf $ \ptr -> liftBool (c_is_FF_in_GF ptr) + +-- | Returns the degree of the Galois field (degree of extension over the prime field) +getGFDegree :: IO Int +getGFDegree = fromIntegral <$> c_get_GF_degree + +-- | Returns the generator of the Galois field +getGFGenerator :: IO CF +getGFGenerator = makeCF =<< c_get_GF_generator + +-------------------------------------------------------------------------------- + diff --git a/src/Math/Singular/Factory/Internal/DList.hs b/src/Math/Singular/Factory/Internal/DList.hs new file mode 100644 --- /dev/null +++ b/src/Math/Singular/Factory/Internal/DList.hs @@ -0,0 +1,29 @@ + +-- | Difference lists, minimalistic edition + +module Math.Singular.Factory.Internal.DList where + +-------------------------------------------------------------------------------- + +import Data.List ( foldr ) + +-------------------------------------------------------------------------------- + +type DList a = [a] -> [a] + +empty :: DList a +empty = id + +singleton :: a -> DList a +singleton = (:) + +toList :: DList a -> [a] +toList dl = dl [] + +append :: DList a -> DList a -> DList a +append = (.) + +concat :: [DList a] -> DList a +concat = foldr append empty + +-------------------------------------------------------------------------------- diff --git a/src/Math/Singular/Factory/Internal/Factory.hs b/src/Math/Singular/Factory/Internal/Factory.hs new file mode 100644 --- /dev/null +++ b/src/Math/Singular/Factory/Internal/Factory.hs @@ -0,0 +1,433 @@ + +-- | Medium-level bindings to singular-factory + +{-# LANGUAGE BangPatterns, DataKinds, TypeSynonymInstances, FlexibleInstances #-} +module Math.Singular.Factory.Internal.Factory where + +-------------------------------------------------------------------------------- + +import Data.List +import Data.Maybe +import Data.Ratio +import Data.Char +import Data.Word + +import Control.Monad + +import GHC.TypeLits +import Data.Proxy + +import System.IO.Unsafe as Unsafe + +import Numeric.GMP.Types +import qualified Numeric.GMP.Utils as GMP + +import Math.Singular.Factory.Internal.CanonicalForm +import Math.Singular.Factory.Internal.DList as DList + +-------------------------------------------------------------------------------- + +-- | The maximum prime characteristic Factory can handle +-- (note: for prime power fields the limit is much smaller) +maxCharacteristic :: Int +maxCharacteristic = 536870909 -- 2^29-3 + +-------------------------------------------------------------------------------- +-- * Version and configuration info + +factoryVersion :: String +factoryVersion = Unsafe.unsafePerformIO getFactoryVersion + +packageVersion :: String +packageVersion = Unsafe.unsafePerformIO getPackageVersion + +haveFLINT, haveNTL, haveGMP :: Bool +haveFLINT = Unsafe.unsafePerformIO (cint2bool <$> c_have_FLINT) +haveNTL = Unsafe.unsafePerformIO (cint2bool <$> c_have_NTL ) +haveGMP = Unsafe.unsafePerformIO (cint2bool <$> c_have_GMP ) + +-------------------------------------------------------------------------------- +-- * Basic operations and instances + +-- | Because the native equality comparison seems to be unreliable, +-- we provide an alternative implementation which subtracts the two +-- polynomials and then tests for the result being zero... +safeEqCF :: CF -> CF -> Bool +safeEqCF x y = isZeroCF (x - y) + +-- | Note: this does not seem to be reliable in practice... +-- Better subtract them and then use isZeroCF... +nativeEqCF :: CF -> CF -> Bool +nativeEqCF x y = Unsafe.unsafePerformIO (isEqualIO x y) + +{- there is already an instance... -} +-- instance Eq CF where (==) = eqCF + +isZeroCF :: CF -> Bool +isZeroCF cf = Unsafe.unsafePerformIO (isZeroIO cf) + +isOneCF :: CF -> Bool +isOneCF cf = Unsafe.unsafePerformIO (isOneIO cf) + +varCF :: Var -> CF +varCF var = Unsafe.unsafePerformIO (varIO var) + +varPowCF :: Var -> Int -> CF +varPowCF var expo = Unsafe.unsafePerformIO (varPowIO var expo) + +instance Num CF where + fromInteger n = Unsafe.unsafePerformIO $ newSmallConstCF $ fromIntegral n -- BIG INTS ARE NOT HANDLED!!! + (+) x y = Unsafe.unsafePerformIO (plusIO x y) + (-) x y = Unsafe.unsafePerformIO (minusIO x y) + (*) x y = Unsafe.unsafePerformIO (timesIO x y) + abs = error "CF: Num/abs is not implemented" + signum = error "CF: Num/signum is not implemented" + +powCF :: CF -> Int -> CF +powCF x n = Unsafe.unsafePerformIO (powIO x n) + +modCF :: CF -> CF -> CF +modCF x y = Unsafe.unsafePerformIO (modIO x y) + +divCF :: CF -> CF -> CF +divCF x y = Unsafe.unsafePerformIO (divIO x y) + +substituteCF :: Var -> CF -> (CF -> CF) +substituteCF var what cf = Unsafe.unsafePerformIO (substituteIO var what cf) + +gcdPolyCF :: CF -> CF -> CF +gcdPolyCF x y = Unsafe.unsafePerformIO (gcdPolyIO x y) + +reduceCF :: CF -> CF -> CF +reduceCF x y = Unsafe.unsafePerformIO (reduceIO x y) + +factorizeCF :: CF -> [(CF,Int)] +factorizeCF x = Unsafe.unsafePerformIO (factorizeIO x) + +{- +-- | Factorization over prime fields. The flag is whether to perform substitution check +fpFactorizeCF :: CF -> Bool -> [(CF,Int)] +fpFactorizeCF x flag = Unsafe.unsafePerformIO (fpFactorizeIO x flag) + +-- | Factorization over GF (prime power) fields. The flag is whether to perform substitution check +gfFactorizeCF :: CF -> Bool -> [(CF,Int)] +gfFactorizeCF x flag = Unsafe.unsafePerformIO (gfFactorizeIO x flag) + +-- | Factorization over algebraic extensions of Q. The flag is whether to perform substitution check +ratFactorizeCF :: CF -> Var -> Bool -> [(CF,Int)] +ratFactorizeCF x alpha flag = Unsafe.unsafePerformIO (ratFactorizeIO x alpha flag) + +-- | Factorization over algebraic extensions of a prime field. The flag is whether to perform substitution check +fqFactorizeCF :: CF -> Var -> Bool -> [(CF,Int)] +fqFactorizeCF x alpha flag = Unsafe.unsafePerformIO (fqFactorizeIO x alpha flag) +-} + +-------------------------------------------------------------------------------- +-- * pretty printing + +showCF :: CF -> String +showCF x = Unsafe.unsafePerformIO (showIO x) + +showCF_with :: (Int -> String) -> CF -> String +showCF_with showVar x = Unsafe.unsafePerformIO (showIO_with showVar x) + +showCF_dense :: CF -> String +showCF_dense x = Unsafe.unsafePerformIO (showIO_dense x) + +printCF :: CF -> IO () +printCF cf = putStrLn =<< showIO cf + +showIO :: CF -> IO String +showIO cf = do + terms <- genericMarshalFromCF makeTerm cf + return $ intercalate " + " (map prettyTerm terms) + +showIO_with :: (Int -> String) -> CF -> IO String +showIO_with showVar cf = do + terms <- genericMarshalFromCF makeTerm cf + return $ intercalate " + " (map (prettyTermWith showVar) terms) + +showIO_dense :: CF -> IO String +showIO_dense cf = do + terms <- genericMarshalFromCF makeTerm cf + return $ intercalate "+" (map prettyTerm terms) + +-------------------------------------------------------------------------------- +-- * Galois fields + +-- | Values in a Galois field. +-- +-- They can either 0, or a power of the canonical generator (of the multiplicative +-- group, which is cyclic). +-- +-- Furthermore, they can be possibly also element of the prime subfield. +-- +data GFValue + = GFZero + | GFSubField { _gfGenExpo :: !Int , _gfModP :: !Int } + | GFGenPow { _gfGenExpo :: !Int } + deriving (Eq,Ord) + +instance Show GFValue where + show = showGFValue1 "#" + +-- | Elements of the prime subfield are shown as numbers, the rest as +-- powers of the generator +showGFValue1 :: String -> GFValue -> String +showGFValue1 gen gfv = case gfv of + GFZero -> "0" + GFSubField expo modp -> show modp + GFGenPow expo -> if expo == 1 then gen else gen ++ "^" ++ show expo + +-- | Elements of the prime subfield are also shown as powers of the generator +showGFValue2 :: String -> GFValue -> String +showGFValue2 gen gfv = case gfv of + GFZero -> "0" + GFSubField expo modp -> if expo == 1 then gen else gen ++ "^" ++ show expo + GFGenPow expo -> if expo == 1 then gen else gen ++ "^" ++ show expo + +-------------------------------------------------------------------------------- +-- * Base domains + +data BaseValue + = ZZ !Integer + | QQ !Rational + | FF !Int + | GF !GFValue + deriving (Eq,Ord,Show) + +baseValueIsOne :: BaseValue -> Bool +baseValueIsOne val = case val of + ZZ n -> n == 1 + QQ q -> q == 1 + GF gfv -> case gfv of { GFSubField _ 1 -> True ; _ -> False } + FF k -> k == 1 + +prettyBaseValue :: BaseValue -> String +prettyBaseValue val = case val of + ZZ n -> show n + QQ q -> case denominator q of + 1 -> show (numerator q) + b -> show (numerator q) ++ "/" ++ show b + GF k -> show k -- "<" ++ show k ++ ">" + FF k -> show k -- "<" ++ show k ++ ">" + +-------------------------------------------------------------------------------- +-- * Monomials + +type Level = Int +type Expo = Int + +newtype Monom + = Monom [(Level,Expo)] + deriving (Eq,Ord,Show) + +monomIsNull :: Monom -> Bool +monomIsNull (Monom list) = null list + +prettyMonom :: Monom -> String +prettyMonom = prettyMonomWith g where + g level = [chr (96 + level)] + +prettyMonomWith :: (Int -> String) -> Monom -> String +prettyMonomWith showVar (Monom [] ) = "1" -- "(1)" +prettyMonomWith showVar (Monom ves) = intercalate "*" (map f ves) where + f (v,0) = "1" + f (v,1) = showVar v + f (v,k) = showVar v ++ "^" ++ show k + +---------------------------------------- +-- * coefficients + +type Coeff = BaseValue -- TEMPORARY HACK + +coeffIsOne :: Coeff -> Bool +coeffIsOne = baseValueIsOne + +prettyCoeff = prettyBaseValue + +---------------------------------------- +-- * Terms + +data Term + = Term !Coeff !Monom + deriving (Eq,Ord,Show) + +prettyTerm :: Term -> String +prettyTerm (Term coeff monom) + | coeffIsOne coeff = prettyMonom monom + | monomIsNull monom = prettyCoeff coeff + | otherwise = prettyCoeff coeff ++ "*" ++ prettyMonom monom + +prettyTermWith :: (Int -> String) -> Term -> String +prettyTermWith showVar (Term coeff monom) + | coeffIsOne coeff = prettyMonomWith showVar monom + | monomIsNull monom = prettyCoeff coeff + | otherwise = prettyCoeff coeff ++ "*" ++ prettyMonomWith showVar monom + +-------------------------------------------------------------------------------- +-- * Marshalling from CF + +getZZ :: CF -> IO (Maybe Integer) +getZZ cf = isInZZ_IO cf >>= \b -> case b of + False -> return Nothing + True -> isImmediateIO cf >>= \b -> case b of + True -> (Just . fromIntegral) <$> getSmallIntValue cf + False -> Just <$> getGmpNumerator cf + -- error "getZZ: bignums are not implemented yet" + +getQQ :: CF -> IO (Maybe Rational) +getQQ cf = isInQQ_IO cf >>= \b -> case b of + False -> return Nothing + True -> do + -- Just numer <- (getZZ =<< getNumer cf) + -- Just denom <- (getZZ =<< getDenom cf) + numer <- getGmpNumerator cf + denom <- getGmpDenominator cf + return $ Just $ (numer % denom) + +getFF :: CF -> IO (Maybe Int) +getFF cf = isInFF_IO cf >>= \b -> case b of + False -> return Nothing + True -> isImmediateIO cf >>= \b -> case b of + True -> (Just . fromIntegral) <$> getSmallIntValue cf + False -> error "getFF: bignums are not implemented yet" + +getGF :: CF -> IO (Maybe GFValue) +getGF cf = isInGF_IO cf >>= \b -> case b of + False -> return Nothing + True -> isFFinGF_IO cf >>= \b -> case b of + True -> do + k <- getSmallIntValue cf + e <- getGFValue cf + return $ Just $ if (k==0) then GFZero else (GFSubField e k) + False -> do + e <- getGFValue cf + return $ Just $ GFGenPow e + +getBaseValue :: CF -> IO (Maybe BaseValue) +getBaseValue cf = isInBaseDomainIO cf >>= \b -> case b of + False -> return Nothing + True -> getBaseValueNotGF cf >>= \mb -> case mb of + Just val -> return $ Just val + Nothing -> getGF cf >>= \mb -> case mb of + Just gf -> return $ Just (GF gf) + Nothing -> return Nothing + +getBaseValueNotGF :: CF -> IO (Maybe BaseValue) +getBaseValueNotGF cf = isInBaseDomainIO cf >>= \b -> case b of + False -> return Nothing + True -> getZZ cf >>= \mb -> case mb of + Just n -> return $ Just (ZZ n) + Nothing -> getQQ cf >>= \mb -> case mb of + Just q -> return $ Just (QQ q) + Nothing -> getFF cf >>= \mb -> case mb of + Just k -> return $ Just (FF k) + Nothing -> return Nothing + +-------------------------------------------------------------------------------- + +isInBaseDomainCF :: CF -> Bool +isInBaseDomainCF cf = Unsafe.unsafePerformIO (isInBaseDomainIO cf) + +isInCoeffDomainCF :: CF -> Bool +isInCoeffDomainCF cf = Unsafe.unsafePerformIO (isInCoeffDomainIO cf) + +isInPolyDomainCF :: CF -> Bool +isInPolyDomainCF cf = Unsafe.unsafePerformIO (isInPolyDomainIO cf) + +isInQuotDomainCF :: CF -> Bool +isInQuotDomainCF cf = Unsafe.unsafePerformIO (isInQuotDomainIO cf) + +isInExtensionCF :: CF -> Bool +isInExtensionCF cf = Unsafe.unsafePerformIO (isInExtensionIO cf) + +-------------------------------------------------------------------------------- + +valueZZ :: CF -> Integer +valueZZ cf = case Unsafe.unsafePerformIO (getZZ cf) of + Just n -> n + Nothing -> error "valueZZ: not an integer" + +valueQQ :: CF -> Rational +valueQQ cf = case Unsafe.unsafePerformIO (getQQ cf) of + Just n -> n + Nothing -> error "valueQQ: not a rational" + +valueGF :: CF -> GFValue +valueGF cf = case Unsafe.unsafePerformIO (getGF cf) of + Just n -> n + Nothing -> error "valueGF: not a GF element" + +valueFF :: CF -> Int +valueFF cf = case Unsafe.unsafePerformIO (getFF cf) of + Just n -> n + Nothing -> error "valueFF: not a FF element" + +-------------------------------------------------------------------------------- + +makeMonom :: [(Level,Expo)] -> Monom +makeMonom ves = Monom $ filter cond ves where cond (lev,expo) = expo > 0 + +makeTerm :: [(Level,Expo)] -> BaseValue -> Term +makeTerm ves cf = Term cf (makeMonom ves) + +genericMarshalFromCF :: ([(Level,Expo)] -> BaseValue -> a) -> CF -> IO [a] +genericMarshalFromCF user cf = DList.toList <$> (genericMarshalFromCF_dlist user cf) + +-------------------------------------------------------------------------------- + +genericMarshalFromCF_dlist :: ([(Level,Expo)] -> BaseValue -> a) -> CF -> IO (DList a) +genericMarshalFromCF_dlist user = worker [] where + worker expos cf = do + getBaseValue cf >>= \mb -> case mb of + Just val -> return (DList.singleton $ user expos val) + Nothing -> do + level <- getLevel cf + deg <- getDegree cf + stuff <- forM [0..deg] $ \d -> do + this <- getCfAtIndex cf d + isZeroIO this >>= \b -> if b + then return (DList.empty) + else if d > 0 + then worker ((level,d):expos) this + else worker expos this + return $ DList.concat stuff + +{- +genericMarshalFromCF_list :: ([(Level,Expo)] -> BaseValue -> a) -> CF -> IO [a] +genericMarshalFromCF_list user = worker [] where + worker expos cf = do + getBaseValue cf >>= \mb -> case mb of + Just val -> return [user expos val] + Nothing -> do + level <- getLevel cf + deg <- getDegree cf + stuff <- forM [0..deg] $ \d -> do + this <- getCfAtIndex cf d + isZero this >>= \b -> if b + then return [] + else if d > 0 + then worker ((level,d):expos) this + else worker expos this + return $ Data.List.concat stuff +-} + +-------------------------------------------------------------------------------- + +marshalUnivariateFromCF :: CF -> IO [(CF,Int)] +marshalUnivariateFromCF = worker where + worker cf = do + deg <- getDegree cf + case deg of + 0 -> return [(cf,0)] + _ -> do + -- level <- getLevel cf + mbs <- forM [0..deg] $ \d -> do + this <- getCfAtIndex cf d + isZeroIO this >>= \b -> if b + then return $ Nothing + else return $ Just (this,d) + return $ catMaybes mbs + +-------------------------------------------------------------------------------- diff --git a/src/Math/Singular/Factory/Parser.hs b/src/Math/Singular/Factory/Parser.hs new file mode 100644 --- /dev/null +++ b/src/Math/Singular/Factory/Parser.hs @@ -0,0 +1,326 @@ + +-- | Parsing polynomials and polynomial expressions + +{-# LANGUAGE BangPatterns #-} +module Math.Singular.Factory.Parser + ( -- * Parsing polynomials + parseExpr , parseGenPoly + , parseStringExpr , parseStringGenPoly + -- * Parser monad + , Parser, runParser + , (<||>) , try + , many , many1 + -- * Parsers + , charP , charP_ , charsP + , spacesP_ , eofP + , signP , natP , integerP , identifierP + , exprP , genPolyP + ) + where + +-------------------------------------------------------------------------------- + +import Data.Char +import Data.List + +import Control.Applicative +import Control.Monad + +import Data.Traversable + +import Data.Proxy + +import Data.Text.Lazy ( Text ) +import qualified Data.Text.Lazy as T +import qualified Data.Text.Lazy.Read as T + +import Math.Singular.Factory.Expr +import Math.Singular.Factory.Variables + +import Math.Singular.Factory.Internal.DList as DList + +-------------------------------------------------------------------------------- + +type Var = String + +-------------------------------------------------------------------------------- +-- * the Parser monad + +newtype Parser a = P { runParser :: Text -> Either String (a,Text) } + +instance Functor Parser where + fmap f (P action) = P $ \text -> case action text of + Right (x, rest) -> Right (f x, rest) + Left msg -> Left msg + +instance Applicative Parser where + pure = return + (<*>) = ap + +instance Monad Parser where + return x = P $ \text -> Right (x,text) + (P action) >>= u = P $ \text -> case action text of + Right (x, rest) -> runParser (u x) rest + Left msg -> Left msg + +-------------------------------------------------------------------------------- + +try :: Parser a -> Parser (Maybe a) +try p = P $ \text -> case runParser p text of + Right (x, rest) -> Right (Just x , rest) + Left msg -> Right (Nothing, text) + +infixr 5 <||> + +(<||>) :: Parser a -> Parser a -> Parser a +(<||>) p q = do + mb <- try p + case mb of + Just y -> return y + Nothing -> q + +instance Alternative Parser where + (<|>) = (<||>) + +many1 :: Parser a -> Parser [a] +many1 p = do + x <- p + xs <- many p + return (x:xs) + +-------------------------------------------------------------------------------- + +charP_ :: (Char -> Bool) -> Parser () +charP_ cond = void $ charP cond + +charP :: (Char -> Bool) -> Parser Char +charP cond = P $ \text -> case T.uncons text of + Nothing -> Left "unexpected end of input" + Just (ch,rest) -> if cond ch + then Right (ch , rest) + else Left "unexpected character" + +charsP :: (Char -> Bool) -> Parser [Char] +charsP cond = P $ \text -> case T.span cond text of + (this,rest) -> Right (T.unpack this, rest) + +eofP :: Parser () +eofP = P $ \text -> if T.null text then Right ((),text) else Left "expected the end of input" + +spacesP_ :: Parser () +spacesP_ = P $ \text -> case T.span isSpace text of + (_,rest) -> Right ((),rest) + +withEof :: Parser a -> Parser a +withEof action = do + y <- action + spacesP_ + eofP + return y + +withSpaces :: Parser a -> Parser a +withSpaces action = do + y <- action + spacesP_ + return y + +-------------------------------------------------------------------------------- +-- * Parsing simple things + +_signP :: Parser Sign +_signP = charP (\c -> (c=='+' || c=='-')) >>= \ch -> return $ if (ch=='+') then Plus else Minus + +signP :: Parser Sign +signP = withSpaces _signP + +optionalSignP :: Parser Sign +optionalSignP = do + mb <- try signP + return $ case mb of + Just pm -> pm + Nothing -> Plus + +inParensP :: Parser a -> Parser a +inParensP action = do + charP_ (=='(') + spacesP_ + y <- withSpaces action + charP_ (==')') + spacesP_ + return y + +natP :: Parser Int +natP = P (T.decimal) + +naturalP :: Parser Integer +naturalP = P (T.decimal) + +integerP :: Parser Integer +integerP = P (T.signed T.decimal) + +identifierP :: Parser String +identifierP = do + x <- charP isAlpha + xs <- charsP $ \c -> isAlpha c || isDigit c || (c == '_') + return (x:xs) + +-------------------------------------------------------------------------------- +-- * parsing polynomials + +varP :: Parser Var +varP = withSpaces identifierP + +kstP :: Parser Integer +kstP = withSpaces integerP + +varPowP :: Parser (Var,Int) +varPowP = do + v <- varP + mb <- try $ do + charP_ (=='^') + spacesP_ + expo <- natP + spacesP_ + return expo + case mb of + Nothing -> return (v,1) + Just e -> return (v,e) + +_monomP :: Parser [(Var,Int)] +_monomP = do + ve <- varPowP + mb <- try $ do + charP_ (=='*') + spacesP_ + ves <- _monomP + return ves + case mb of + Nothing -> return (ve:[] ) + Just ves -> return (ve:ves) + +monomP :: Parser (Monom Var) +monomP = Monom <$> _monomP + +_termP :: Parser (Integer, Monom Var) +_termP = do + pm <- withSpaces optionalSignP + mbcf <- try (withSpaces integerP) + void $ try $ withSpaces (charP_ (=='*')) + monom <- case mbcf of + Nothing -> monomP + Just _ -> maybe (Monom []) id <$> try monomP + let cf = maybe 1 id mbcf + return (negateIfMinus pm cf , monom) + +termP :: Parser (Term Integer Var) +termP = (uncurry Term) <$> _termP + +genPolyP :: Parser (GenPoly Integer Var) +genPolyP = (GenPoly . filter isNotZero) <$> (spacesP_ >> many1 termP) where + isNotZero (Term cf _) = cf /= 0 + +-------------------------------------------------------------------------------- +-- * Parsing polynomial expressions + +exprP :: Parser (Expr Var) +exprP = do + spacesP_ + level3 + +atomicP = (VarE <$> varP) <||> (KstE <$> kstP) <||> inParensP exprP + +level0 = atomicP +level1 = powP <||> atomicP +level2 = productP <||> powP <||> atomicP +level3 = sumP <||> productP <||> powP <||> atomicP + +powP :: Parser (Expr Var) +powP = do + (e,n) <- _powP + return $ case n of + 0 -> KstE 1 + 1 -> e + _ -> PowE e n + +_powP :: Parser (Expr Var,Int) +_powP = do + e <- level0 + spacesP_ + charP_ (=='^') + spacesP_ + n <- natP + spacesP_ + return (e,n) + +productP :: Parser (Expr Var) +productP = do + es <- _productP + return $ case es of + [x] -> x + _ -> MulE es + +_productP :: Parser [Expr Var] +_productP = do + e1 <- level1 + spacesP_ + mb <- try $ do + charP_ (=='*') + spacesP_ + es <- _productP + return es + case mb of + Nothing -> return (e1:[]) + Just es -> return (e1:es) + +sumP :: Parser (Expr Var) +sumP = do + es <- _sumP + return $ case es of + [(Plus ,x)] -> x + [(Minus,x)] -> NegE x + _ -> LinE es + +_sumP :: Parser [(Sign,Expr Var)] +_sumP = do + pm <- optionalSignP + (this,rest) <- __sumP + return ((pm,this):rest) + +__sumP :: Parser ( Expr Var , [(Sign,Expr Var)] ) +__sumP = do + e1 <- level2 + spacesP_ + mb <- try $ do + pm <- signP + spacesP_ + (e,lin) <- __sumP + return ((pm,e):lin) + case mb of + Nothing -> return (e1,[] ) + Just lin -> return (e1,lin) + +-------------------------------------------------------------------------------- + +parseStringExpr :: Text -> Either String (Expr Var) +parseStringExpr text = case runParser (withEof exprP) text of + Right (y,_) -> Right y + Left msg -> Left msg + +parseStringGenPoly :: Text -> Either String (GenPoly Integer Var) +parseStringGenPoly text = case runParser (withEof genPolyP) text of + Right (y,_) -> Right y + Left msg -> Left msg + +-------------------------------------------------------------------------------- + +parseGenPoly :: forall vars. VariableSet vars => Proxy vars -> Text -> Maybe (GenPoly Integer VarIdx) +parseGenPoly pxy text = case parseStringGenPoly text of + Left {} -> Nothing + Right gp -> sequence $ fmap (recogVarName pxy) gp + +parseExpr :: forall vars. VariableSet vars => Proxy vars -> Text -> Maybe (Expr VarIdx) +parseExpr pxy text = case parseStringExpr text of + Left {} -> Nothing + Right e -> sequence $ fmap (recogVarName pxy) e + +-------------------------------------------------------------------------------- diff --git a/src/Math/Singular/Factory/Polynomial.hs b/src/Math/Singular/Factory/Polynomial.hs new file mode 100644 --- /dev/null +++ b/src/Math/Singular/Factory/Polynomial.hs @@ -0,0 +1,158 @@ + +-- | High-level bindings to singular-factory + +{-# LANGUAGE + BangPatterns, PatternSynonyms, KindSignatures, DataKinds, + FlexibleInstances, TypeSynonymInstances, ScopedTypeVariables, + EmptyDataDecls + #-} +module Math.Singular.Factory.Polynomial where + +-------------------------------------------------------------------------------- + +import GHC.TypeLits +import Data.Proxy + +import Data.Text.Lazy ( Text , pack ) + +import System.IO.Unsafe as Unsafe + +import Math.Singular.Factory.Internal.CanonicalForm +import Math.Singular.Factory.Internal.Factory + +import Math.Singular.Factory.Variables +import Math.Singular.Factory.Domains + +import Math.Singular.Factory.Expr +import Math.Singular.Factory.Parser + +-------------------------------------------------------------------------------- +-- * Polynomials + +-- | A multivariate polynomial over a base domain. +-- +-- Typically, you want to fix your variable set (see the module "Math.Singular.Factory.Variables"), +-- make a type synonym, and use that; for example: +-- +-- > type Poly domain = Polynomial (VarN "x") domain +-- +newtype Polynomial varset domain + = Poly { unPoly :: CF } + +instance Eq (Polynomial vars domain) where + (==) (Poly cf1) (Poly cf2) = safeEqCF cf1 cf2 + +instance forall vars domain. VariableSet vars => Show (Polynomial vars domain) where + show (Poly cf) = showCF_with (varIdxName (Proxy :: Proxy vars)) cf + +polyIsZero :: Polynomial vars domain -> Bool +polyIsZero (Poly cf) = isZeroCF cf + +polyIsOne :: Polynomial vars domain -> Bool +polyIsOne (Poly cf) = isOneCF cf + +-- | Returns true if the polynomial is a constant +inBaseDomain :: Polynomial vars domain -> Bool +inBaseDomain (Poly cf) = isInBaseDomainCF cf + +-- | If it is a constant, returns the value +mbConstant :: BaseDomain domain => Polynomial vars domain -> Maybe domain +mbConstant (Poly cf) = if isInBaseDomainCF cf + then Just (unsafeCfToBase cf) + else Nothing + +-- | A constant polynomial +konst :: BaseDomain domain => domain -> Polynomial vars domain +konst = Poly . baseToCF + +-- | A variable as a polynomial +var :: VarIdx -> Polynomial vars domain +var idx = Poly $ varCF (theNthVar idx) + +-- | A power of a variable +varPow :: VarIdx -> Int -> Polynomial vars domain +varPow idx expo = Poly $ varPowCF (theNthVar idx) expo + +-------------------------------------------------------------------------------- +-- * Operations on polynomials + +mapIntoDomain + :: forall domain1 domain2 vars. (BaseDomain domain1, BaseDomain domain2) + => Polynomial vars domain1 -> Polynomial vars domain2 +mapIntoDomain (Poly cf) = result where + result = Poly (mapIntoCF (factoryChar pxy) cf) + pxy = Proxy :: Proxy domain2 + +instance forall vars domain. BaseDomain domain => Num (Polynomial vars domain) where + fromInteger = Poly . (baseToCF :: domain -> CF) . fromInteger + negate (Poly cf) = Poly (negate cf) + (Poly cf1) + (Poly cf2) = Poly (cf1 + cf2) + (Poly cf1) - (Poly cf2) = Poly (cf1 - cf2) + (Poly cf1) * (Poly cf2) = Poly (cf1 * cf2) + abs = id + signum = const 1 + +pow :: BaseDomain domain => Polynomial vars domain -> Int -> Polynomial vars domain +pow (Poly cf) expo = Poly $ powCF cf expo + +-- | Polynomial GCD +polyGCD :: BaseDomain domain => Polynomial vars domain -> Polynomial vars domain -> Polynomial vars domain +polyGCD (Poly cf1) (Poly cf2) = Poly $ gcdPolyCF cf1 cf2 + +-- | Polynomial reduction +polyReduce :: BaseDomain domain => Polynomial vars domain -> Polynomial vars domain -> Polynomial vars domain +polyReduce (Poly cf1) (Poly cf2) = Poly $ reduceCF cf1 cf2 + +-- | Polynomial factorization +factorize :: BaseDomain domain => Polynomial vars domain -> [(Polynomial vars domain, Int)] +factorize (Poly cf) = map f (factorizeCF cf) where + f (p, expo) = (Poly p, expo) + +{- +-- | Polynomial factorization (the new algorithms by Martin Lee) +factorizeNew :: forall vars domain. BaseDomain domain => Polynomial vars domain -> [(Polynomial vars domain, Int)] +factorizeNew (Poly cf) = + case factoryChar (Proxy :: Proxy domain) of + CharZero -> map f (ratFactorizeCF cf alphaOne True) + CharFp p -> map f (fpFactorizeCF cf True) + CharGF p n -> map f (gfFactorizeCF cf True) + where + f (p, expo) = (Poly p, expo) + alphaOne = theNthVar 1 -- const Variable& v= Variable (1) + +-- | Polynomial factorization (old way) +factorizeOld :: BaseDomain domain => Polynomial vars domain -> [(Polynomial vars domain, Int)] +factorizeOld (Poly cf) = map f (oldFactorizeCF cf) where + f (p, expo) = (Poly p, expo) +-} + +-- | Substitution +substitute1 :: BaseDomain domain => VarIdx -> Polynomial vars domain -> Polynomial vars domain -> Polynomial vars domain +substitute1 idx (Poly what) (Poly cf) = Poly (substituteCF (theNthVar idx) what cf) + +-- | Evaluate a polynomial at the given point +evaluate :: BaseDomain domain => (VarIdx -> domain) -> Polynomial vars domain -> domain +evaluate fun (Poly cf0) = unsafeCfToBase (go 1 cf0) where + go :: Int -> CF -> CF + go !idx !cf = if isInBaseDomainCF cf + then cf + else go (idx+1) (substituteCF (theNthVar idx) (baseToCF $ fun idx) cf) + +-------------------------------------------------------------------------------- +-- * Parsing + +-- | Parse a polynomial in expanded form +parsePolynomial :: forall vars. VariableSet vars => Text -> Maybe (Polynomial vars Integer) +parsePolynomial text = evalGenPoly konst var <$> parseGenPoly (Proxy :: Proxy vars) text + +-- | Parse a polynomial expression (for example the product of two polynomials) +parsePolyExpr :: forall vars. VariableSet vars => Text -> Maybe (Polynomial vars Integer) +parsePolyExpr text = evalExpr {- konst -} var <$> parseExpr (Proxy :: Proxy vars) text + +parsePolynomialStr :: forall vars. VariableSet vars => String -> Maybe (Polynomial vars Integer) +parsePolynomialStr = parsePolynomial . pack + +parsePolyExprStr :: forall vars. VariableSet vars => String -> Maybe (Polynomial vars Integer) +parsePolyExprStr = parsePolyExpr . pack + +-------------------------------------------------------------------------------- diff --git a/src/Math/Singular/Factory/Variables.hs b/src/Math/Singular/Factory/Variables.hs new file mode 100644 --- /dev/null +++ b/src/Math/Singular/Factory/Variables.hs @@ -0,0 +1,208 @@ + +-- | Variables. +-- +-- Singular-factory supports only a single linear sequence of variables, +-- indexed starting from 1, optionally having single-character names. +-- +-- So what we do is to use (phantom) types to encode naming conventions. +-- + +{-# LANGUAGE BangPatterns, DataKinds, ScopedTypeVariables, KindSignatures #-} +module Math.Singular.Factory.Variables where + +-------------------------------------------------------------------------------- + +import Data.Char +import Data.List ( findIndex ) +import Text.Read + +import Data.Proxy +import GHC.TypeLits + +import System.IO.Unsafe as Unsafe + +import Math.Singular.Factory.Internal.CanonicalForm -- Factory + +-------------------------------------------------------------------------------- +-- * Raw factory variables + +-- | A variable index. +-- +-- In factory, there is a single linear sequence of variables. +-- Variables are indexed starting from 1. +-- +type VarIdx = Int + +-- | In factory, there is a single linear sequence of variables. +-- We \"precalculate\" these (lazily). +theFactoryVars :: [Var] +theFactoryVars = map mk [1..] where + mk i = Unsafe.unsafePerformIO $ newVarL i + +theNthVar :: VarIdx -> Var +theNthVar idx = theFactoryVars !! (idx-1) + +-------------------------------------------------------------------------------- +-- * Variable sets + +-- | The class of variable sets. Since Factory only supports a single linear +-- variable set, these differ only by naming conventions. +-- +class VariableSet v where + varIdxName :: Proxy v -> VarIdx -> String + recogVarName :: Proxy v -> String -> Maybe VarIdx + +-- | The variable set @x1, x2, x3, x4...@ (where \"x\" can be any string) +data VarN (s :: Symbol) + +-- | The variable set @x_1, x_2, x_3, x_4...@ (where \"x\" can be any string) +data Var_N (s :: Symbol) + +-- | The variable set @x[1], x[2], x[3], x[4]...@ (where \"x\" can be any string) +data VarBracketN (s :: Symbol) + +-- | The variable set @a, b, c, d...@ +data VarAbc + +-- | The variable set @A, B, C, D...@ +data VarABC + +-- | The variable set @x, y, z, u, v, w, a, b, c...@ +data VarXyz + +-- | The variable set @X, Y, Z, U, V, W, A, B, C...@ +data VarXYZ + +instance forall s. KnownSymbol s => VariableSet (VarN s) where + varIdxName _ = indexedVars (symbolVal (Proxy :: Proxy s)) + recogVarName _ = recogIndexed (symbolVal (Proxy :: Proxy s)) + +instance forall s. KnownSymbol s => VariableSet (Var_N s) where + varIdxName _ = indexedVarsUnderscore (symbolVal (Proxy :: Proxy s)) + recogVarName _ = recogIndexedUnderscore (symbolVal (Proxy :: Proxy s)) + +instance forall s. KnownSymbol s => VariableSet (VarBracketN s) where + varIdxName _ = indexedVarsBracket (symbolVal (Proxy :: Proxy s)) + recogVarName _ = recogIndexedBracket (symbolVal (Proxy :: Proxy s)) + +instance VariableSet VarAbc where + varIdxName _ = abcVars + recogVarName _ = recogAbc + +instance VariableSet VarABC where + varIdxName _ = capitalAbcVars + recogVarName _ = recogABC + +instance VariableSet VarXyz where + varIdxName _ = xyzVars + recogVarName _ = recogXyz + +instance VariableSet VarXYZ where + varIdxName _ = capitalXyzVars + recogVarName _ = recogXYZ + +-------------------------------------------------------------------------------- +-- * Standard naming conventions of variables + +-- | Eg. @x1, x2, x3...@ +indexedVars :: String -> VarIdx -> String +indexedVars prefix = \i -> prefix ++ show i + +-- | Eg. @x_1, x_2, x_3...@ +indexedVarsUnderscore :: String -> VarIdx -> String +indexedVarsUnderscore prefix = \i -> prefix ++ "_" ++ show i + +-- | Eg. @x[1], x[2], x[3]...@ +indexedVarsBracket :: String -> VarIdx -> String +indexedVarsBracket prefix = \i -> prefix ++ "[" ++ show i ++ "]" + +-- | That is, @a, b, c...@ +abcVars :: VarIdx -> String +abcVars idx = lowerVarList !! (idx-1) + +-- | That is, @A, B, C...@ +capitalAbcVars :: VarIdx -> String +capitalAbcVars idx = lowerVarList !! (idx-1) + +-- | @x, y, z, u, v, w, a, b, c ... , t@ +xyzVars :: VarIdx -> String +xyzVars idx = [ varListXyz !! (idx-1) ] + +capitalXyzVars :: VarIdx -> String +capitalXyzVars idx = [ varListCapitalXYZ !! (idx-1) ] + +-------------------------------------------------------------------------------- +-- * Variable lists + +varListXyz :: [Char] +varListXyz = "xyzuvwabcdefghijklmnopqrst" + +varListCapitalXYZ :: [Char] +varListCapitalXYZ = "XYZUVWABCDEFGHIJKLMNOPQRST" + +-- | The infinite list of variables @a, b ..., z, aa, ab, ac, ...@ +lowerVarList :: [String] +lowerVarList = map (:[]) abc ++ [ ys ++ [y] | ys<-lowerVarList , y<-abc ] where + abc = ['a'..'z'] + +upperVarList :: [String] +upperVarList = map (:[]) abc ++ [ ys ++ [y] | ys<-lowerVarList , y<-abc ] where + abc = ['A'..'Z'] + +-------------------------------------------------------------------------------- +-- * Parsing standard variable names + +readPosIdxMaybe :: String -> Maybe VarIdx +readPosIdxMaybe s = case readMaybe s :: Maybe Word of + Nothing -> Nothing + Just w -> let j = (fromIntegral w :: Int) + in if j >= 1 then Just j else Nothing + +recogIndexed :: String -> String -> Maybe VarIdx +recogIndexed !prefix = recog where + !n = length prefix + recog s = case splitAt n s of + (p,q) -> if p /= prefix + then Nothing + else readPosIdxMaybe q + +recogIndexedUnderscore :: String -> String -> Maybe VarIdx +recogIndexedUnderscore !prefix = recog where + !n = length prefix + recog s = case splitAt n s of + (p,q) -> if p /= prefix + then Nothing + else case q of + ('_':r ) -> readPosIdxMaybe r + _ -> Nothing + +recogIndexedBracket :: String -> String -> Maybe VarIdx +recogIndexedBracket !prefix = recog where + !n = length prefix + recog s = case splitAt n s of + (p,q) -> if p /= prefix + then Nothing + else if length q >= 3 && head q == '[' && last q == ']' + then readPosIdxMaybe (init $ tail q) + else Nothing + +recogAbc :: String -> Maybe VarIdx +recogAbc [c] = if c >= 'a' && c <= 'z' then Just (ord c - 96) else Nothing +recogAbc [] = Nothing +recogAbc _ = Nothing + +recogABC :: String -> Maybe VarIdx +recogABC [c] = if c >= 'A' && c <= 'Z' then Just (ord c - 64) else Nothing +recogABC [] = Nothing +recogABC _ = Nothing + +recogXyz :: String -> Maybe VarIdx +recogXyz [c] = if c >= 'a' && c <= 'z' then (+1) <$> findIndex (==c) varListXyz else Nothing +recogXyz _ = Nothing + +recogXYZ :: String -> Maybe VarIdx +recogXYZ [c] = if c >= 'A' && c <= 'Z' then (+1) <$> findIndex (==c) varListCapitalXYZ else Nothing +recogXYZ _ = Nothing + +-------------------------------------------------------------------------------- + diff --git a/test/TestSuite.hs b/test/TestSuite.hs new file mode 100644 --- /dev/null +++ b/test/TestSuite.hs @@ -0,0 +1,256 @@ + +-- | The test-suite + +{-# LANGUAGE BangPatterns, DataKinds, DeriveFunctor #-} +module Main where + +-------------------------------------------------------------------------------- + +import Data.List + +import Control.Monad +import System.Random + +import Test.Tasty +import Test.Tasty.HUnit + +import Math.Singular.Factory + +import Math.Singular.Factory.Domains +import Math.Singular.Factory.Polynomial +import Math.Singular.Factory.Counting +import Math.Singular.Factory.Variables + +-------------------------------------------------------------------------------- + +main = do + initialize + printConfig + defaultMain tests + +tests :: TestTree +tests = testGroup "Tests" + [ unit_tests +-- , randomized_tests + ] + +unit_tests :: TestTree +unit_tests = testGroup "Unit tests" + [ testCase "reconstruction for some polys over Z" (assertBool "failed" $ reconstr_some_polys some_polys_ZZ) + , testGroup "reconstr. of some more polys over Z" + [ testCase ("poly #" ++ show i) (assertBool "failed" $ prop_reconstr_from_factors p) + | (i,p) <- zip [0..] some_more_polys_ZZ + ] + , testCase "reconstruction for some polys over F_2" (assertBool "failed" $ reconstr_some_polys some_polys_F2) + , testCase "reconstruction for some polys over F_3" (assertBool "failed" $ reconstr_some_polys some_polys_F3) + , testCase "reconstruction for some polys over F_5" (assertBool "failed" $ reconstr_some_polys some_polys_F5) + , testCase "reconstruction for some polys over GF(4)" (assertBool "failed" $ reconstr_some_polys some_polys_GF4) + , testCase "reconstruction for some polys over GF(8)" (assertBool "failed" $ reconstr_some_polys some_polys_GF8) + , testCase "reconstruction for some polys over GF(9)" (assertBool "failed" $ reconstr_some_polys some_polys_GF9) + ] + +-------------------------------------------------------------------------------- + +type Poly domain = Polynomial VarAbc domain + +some_polys_ZZ :: [Poly Integer] +some_polys_ZZ = + [ (x^2 + 1) * (x^5 + 1) * (x^7 + 1) + + , x^2 - y^2 + , y^2 - z^2 + , x^5 - y^5 + , x^7 - y^7 + , x^10 - y^10 + + , (1 + x + y + z)^2 + , (1 + x + y + z)^4 + , (1 + x + y + z)^8 + + , (1 + x + y + z)^2 - 1 + , (1 + x + y + z)^4 - 1 + , (1 + x + y + z)^8 - 1 + + , (1 + x^2 + y^2 + z^2)^2 - 1 + , (1 + x^2 + y^2 + z^2)^4 - 1 + , (1 + x^2 + y^2 + z^2)^8 - 1 + + , ((1 - x - y - z)^2 - 1) * ((1 - x - y - z )^2 + 1) + , ((1 + x + y + z)^2 - 1) * ((1 - x - y - z )^2 + 1) + , ((1 + x + y + z)^2 - 1) * ((1 + x + y + z )^2 + 1) + ] + where + myvars@[x,y,z,u,v,w] = map var [1..6] + +-- | These seem to be too big for finite fields except GF(2) ??? +some_more_polys_ZZ :: [Poly Integer] +some_more_polys_ZZ= + [ (x-13) * (x+27) * (x-42)^2 * (x+7)^3 * (x-11)^5 + , ( (1+x+y+z)^4 + 1 ) * ( (1+x+y+z)^4 + 2 ) + , ( (1+x+y+z)^5 + 1 ) * ( (1+x+y+z)^5 + 2 ) + , ( (1+x+y+z)^7 + 1 ) * ( (1+x+y+z)^7 + 2 ) + , ( (1+x+y+z)^10 + 1 ) * ( (1+x+y+z)^10 + 2 ) + , ( (1+x+y+z)^20 + 1 ) * ( (1+x+y+z)^20 + 2 ) + ] ++ + [ sparse1 hk | hk <- [5,8..25] ] ++ + [ sparse2 k | k <- [5..15] ] + where + myvars@[x,y,z,u,v,w] = map var [1..6] + +-- from: +-- Martin Mok-Don Lee: Factorization of multivariate polynomials +sparse1 :: Int -> Poly Integer +sparse1 halfk + = (x^2*y^2*z + x*(y^k+z^k) + 3*y + 3*z - 3*z^2 - 2*y^halfk*z^halfk ) + * (x^2*y^2*z^2 + x*(y^k+z^k) - 2*y - 5*z + 4*y^2 + 3*y^halfk*z^halfk ) + where + k = 2*halfk + [x,y,z] = map var [1..3] + +-- from: +-- Martin Mok-Don Lee: Factorization of multivariate polynomials +sparse2 :: Int -> Poly Integer +sparse2 k + = ( (x*( y^3+2*z^3) + 5*y*z)*(x*(y+4*z)+2) + (2*x-7)*(y^k*z^k-y^(k-1)*z^(k-1)) ) + * ( (x*(3*y^3+4*z^3) + 3*y*z)*(x*(y+3*z)+7) - (3*x+5)*(y^k*z^k-y^(k-1)*z^(k-1)) ) + where + [x,y,z] = map var [1..3] + +-------------------------------------------------------------------------------- + + +some_polys_F2 :: [Poly (FF 2)] +some_polys_F2 = map mapIntoDomain some_polys_ZZ + +some_polys_F3 :: [Poly (FF 3)] +some_polys_F3 = map mapIntoDomain some_polys_ZZ + +some_polys_F5 :: [Poly (FF 5)] +some_polys_F5 = map mapIntoDomain some_polys_ZZ + +some_polys_GF4 :: [Poly (GF 2 2 "q")] +some_polys_GF4 = map mapIntoDomain some_polys_ZZ + +some_polys_GF8 :: [Poly (GF 2 3 "q")] +some_polys_GF8 = map mapIntoDomain some_polys_ZZ + +some_polys_GF9 :: [Poly (GF 3 2 "q")] +some_polys_GF9 = map mapIntoDomain some_polys_ZZ + +-------------------------------------------------------------------------------- +-- TODO: collect examples where the factorization is known + +data KnownFactors poly = KnownFactors + { _poly :: poly + , _facs :: [(poly,Int)] + } + deriving (Show,Functor) + +-- A bivariate polynomial over F2 +-- taken from +-- F. K. Abu Salem: An efficient sparse adaptation of the polytope method over Fp +-- and a record-high binary bivariate factorisation +abuSalemBenchmark = KnownFactors + { _poly = "x^4120 + x^4118*y^2 + x^3708*y^400 + x^3706*y^402+ x^2781*y^1300 + x^2779*y^1302 + x^1339*y^2700+ x^927*y^3100 + y^4000 + x^7172*y^4167 + x^8349*y^4432+ x^8347*y^4434 + x^6760*y^4567 + x^5833*y^5467+ x^5568*y^7132 + x^11401*y^8599" + , _facs = [ ("x^5568*y^4432 + x^1339 + x^927*y^400 + y^1300" , 1) + , ("x^5833*y^4167 + x^2781 + x^2779*y^2 + y^2700" , 1) + ] + } + +-------------------------------------------------------------------------------- + +reconstructFromFactors :: BaseDomain dom => [(Poly dom, Int)] -> Poly dom +reconstructFromFactors = product . map (uncurry pow) + +prop_reconstr_from_factors :: BaseDomain dom => Poly dom -> Bool +prop_reconstr_from_factors p = (p == reconstructFromFactors (factorize p)) + +reconstr_some_polys :: BaseDomain dom => [Poly dom] -> Bool +reconstr_some_polys list = and (map prop_reconstr_from_factors list) + +-------------------------------------------------------------------------------- + +{- +counting_main = do + + let [x,y,z] = map var [1..3] + let f = - y^2*z + x^3 + 2*z^3 :: Poly V Integer + print f + + let [x,y,z,u,v,w] = map var [1..6] + let g1 = x*v - y*u + g2 = x*w - z*u + g3 = y*w - z*v + gs = [g1,g2,g3] :: [Poly V Integer] +-} + +{- + cnts @(FF 5) (map mapIntoDomain gs) + cntsP @(FF 5) (map mapIntoDomain gs) +-} + +{- + putStrLn "==========================\np=2" + cnt @(FF 2) (mapIntoDomain f) + cnt @(GF 2 2 "x") (mapIntoDomain f) + cnt @(GF 2 3 "x") (mapIntoDomain f) + cnt @(GF 2 4 "x") (mapIntoDomain f) + cnt @(GF 2 5 "x") (mapIntoDomain f) + cnt @(GF 2 6 "x") (mapIntoDomain f) + cnt @(GF 2 7 "x") (mapIntoDomain f) + + + putStrLn "==========================\np=3" + cnt @(FF 3) (mapIntoDomain f) + cnt @(GF 3 2 "x") (mapIntoDomain f) + cnt @(GF 3 3 "x") (mapIntoDomain f) + cnt @(GF 3 4 "x") (mapIntoDomain f) + cnt @(GF 3 5 "x") (mapIntoDomain f) + cnt @(GF 3 6 "x") (mapIntoDomain f) + cnt @(GF 3 7 "x") (mapIntoDomain f) +-} + +{- + putStrLn "==========================\np=5" + cntP @(FF 5) (mapIntoDomain f) + cntP @(GF 5 2 "x") (mapIntoDomain f) + cntP @(GF 5 3 "x") (mapIntoDomain f) + cntP @(GF 5 4 "x") (mapIntoDomain f) + cntP @(GF 5 5 "x") (mapIntoDomain f) +-} + +{- + putStrLn "==========================\np=7" + cntP @(FF 7) (mapIntoDomain f) + cntP @(GF 7 2 "x") (mapIntoDomain f) + cntP @(GF 7 3 "x") (mapIntoDomain f) + cntP @(GF 7 4 "x") (mapIntoDomain f) + cntP @(GF 7 5 "x") (mapIntoDomain f) + +cnt :: forall domain. FiniteDomain domain => Poly V domain -> IO () +cnt poly = do + let cnt = countAffineHypersurface 3 poly + let q = domainSize (Proxy :: Proxy domain) + print (q,cnt,divMod (cnt-1) (q-1)) + +cntP :: forall domain. FiniteDomain domain => Poly V domain -> IO () +cntP poly = do + let cnt = countProjectiveHypersurface 3 poly + let q = domainSize (Proxy :: Proxy domain) + print (q,cnt) + +cnts :: forall domain. FiniteDomain domain => [Poly V domain] -> IO () +cnts polys = do + let cnt = countAffineSolutions 6 polys + let q = domainSize (Proxy :: Proxy domain) + print (q,cnt,divMod (cnt-1) (q-1)) + +cntsP :: forall domain. FiniteDomain domain => [Poly V domain] -> IO () +cntsP polys = do + let cnt = countProjectiveSolutions 6 polys + let q = domainSize (Proxy :: Proxy domain) + print (q,cnt) + +-} + +--------------------------------------------------------------------------------