diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,675 @@ + GNU GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The GNU General Public License is a free, copyleft license for +software and other kinds of works. + + The licenses for most software and other practical works are designed +to take away your freedom to share and change the works. By contrast, +the GNU General Public License is intended to guarantee your freedom to +share and change all versions of a program--to make sure it remains free +software for all its users. We, the Free Software Foundation, use the +GNU General Public License for most of our software; it applies also to +any other work released this way by its authors. You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +them if you wish), that you receive source code or can get it if you +want it, that you can change the software or use pieces of it in new +free programs, and that you know you can do these things. + + To protect your rights, we need to prevent others from denying you +these rights or asking you to surrender the rights. Therefore, you have +certain responsibilities if you distribute copies of the software, or if +you modify it: responsibilities to respect the freedom of others. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must pass on to the recipients the same +freedoms that you received. You must make sure that they, too, receive +or can get the source code. And you must show them these terms so they +know their rights. + + Developers that use the GNU GPL protect your rights with two steps: +(1) assert copyright on the software, and (2) offer you this License +giving you legal permission to copy, distribute and/or modify it. + + For the developers' and authors' protection, the GPL clearly explains +that there is no warranty for this free software. For both users' and +authors' sake, the GPL requires that modified versions be marked as +changed, so that their problems will not be attributed erroneously to +authors of previous versions. + + Some devices are designed to deny users access to install or run +modified versions of the software inside them, although the manufacturer +can do so. This is fundamentally incompatible with the aim of +protecting users' freedom to change the software. The systematic +pattern of such abuse occurs in the area of products for individuals to +use, which is precisely where it is most unacceptable. Therefore, we +have designed this version of the GPL to prohibit the practice for those +products. If such problems arise substantially in other domains, we +stand ready to extend this provision to those domains in future versions +of the GPL, as needed to protect the freedom of users. + + Finally, every program is threatened constantly by software patents. +States should not allow patents to restrict development and use of +software on general-purpose computers, but in those that do, we wish to +avoid the special danger that patents applied to a free program could +make it effectively proprietary. To prevent this, the GPL assures that +patents cannot be used to render the program non-free. + + The precise terms and conditions for copying, distribution and +modification follow. + + TERMS AND CONDITIONS + + 0. Definitions. + + "This License" refers to version 3 of the GNU General Public License. + + "Copyright" also means copyright-like laws that apply to other kinds of +works, such as semiconductor masks. + + "The Program" refers to any copyrightable work licensed under this +License. Each licensee is addressed as "you". "Licensees" and +"recipients" may be individuals or organizations. + + To "modify" a work means to copy from or adapt all or part of the work +in a fashion requiring copyright permission, other than the making of an +exact copy. The resulting work is called a "modified version" of the +earlier work or a work "based on" the earlier work. + + A "covered work" means either the unmodified Program or a work based +on the Program. + + To "propagate" a work means to do anything with it that, without +permission, would make you directly or secondarily liable for +infringement under applicable copyright law, except executing it on a +computer or modifying a private copy. Propagation includes copying, +distribution (with or without modification), making available to the +public, and in some countries other activities as well. + + To "convey" a work means any kind of propagation that enables other +parties to make or receive copies. Mere interaction with a user through +a computer network, with no transfer of a copy, is not conveying. + + An interactive user interface displays "Appropriate Legal Notices" +to the extent that it includes a convenient and prominently visible +feature that (1) displays an appropriate copyright notice, and (2) +tells the user that there is no warranty for the work (except to the +extent that warranties are provided), that licensees may convey the +work under this License, and how to view a copy of this License. If +the interface presents a list of user commands or options, such as a +menu, a prominent item in the list meets this criterion. + + 1. Source Code. + + The "source code" for a work means the preferred form of the work +for making modifications to it. "Object code" means any non-source +form of a work. + + A "Standard Interface" means an interface that either is an official +standard defined by a recognized standards body, or, in the case of +interfaces specified for a particular programming language, one that +is widely used among developers working in that language. + + The "System Libraries" of an executable work include anything, other +than the work as a whole, that (a) is included in the normal form of +packaging a Major Component, but which is not part of that Major +Component, and (b) serves only to enable use of the work with that +Major Component, or to implement a Standard Interface for which an +implementation is available to the public in source code form. A +"Major Component", in this context, means a major essential component +(kernel, window system, and so on) of the specific operating system +(if any) on which the executable work runs, or a compiler used to +produce the work, or an object code interpreter used to run it. + + The "Corresponding Source" for a work in object code form means all +the source code needed to generate, install, and (for an executable +work) run the object code and to modify the work, including scripts to +control those activities. However, it does not include the work's +System Libraries, or general-purpose tools or generally available free +programs which are used unmodified in performing those activities but +which are not part of the work. For example, Corresponding Source +includes interface definition files associated with source files for +the work, and the source code for shared libraries and dynamically +linked subprograms that the work is specifically designed to require, +such as by intimate data communication or control flow between those +subprograms and other parts of the work. + + The Corresponding Source need not include anything that users +can regenerate automatically from other parts of the Corresponding +Source. + + The Corresponding Source for a work in source code form is that +same work. + + 2. Basic Permissions. + + All rights granted under this License are granted for the term of +copyright on the Program, and are irrevocable provided the stated +conditions are met. This License explicitly affirms your unlimited +permission to run the unmodified Program. The output from running a +covered work is covered by this License only if the output, given its +content, constitutes a covered work. This License acknowledges your +rights of fair use or other equivalent, as provided by copyright law. + + You may make, run and propagate covered works that you do not +convey, without conditions so long as your license otherwise remains +in force. You may convey covered works to others for the sole purpose +of having them make modifications exclusively for you, or provide you +with facilities for running those works, provided that you comply with +the terms of this License in conveying all material for which you do +not control copyright. Those thus making or running the covered works +for you must do so exclusively on your behalf, under your direction +and control, on terms that prohibit them from making any copies of +your copyrighted material outside their relationship with you. + + Conveying under any other circumstances is permitted solely under +the conditions stated below. Sublicensing is not allowed; section 10 +makes it unnecessary. + + 3. Protecting Users' Legal Rights From Anti-Circumvention Law. + + No covered work shall be deemed part of an effective technological +measure under any applicable law fulfilling obligations under article +11 of the WIPO copyright treaty adopted on 20 December 1996, or +similar laws prohibiting or restricting circumvention of such +measures. + + When you convey a covered work, you waive any legal power to forbid +circumvention of technological measures to the extent such circumvention +is effected by exercising rights under this License with respect to +the covered work, and you disclaim any intention to limit operation or +modification of the work as a means of enforcing, against the work's +users, your or third parties' legal rights to forbid circumvention of +technological measures. + + 4. Conveying Verbatim Copies. + + You may convey verbatim copies of the Program's source code as you +receive it, in any medium, provided that you conspicuously and +appropriately publish on each copy an appropriate copyright notice; +keep intact all notices stating that this License and any +non-permissive terms added in accord with section 7 apply to the code; +keep intact all notices of the absence of any warranty; and give all +recipients a copy of this License along with the Program. + + You may charge any price or no price for each copy that you convey, +and you may offer support or warranty protection for a fee. + + 5. Conveying Modified Source Versions. + + You may convey a work based on the Program, or the modifications to +produce it from the Program, in the form of source code under the +terms of section 4, provided that you also meet all of these conditions: + + a) The work must carry prominent notices stating that you modified + it, and giving a relevant date. + + b) The work must carry prominent notices stating that it is + released under this License and any conditions added under section + 7. This requirement modifies the requirement in section 4 to + "keep intact all notices". + + c) You must license the entire work, as a whole, under this + License to anyone who comes into possession of a copy. This + License will therefore apply, along with any applicable section 7 + additional terms, to the whole of the work, and all its parts, + regardless of how they are packaged. This License gives no + permission to license the work in any other way, but it does not + invalidate such permission if you have separately received it. + + d) If the work has interactive user interfaces, each must display + Appropriate Legal Notices; however, if the Program has interactive + interfaces that do not display Appropriate Legal Notices, your + work need not make them do so. + + A compilation of a covered work with other separate and independent +works, which are not by their nature extensions of the covered work, +and which are not combined with it such as to form a larger program, +in or on a volume of a storage or distribution medium, is called an +"aggregate" if the compilation and its resulting copyright are not +used to limit the access or legal rights of the compilation's users +beyond what the individual works permit. Inclusion of a covered work +in an aggregate does not cause this License to apply to the other +parts of the aggregate. + + 6. Conveying Non-Source Forms. + + You may convey a covered work in object code form under the terms +of sections 4 and 5, provided that you also convey the +machine-readable Corresponding Source under the terms of this License, +in one of these ways: + + a) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by the + Corresponding Source fixed on a durable physical medium + customarily used for software interchange. + + b) Convey the object code in, or embodied in, a physical product + (including a physical distribution medium), accompanied by a + written offer, valid for at least three years and valid for as + long as you offer spare parts or customer support for that product + model, to give anyone who possesses the object code either (1) a + copy of the Corresponding Source for all the software in the + product that is covered by this License, on a durable physical + medium customarily used for software interchange, for a price no + more than your reasonable cost of physically performing this + conveying of source, or (2) access to copy the + Corresponding Source from a network server at no charge. + + c) Convey individual copies of the object code with a copy of the + written offer to provide the Corresponding Source. This + alternative is allowed only occasionally and noncommercially, and + only if you received the object code with such an offer, in accord + with subsection 6b. + + d) Convey the object code by offering access from a designated + place (gratis or for a charge), and offer equivalent access to the + Corresponding Source in the same way through the same place at no + further charge. You need not require recipients to copy the + Corresponding Source along with the object code. If the place to + copy the object code is a network server, the Corresponding Source + may be on a different server (operated by you or a third party) + that supports equivalent copying facilities, provided you maintain + clear directions next to the object code saying where to find the + Corresponding Source. Regardless of what server hosts the + Corresponding Source, you remain obligated to ensure that it is + available for as long as needed to satisfy these requirements. + + e) Convey the object code using peer-to-peer transmission, provided + you inform other peers where the object code and Corresponding + Source of the work are being offered to the general public at no + charge under subsection 6d. + + A separable portion of the object code, whose source code is excluded +from the Corresponding Source as a System Library, need not be +included in conveying the object code work. + + A "User Product" is either (1) a "consumer product", which means any +tangible personal property which is normally used for personal, family, +or household purposes, or (2) anything designed or sold for incorporation +into a dwelling. In determining whether a product is a consumer product, +doubtful cases shall be resolved in favor of coverage. For a particular +product received by a particular user, "normally used" refers to a +typical or common use of that class of product, regardless of the status +of the particular user or of the way in which the particular user +actually uses, or expects or is expected to use, the product. A product +is a consumer product regardless of whether the product has substantial +commercial, industrial or non-consumer uses, unless such uses represent +the only significant mode of use of the product. + + "Installation Information" for a User Product means any methods, +procedures, authorization keys, or other information required to install +and execute modified versions of a covered work in that User Product from +a modified version of its Corresponding Source. The information must +suffice to ensure that the continued functioning of the modified object +code is in no case prevented or interfered with solely because +modification has been made. + + If you convey an object code work under this section in, or with, or +specifically for use in, a User Product, and the conveying occurs as +part of a transaction in which the right of possession and use of the +User Product is transferred to the recipient in perpetuity or for a +fixed term (regardless of how the transaction is characterized), the +Corresponding Source conveyed under this section must be accompanied +by the Installation Information. But this requirement does not apply +if neither you nor any third party retains the ability to install +modified object code on the User Product (for example, the work has +been installed in ROM). + + The requirement to provide Installation Information does not include a +requirement to continue to provide support service, warranty, or updates +for a work that has been modified or installed by the recipient, or for +the User Product in which it has been modified or installed. Access to a +network may be denied when the modification itself materially and +adversely affects the operation of the network or violates the rules and +protocols for communication across the network. + + Corresponding Source conveyed, and Installation Information provided, +in accord with this section must be in a format that is publicly +documented (and with an implementation available to the public in +source code form), and must require no special password or key for +unpacking, reading or copying. + + 7. Additional Terms. + + "Additional permissions" are terms that supplement the terms of this +License by making exceptions from one or more of its conditions. +Additional permissions that are applicable to the entire Program shall +be treated as though they were included in this License, to the extent +that they are valid under applicable law. If additional permissions +apply only to part of the Program, that part may be used separately +under those permissions, but the entire Program remains governed by +this License without regard to the additional permissions. + + When you convey a copy of a covered work, you may at your option +remove any additional permissions from that copy, or from any part of +it. (Additional permissions may be written to require their own +removal in certain cases when you modify the work.) You may place +additional permissions on material, added by you to a covered work, +for which you have or can give appropriate copyright permission. + + Notwithstanding any other provision of this License, for material you +add to a covered work, you may (if authorized by the copyright holders of +that material) supplement the terms of this License with terms: + + a) Disclaiming warranty or limiting liability differently from the + terms of sections 15 and 16 of this License; or + + b) Requiring preservation of specified reasonable legal notices or + author attributions in that material or in the Appropriate Legal + Notices displayed by works containing it; or + + c) Prohibiting misrepresentation of the origin of that material, or + requiring that modified versions of such material be marked in + reasonable ways as different from the original version; or + + d) Limiting the use for publicity purposes of names of licensors or + authors of the material; or + + e) Declining to grant rights under trademark law for use of some + trade names, trademarks, or service marks; or + + f) Requiring indemnification of licensors and authors of that + material by anyone who conveys the material (or modified versions of + it) with contractual assumptions of liability to the recipient, for + any liability that these contractual assumptions directly impose on + those licensors and authors. + + All other non-permissive additional terms are considered "further +restrictions" within the meaning of section 10. If the Program as you +received it, or any part of it, contains a notice stating that it is +governed by this License along with a term that is a further +restriction, you may remove that term. If a license document contains +a further restriction but permits relicensing or conveying under this +License, you may add to a covered work material governed by the terms +of that license document, provided that the further restriction does +not survive such relicensing or conveying. + + If you add terms to a covered work in accord with this section, you +must place, in the relevant source files, a statement of the +additional terms that apply to those files, or a notice indicating +where to find the applicable terms. + + Additional terms, permissive or non-permissive, may be stated in the +form of a separately written license, or stated as exceptions; +the above requirements apply either way. + + 8. Termination. + + You may not propagate or modify a covered work except as expressly +provided under this License. Any attempt otherwise to propagate or +modify it is void, and will automatically terminate your rights under +this License (including any patent licenses granted under the third +paragraph of section 11). + + However, if you cease all violation of this License, then your +license from a particular copyright holder is reinstated (a) +provisionally, unless and until the copyright holder explicitly and +finally terminates your license, and (b) permanently, if the copyright +holder fails to notify you of the violation by some reasonable means +prior to 60 days after the cessation. + + Moreover, your license from a particular copyright holder is +reinstated permanently if the copyright holder notifies you of the +violation by some reasonable means, this is the first time you have +received notice of violation of this License (for any work) from that +copyright holder, and you cure the violation prior to 30 days after +your receipt of the notice. + + Termination of your rights under this section does not terminate the +licenses of parties who have received copies or rights from you under +this License. If your rights have been terminated and not permanently +reinstated, you do not qualify to receive new licenses for the same +material under section 10. + + 9. Acceptance Not Required for Having Copies. + + You are not required to accept this License in order to receive or +run a copy of the Program. Ancillary propagation of a covered work +occurring solely as a consequence of using peer-to-peer transmission +to receive a copy likewise does not require acceptance. However, +nothing other than this License grants you permission to propagate or +modify any covered work. These actions infringe copyright if you do +not accept this License. Therefore, by modifying or propagating a +covered work, you indicate your acceptance of this License to do so. + + 10. Automatic Licensing of Downstream Recipients. + + Each time you convey a covered work, the recipient automatically +receives a license from the original licensors, to run, modify and +propagate that work, subject to this License. You are not responsible +for enforcing compliance by third parties with this License. + + An "entity transaction" is a transaction transferring control of an +organization, or substantially all assets of one, or subdividing an +organization, or merging organizations. If propagation of a covered +work results from an entity transaction, each party to that +transaction who receives a copy of the work also receives whatever +licenses to the work the party's predecessor in interest had or could +give under the previous paragraph, plus a right to possession of the +Corresponding Source of the work from the predecessor in interest, if +the predecessor has it or can get it with reasonable efforts. + + You may not impose any further restrictions on the exercise of the +rights granted or affirmed under this License. For example, you may +not impose a license fee, royalty, or other charge for exercise of +rights granted under this License, and you may not initiate litigation +(including a cross-claim or counterclaim in a lawsuit) alleging that +any patent claim is infringed by making, using, selling, offering for +sale, or importing the Program or any portion of it. + + 11. Patents. + + A "contributor" is a copyright holder who authorizes use under this +License of the Program or a work on which the Program is based. The +work thus licensed is called the contributor's "contributor version". + + A contributor's "essential patent claims" are all patent claims +owned or controlled by the contributor, whether already acquired or +hereafter acquired, that would be infringed by some manner, permitted +by this License, of making, using, or selling its contributor version, +but do not include claims that would be infringed only as a +consequence of further modification of the contributor version. For +purposes of this definition, "control" includes the right to grant +patent sublicenses in a manner consistent with the requirements of +this License. + + Each contributor grants you a non-exclusive, worldwide, royalty-free +patent license under the contributor's essential patent claims, to +make, use, sell, offer for sale, import and otherwise run, modify and +propagate the contents of its contributor version. + + In the following three paragraphs, a "patent license" is any express +agreement or commitment, however denominated, not to enforce a patent +(such as an express permission to practice a patent or covenant not to +sue for patent infringement). To "grant" such a patent license to a +party means to make such an agreement or commitment not to enforce a +patent against the party. + + If you convey a covered work, knowingly relying on a patent license, +and the Corresponding Source of the work is not available for anyone +to copy, free of charge and under the terms of this License, through a +publicly available network server or other readily accessible means, +then you must either (1) cause the Corresponding Source to be so +available, or (2) arrange to deprive yourself of the benefit of the +patent license for this particular work, or (3) arrange, in a manner +consistent with the requirements of this License, to extend the patent +license to downstream recipients. "Knowingly relying" means you have +actual knowledge that, but for the patent license, your conveying the +covered work in a country, or your recipient's use of the covered work +in a country, would infringe one or more identifiable patents in that +country that you have reason to believe are valid. + + If, pursuant to or in connection with a single transaction or +arrangement, you convey, or propagate by procuring conveyance of, a +covered work, and grant a patent license to some of the parties +receiving the covered work authorizing them to use, propagate, modify +or convey a specific copy of the covered work, then the patent license +you grant is automatically extended to all recipients of the covered +work and works based on it. + + A patent license is "discriminatory" if it does not include within +the scope of its coverage, prohibits the exercise of, or is +conditioned on the non-exercise of one or more of the rights that are +specifically granted under this License. You may not convey a covered +work if you are a party to an arrangement with a third party that is +in the business of distributing software, under which you make payment +to the third party based on the extent of your activity of conveying +the work, and under which the third party grants, to any of the +parties who would receive the covered work from you, a discriminatory +patent license (a) in connection with copies of the covered work +conveyed by you (or copies made from those copies), or (b) primarily +for and in connection with specific products or compilations that +contain the covered work, unless you entered into that arrangement, +or that patent license was granted, prior to 28 March 2007. + + Nothing in this License shall be construed as excluding or limiting +any implied license or other defenses to infringement that may +otherwise be available to you under applicable patent law. + + 12. No Surrender of Others' Freedom. + + If conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot convey a +covered work so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you may +not convey it at all. For example, if you agree to terms that obligate you +to collect a royalty for further conveying from those to whom you convey +the Program, the only way you could satisfy both those terms and this +License would be to refrain entirely from conveying the Program. + + 13. Use with the GNU Affero General Public License. + + Notwithstanding any other provision of this License, you have +permission to link or combine any covered work with a work licensed +under version 3 of the GNU Affero General Public License into a single +combined work, and to convey the resulting work. The terms of this +License will continue to apply to the part which is the covered work, +but the special requirements of the GNU Affero General Public License, +section 13, concerning interaction through a network will apply to the +combination as such. + + 14. Revised Versions of this License. + + The Free Software Foundation may publish revised and/or new versions of +the GNU General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + + Each version is given a distinguishing version number. If the +Program specifies that a certain numbered version of the GNU General +Public License "or any later version" applies to it, you have the +option of following the terms and conditions either of that numbered +version or of any later version published by the Free Software +Foundation. If the Program does not specify a version number of the +GNU General Public License, you may choose any version ever published +by the Free Software Foundation. + + If the Program specifies that a proxy can decide which future +versions of the GNU General Public License can be used, that proxy's +public statement of acceptance of a version permanently authorizes you +to choose that version for the Program. + + Later license versions may give you additional or different +permissions. However, no additional obligations are imposed on any +author or copyright holder as a result of your choosing to follow a +later version. + + 15. Disclaimer of Warranty. + + THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY +APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT +HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY +OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, +THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM +IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF +ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. Limitation of Liability. + + IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS +THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY +GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE +USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF +DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD +PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), +EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF +SUCH DAMAGES. + + 17. Interpretation of Sections 15 and 16. + + If the disclaimer of warranty and limitation of liability provided +above cannot be given local legal effect according to their terms, +reviewing courts shall apply local law that most closely approximates +an absolute waiver of all civil liability in connection with the +Program, unless a warranty or assumption of liability accompanies a +copy of the Program in return for a fee. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +state the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + +Also add information on how to contact you by electronic and paper mail. + + If the program does terminal interaction, make it output a short +notice like this when it starts in an interactive mode: + + Copyright (C) + This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, your program's commands +might be different; for a GUI interface, you would use an "about box". + + You should also get your employer (if you work as a programmer) or school, +if any, to sign a "copyright disclaimer" for the program, if necessary. +For more information on this, and how to apply and follow the GNU GPL, see +. + + The GNU General Public License does not permit incorporating your program +into proprietary programs. If your program is a subroutine library, you +may consider it more useful to permit linking proprietary applications with +the library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. But first, please read +. + diff --git a/Numeric/FFT.hs b/Numeric/FFT.hs new file mode 100644 --- /dev/null +++ b/Numeric/FFT.hs @@ -0,0 +1,40 @@ +-- | Mixed-radix FFT calculation. +-- +-- Arbitrary input vector lengths are handled using a mixed-radix +-- Cooley-Tukey decimation in time algorithm with residual prime +-- length vectors being treated using Rader's algorithm or hand-coded +-- codelets for small primes. +module Numeric.FFT + ( fft, ifft, fftWith, ifftWith + , plan, planFromFactors, execute + , Plan (..), Direction (..), BaseTransform (..) + ) where + +import Prelude hiding (length, map, sum, zipWith) +import Data.Vector +import Data.Complex + +import Numeric.FFT.Types +import Numeric.FFT.Plan +import Numeric.FFT.Execute + + +-- | Forward FFT with embedded plan calculation. +fft :: Vector (Complex Double) -> IO (Vector (Complex Double)) +fft xs = do + p <- plan $ length xs + return $ fftWith p xs + +-- | Inverse FFT with embedded plan calculation. +ifft :: Vector (Complex Double) -> IO (Vector (Complex Double)) +ifft xs = do + p <- plan $ length xs + return $ ifftWith p xs + +-- | Forward FFT with pre-computed plan. +fftWith :: Plan -> Vector (Complex Double) -> Vector (Complex Double) +fftWith p = convert . execute p Forward . convert + +-- | Inverse FFT with pre-computed plan. +ifftWith :: Plan -> Vector (Complex Double) -> Vector (Complex Double) +ifftWith p = convert . execute p Inverse . convert diff --git a/Numeric/FFT/Execute.hs b/Numeric/FFT/Execute.hs new file mode 100644 --- /dev/null +++ b/Numeric/FFT/Execute.hs @@ -0,0 +1,228 @@ +module Numeric.FFT.Execute ( execute ) where + +import Prelude hiding (concatMap, foldr, length, map, mapM_, + null, reverse, sum, zip, zipWith) +import qualified Prelude as P +import Control.Monad (when) +import qualified Control.Monad as CM +import Control.Monad.ST +import Control.Monad.Primitive (PrimMonad) +import Data.Complex +import Data.STRef +import qualified Data.Vector as V +import Data.Vector.Unboxed +import qualified Data.Vector.Unboxed.Mutable as MV +import qualified Data.IntMap.Strict as IM +import qualified Data.Map as M + +import Numeric.FFT.Types +import Numeric.FFT.Utils +import Numeric.FFT.Special + + +-- | Main FFT plan execution driver. +execute :: Plan -> Direction -> VCD -> VCD +execute (Plan dlinfo perm base) dir h = + if n == 1 then h else if V.null dlinfo + then runST $ do + mhin <- case perm of + Nothing -> thaw h + Just p -> unsafeThaw $ backpermute h p + mhout <- MV.replicate n 0 + applyBase base sign mhin mhout + when (dir == Inverse) $ do + let s = 1.0 / fromIntegral n :+ 0 + CM.forM_ [0..n-1] $ \i -> do + x <- MV.unsafeRead mhout i + MV.unsafeWrite mhout i $ s * x + unsafeFreeze mhout + else fullfft + where + n = length h -- Input vector length. + bsize = baseSize base -- Size of base transform. + + -- Root of unity sign. + sign = case dir of + Forward -> 1 + Inverse -> -1 + + -- Apply Danielson-Lanczos steps and base transform to digit + -- reversal ordered input vector. + fullfft = runST $ do + mhin <- case perm of + Nothing -> thaw h + Just p -> unsafeThaw $ backpermute h p + mhtmp <- MV.replicate n 0 + multBase mhin mhtmp + mhr <- newSTRef (mhtmp, mhin) + V.forM_ dlinfo $ \dlstep -> do + (mh0, mh1) <- readSTRef mhr + dl sign dlstep mh0 mh1 + writeSTRef mhr (mh1, mh0) + mhs <- readSTRef mhr + let vout = fst mhs + when (dir == Inverse) $ do + let s = 1.0 / fromIntegral n :+ 0 + CM.forM_ [0..n-1] $ \i -> do + x <- MV.unsafeRead vout i + MV.unsafeWrite vout i $ s * x + unsafeFreeze vout + + -- Multiple base transform application for "bottom" of algorithm. + multBase :: MVCD s -> MVCD s -> ST s () + multBase xmin xmout = + V.zipWithM_ (applyBase base sign) + (slicemvecs bsize xmin) (slicemvecs bsize xmout) + + +-- | Monadic FFT plan execution driver -- used by Rader's algorithm +-- for convolutions. +executeM :: Plan -> Direction -> MVCD s -> MVCD s -> ST s () +executeM (Plan dlinfo perm base) dir hin hout = + if n == 1 + then MV.copy hout hin + else do + htmp <- MV.replicate n 0 + + -- Input permutation. + case perm of + Nothing -> MV.copy htmp hin + Just p -> backpermuteM n p hin htmp + + -- Apply Danielson-Lanczos steps and base transform to digit + -- reversal ordered input vector. + multBase htmp hout + mhr <- newSTRef (hout, htmp) + V.forM_ dlinfo $ \dlstep -> do + (mh0, mh1) <- readSTRef mhr + dl sign dlstep mh0 mh1 + writeSTRef mhr (mh1, mh0) + when (odd $ V.length dlinfo) $ MV.copy hout htmp + + -- Output scaling for inverse transform. + when (dir == Inverse) $ do + let s = 1.0 / fromIntegral n :+ 0 + forM_ (enumFromN 0 n) $ \i -> do + x <- MV.unsafeRead hout i + MV.unsafeWrite hout i $ s * x + where + n = MV.length hin -- Input vector length. + bsize = baseSize base -- Size of base transform. + + -- Root of unity sign. + sign = case dir of + Forward -> 1 + Inverse -> -1 + + -- Multiple base transform application for "bottom" of algorithm. + multBase :: MVCD s -> MVCD s -> ST s () + multBase xmin xmout = + V.zipWithM_ (applyBase base sign) + (slicemvecs bsize xmin) (slicemvecs bsize xmout) + + +-- | Single Danielson-Lanczos step: process all duplicates and +-- concatenate into a single vector. +dl :: Int -> (Int, Int, VVVCD, VVVCD) -> MVCD s -> MVCD s -> ST s () +dl sign (wfac, split, dmatp, dmatm) mhin mhout = + V.zipWithM_ doone (slicemvecs wfac mhin) (slicemvecs wfac mhout) + where + -- Twiddled diagonal entries in row r, column c (both + -- zero-indexed), where each row and column if a wfac x wfac + -- matrix. + dmat = if sign == 1 then dmatp else dmatm + d r c = (dmat V.! r) V.! c + + -- Size of each diagonal sub-matrix. + ns = wfac `div` split + + -- Index vectors. + nsidxs = enumFromN 0 ns + splitidxs = enumFromN 1 (split-1) + + -- Process one duplicate by processing all rows and writing the + -- results into a single output vector. + doone :: MVCD s -> MVCD s -> ST s () + doone vin vout = do + let vs = (slicemvecs ns vin, slicemvecs ns vout) + mapM_ (single vs) $ enumFromN 0 split + where + -- Multiply a single block by its appropriate diagonal + -- elements and accumulate the result. + mult :: VMVCD s -> MVCD s -> Int -> Bool -> Int -> ST s () + mult vins vo r first c = do + let vi = vins V.! c + dvals = d r c + forM_ nsidxs $ \i -> do + xi <- MV.unsafeRead vi i + xo <- if first then return 0 else MV.unsafeRead vo i + MV.unsafeWrite vo i (xo + xi * dvals ! i) + -- Multiply all blocks by the corresponding diagonal + -- elements in a single row. + single :: (VMVCD s, VMVCD s) -> Int -> ST s () + single (vis, vos) r = do + mult vis (vos V.! r) r True 0 + mapM_ (mult vis (vos V.! r) r False) splitidxs + -- single (vis, vos) r = + -- let m = mult vis (vos V.! r) r + -- in do + -- m True 0 + -- mapM_ (m False) splitidxs + + +-- | Apply a base transform to a single vector. +applyBase :: BaseTransform -> Int -> MVCD s -> MVCD s -> ST s () + +-- Simple DFT algorithm. +applyBase (DFTBase sz wsfwd wsinv) sign mhin mhout = do + h <- freeze mhin + forM_ (enumFromN 0 sz) $ \i -> MV.unsafeWrite mhout i (doone h i) + where ws = if sign == 1 then wsfwd else wsinv + doone h i = sum $ zipWith (*) h $ + generate sz (\k -> ws ! (i * k `mod` sz)) + +-- Special hard-coded cases. +applyBase (SpecialBase sz) sign mhin mhout = + case IM.lookup sz specialBases of + Just f -> f sign mhin mhout + Nothing -> error "invalid problem size for SpecialBase" + +-- Rader prime-length FFT. +applyBase (RaderBase sz outperm bfwd binv csz cplan) sign mhin mhout = do + -- Padding size. + let pad = csz - (sz - 1) + + -- Permuted input vector padded to next greater power of two size + -- for fast convolution. + apad <- MV.replicate csz 0 + forM_ (enumFromN 0 csz) $ \i -> do + val <- if i == 0 then MV.unsafeRead mhin 1 + else if i > pad + then MV.unsafeRead mhin $ i - pad + 1 + else return 0 + MV.unsafeWrite apad i val + + -- FFT-based convolution calculation. + convtmp <- MV.replicate csz 0 + executeM cplan Forward apad convtmp + let bmult = if sign == 1 then bfwd else binv + forM_ (enumFromN 0 csz) $ \i -> do + x <- MV.unsafeRead convtmp i + MV.unsafeWrite convtmp i $ x * (bmult ! i) + executeM cplan Inverse convtmp apad + conv <- unsafeFreeze apad + + -- Input vector sum. + sumhref <- newSTRef 0 + forM_ (enumFromN 0 sz) $ \i -> do + val <- MV.unsafeRead mhin i + modifySTRef sumhref (+ val) + sumh <- readSTRef sumhref + + -- Write output based on output generator index ordering. + h0 <- MV.unsafeRead mhin 0 + forM_ (enumFromN 0 sz) $ \i -> do + let (idx, val) = case i of + 0 -> (0, sumh) + _ -> (outperm ! (i - 1), h0 + conv ! (i - 1)) + MV.unsafeWrite mhout idx val diff --git a/Numeric/FFT/Plan.hs b/Numeric/FFT/Plan.hs new file mode 100644 --- /dev/null +++ b/Numeric/FFT/Plan.hs @@ -0,0 +1,312 @@ +module Numeric.FFT.Plan ( plan, planFromFactors ) where + +import Prelude hiding ((++), any, concatMap, enumFromTo, filter, length, map, + maximum, null, reverse, scanl, sum, zip, zipWith) +import qualified Prelude as P +import Control.Applicative ((<$>)) +import qualified Control.Monad as CM +import Data.Complex +import Data.Function (on) +import Data.IORef +import qualified Data.IntMap.Strict as IM +import Data.List (nub, (\\)) +import qualified Data.List as L +import Data.Ord +import qualified Data.Set as S +import qualified Data.Vector as V +import Data.Vector.Unboxed +import Control.Monad.IO.Class +import System.Directory +import System.Environment +import System.FilePath +import System.IO +import System.IO.Unsafe (unsafePerformIO) +import Criterion +import Criterion.Config +import Criterion.Monad +import Criterion.Environment + +import Numeric.FFT.Types +import Numeric.FFT.Execute +import Numeric.FFT.Utils +import Numeric.FFT.Special + + +-- | Number of plans to test empirically. +nTestPlans :: Int +nTestPlans = 50 + +-- | Globally shared timing environment. (Not thread-safe...) +timingEnv :: IORef (Maybe Environment) +timingEnv = unsafePerformIO (newIORef Nothing) +{-# NOINLINE timingEnv #-} + +-- | Plan calculation for a given problem size. +plan :: Int -> IO Plan +plan 1 = return $ Plan V.empty Nothing (SpecialBase 1) +plan n = do + wis <- readWisdom n + let fixRader p = case plBase p of + bpl@(RaderBase _ _ _ _ csz _) -> do + cplan <- liftIO $ plan csz + return $ p { plBase = bpl { raderConvPlan = cplan } } + _ -> return p + pret <- case wis of + Just (p, t) -> planFromFactors n p + Nothing -> do + let ps = testPlans n nTestPlans + withConfig (defaultConfig { cfgVerbosity = ljust Quiet + , cfgSamples = ljust 1 }) $ do + menv <- liftIO $ readIORef timingEnv + env <- case menv of + Just e -> return e + Nothing -> do + meas <- measureEnvironment + liftIO $ writeIORef timingEnv $ Just meas + return meas + let v = generate n (\i -> sin (2 * pi * fromIntegral i / 511) :+ 0) + tps <- CM.forM ps $ \p -> do + ptest <- liftIO $ planFromFactors n p >>= fixRader + ts <- runBenchmark env $ nf (execute ptest Forward) v + return (sum ts / fromIntegral (length ts), p) + let (rest, resp) = L.minimumBy (compare `on` fst) tps + liftIO $ writeWisdom n resp rest + liftIO $ planFromFactors n resp + fixRader pret + +-- | Get execution time for plan for a given problem size. +planTime :: Int -> IO Double +planTime n = do + wis <- readWisdom n + case wis of + Just (_, t) -> return t + Nothing -> do + -- Force generation of wisdom if it's not already there. + _ <- plan n + Just (_, t) <- readWisdom n + return t + +-- | Plan calculation for a given problem factorisation. +planFromFactors :: Int -> (Int, Vector Int) -> IO Plan +planFromFactors n (lastf, fs) = do + -- Base transform. + (base, mextraperm) <- makeBase lastf + + -- Include permutation of base transform if needed. + let perm = case (digperm, mextraperm) of + (Just dp, Just ep) -> Just $ dupperm n ep %.% dp + (Nothing, Just ep) -> Just $ dupperm n ep + (Just dp, Nothing) -> Just dp + (Nothing, Nothing) -> Nothing + + return $ Plan dlinfo perm base + where + -- Input data "digit reversal" permutation. + digperm = digrev n fs + + -- Size information for Danielson-Lanczos steps. + wfacs = map (n `div`) $ scanl (*) 1 fs + vwfacs = convert wfacs + vfs = convert fs + dmatps = V.zipWith (dmat 1) vwfacs vfs + dmatms = V.zipWith (dmat (-1)) vwfacs vfs + dlinfo = V.reverse $ V.zip4 vwfacs vfs dmatps dmatms + + -- Calculate diagonal matrix entries used in Danielson-Lanczos steps. + dmat :: Int -> Int -> Int -> VVVCD + dmat sign wfac split = + let ns = wfac `div` split + w = omega $ sign * wfac + in V.generate split $ + \r -> V.generate split $ + \c -> map (w^(ns*r*c) *) $ map ((w^^) . (c *)) $ enumFromN 0 ns + +-- | Read from wisdom for a given problem size. +readWisdom :: Int -> IO (Maybe ((Int, Vector Int), Double)) +readWisdom n = do + home <- getEnv "HOME" + let wisf = home ".fft-plan" show n + ex <- doesFileExist wisf + case ex of + False -> return Nothing + True -> do + wist <- readFile wisf + let ((wisb, wisfs), wistim) = read wist :: ((Int, [Int]), Double) + return $ Just ((wisb, fromList wisfs), wistim) + +-- | Write wisdom for a given problem size. +writeWisdom :: Int -> (Int, Vector Int) -> Double -> IO () +writeWisdom n (b, fs) tim = do + home <- getEnv "HOME" + let wisd = home ".fft-plan" + wisf = wisd show n + createDirectoryIfMissing True wisd + writeFile wisf $ show ((b, toList fs), tim) P.++ "\n" + +-- | Make base transform for a given sub-problem size. +makeBase :: Int -> IO (BaseTransform, Maybe VI) +makeBase sz + | sz `IM.member` specialBases = return (SpecialBase sz, Nothing) + | isPrime sz = makeRaderBase sz + | otherwise = return (makeDFTBase sz, Nothing) + +-- | Generate digit reversal permutation using elementary "modulo" +-- permutations: last digit is not permuted to match with using a +-- simple DFT or instance of Rader's algorithm at the "bottom" of the +-- overall algorithm. +digrev :: Int -> VI -> Maybe VI +digrev n fs + | null fs = Nothing + | otherwise = Just $ V.foldl1' (%.%) $ V.map (dupperm n) subperms + where + vfs = convert fs + + -- Sizes of the individual permutations that we need, one per + -- factor. + sizes = V.scanl div n vfs + + -- Partial sub-permutations, one per factor. + subperms = V.reverse $ V.zipWith perm sizes vfs + + -- Generate a single "modulo" permutation. + perm sz fac = concatMap doone $ enumFromN 0 fac + where doone i = generate (sz `div` fac) (\j -> j * fac + i) + +-- | Pre-computation plan for basic DFT transform. +makeDFTBase :: Int -> BaseTransform +makeDFTBase sz = DFTBase sz wsfwd wsinv + where w = omega sz + wsfwd = generate sz (w ^) + wsinv = map (1 /) wsfwd + +-- | Pre-compute plan for prime-length Rader FFT transform. +makeRaderBase :: Int -> IO (BaseTransform, Maybe VI) +makeRaderBase sz = do + -- Forward FFT with embedded plan calculation. + let fft p xs = execute p Forward xs + + -- Times for unpadded and padded convolution transforms. + unpadtime <- planTime sz1 + padtime <- planTime pow2sz + + -- Should we use a padded or an unpadded transform for the + -- convolution? + let pad = padtime < unpadtime + csz = if pad then pow2sz else sz1 + + -- Plan for convolution transforms. + cplan <- plan csz + + -- FFT transforms of b sequences for use in convolution, one for + -- forward transform and one for inverse transform. Either just the + -- basic root of unity powers for the convolution (if we're not + -- padding), or the powers cyclically repeated to make a vector of + -- the next power-of-two length. + let convb = + fft cplan $ if pad + then generate csz (\idx -> bs ! (idx `mod` sz1)) + else bs + convbinv = + fft cplan $ if pad + then generate csz (\idx -> bsinv ! (idx `mod` sz1)) + else bsinv + + return (RaderBase sz outperm convb convbinv csz cplan, Just inperm) + where + -- Convolution length. + sz1 = sz - 1 + + -- Convolution length padded to next greater power of two. + pow2sz = if sz1 == 2^(log2 sz1) + then sz1 + else 2 ^ (1 + log2 (2 * sz1 - 3)) + + -- Group generator and inverse group generator. + g = primitiveRoot sz + ig = invModN sz g + + -- Input value permutation according to group generator indexing. + inperm = 0 `cons` iterateN sz1 (\n -> (g * n) `mod` sz) 1 + + -- Index vector based on inverse group generator ordering. + outperm = iterateN sz1 (\n -> (ig * n) `mod` sz) 1 + + -- Root of unity powers based on inverse group generator indexing, + -- for forward and inverse transform. + w = omega sz + bs = backpermute (map (w ^^) $ enumFromTo 0 sz1) outperm + bsinv = backpermute (map ((w ^^) . negate) $ enumFromTo 0 sz1) outperm + +-- | Base transform type with heuristic ordering. +data BaseType = Special Int | Rader Int deriving (Eq, Show) + +-- | Newtype wrapper for custom sorting. +newtype SPlan = SPlan (BaseType, Vector Int) deriving (Eq, Show) + +-- | Base transform size. +bSize :: BaseType -> Int +bSize (Special b) = b +bSize (Rader b) = b + +-- | Heuristic ordering for base transform types: special bases come +-- first, then prime bases using Rader's algorithm, ordered according +-- to size compensating for padding needed in the Rader's algorithm +-- convolution. +instance Ord BaseType where + compare (Special _) (Rader _) = LT + compare (Rader _) (Special _) = GT + compare (Special s1) (Special s2) = compare s2 s1 + compare (Rader r1) (Rader r2) = case (isPow2 $ r1 - 1, isPow2 $ r2 - 1) of + (True, True) -> compare r1 r2 + (True, False) -> compare r1 (2 * r2) + (False, True) -> compare (2 * r1) r2 + (False, False) -> compare r1 r2 + +-- | Heuristic ordering for full plans, based first on base type, then +-- on the maximum size of Danielson-Lanczos step. +instance Ord SPlan where + compare (SPlan (b1, fs1)) (SPlan (b2, fs2)) = case compare b1 b2 of + LT -> LT + EQ -> compare (maximum fs2) (maximum fs1) + GT -> GT + +-- | Generate test plans for a given input size, sorted in heuristic +-- order. +testPlans :: Int -> Int -> [(Int, Vector Int)] +testPlans n nplans = L.take nplans $ L.map clean $ L.sort okplans + where vfs = allFactors n + bs = usableBases n vfs + doone b = basePlans n vfs b + clean (SPlan (b, fs)) = (bSize b, fs) + allplans = P.concatMap doone bs + okplans = case L.filter (not . ridiculous) allplans of + [] -> L.filter (not . reallyRidiculous) allplans + oks -> oks + ridiculous (SPlan (_, fs)) = any (> 128) fs + reallyRidiculous (SPlan (_, fs)) = + any (> 128) $ filter (not . isPrime) fs + +-- | List plans from a single base. +basePlans :: Int -> Vector Int -> BaseType -> [SPlan] +basePlans n vfs bt = if null lfs + then [SPlan (bt, empty)] + else P.map (\v -> SPlan (bt, v)) $ leftOvers lfs + where lfs = fromList $ (toList vfs) \\ (toList $ allFactors b) + b = bSize bt + +-- | Produce all distinct permutations and compositions constructable +-- from a given list of factors. +leftOvers :: Vector Int -> [Vector Int] +leftOvers fs = + if null fs + then [] + else S.toList $ L.foldl' go S.empty (multisetPerms fs) + where n = length fs + go fset perm = foldl' doone fset (enumFromN 0 (2^(n - 1))) + where doone s i = S.insert (makeComp perm i) s + +-- | Usable base transform sizes. +usableBases :: Int -> Vector Int -> [BaseType] +usableBases n fs = P.map Special bs P.++ P.map Rader ps + where bs = toList $ filter ((== 0) . (n `mod`)) specialBaseSizes + ps = toList $ filter isPrime $ filter (> maxPrimeSpecialBaseSize) fs diff --git a/Numeric/FFT/Special.hs b/Numeric/FFT/Special.hs new file mode 100644 --- /dev/null +++ b/Numeric/FFT/Special.hs @@ -0,0 +1,54 @@ +module Numeric.FFT.Special + ( specialBases + , specialBaseSizes + , maxPrimeSpecialBaseSize ) where + +import Prelude hiding (filter, maximum) +import Control.Monad.ST +import Data.IntMap.Strict (IntMap) +import qualified Data.IntMap.Strict as IM +import Data.Complex +import Data.Vector.Unboxed +import qualified Data.Vector.Unboxed.Mutable as MV + +import Numeric.FFT.Types +import Numeric.FFT.Utils +import Numeric.FFT.Special.Primes +import Numeric.FFT.Special.PowersOfTwo +import Numeric.FFT.Special.Miscellaneous + + +-- | Map from input vector lengths to hard-coded FFT transforms for +-- small problem sizes. Each function in the map takes a transform +-- direction (+1 for forward, -1 for inverse) and returns the +-- /unscaled/ transform (scaling for inverse transforms is applied at +-- the top-level). +specialBases :: IntMap (Int -> MVCD s -> MVCD s -> ST s ()) +specialBases = IM.fromList [ ( 2, special2) + , ( 4, special4) + , ( 8, special8) + , (16, special16) + , (32, special32) + , (64, special64) + , ( 3, special3) + , ( 5, special5) + , ( 7, special7) + , (11, special11) + , (13, special13) + , ( 6, special6) + , ( 9, special9) + , (10, special10) + , (12, special12) + , (14, special14) + , (15, special15) + , (20, special20) + , (25, special25) + ] + +-- | Available base transform sizes. +specialBaseSizes :: Vector Int +specialBaseSizes = fromList $ IM.keys specialBases + +-- | Largest prime for which we have a specialised base transform. +maxPrimeSpecialBaseSize :: Int +maxPrimeSpecialBaseSize = maximum $ filter isPrime specialBaseSizes diff --git a/Numeric/FFT/Special/Miscellaneous.hs b/Numeric/FFT/Special/Miscellaneous.hs new file mode 100644 --- /dev/null +++ b/Numeric/FFT/Special/Miscellaneous.hs @@ -0,0 +1,779 @@ +module Numeric.FFT.Special.Miscellaneous + ( special6, special9, special10, special12, special14 + , special15, special20, special25 + ) where + +import Control.Monad.ST +import Data.IntMap.Strict (IntMap) +import qualified Data.IntMap.Strict as IM +import Data.Complex +import Data.Vector.Unboxed +import qualified Data.Vector.Unboxed.Mutable as MV + +import Numeric.FFT.Types +import Numeric.FFT.Utils + + +-- | Length 6 hard-coded FFT. +kp866025403, kp500000000 :: Double +kp866025403 = 0.866025403784438646763723170752936183471402627 +kp500000000 = 0.500000000000000000000000000000000000000000000 +special6 :: Int -> MVCD s -> MVCD s -> ST s () +special6 sign xsin xsout = do + xr0 :+ xi0 <- MV.unsafeRead xsin 0 ; xr1 :+ xi1 <- MV.unsafeRead xsin 1 + xr2 :+ xi2 <- MV.unsafeRead xsin 2 ; xr3 :+ xi3 <- MV.unsafeRead xsin 3 + xr4 :+ xi4 <- MV.unsafeRead xsin 4 ; xr5 :+ xi5 <- MV.unsafeRead xsin 5 + let tb = xr0 + xr3 ; t3 = xr0 - xr3 ; tx = xi0 + xi3 ; tp = xi0 - xi3 + tc = xr2 + xr5 ; t6 = xr2 - xr5 ; td = xr4 + xr1 ; t9 = xr4 - xr1 + te = tc + td ; tA = td - tc ; ts = t9 - t6 ; ta = t6 + t9 + tu = xi2 + xi5 ; ti = xi2 - xi5 ; tf = t3 - kp500000000 * ta + tv = xi4 + xi1 ; tl = xi4 - xi1 ; tt = tb - kp500000000 * te + ty = tu + tv ; tw = tu - tv ; tq = ti + tl ; tm = ti - tl + tr = tp - kp500000000 * tq ; tz = tx - kp500000000 * ty + r5 = (tf + kp866025403 * tm) :+ (tr + kp866025403 * ts) + r4 = (tt - kp866025403 * tw) :+ (tz - kp866025403 * tA) + r3 = (t3 + ta) :+ (tp + tq) + r2 = (tt + kp866025403 * tw) :+ (tz + kp866025403 * tA) + r1 = (tf - kp866025403 * tm) :+ (tr - kp866025403 * ts) + MV.unsafeWrite xsout 0 $ (tb + te) :+ (tx + ty) + MV.unsafeWrite xsout 1 $ if sign == 1 then r1 else r5 + MV.unsafeWrite xsout 2 $ if sign == 1 then r2 else r4 + MV.unsafeWrite xsout 3 $ if sign == 1 then r3 else r3 + MV.unsafeWrite xsout 4 $ if sign == 1 then r4 else r2 + MV.unsafeWrite xsout 5 $ if sign == 1 then r5 else r1 + +-- | Length 9 hard-coded FFT. +kp954188894, kp363970234, kp852868531, kp984807753 :: Double +kp492403876, kp777861913, kp839099631, kp176326980 :: Double +--kp866025403, kp500000000 :: Double +kp954188894 = 0.954188894138671133499268364187245676532219158 +kp363970234 = 0.363970234266202361351047882776834043890471784 +kp852868531 = 0.852868531952443209628250963940074071936020296 +kp984807753 = 0.984807753012208059366743024589523013670643252 +kp492403876 = 0.492403876506104029683371512294761506835321626 +kp777861913 = 0.777861913430206160028177977318626690410586096 +kp839099631 = 0.839099631177280011763127298123181364687434283 +kp176326980 = 0.176326980708464973471090386868618986121633062 +--kp866025403 = 0.866025403784438646763723170752936183471402627 +--kp500000000 = 0.500000000000000000000000000000000000000000000 +special9 :: Int -> MVCD s -> MVCD s -> ST s () +special9 sign xsin xsout = do + xr0 :+ xi0 <- MV.unsafeRead xsin 0 ; xr1 :+ xi1 <- MV.unsafeRead xsin 1 + xr2 :+ xi2 <- MV.unsafeRead xsin 2 ; xr3 :+ xi3 <- MV.unsafeRead xsin 3 + xr4 :+ xi4 <- MV.unsafeRead xsin 4 ; xr5 :+ xi5 <- MV.unsafeRead xsin 5 + xr6 :+ xi6 <- MV.unsafeRead xsin 6 ; xr7 :+ xi7 <- MV.unsafeRead xsin 7 + xr8 :+ xi8 <- MV.unsafeRead xsin 8 + let t4 = xr3 + xr6 ; tm = xr6 - xr3 ; tM = xi3 - xi6 ; tk = xi3 + xi6 + tL = xr0 - kp500000000 * t4 ; t5 = xr0 + t4 + tl = xi0 - kp500000000 * tk ; t1f = xi0 + tk + tE = xr4 - xr7 ; t9 = xr4 + xr7 ; tH = xi7 - xi4 ; tC = xi4 + xi7 + ta = xr1 + t9 ; tG = xr1 - kp500000000 * t9 + t1c = xi1 + tC ; tD = xi1 - kp500000000 * tC + tI = tG - kp866025403 * tH ; tX = tG + kp866025403 * tH + tF = tD - kp866025403 * tE ; tW = tD + kp866025403 * tE + t17 = tl - kp866025403 * tm ; tn = tl + kp866025403 * tm + tw = xr8 - xr5 ; te = xr5 + xr8 ; tu = xi5 + xi8 ; tr = xi5 - xi8 + tN = tL + kp866025403 * tM ; tV = tL - kp866025403 * tM + tf = xr2 + te ; to = xr2 - kp500000000 * te + t1d = xi2 + tu ; tv = xi2 - kp500000000 * tu + ts = to + kp866025403 * tr ; tZ = to - kp866025403 * tr + tg = ta + tf ; t1i = tf - ta + tx = tv + kp866025403 * tw ; t10 = tv - kp866025403 * tw + t1e = t1c - t1d ; t1g = t1c + t1d + t1b = t5 - kp500000000 * tg ; t1h = t1f - kp500000000 * t1g + tO = tx + kp176326980 * ts ; ty = ts - kp176326980 * tx + tJ = tF - kp839099631 * tI ; tP = tI + kp839099631 * tF + tS = ty + kp777861913 * tJ ; tK = ty - kp777861913 * tJ + tU = tO - kp777861913 * tP ; tQ = tO + kp777861913 * tP + tT = tn + kp492403876 * tK ; tR = tN - kp492403876 * tQ + t14 = tX - kp176326980 * tW ; tY = tW + kp176326980 * tX + t11 = tZ - kp363970234 * t10 ; t15 = t10 + kp363970234 * tZ + t12 = tY - kp954188894 * t11 ; t1a = tY + kp954188894 * t11 + t16 = t14 - kp954188894 * t15 ; t18 = t14 + kp954188894 * t15 + t13 = tV - kp492403876 * t12 ; t19 = t17 + kp492403876 * t18 + r8 = (tN + kp984807753 * tQ) :+ (tn - kp984807753 * tK) + r7 = (tV + kp984807753 * t12) :+ (t17 - kp984807753 * t18) + r6 = (t1b + kp866025403 * t1e) :+ (t1h + kp866025403 * t1i) + r5 = (tR + kp852868531 * tS) :+ (tT + kp852868531 * tU) + r4 = (t13 - kp852868531 * t16) :+ (t19 - kp852868531 * t1a) + r3 = (t1b - kp866025403 * t1e) :+ (t1h - kp866025403 * t1i) + r2 = (tR - kp852868531 * tS) :+ (tT - kp852868531 * tU) + r1 = (t13 + kp852868531 * t16) :+ (t19 + kp852868531 * t1a) + MV.unsafeWrite xsout 0 $ (t5 + tg) :+ (t1f + t1g) + MV.unsafeWrite xsout 1 $ if sign == 1 then r1 else r8 + MV.unsafeWrite xsout 2 $ if sign == 1 then r2 else r7 + MV.unsafeWrite xsout 3 $ if sign == 1 then r3 else r6 + MV.unsafeWrite xsout 4 $ if sign == 1 then r4 else r5 + MV.unsafeWrite xsout 5 $ if sign == 1 then r5 else r4 + MV.unsafeWrite xsout 6 $ if sign == 1 then r6 else r3 + MV.unsafeWrite xsout 7 $ if sign == 1 then r7 else r2 + MV.unsafeWrite xsout 8 $ if sign == 1 then r8 else r1 + +-- | Length 10 hard-coded FFT. +kp951056516, kp559016994, kp250000000, kp618033988 :: Double +kp951056516 = 0.951056516295153572116439333379382143405698634 +kp559016994 = 0.559016994374947424102293417182819058860154590 +kp250000000 = 0.250000000000000000000000000000000000000000000 +kp618033988 = 0.618033988749894848204586834365638117720309180 +special10 :: Int -> MVCD s -> MVCD s -> ST s () +special10 sign xsin xsout = do + xr0 :+ xi0 <- MV.unsafeRead xsin 0 ; xr1 :+ xi1 <- MV.unsafeRead xsin 1 + xr2 :+ xi2 <- MV.unsafeRead xsin 2 ; xr3 :+ xi3 <- MV.unsafeRead xsin 3 + xr4 :+ xi4 <- MV.unsafeRead xsin 4 ; xr5 :+ xi5 <- MV.unsafeRead xsin 5 + xr6 :+ xi6 <- MV.unsafeRead xsin 6 ; xr7 :+ xi7 <- MV.unsafeRead xsin 7 + xr8 :+ xi8 <- MV.unsafeRead xsin 8 ; xr9 :+ xi9 <- MV.unsafeRead xsin 9 + let tj = xr0 + xr5 ; t3 = xr0 - xr5 ; t1b = xi0 + xi5 ; tN = xi0 - xi5 + tk = xr2 + xr7 ; t6 = xr2 - xr7 ; to = xr6 + xr1 ; tg = xr6 - xr1 + tl = xr8 + xr3 ; t9 = xr8 - xr3 ; tn = xr4 + xr9 ; td = xr4 - xr9 + tm = tk + tl ; t1j = tk - tl ; ta = t6 + t9 ; tU = t6 - t9 + tp = tn + to ; t1i = tn - to ; th = td + tg ; tV = td - tg + tq = tm + tp ; t10 = tm - tp ; ti = ta + th ; ts = ta - th + tw = xi2 - xi7 ; t15 = xi2 + xi7 ; t13 = xi6 + xi1 ; tG = xi6 - xi1 + t16 = xi8 + xi3 ; tz = xi8 - xi3 ; t12 = xi4 + xi9 ; tD = xi4 - xi9 + t1c = t15 + t16 ; t17 = t15 - t16 ; tO = tw + tz ; tA = tw - tz + t1d = t12 + t13 ; t14 = t12 - t13 ; tP = tD + tG ; tH = tD - tG + t1e = t1c + t1d ; t1g = t1c - t1d ; tQ = tO + tP ; tS = tO - tP + tK = tH - kp618033988 * tA ; tI = tA + kp618033988 * tH + tr = t3 - kp250000000 * ti ; tY = tV - kp618033988 * tU + tW = tU + kp618033988 * tV ; tR = tN - kp250000000 * tQ + tJ = tr - kp559016994 * ts ; tt = tr + kp559016994 * ts + t1a = t17 + kp618033988 * t14 ; t18 = t14 - kp618033988 * t17 + tX = tR - kp559016994 * tS ; tT = tR + kp559016994 * tS + tZ = tj - kp250000000 * tq ; t1m = t1j + kp618033988 * t1i + t1k = t1i - kp618033988 * t1j ; t1f = t1b - kp250000000 * t1e + t19 = tZ + kp559016994 * t10 ; t11 = tZ - kp559016994 * t10 + t1h = t1f - kp559016994 * t1g ; t1l = t1f + kp559016994 * t1g + r9 = (tt + kp951056516 * tI) :+ (tT - kp951056516 * tW) + r8 = (t11 - kp951056516 * t18) :+ (t1h + kp951056516 * t1k) + r7 = (tJ + kp951056516 * tK) :+ (tX - kp951056516 * tY) + r6 = (t19 - kp951056516 * t1a) :+ (t1l + kp951056516 * t1m) + r5 = (t3 + ti) :+ (tN + tQ) + r4 = (t19 + kp951056516 * t1a) :+ (t1l - kp951056516 * t1m) + r3 = (tJ - kp951056516 * tK) :+ (tX + kp951056516 * tY) + r2 = (t11 + kp951056516 * t18) :+ (t1h - kp951056516 * t1k) + r1 = (tt - kp951056516 * tI) :+ (tT + kp951056516 * tW) + MV.unsafeWrite xsout 0 $ (tj + tq) :+ (t1b + t1e) + MV.unsafeWrite xsout 1 $ if sign == 1 then r1 else r9 + MV.unsafeWrite xsout 2 $ if sign == 1 then r2 else r8 + MV.unsafeWrite xsout 3 $ if sign == 1 then r3 else r7 + MV.unsafeWrite xsout 4 $ if sign == 1 then r4 else r6 + MV.unsafeWrite xsout 5 $ if sign == 1 then r5 else r5 + MV.unsafeWrite xsout 6 $ if sign == 1 then r6 else r4 + MV.unsafeWrite xsout 7 $ if sign == 1 then r7 else r3 + MV.unsafeWrite xsout 8 $ if sign == 1 then r8 else r2 + MV.unsafeWrite xsout 9 $ if sign == 1 then r9 else r1 + +-- | Length 12 hard-coded FFT. +--kp866025403, kp500000000 :: Double +--kp866025403 = 0.866025403784438646763723170752936183471402627 +--kp500000000 = 0.500000000000000000000000000000000000000000000 +special12 :: Int -> MVCD s -> MVCD s -> ST s () +special12 sign xsin xsout = do + xr0 :+ xi0 <- MV.unsafeRead xsin 0 ; xr1 :+ xi1 <- MV.unsafeRead xsin 1 + xr2 :+ xi2 <- MV.unsafeRead xsin 2 ; xr3 :+ xi3 <- MV.unsafeRead xsin 3 + xr4 :+ xi4 <- MV.unsafeRead xsin 4 ; xr5 :+ xi5 <- MV.unsafeRead xsin 5 + xr6 :+ xi6 <- MV.unsafeRead xsin 6 ; xr7 :+ xi7 <- MV.unsafeRead xsin 7 + xr8 :+ xi8 <- MV.unsafeRead xsin 8 ; xr9 :+ xi9 <- MV.unsafeRead xsin 9 + xr10 :+ xi10 <- MV.unsafeRead xsin 10 ; xr11 :+ xi11 <- MV.unsafeRead xsin 11 + let t4 = xr4 + xr8 ; tA = xr8 - xr4 ; tS = xi4 - xi8 ; tr = xi4 + xi8 + tR = xr0 - kp500000000 * t4 ; t5 = xr0 + t4 ; ts = xi0 + tr + tz = xi0 - kp500000000 * tr ; t9 = xr10 + xr2 ; tD = xr2 - xr10 + tV = xi10 - xi2 ; tw = xi10 + xi2 ; tU = xr6 - kp500000000 * t9 + ta = xr6 + t9 ; tx = xi6 + tw ; tC = xi6 - kp500000000 * tw + tf = xr7 + xr11 ; t1d = xr11 - xr7 ; tJ = xi7 - xi11 ; t1b = xi7 + xi11 + tG = xr3 - kp500000000 * tf ; tg = xr3 + tf ; t1u = xi3 + t1b + t1c = xi3 - kp500000000 * t1b ; tk = xr1 + xr5 ; t1i = xr5 - xr1 + t1t = t5 - ta ; tb = t5 + ta ; tO = xi1 - xi5 ; t1g = xi1 + xi5 + tL = xr9 - kp500000000 * tk ; tl = xr9 + tk ; t1x = ts + tx + ty = ts - tx ; t1v = xi9 + t1g ; t1h = xi9 - kp500000000 * t1g + tn = tg - tl ; tm = tg + tl ; t1y = t1u + t1v ; t1w = t1u - t1v + tB = tz - kp866025403 * tA ; tZ = tz + kp866025403 * tA + t10 = tC + kp866025403 * tD ; tE = tC - kp866025403 * tD + t1o = t1c - kp866025403 * t1d ; t1e = t1c + kp866025403 * t1d + t1l = tZ + t10 ; t11 = tZ - t10 ; t1j = t1h + kp866025403 * t1i + t1p = t1h - kp866025403 * t1i ; tK = tG - kp866025403 * tJ + t12 = tG + kp866025403 * tJ ; t13 = tL + kp866025403 * tO + tP = tL - kp866025403 * tO ; tT = tR - kp866025403 * tS + t15 = tR + kp866025403 * tS ; t1m = t1e + t1j ; t1k = t1e - t1j + t18 = t12 + t13 ; t14 = t12 - t13 ; t16 = tU + kp866025403 * tV + tW = tU - kp866025403 * tV ; t17 = t15 + t16 ; t19 = t15 - t16 + t1r = tB + tE ; tF = tB - tE ; t1s = t1o + t1p ; t1q = t1o - t1p + tY = tK + tP ; tQ = tK - tP ; tX = tT + tW ; t1n = tT - tW + r11 = (t19 + t1k) :+ (t11 - t14) + r10 = (tX - tY) :+ (t1r - t1s) + r9 = (t1t - t1w) :+ (tn + ty) + r8 = (t17 + t18) :+ (t1l + t1m) + r7 = (t1n + t1q) :+ (tF - tQ) + r6 = (tb - tm) :+ (t1x - t1y) + r5 = (t19 - t1k) :+ (t11 + t14) + r4 = (tX + tY) :+ (t1r + t1s) + r3 = (t1t + t1w) :+ (ty - tn) + r2 = (t17 - t18) :+ (t1l - t1m) + r1 = (t1n - t1q) :+ (tF + tQ) + MV.unsafeWrite xsout 0 $ (tb + tm) :+ (t1x + t1y) + MV.unsafeWrite xsout 1 $ if sign == 1 then r1 else r11 + MV.unsafeWrite xsout 2 $ if sign == 1 then r2 else r10 + MV.unsafeWrite xsout 3 $ if sign == 1 then r3 else r9 + MV.unsafeWrite xsout 4 $ if sign == 1 then r4 else r8 + MV.unsafeWrite xsout 5 $ if sign == 1 then r5 else r7 + MV.unsafeWrite xsout 6 $ if sign == 1 then r6 else r6 + MV.unsafeWrite xsout 7 $ if sign == 1 then r7 else r5 + MV.unsafeWrite xsout 8 $ if sign == 1 then r8 else r4 + MV.unsafeWrite xsout 9 $ if sign == 1 then r9 else r3 + MV.unsafeWrite xsout 10 $ if sign == 1 then r10 else r2 + MV.unsafeWrite xsout 11 $ if sign == 1 then r11 else r1 + +-- | Length 14 hard-coded FFT. +kp974927912, kp801937735, kp900968867 :: Double +kp554958132, kp692021471, kp356895867 :: Double +kp974927912 = 0.974927912181823607018131682993931217232785801 +kp801937735 = 0.801937735804838252472204639014890102331838324 +kp900968867 = 0.900968867902419126236102319507445051165919162 +kp554958132 = 0.554958132087371191422194871006410481067288862 +kp692021471 = 0.692021471630095869627814897002069140197260599 +kp356895867 = 0.356895867892209443894399510021300583399127187 +special14 :: Int -> MVCD s -> MVCD s -> ST s () +special14 sign xsin xsout = do + xr0 :+ xi0 <- MV.unsafeRead xsin 0 ; xr1 :+ xi1 <- MV.unsafeRead xsin 1 + xr2 :+ xi2 <- MV.unsafeRead xsin 2 ; xr3 :+ xi3 <- MV.unsafeRead xsin 3 + xr4 :+ xi4 <- MV.unsafeRead xsin 4 ; xr5 :+ xi5 <- MV.unsafeRead xsin 5 + xr6 :+ xi6 <- MV.unsafeRead xsin 6 ; xr7 :+ xi7 <- MV.unsafeRead xsin 7 + xr8 :+ xi8 <- MV.unsafeRead xsin 8 ; xr9 :+ xi9 <- MV.unsafeRead xsin 9 + xr10 :+ xi10 <- MV.unsafeRead xsin 10 ; xr11 :+ xi11 <- MV.unsafeRead xsin 11 + xr12 :+ xi12 <- MV.unsafeRead xsin 12 ; xr13 :+ xi13 <- MV.unsafeRead xsin 13 + let tp = xr0 + xr7 ; t3 = xr0 - xr7 ; t1x = xi0 + xi7 ; t1b = xi0 - xi7 + tq = xr2 + xr9 ; t6 = xr2 - xr9 ; tr = xr12 + xr5 ; t9 = xr12 - xr5 + tx = xr8 + xr1 ; tn = xr8 - xr1 ; tw = xr6 + xr13 ; tk = xr6 - xr13 + to = tk + tn ; t1i = tn - tk ; tu = xr10 + xr3 ; tg = xr10 - xr3 + tt = xr4 + xr11 ; td = xr4 - xr11 ; t1M = tr - tq ; ts = tq + tr + ta = t6 + t9 ; t1k = t9 - t6 ; t1L = tt - tu ; tv = tt + tu + th = td + tg ; t1j = tg - td ; t1K = tw - tx ; ty = tw + tx + tZ = to - kp356895867 * ta ; t14 = th - kp356895867 * to + tz = ta - kp356895867 * th ; t1Z = ty - kp356895867 * ts + t27 = ts - kp356895867 * tv ; t2c = tv - kp356895867 * ty + t1B = xi4 + xi11 ; tE = xi4 - xi11 ; t1C = xi10 + xi3 ; tH = xi10 - xi3 + t1F = xi8 + xi1 ; tV = xi8 - xi1 ; t1E = xi6 + xi13 ; tS = xi6 - xi13 + t1z = xi12 + xi5 ; tO = xi12 - xi5 ; t1d = tE + tH ; tI = tE - tH + t23 = t1F - t1E ; t1G = t1E + t1F ; t1D = t1B + t1C ; t24 = t1C - t1B + t1y = xi2 + xi9 ; tL = xi2 - xi9 ; tW = tS - tV ; t1e = tS + tV + t22 = t1y - t1z ; t1A = t1y + t1z ; tP = tL - tO ; t1c = tL + tO + t1n = t1e - kp356895867 * t1c ; t1s = t1c - kp356895867 * t1d + t1f = t1d - kp356895867 * t1e ; t1P = t1G - kp356895867 * t1A + t1U = t1A - kp356895867 * t1D ; t1H = t1D - kp356895867 * t1G + tA = to - kp692021471 * tz ; tX = tP + kp554958132 * tW + t1t = t1e - kp692021471 * t1s ; t1v = t1k + kp554958132 * t1i + tB = t3 - kp900968867 * tA ; tY = tI + kp801937735 * tX + t1u = t1b - kp900968867 * t1t ; t1w = t1j + kp801937735 * t1v + t10 = th - kp692021471 * tZ ; t11 = t3 - kp900968867 * t10 + t12 = tW + kp554958132 * tI ; t1o = t1d - kp692021471 * t1n + t1q = t1i + kp554958132 * t1j ; t15 = ta - kp692021471 * t14 + t13 = tP - kp801937735 * t12 ; t1p = t1b - kp900968867 * t1o + t1r = t1k - kp801937735 * t1q ; t16 = t3 - kp900968867 * t15 + t17 = tI - kp554958132 * tP ; t1g = t1c - kp692021471 * t1f + t1l = t1j - kp554958132 * t1k ; t1I = t1A - kp692021471 * t1H + t18 = tW - kp801937735 * t17 ; t1h = t1b - kp900968867 * t1g + t1m = t1i - kp801937735 * t1l ; t1J = t1x - kp900968867 * t1I + t1N = t1L + kp554958132 * t1M ; t2d = ts - kp692021471 * t2c + t2f = t24 + kp554958132 * t22 ; t1Q = t1D - kp692021471 * t1P + t1O = t1K - kp801937735 * t1N ; t2e = tp - kp900968867 * t2d + t2g = t23 - kp801937735 * t2f ; t1R = t1x - kp900968867 * t1Q + t1S = t1K + kp554958132 * t1L ; t20 = tv - kp692021471 * t1Z + t25 = t23 + kp554958132 * t24 ; t1V = t1G - kp692021471 * t1U + t1T = t1M + kp801937735 * t1S ; t21 = tp - kp900968867 * t20 + t26 = t22 + kp801937735 * t25 ; t1W = t1x - kp900968867 * t1V + t1X = t1M - kp554958132 * t1K ; t28 = ty - kp692021471 * t27 + t2a = t22 - kp554958132 * t23 ; t1Y = t1L - kp801937735 * t1X + t29 = tp - kp900968867 * t28 ; t2b = t24 - kp801937735 * t2a + r13 = (tB + kp974927912 * tY) :+ (t1u + kp974927912 * t1w) + r12 = (t21 + kp974927912 * t26) :+ (t1R + kp974927912 * t1T) + r11 = (t16 + kp974927912 * t18) :+ (t1h + kp974927912 * t1m) + r10 = (t2e + kp974927912 * t2g) :+ (t1J + kp974927912 * t1O) + r9 = (t11 - kp974927912 * t13) :+ (t1p - kp974927912 * t1r) + r8 = (t29 + kp974927912 * t2b) :+ (t1W + kp974927912 * t1Y) + r7 = (t3 + ta + th + to) :+ (t1b + t1c + t1d + t1e) + r6 = (t29 - kp974927912 * t2b) :+ (t1W - kp974927912 * t1Y) + r5 = (t11 + kp974927912 * t13) :+ (t1p + kp974927912 * t1r) + r4 = (t2e - kp974927912 * t2g) :+ (t1J - kp974927912 * t1O) + r3 = (t16 - kp974927912 * t18) :+ (t1h - kp974927912 * t1m) + r2 = (t21 - kp974927912 * t26) :+ (t1R - kp974927912 * t1T) + r1 = (tB - kp974927912 * tY) :+ (t1u - kp974927912 * t1w) + MV.unsafeWrite xsout 0 $ (tp + ts + tv + ty) :+ (t1x + t1A + t1D + t1G) + MV.unsafeWrite xsout 1 $ if sign == 1 then r1 else r13 + MV.unsafeWrite xsout 2 $ if sign == 1 then r2 else r12 + MV.unsafeWrite xsout 3 $ if sign == 1 then r3 else r11 + MV.unsafeWrite xsout 4 $ if sign == 1 then r4 else r10 + MV.unsafeWrite xsout 5 $ if sign == 1 then r5 else r9 + MV.unsafeWrite xsout 6 $ if sign == 1 then r6 else r8 + MV.unsafeWrite xsout 7 $ if sign == 1 then r7 else r7 + MV.unsafeWrite xsout 8 $ if sign == 1 then r8 else r6 + MV.unsafeWrite xsout 9 $ if sign == 1 then r9 else r5 + MV.unsafeWrite xsout 10 $ if sign == 1 then r10 else r4 + MV.unsafeWrite xsout 11 $ if sign == 1 then r11 else r3 + MV.unsafeWrite xsout 12 $ if sign == 1 then r12 else r2 + MV.unsafeWrite xsout 13 $ if sign == 1 then r13 else r1 + +-- | Length 15 hard-coded FFT. +--kp951056516, kp559016994, kp618033988 :: Double +--kp250000000, kp866025403, kp500000000 :: Double +--kp951056516 = 0.951056516295153572116439333379382143405698634 +--kp559016994 = 0.559016994374947424102293417182819058860154590 +--kp618033988 = 0.618033988749894848204586834365638117720309180 +--kp250000000 = 0.250000000000000000000000000000000000000000000 +--kp866025403 = 0.866025403784438646763723170752936183471402627 +--kp500000000 = 0.500000000000000000000000000000000000000000000 +special15 :: Int -> MVCD s -> MVCD s -> ST s () +special15 sign xsin xsout = do + xr0 :+ xi0 <- MV.unsafeRead xsin 0 ; xr1 :+ xi1 <- MV.unsafeRead xsin 1 + xr2 :+ xi2 <- MV.unsafeRead xsin 2 ; xr3 :+ xi3 <- MV.unsafeRead xsin 3 + xr4 :+ xi4 <- MV.unsafeRead xsin 4 ; xr5 :+ xi5 <- MV.unsafeRead xsin 5 + xr6 :+ xi6 <- MV.unsafeRead xsin 6 ; xr7 :+ xi7 <- MV.unsafeRead xsin 7 + xr8 :+ xi8 <- MV.unsafeRead xsin 8 ; xr9 :+ xi9 <- MV.unsafeRead xsin 9 + xr10 :+ xi10 <- MV.unsafeRead xsin 10 ; xr11 :+ xi11 <- MV.unsafeRead xsin 11 + xr12 :+ xi12 <- MV.unsafeRead xsin 12 ; xr13 :+ xi13 <- MV.unsafeRead xsin 13 + xr14 :+ xi14 <- MV.unsafeRead xsin 14 + let t1y = xr10 - xr5 ; t4 = xr5 + xr10 ; t1w = xi5 + xi10 ; tw = xi5 - xi10 + t5 = xr0 + t4 ; tt = xr0 - kp500000000 * t4 + t2l = xi0 + t1w ; t1x = xi0 - kp500000000 * t1w + tx = tt - kp866025403 * tw ; tV = tt + kp866025403 * tw + t1z = t1x + kp866025403 * t1y ; t1X = t1x - kp866025403 * t1y + tk = xr11 + xr1 ; t1k = xr1 - xr11 ; tM = xi11 - xi1 ; t1i = xi11 + xi1 + tJ = xr6 - kp500000000 * tk ; tl = xr6 + tk ; t2c = xi6 + t1i + t1j = xi6 - kp500000000 * t1i ; t1p = xr4 - xr14 ; tp = xr14 + xr4 + tN = tJ - kp866025403 * tM ; tZ = tJ + kp866025403 * tM + tO = xr9 - kp500000000 * tp ; tq = xr9 + tp ; t1n = xi14 + xi4 + tR = xi14 - xi4 ; t2s = tl - tq ; tr = tl + tq + t10 = tO + kp866025403 * tR ; tS = tO - kp866025403 * tR + t1o = xi9 - kp500000000 * t1n ; t2d = xi9 + t1n + t1O = t1j - kp866025403 * t1k ; t1l = t1j + kp866025403 * t1k + t24 = tN - tS ; tT = tN + tS ; t1P = t1o - kp866025403 * t1p + t1q = t1o + kp866025403 * t1p ; t2e = t2c - t2d ; t2n = t2c + t2d + t1Z = t1O + t1P ; t1Q = t1O - t1P ; t1r = t1l - t1q ; t1B = t1l + t1q + t11 = tZ + t10 ; t1H = tZ - t10 ; t9 = xr8 + xr13 ; t19 = xr13 - xr8 + tB = xi8 - xi13 ; t17 = xi8 + xi13 ; ty = xr3 - kp500000000 * t9 + ta = xr3 + t9 ; t2f = xi3 + t17 ; t18 = xi3 - kp500000000 * t17 + t1e = xr7 - xr2 ; te = xr2 + xr7 ; tC = ty - kp866025403 * tB + tW = ty + kp866025403 * tB ; tD = xr12 - kp500000000 * te + tf = xr12 + te ; t1c = xi2 + xi7 ; tG = xi2 - xi7 ; t2t = ta - tf + tg = ta + tf ; tX = tD + kp866025403 * tG + tH = tD - kp866025403 * tG ; t1d = xi12 - kp500000000 * t1c + t2g = xi12 + t1c ; t1R = t18 - kp866025403 * t19 + t1a = t18 + kp866025403 * t19 ; t25 = tC - tH ; tI = tC + tH + t1S = t1d - kp866025403 * t1e ; t1f = t1d + kp866025403 * t1e + t2h = t2f - t2g ; t2m = t2f + t2g ; t1Y = t1R + t1S ; t1T = t1R - t1S + t1g = t1a - t1f ; t1A = t1a + t1f ; t2a = tg - tr ; ts = tg + tr + tY = tW + tX ; t1G = tW - tX ; t29 = t5 - kp250000000 * ts + t2o = t2m + t2n ; t2q = t2m - t2n ; t2k = t2h + kp618033988 * t2e + t2i = t2e - kp618033988 * t2h ; t2b = t29 - kp559016994 * t2a + t2j = t29 + kp559016994 * t2a ; t2p = t2l - kp250000000 * t2o + tU = tI + tT ; t1M = tI - tT ; t2r = t2p - kp559016994 * t2q + t2v = t2p + kp559016994 * t2q ; t2w = t2t + kp618033988 * t2s + t2u = t2s - kp618033988 * t2t ; t1L = tx - kp250000000 * tU + t20 = t1Y + t1Z ; t22 = t1Y - t1Z ; t1N = t1L - kp559016994 * t1M + t1V = t1L + kp559016994 * t1M ; t1W = t1T + kp618033988 * t1Q + t1U = t1Q - kp618033988 * t1T ; t21 = t1X - kp250000000 * t20 + t1C = t1A + t1B ; t1E = t1A - t1B ; t23 = t21 - kp559016994 * t22 + t27 = t21 + kp559016994 * t22 ; t28 = t25 + kp618033988 * t24 + t26 = t24 - kp618033988 * t25 ; t1D = t1z - kp250000000 * t1C + t12 = tY + t11 ; t14 = tY - t11 ; t1F = t1D + kp559016994 * t1E + t1J = t1D - kp559016994 * t1E ; t1K = t1H - kp618033988 * t1G + t1I = t1G + kp618033988 * t1H ; t13 = tV - kp250000000 * t12 + t1t = t13 - kp559016994 * t14 ; t15 = t13 + kp559016994 * t14 + t1s = t1g + kp618033988 * t1r ; t1u = t1r - kp618033988 * t1g + r14 = (t15 + kp951056516 * t1s) :+ (t1F - kp951056516 * t1I) + r13 = (t1N - kp951056516 * t1U) :+ (t23 + kp951056516 * t26) + r12 = (t2b + kp951056516 * t2i) :+ (t2r - kp951056516 * t2u) + r11 = (t15 - kp951056516 * t1s) :+ (t1F + kp951056516 * t1I) + r10 = (tx + tU) :+ (t1X + t20) + r9 = (t2j + kp951056516 * t2k) :+ (t2v - kp951056516 * t2w) + r8 = (t1t - kp951056516 * t1u) :+ (t1J + kp951056516 * t1K) + r7 = (t1N + kp951056516 * t1U) :+ (t23 - kp951056516 * t26) + r6 = (t2j - kp951056516 * t2k) :+ (t2v + kp951056516 * t2w) + r5 = (tV + t12) :+ (t1z + t1C) + r4 = (t1V + kp951056516 * t1W) :+ (t27 - kp951056516 * t28) + r3 = (t2b - kp951056516 * t2i) :+ (t2r + kp951056516 * t2u) + r2 = (t1t + kp951056516 * t1u) :+ (t1J - kp951056516 * t1K) + r1 = (t1V - kp951056516 * t1W) :+ (t27 + kp951056516 * t28) + MV.unsafeWrite xsout 0 $ (t5 + ts) :+ (t2l + t2o) + MV.unsafeWrite xsout 1 $ if sign == 1 then r1 else r14 + MV.unsafeWrite xsout 2 $ if sign == 1 then r2 else r13 + MV.unsafeWrite xsout 3 $ if sign == 1 then r3 else r12 + MV.unsafeWrite xsout 4 $ if sign == 1 then r4 else r11 + MV.unsafeWrite xsout 5 $ if sign == 1 then r5 else r10 + MV.unsafeWrite xsout 6 $ if sign == 1 then r6 else r9 + MV.unsafeWrite xsout 7 $ if sign == 1 then r7 else r8 + MV.unsafeWrite xsout 8 $ if sign == 1 then r8 else r7 + MV.unsafeWrite xsout 9 $ if sign == 1 then r9 else r6 + MV.unsafeWrite xsout 10 $ if sign == 1 then r10 else r5 + MV.unsafeWrite xsout 11 $ if sign == 1 then r11 else r4 + MV.unsafeWrite xsout 12 $ if sign == 1 then r12 else r3 + MV.unsafeWrite xsout 13 $ if sign == 1 then r13 else r2 + MV.unsafeWrite xsout 14 $ if sign == 1 then r14 else r1 + +-- | Length 20 hard-coded FFT. +--kp951056516, kp559016994, kp618033988, kp250000000 :: Double +--kp951056516 = 0.951056516295153572116439333379382143405698634 +--kp559016994 = 0.559016994374947424102293417182819058860154590 +--kp618033988 = 0.618033988749894848204586834365638117720309180 +--kp250000000 = 0.250000000000000000000000000000000000000000000 +special20 :: Int -> MVCD s -> MVCD s -> ST s () +special20 sign xsin xsout = do + xr0 :+ xi0 <- MV.unsafeRead xsin 0 ; xr1 :+ xi1 <- MV.unsafeRead xsin 1 + xr2 :+ xi2 <- MV.unsafeRead xsin 2 ; xr3 :+ xi3 <- MV.unsafeRead xsin 3 + xr4 :+ xi4 <- MV.unsafeRead xsin 4 ; xr5 :+ xi5 <- MV.unsafeRead xsin 5 + xr6 :+ xi6 <- MV.unsafeRead xsin 6 ; xr7 :+ xi7 <- MV.unsafeRead xsin 7 + xr8 :+ xi8 <- MV.unsafeRead xsin 8 ; xr9 :+ xi9 <- MV.unsafeRead xsin 9 + xr10 :+ xi10 <- MV.unsafeRead xsin 10 ; xr11 :+ xi11 <- MV.unsafeRead xsin 11 + xr12 :+ xi12 <- MV.unsafeRead xsin 12 ; xr13 :+ xi13 <- MV.unsafeRead xsin 13 + xr14 :+ xi14 <- MV.unsafeRead xsin 14 ; xr15 :+ xi15 <- MV.unsafeRead xsin 15 + xr16 :+ xi16 <- MV.unsafeRead xsin 16 ; xr17 :+ xi17 <- MV.unsafeRead xsin 17 + xr18 :+ xi18 <- MV.unsafeRead xsin 18 ; xr19 :+ xi19 <- MV.unsafeRead xsin 19 + let t1N = xr0 - xr10 ; t3 = xr0 + xr10 ; t2L = xi0 + xi10 ; tN = xi0 - xi10 + tO = xr5 - xr15 ; t6 = xr5 + xr15 ; t2M = xi5 + xi15 ; t1Q = xi5 - xi15 + t1d = tO + tN ; tP = tN - tO ; tD = t3 + t6 ; t7 = t3 - t6 + t3b = t2L + t2M ; t2N = t2L - t2M ; t2f = t1N + t1Q ; t1R = t1N - t1Q + t1o = xr8 - xr18 ; tp = xr8 + xr18 ; t2u = xi8 + xi18 ; t13 = xi8 - xi18 + t14 = xr13 - xr3 ; ts = xr13 + xr3 ; t2v = xi13 + xi3 ; t1r = xi13 - xi3 + t1t = xr12 - xr2 ; tw = xr12 + xr2 ; t2x = xi12 + xi2 ; t18 = xi12 - xi2 + tH = tp + ts ; tt = tp - ts ; t19 = xr17 - xr7 ; tz = xr17 + xr7 + t2y = xi17 + xi7 ; t1w = xi17 - xi7 ; t2w = t2u - t2v ; t35 = t2u + t2v + tI = tw + tz ; tA = tw - tz ; t2z = t2x - t2y ; t36 = t2x + t2y + t2U = tt - tA ; tB = tt + tA ; t2P = t2w + t2z ; t2A = t2w - t2z + t3d = t35 + t36 ; t37 = t35 - t36 ; t15 = t13 - t14 ; t1h = t14 + t13 + t1i = t19 + t18 ; t1a = t18 - t19 ; t1s = t1o - t1r ; t29 = t1o + t1r + t3j = tH - tI ; tJ = tH + tI ; t1x = t1t - t1w ; t2a = t1t + t1w + t2n = t15 - t1a ; t1b = t15 + t1a ; t1T = t1s + t1x ; t1y = t1s - t1x + t2b = t29 - t2a ; t2h = t29 + t2a ; t1j = t1h + t1i ; t1Y = t1h - t1i + ta = xr4 + xr14 ; t1z = xr4 - xr14 ; t2B = xi4 + xi14 ; tS = xi4 - xi14 + tT = xr9 - xr19 ; td = xr9 + xr19 ; t2C = xi9 + xi19 ; t1C = xi9 - xi19 + t1E = xr16 - xr6 ; th = xr16 + xr6 ; t2E = xi16 + xi6 ; tX = xi16 - xi6 + tE = ta + td ; te = ta - td ; tY = xr1 - xr11 ; tk = xr1 + xr11 + t2F = xi1 + xi11 ; t1H = xi1 - xi11 ; t2D = t2B - t2C ; t32 = t2B + t2C + tF = th + tk ; tl = th - tk ; t2G = t2E - t2F ; t33 = t2E + t2F + t2V = te - tl ; tm = te + tl ; t2O = t2D + t2G ; t2H = t2D - t2G + t3c = t32 + t33 ; t34 = t32 - t33 ; tU = tS - tT ; t1e = tT + tS + t1f = tY + tX ; tZ = tX - tY ; t1D = t1z - t1C ; t26 = t1z + t1C + t3i = tE - tF ; tG = tE + tF ; t1I = t1E - t1H ; t27 = t1E + t1H + t2m = tU - tZ ; t10 = tU + tZ ; t1S = t1D + t1I ; t1J = t1D - t1I + t28 = t26 - t27 ; t2g = t26 + t27 ; t2s = tm - tB ; tC = tm + tB + t1g = t1e + t1f ; t1Z = t1e - t1f ; t2r = t7 - kp250000000 * tC + t2Q = t2O + t2P ; t2S = t2O - t2P ; t2K = t2H + kp618033988 * t2A + t2I = t2A - kp618033988 * t2H ; t2t = t2r - kp559016994 * t2s + t2J = t2r + kp559016994 * t2s ; t2R = t2N - kp250000000 * t2Q + tK = tG + tJ ; t30 = tG - tJ ; t2T = t2R - kp559016994 * t2S + t2X = t2R + kp559016994 * t2S ; t2Y = t2V + kp618033988 * t2U + t2W = t2U - kp618033988 * t2V ; t2Z = tD - kp250000000 * tK + t3e = t3c + t3d ; t3g = t3c - t3d ; t31 = t2Z + kp559016994 * t30 + t39 = t2Z - kp559016994 * t30 ; t3a = t37 - kp618033988 * t34 + t38 = t34 + kp618033988 * t37 ; t3f = t3b - kp250000000 * t3e + t1c = t10 + t1b ; t24 = t10 - t1b ; t3h = t3f + kp559016994 * t3g + t3l = t3f - kp559016994 * t3g ; t3m = t3j - kp618033988 * t3i + t3k = t3i + kp618033988 * t3j ; t23 = tP - kp250000000 * t1c + t2i = t2g + t2h ; t2k = t2g - t2h ; t25 = t23 + kp559016994 * t24 + t2d = t23 - kp559016994 * t24 ; t2e = t2b - kp618033988 * t28 + t2c = t28 + kp618033988 * t2b ; t2j = t2f - kp250000000 * t2i + t1k = t1g + t1j ; t1m = t1g - t1j ; t2l = t2j + kp559016994 * t2k + t2p = t2j - kp559016994 * t2k ; t2q = t2n - kp618033988 * t2m + t2o = t2m + kp618033988 * t2n ; t1l = t1d - kp250000000 * t1k + t1U = t1S + t1T ; t1W = t1S - t1T ; t1n = t1l - kp559016994 * t1m + t1L = t1l + kp559016994 * t1m ; t1M = t1J + kp618033988 * t1y + t1K = t1y - kp618033988 * t1J ; t1V = t1R - kp250000000 * t1U + t21 = t1V + kp559016994 * t1W ; t1X = t1V - kp559016994 * t1W + t20 = t1Y - kp618033988 * t1Z ; t22 = t1Z + kp618033988 * t1Y + r19 = (t2l + kp951056516 * t2o) :+ (t25 - kp951056516 * t2c) + r18 = (t2t - kp951056516 * t2I) :+ (t2T + kp951056516 * t2W) + r17 = (t1X + kp951056516 * t20) :+ (t1n - kp951056516 * t1K) + r16 = (t31 - kp951056516 * t38) :+ (t3h + kp951056516 * t3k) + r15 = (t2f + t2i) :+ (tP + t1c) + r14 = (t2J + kp951056516 * t2K) :+ (t2X - kp951056516 * t2Y) + r13 = (t1X - kp951056516 * t20) :+ (t1n + kp951056516 * t1K) + r12 = (t39 + kp951056516 * t3a) :+ (t3l - kp951056516 * t3m) + r11 = (t2l - kp951056516 * t2o) :+ (t25 + kp951056516 * t2c) + r10 = (t7 + tC) :+ (t2N + t2Q) + r9 = (t21 + kp951056516 * t22) :+ (t1L - kp951056516 * t1M) + r8 = (t39 - kp951056516 * t3a) :+ (t3l + kp951056516 * t3m) + r7 = (t2p + kp951056516 * t2q) :+ (t2d - kp951056516 * t2e) + r6 = (t2J - kp951056516 * t2K) :+ (t2X + kp951056516 * t2Y) + r5 = (t1R + t1U) :+ (t1d + t1k) + r4 = (t31 + kp951056516 * t38) :+ (t3h - kp951056516 * t3k) + r3 = (t2p - kp951056516 * t2q) :+ (t2d + kp951056516 * t2e) + r2 = (t2t + kp951056516 * t2I) :+ (t2T - kp951056516 * t2W) + r1 = (t21 - kp951056516 * t22) :+ (t1L + kp951056516 * t1M) + MV.unsafeWrite xsout 0 $ (tD + tK) :+ (t3b + t3e) + MV.unsafeWrite xsout 1 $ if sign == 1 then r1 else r19 + MV.unsafeWrite xsout 2 $ if sign == 1 then r2 else r18 + MV.unsafeWrite xsout 3 $ if sign == 1 then r3 else r17 + MV.unsafeWrite xsout 4 $ if sign == 1 then r4 else r16 + MV.unsafeWrite xsout 5 $ if sign == 1 then r5 else r15 + MV.unsafeWrite xsout 6 $ if sign == 1 then r6 else r14 + MV.unsafeWrite xsout 7 $ if sign == 1 then r7 else r13 + MV.unsafeWrite xsout 8 $ if sign == 1 then r8 else r12 + MV.unsafeWrite xsout 9 $ if sign == 1 then r9 else r11 + MV.unsafeWrite xsout 10 $ if sign == 1 then r10 else r10 + MV.unsafeWrite xsout 11 $ if sign == 1 then r11 else r9 + MV.unsafeWrite xsout 12 $ if sign == 1 then r12 else r8 + MV.unsafeWrite xsout 13 $ if sign == 1 then r13 else r7 + MV.unsafeWrite xsout 14 $ if sign == 1 then r14 else r6 + MV.unsafeWrite xsout 15 $ if sign == 1 then r15 else r5 + MV.unsafeWrite xsout 16 $ if sign == 1 then r16 else r4 + MV.unsafeWrite xsout 17 $ if sign == 1 then r17 else r3 + MV.unsafeWrite xsout 18 $ if sign == 1 then r18 else r2 + MV.unsafeWrite xsout 19 $ if sign == 1 then r19 else r1 + +-- | Length 25 hard-coded FFT. +kp803003575, kp554608978, kp248028675, kp726211448 :: Double +kp525970792, kp992114701, kp851038619, kp912575812 :: Double +kp912018591, kp943557151, kp614372930, kp621716863 :: Double +kp994076283, kp734762448, kp772036680, kp126329378 :: Double +kp827271945, kp949179823, kp860541664, kp557913902 :: Double +kp249506682, kp681693190, kp560319534, kp998026728 :: Double +kp906616052, kp968479752, kp845997307, kp470564281 :: Double +kp062914667, kp921177326, kp833417178, kp541454447 :: Double +kp242145790, kp683113946, kp559154169, kp968583161 :: Double +kp904730450, kp831864738, kp871714437, kp939062505 :: Double +kp549754652, kp634619297, kp256756360 :: Double +--kp951056516, kp559016994, kp250000000, kp618033988 :: Double +kp803003575 = 0.803003575438660414833440593570376004635464850 +kp554608978 = 0.554608978404018097464974850792216217022558774 +kp248028675 = 0.248028675328619457762448260696444630363259177 +kp726211448 = 0.726211448929902658173535992263577167607493062 +kp525970792 = 0.525970792408939708442463226536226366643874659 +kp992114701 = 0.992114701314477831049793042785778521453036709 +kp851038619 = 0.851038619207379630836264138867114231259902550 +kp912575812 = 0.912575812670962425556968549836277086778922727 +kp912018591 = 0.912018591466481957908415381764119056233607330 +kp943557151 = 0.943557151597354104399655195398983005179443399 +kp614372930 = 0.614372930789563808870829930444362096004872855 +kp621716863 = 0.621716863012209892444754556304102309693593202 +kp994076283 = 0.994076283785401014123185814696322018529298887 +kp734762448 = 0.734762448793050413546343770063151342619912334 +kp772036680 = 0.772036680810363904029489473607579825330539880 +kp126329378 = 0.126329378446108174786050455341811215027378105 +kp827271945 = 0.827271945972475634034355757144307982555673741 +kp949179823 = 0.949179823508441261575555465843363271711583843 +kp860541664 = 0.860541664367944677098261680920518816412804187 +kp557913902 = 0.557913902031834264187699648465567037992437152 +kp249506682 = 0.249506682107067890488084201715862638334226305 +kp681693190 = 0.681693190061530575150324149145440022633095390 +kp560319534 = 0.560319534973832390111614715371676131169633784 +kp998026728 = 0.998026728428271561952336806863450553336905220 +kp906616052 = 0.906616052148196230441134447086066874408359177 +kp968479752 = 0.968479752739016373193524836781420152702090879 +kp845997307 = 0.845997307939530944175097360758058292389769300 +kp470564281 = 0.470564281212251493087595091036643380879947982 +kp062914667 = 0.062914667253649757225485955897349402364686947 +kp921177326 = 0.921177326965143320250447435415066029359282231 +kp833417178 = 0.833417178328688677408962550243238843138996060 +kp541454447 = 0.541454447536312777046285590082819509052033189 +kp242145790 = 0.242145790282157779872542093866183953459003101 +kp683113946 = 0.683113946453479238701949862233725244439656928 +kp559154169 = 0.559154169276087864842202529084232643714075927 +kp968583161 = 0.968583161128631119490168375464735813836012403 +kp904730450 = 0.904730450839922351881287709692877908104763647 +kp831864738 = 0.831864738706457140726048799369896829771167132 +kp871714437 = 0.871714437527667770979999223229522602943903653 +kp939062505 = 0.939062505817492352556001843133229685779824606 +kp549754652 = 0.549754652192770074288023275540779861653779767 +kp634619297 = 0.634619297544148100711287640319130485732531031 +kp256756360 = 0.256756360367726783319498520922669048172391148 +--kp951056516 = 0.951056516295153572116439333379382143405698634 +--kp559016994 = 0.559016994374947424102293417182819058860154590 +--kp250000000 = 0.250000000000000000000000000000000000000000000 +--kp618033988 = 0.618033988749894848204586834365638117720309180 +special25 :: Int -> MVCD s -> MVCD s -> ST s () +special25 sign xsin xsout = do + xr0 :+ xi0 <- MV.unsafeRead xsin 0 ; xr1 :+ xi1 <- MV.unsafeRead xsin 1 + xr2 :+ xi2 <- MV.unsafeRead xsin 2 ; xr3 :+ xi3 <- MV.unsafeRead xsin 3 + xr4 :+ xi4 <- MV.unsafeRead xsin 4 ; xr5 :+ xi5 <- MV.unsafeRead xsin 5 + xr6 :+ xi6 <- MV.unsafeRead xsin 6 ; xr7 :+ xi7 <- MV.unsafeRead xsin 7 + xr8 :+ xi8 <- MV.unsafeRead xsin 8 ; xr9 :+ xi9 <- MV.unsafeRead xsin 9 + xr10 :+ xi10 <- MV.unsafeRead xsin 10 ; xr11 :+ xi11 <- MV.unsafeRead xsin 11 + xr12 :+ xi12 <- MV.unsafeRead xsin 12 ; xr13 :+ xi13 <- MV.unsafeRead xsin 13 + xr14 :+ xi14 <- MV.unsafeRead xsin 14 ; xr15 :+ xi15 <- MV.unsafeRead xsin 15 + xr16 :+ xi16 <- MV.unsafeRead xsin 16 ; xr17 :+ xi17 <- MV.unsafeRead xsin 17 + xr18 :+ xi18 <- MV.unsafeRead xsin 18 ; xr19 :+ xi19 <- MV.unsafeRead xsin 19 + xr20 :+ xi20 <- MV.unsafeRead xsin 20 ; xr21 :+ xi21 <- MV.unsafeRead xsin 21 + xr22 :+ xi22 <- MV.unsafeRead xsin 22 ; xr23 :+ xi23 <- MV.unsafeRead xsin 23 + xr24 :+ xi24 <- MV.unsafeRead xsin 24 + let t4 = xr5+xr20 ; t1S = xr5-xr20 ; t7 = xr10+xr15 ; t1T = xr10-xr15 + t4Q = t1T-kp618033988*t1S ; t1U = t1S+kp618033988*t1T + t8 = t4+t7 ; t3a = t4-t7 ; t3c = xi5-xi20 ; t1y = xi5+xi20 + t39 = xr0-kp250000000*t8 ; t9 = xr0+t8 + t1B = xi10+xi15 ; t3d = xi10-xi15 + t3b = t39+kp559016994*t3a ; t45 = t39-kp559016994*t3a + t3e = t3c+kp618033988*t3d ; t46 = t3d-kp618033988*t3c + t1C = t1y+t1B ; t1Q = t1y-t1B ; t1P = xi0-kp250000000*t1C + t1D = xi0+t1C ; t4P = t1P-kp559016994*t1Q + t1R = t1P+kp559016994*t1Q ; t1Z = xr21-xr6 ; td = xr6+xr21 + t20 = xr16-xr11 ; tg = xr11+xr16 ; th = td+tg ; t24 = td-tg + t26 = xi6-xi21 ; tT = xi6+xi21 ; tW = xi11+xi16 ; t27 = xi16-xi11 + t1X = tT-tW ; tX = tT+tW ; t2l = xr24-xr9 ; tm = xr9+xr24 + t2m = xr19-xr14 ; tp = xr14+xr19 ; tq = tm+tp ; t2c = tm-tp + t2e = xi24-xi9 ; t12 = xi9+xi24 ; t15 = xi14+xi19 ; t2f = xi19-xi14 + t23 = xr1-kp250000000*th ; ti = xr1+th ; t2j = t15-t12 + t16 = t12+t15 ; tr = xr4+tq ; t2b = kp250000000 * tq - xr4 + t1W = xi1-kp250000000*tX ; tY = xi1+tX + t21 = t1Z+kp618033988*t20 ; t4y = t20-kp618033988*t1Z + t2i = xi4-kp250000000*t16 ; t17 = xi4+t16 ; ts = ti+tr + t1K = ti-tr ; t18 = tY-t17 ; t1E = tY+t17 + t2n = t2l+kp618033988*t2m ; t4r = t2m-kp618033988*t2l + t4x = t1W-kp559016994*t1X ; t1Y = t1W+kp559016994*t1X + t4o = t2f-kp618033988*t2e ; t2g = t2e+kp618033988*t2f + t4z = t4x+kp951056516*t4y ; t5f = t4x-kp951056516*t4y + t3z = t1Y-kp951056516*t21 ; t22 = t1Y+kp951056516*t21 + t4q = t2i+kp559016994*t2j ; t2k = t2i-kp559016994*t2j + t4s = t4q+kp951056516*t4r ; t5b = t4q-kp951056516*t4r + t3C = t2k-kp951056516*t2n ; t2o = t2k+kp951056516*t2n + t2d = t2b-kp559016994*t2c ; t4n = t2b+kp559016994*t2c + t28 = t26-kp618033988*t27 ; t4v = t27+kp618033988*t26 + t3D = t2d-kp951056516*t2g ; t2h = t2d+kp951056516*t2g + t4p = t4n+kp951056516*t4o ; t5c = t4n-kp951056516*t4o + t4u = t23-kp559016994*t24 ; t25 = t23+kp559016994*t24 + t4w = t4u-kp951056516*t4v ; t5e = t4u+kp951056516*t4v + t3A = t25-kp951056516*t28 ; t29 = t25+kp951056516*t28 + t2u = xr22-xr7 ; tw = xr7+xr22 ; t2v = xr17-xr12 ; tz = xr12+xr17 + tA = tw+tz ; t2z = tz-tw ; t2B = xi22-xi7 ; t1c = xi7+xi22 + t1f = xi12+xi17 ; t2C = xi12-xi17 ; t2s = t1f-t1c ; t1g = t1c+t1f + t2J = xr8-xr23 ; tF = xr8+xr23 ; t2K = xr13-xr18 ; tI = xr13+xr18 + tJ = tF+tI ; t2O = tI-tF ; t2Q = xi23-xi8 ; t1l = xi8+xi23 + t1o = xi13+xi18 ; t2R = xi18-xi13 ; t2y = xr2-kp250000000*tA + tB = xr2+tA ; t2H = t1o-t1l ; t1p = t1l+t1o ; tK = xr3+tJ + t2N = xr3-kp250000000*tJ ; t2r = xi2-kp250000000*t1g + t1h = xi2+t1g ; t2w = t2u+kp618033988*t2v + t49 = t2v-kp618033988*t2u ; t2G = xi3-kp250000000*t1p + t1q = xi3+t1p ; tL = tB+tK ; t1L = tB-tK ; t1r = t1h-t1q + t1F = t1h+t1q ; t2S = t2Q+kp618033988*t2R + t4j = t2R-kp618033988*t2Q ; t48 = t2r+kp559016994*t2s + t2t = t2r-kp559016994*t2s ; t4g = t2K-kp618033988*t2J + t2L = t2J+kp618033988*t2K ; t4a = t48+kp951056516*t49 + t57 = t48-kp951056516*t49 ; t3v = t2t-kp951056516*t2w + t2x = t2t+kp951056516*t2w ; t4i = t2N+kp559016994*t2O + t2P = t2N-kp559016994*t2O ; t4k = t4i-kp951056516*t4j + t55 = t4i+kp951056516*t4j ; t3s = t2P+kp951056516*t2S + t2T = t2P-kp951056516*t2S ; t2I = t2G-kp559016994*t2H + t4f = t2G+kp559016994*t2H ; t2D = t2B-kp618033988*t2C + t4c = t2C+kp618033988*t2B ; t3t = t2I+kp951056516*t2L + t2M = t2I-kp951056516*t2L ; t4h = t4f-kp951056516*t4g + t54 = t4f+kp951056516*t4g ; tM = ts+tL ; tO = ts-tL + t4b = t2y+kp559016994*t2z ; t2A = t2y-kp559016994*t2z + tN = t9-kp250000000*tM ; t4d = t4b+kp951056516*t4c + t58 = t4b-kp951056516*t4c ; t3w = t2A+kp951056516*t2D + t2E = t2A-kp951056516*t2D ; t1s = t18+kp618033988*t1r + t1u = t1r-kp618033988*t18 ; tP = tN+kp559016994*tO + t1t = tN-kp559016994*tO ; t1G = t1E+t1F ; t1I = t1E-t1F + t1H = t1D-kp250000000*t1G ; t1J = t1H+kp559016994*t1I + t1N = t1H-kp559016994*t1I ; t1M = t1K+kp618033988*t1L + t1O = t1L-kp618033988*t1K ; t3H = t1R+kp951056516*t1U + t1V = t1R-kp951056516*t1U ; t3f = t3b+kp951056516*t3e + t3r = t3b-kp951056516*t3e ; t30 = t29+kp256756360*t22 + t2a = t22-kp256756360*t29 ; t2p = t2h+kp634619297*t2o + t31 = t2o-kp634619297*t2h ; t33 = t2E+kp549754652*t2x + t2F = t2x-kp549754652*t2E ; t2U = t2M-kp939062505*t2T + t34 = t2T+kp939062505*t2M ; t3m = t2a-kp871714437*t2p + t2q = t2a+kp871714437*t2p ; t3n = t2F-kp831864738*t2U + t2V = t2F+kp831864738*t2U ; t2W = t2q+kp904730450*t2V + t2Y = t2q-kp904730450*t2V ; t32 = t30-kp871714437*t31 + t3g = t30+kp871714437*t31 ; t3h = t33+kp831864738*t34 + t35 = t33-kp831864738*t34 ; t3i = t3g+kp904730450*t3h + t3k = t3g-kp904730450*t3h ; t36 = t32+kp559154169*t35 + t38 = t35-kp683113946*t32 ; t2X = t1V-kp242145790*t2W + t3o = t3m+kp559154169*t3n ; t3q = t3n-kp683113946*t3m + t3j = t3f-kp242145790*t3i ; t2Z = t2X+kp541454447*t2Y + t37 = t2X-kp541454447*t2Y ; t47 = t45+kp951056516*t46 + t53 = t45-kp951056516*t46 ; t3p = t3j-kp541454447*t3k + t3l = t3j+kp541454447*t3k ; t5j = t4P+kp951056516*t4Q + t4R = t4P-kp951056516*t4Q ; t5k = t55-kp062914667*t54 + t56 = t54+kp062914667*t55 ; t59 = t57+kp634619297*t58 + t5l = t58-kp634619297*t57 ; t5n = t5c-kp470564281*t5b + t5d = t5b+kp470564281*t5c ; t5g = t5e+kp549754652*t5f + t5o = t5f-kp549754652*t5e ; t5u = t56-kp845997307*t59 + t5a = t56+kp845997307*t59 ; t5v = t5d-kp968479752*t5g + t5h = t5d+kp968479752*t5g ; t5i = t5a+kp906616052*t5h + t5A = t5a-kp906616052*t5h ; t5D = t5k-kp845997307*t5l + t5m = t5k+kp845997307*t5l ; t5p = t5n+kp968479752*t5o + t5C = t5n-kp968479752*t5o ; t5s = t5m+kp906616052*t5p + t5q = t5m-kp906616052*t5p ; t5w = t5u-kp560319534*t5v + t5y = t5v+kp681693190*t5u ; t5E = t5C-kp681693190*t5D + t5G = t5D+kp560319534*t5C ; t5r = t5j+kp249506682*t5q + t5z = t53-kp249506682*t5i ; t5t = t5r-kp557913902*t5s + t5x = t5r+kp557913902*t5s ; t5F = t5z+kp557913902*t5A + t5B = t5z-kp557913902*t5A ; t4J = t4d-kp062914667*t4a + t4e = t4a+kp062914667*t4d ; t4l = t4h-kp827271945*t4k + t4K = t4k+kp827271945*t4h ; t4G = t4s-kp126329378*t4p + t4t = t4p+kp126329378*t4s ; t4A = t4w+kp939062505*t4z + t4H = t4z-kp939062505*t4w ; t4Y = t4e-kp772036680*t4l + t4m = t4e+kp772036680*t4l ; t4Z = t4t-kp734762448*t4A + t4B = t4t+kp734762448*t4A ; t4C = t4m+kp994076283*t4B + t4E = t4m-kp994076283*t4B ; t4I = t4G+kp734762448*t4H + t4T = t4G-kp734762448*t4H ; t4S = t4J+kp772036680*t4K + t4L = t4J-kp772036680*t4K ; t4U = t4S+kp994076283*t4T + t4W = t4S-kp994076283*t4T ; t4M = t4I-kp621716863*t4L + t4O = t4L+kp614372930*t4I ; t4D = t47-kp249506682*t4C + t50 = t4Y+kp614372930*t4Z ; t52 = t4Z-kp621716863*t4Y + t4V = t4R+kp249506682*t4U ; t4F = t4D-kp557913902*t4E + t4N = t4D+kp557913902*t4E ; t51 = t4V+kp557913902*t4W + t4X = t4V-kp557913902*t4W ; t3I = t3t+kp126329378*t3s + t3u = t3s-kp126329378*t3t ; t3x = t3v-kp470564281*t3w + t3J = t3w+kp470564281*t3v ; t3L = t3A-kp634619297*t3z + t3B = t3z+kp634619297*t3A ; t3E = t3C-kp827271945*t3D + t3M = t3D+kp827271945*t3C ; t3S = t3u+kp912018591*t3x + t3y = t3u-kp912018591*t3x ; t3T = t3B+kp912575812*t3E + t3F = t3B-kp912575812*t3E ; t3G = t3y-kp851038619*t3F + t3Y = t3y+kp851038619*t3F ; t41 = t3I-kp912018591*t3J + t3K = t3I+kp912018591*t3J ; t3N = t3L+kp912575812*t3M + t40 = t3L-kp912575812*t3M ; t3Q = t3K-kp851038619*t3N + t3O = t3K+kp851038619*t3N ; t3U = t3S-kp525970792*t3T + t3W = t3T+kp726211448*t3S ; t42 = t40-kp726211448*t41 + t44 = t41+kp525970792*t40 ; t3P = t3H+kp248028675*t3O + t3X = t3r+kp248028675*t3G ; t3R = t3P-kp554608978*t3Q + t3V = t3P+kp554608978*t3Q ; t3Z = t3X+kp554608978*t3Y + t43 = t3X-kp554608978*t3Y + r24 = (t3f+kp968583161*t3i) :+ (t1V+kp968583161*t2W) + r23 = (t53+kp998026728*t5i) :+ (t5j-kp998026728*t5q) + r22 = (t47+kp998026728*t4C) :+ (t4R-kp998026728*t4U) + r21 = (t3r-kp992114701*t3G) :+ (t3H-kp992114701*t3O) + r20 = (tP+kp951056516*t1s) :+ (t1J-kp951056516*t1M) + r19 = (t3l+kp921177326*t3o) :+ (t2Z-kp921177326*t36) + r18 = (t5B-kp860541664*t5E) :+ (t5x+kp860541664*t5y) + r17 = (t4F+kp943557151*t4M) :+ (t51+kp943557151*t52) + r16 = (t3Z-kp803003575*t42) :+ (t3V-kp803003575*t3W) + r15 = (t1t-kp951056516*t1u) :+ (t1N+kp951056516*t1O) + r14 = (t3p-kp833417178*t3q) :+ (t37+kp833417178*t38) + r13 = (t5F-kp949179823*t5G) :+ (t5t-kp949179823*t5w) + r12 = (t4N+kp949179823*t4O) :+ (t4X+kp949179823*t50) + r11 = (t43-kp943557151*t44) :+ (t3R+kp943557151*t3U) + r10 = (t1t+kp951056516*t1u) :+ (t1N-kp951056516*t1O) + r9 = (t3p+kp833417178*t3q) :+ (t37-kp833417178*t38) + r8 = (t5F+kp949179823*t5G) :+ (t5t+kp949179823*t5w) + r7 = (t4N-kp949179823*t4O) :+ (t4X-kp949179823*t50) + r6 = (t43+kp943557151*t44) :+ (t3R-kp943557151*t3U) + r5 = (tP-kp951056516*t1s) :+ (t1J+kp951056516*t1M) + r4 = (t3l-kp921177326*t3o) :+ (t2Z+kp921177326*t36) + r3 = (t5B+kp860541664*t5E) :+ (t5x-kp860541664*t5y) + r2 = (t4F-kp943557151*t4M) :+ (t51-kp943557151*t52) + r1 = (t3Z+kp803003575*t42) :+ (t3V+kp803003575*t3W) + MV.unsafeWrite xsout 0 $ (t9+tM) :+ (t1D+t1G) + MV.unsafeWrite xsout 1 $ if sign == 1 then r1 else r24 + MV.unsafeWrite xsout 2 $ if sign == 1 then r2 else r23 + MV.unsafeWrite xsout 3 $ if sign == 1 then r3 else r22 + MV.unsafeWrite xsout 4 $ if sign == 1 then r4 else r21 + MV.unsafeWrite xsout 5 $ if sign == 1 then r5 else r20 + MV.unsafeWrite xsout 6 $ if sign == 1 then r6 else r19 + MV.unsafeWrite xsout 7 $ if sign == 1 then r7 else r18 + MV.unsafeWrite xsout 8 $ if sign == 1 then r8 else r17 + MV.unsafeWrite xsout 9 $ if sign == 1 then r9 else r16 + MV.unsafeWrite xsout 10 $ if sign == 1 then r10 else r15 + MV.unsafeWrite xsout 11 $ if sign == 1 then r11 else r14 + MV.unsafeWrite xsout 12 $ if sign == 1 then r12 else r13 + MV.unsafeWrite xsout 13 $ if sign == 1 then r13 else r12 + MV.unsafeWrite xsout 14 $ if sign == 1 then r14 else r11 + MV.unsafeWrite xsout 15 $ if sign == 1 then r15 else r10 + MV.unsafeWrite xsout 16 $ if sign == 1 then r16 else r9 + MV.unsafeWrite xsout 17 $ if sign == 1 then r17 else r8 + MV.unsafeWrite xsout 18 $ if sign == 1 then r18 else r7 + MV.unsafeWrite xsout 19 $ if sign == 1 then r19 else r6 + MV.unsafeWrite xsout 20 $ if sign == 1 then r20 else r5 + MV.unsafeWrite xsout 21 $ if sign == 1 then r21 else r4 + MV.unsafeWrite xsout 22 $ if sign == 1 then r22 else r3 + MV.unsafeWrite xsout 23 $ if sign == 1 then r23 else r2 + MV.unsafeWrite xsout 24 $ if sign == 1 then r24 else r1 + diff --git a/Numeric/FFT/Special/PowersOfTwo.hs b/Numeric/FFT/Special/PowersOfTwo.hs new file mode 100644 --- /dev/null +++ b/Numeric/FFT/Special/PowersOfTwo.hs @@ -0,0 +1,808 @@ +module Numeric.FFT.Special.PowersOfTwo + ( special2, special4, special8, special16, special32, special64 + ) where + +import Control.Monad.ST +import Data.IntMap.Strict (IntMap) +import qualified Data.IntMap.Strict as IM +import Data.Complex +import Data.Vector.Unboxed +import qualified Data.Vector.Unboxed.Mutable as MV + +import Numeric.FFT.Types +import Numeric.FFT.Utils + + +-- | Length 2 hard-coded FFT. +special2 :: Int -> MVCD s -> MVCD s -> ST s () +special2 _ xsin xsout = do + a <- MV.unsafeRead xsin 0 + b <- MV.unsafeRead xsin 1 + MV.unsafeWrite xsout 0 $ a + b + MV.unsafeWrite xsout 1 $ a - b + +-- | Length 4 hard-coded FFT. +special4 :: Int -> MVCD s -> MVCD s -> ST s () +special4 sign xsin xsout = do + xr0 :+ xi0 <- MV.unsafeRead xsin 0 ; xr1 :+ xi1 <- MV.unsafeRead xsin 1 + xr2 :+ xi2 <- MV.unsafeRead xsin 2 ; xr3 :+ xi3 <- MV.unsafeRead xsin 3 + let tb = xr0 - xr2 ; t3 = xr0 + xr2 ; tf = xi0 + xi2 ; t9 = xi0 - xi2 + t6 = xr1 + xr3 ; ta = xr1 - xr3 ; te = xi1 - xi3 ; tg = xi1 + xi3 + r3 = (tb + te) :+ (t9 - ta) + r2 = (t3 - t6) :+ (tf - tg) + r1 = (tb - te) :+ (ta + t9) + MV.unsafeWrite xsout 0 $ (t3 + t6) :+ (tf + tg) + MV.unsafeWrite xsout 1 $ if sign == 1 then r1 else r3 + MV.unsafeWrite xsout 2 $ if sign == 1 then r2 else r2 + MV.unsafeWrite xsout 3 $ if sign == 1 then r3 else r1 + +-- | Length 8 hard-coded FFT. +kp707106781 :: Double +kp707106781 = 0.707106781186547524400844362104849039284835938 +special8 :: Int -> MVCD s -> MVCD s -> ST s () +special8 sign xsin xsout = do + xr0 :+ xi0 <- MV.unsafeRead xsin 0 ; xr1 :+ xi1 <- MV.unsafeRead xsin 1 + xr2 :+ xi2 <- MV.unsafeRead xsin 2 ; xr3 :+ xi3 <- MV.unsafeRead xsin 3 + xr4 :+ xi4 <- MV.unsafeRead xsin 4 ; xr5 :+ xi5 <- MV.unsafeRead xsin 5 + xr6 :+ xi6 <- MV.unsafeRead xsin 6 ; xr7 :+ xi7 <- MV.unsafeRead xsin 7 + let tn = xr0 - xr4 ; t3 = xr0 + xr4 ; tC = xi0 - xi4 ; ti = xi0 + xi4 + tB = xr2 - xr6 ; t6 = xr2 + xr6 ; to = xi2 - xi6 ; tl = xi2 + xi6 + td = xr7 + xr3 ; tv = xr7 - xr3 ; tN = xi7 + xi3 ; ty = xi7 - xi3 + tz = tv - ty ; tH = tv + ty ; ta = xr1 + xr5 ; tq = xr1 - xr5 + tt = xi1 - xi5 ; tM = xi1 + xi5 ; tL = t3 - t6 ; t7 = t3 + t6 + tG = tt - tq ; tu = tq + tt ; te = ta + td ; tf = td - ta + tm = ti - tl ; tP = ti + tl ; tQ = tM + tN ; tO = tM - tN + tF = tn - to ; tp = tn + to ; tA = tu + tz ; tE = tz - tu + tD = tB + tC ; tJ = tC - tB ; tK = tG + tH ; tI = tG - tH + r7 = (tp + kp707106781 * tA) :+ (tJ + kp707106781 * tK) + r6 = (tL + tO) :+ (tf + tm) + r5 = (tF + kp707106781 * tI) :+ (tD + kp707106781 * tE) + r4 = (t7 - te) :+ (tP - tQ) + r3 = (tp - kp707106781 * tA) :+ (tJ - kp707106781 * tK) + r2 = (tL - tO) :+ (tm - tf) + r1 = (tF - kp707106781 * tI) :+ (tD - kp707106781 * tE) + MV.unsafeWrite xsout 0 $ (t7 + te) :+ (tP + tQ) + MV.unsafeWrite xsout 1 $ if sign == 1 then r1 else r7 + MV.unsafeWrite xsout 2 $ if sign == 1 then r2 else r6 + MV.unsafeWrite xsout 3 $ if sign == 1 then r3 else r5 + MV.unsafeWrite xsout 4 $ if sign == 1 then r4 else r4 + MV.unsafeWrite xsout 5 $ if sign == 1 then r5 else r3 + MV.unsafeWrite xsout 6 $ if sign == 1 then r6 else r2 + MV.unsafeWrite xsout 7 $ if sign == 1 then r7 else r1 + +-- | Length 16 hard-coded FFT. +kp923879532, kp414213562 :: Double +--kp707106781 :: Double +kp923879532 = 0.923879532511286756128183189396788286822416626 +kp414213562 = 0.414213562373095048801688724209698078569671875 +--kp707106781 = 0.707106781186547524400844362104849039284835938 +special16 :: Int -> MVCD s -> MVCD s -> ST s () +special16 sign xsin xsout = do + xr0 :+ xi0 <- MV.unsafeRead xsin 0 ; xr1 :+ xi1 <- MV.unsafeRead xsin 1 + xr2 :+ xi2 <- MV.unsafeRead xsin 2 ; xr3 :+ xi3 <- MV.unsafeRead xsin 3 + xr4 :+ xi4 <- MV.unsafeRead xsin 4 ; xr5 :+ xi5 <- MV.unsafeRead xsin 5 + xr6 :+ xi6 <- MV.unsafeRead xsin 6 ; xr7 :+ xi7 <- MV.unsafeRead xsin 7 + xr8 :+ xi8 <- MV.unsafeRead xsin 8 ; xr9 :+ xi9 <- MV.unsafeRead xsin 9 + xr10 :+ xi10 <- MV.unsafeRead xsin 10 ; xr11 :+ xi11 <- MV.unsafeRead xsin 11 + xr12 :+ xi12 <- MV.unsafeRead xsin 12 ; xr13 :+ xi13 <- MV.unsafeRead xsin 13 + xr14 :+ xi14 <- MV.unsafeRead xsin 14 ; xr15 :+ xi15 <- MV.unsafeRead xsin 15 + let tL = xr0 - xr8 ; t3 = xr0 + xr8 ; t1k = xi0 - xi8 ; ty = xi0 + xi8 + t1j = xr4 - xr12 ; t6 = xr4 + xr12 ; tM = xi4 - xi12 ; tB = xi4 + xi12 + t1l = t1j + t1k ; t1H = t1k - t1j ; t1R = t3 - t6 ; t7 = t3 + t6 + t1x = tL + tM ; tN = tL - tM ; tC = ty + tB ; t25 = ty - tB + t1c = xr15 - xr7 ; tp = xr15 + xr7 ; t20 = xi15 + xi7 ; t1a = xi15 - xi7 + t17 = xr3 - xr11 ; ts = xr3 + xr11 ; t21 = xi3 + xi11 ; t1f = xi3 - xi11 + t1E = t1a - t17 ; t1b = t17 + t1a ; t1Z = tp - ts ; tt = tp + ts + t2h = t20 + t21 ; t22 = t20 - t21 ; t1D = t1c + t1f ; t1g = t1c - t1f + tP = xr2 - xr10 ; ta = xr2 + xr10 ; tO = xi2 - xi10 ; tF = xi2 + xi10 + t1n = tP + tO ; tQ = tO - tP ; tR = xr14 - xr6 ; td = xr14 + xr6 + tS = xi14 - xi6 ; tI = xi14 + xi6 ; te = ta + td ; t26 = td - ta + tT = tR + tS ; t1m = tR - tS ; tJ = tF + tI ; t1S = tF - tI + t11 = xr1 - xr9 ; ti = xr1 + xr9 ; t1V = xi1 + xi9 ; tZ = xi1 - xi9 + t2f = t7 - te ; tf = t7 + te ; tW = xr5 - xr13 ; tl = xr5 + xr13 + t1W = xi5 + xi13 ; t14 = xi5 - xi13 ; t1B = tZ - tW ; t10 = tW + tZ + t1U = ti - tl ; tm = ti + tl ; t2g = t1V + t1W ; t1X = t1V - t1W + t1A = t11 + t14 ; t15 = t11 - t14 ; tu = tm + tt ; tv = tt - tm + tK = tC - tJ ; t2j = tC + tJ ; t2k = t2g + t2h ; t2i = t2g - t2h + t29 = t1R - t1S ; t1T = t1R + t1S ; t27 = t25 - t26 ; t2d = t26 + t25 + t2a = t1X - t1U ; t1Y = t1U + t1X ; t23 = t1Z - t22 ; t2b = t1Z + t22 + t28 = t23 - t1Y ; t24 = t1Y + t23 ; t1I = tQ + tT ; tU = tQ - tT + t2e = t2a + t2b ; t2c = t2a - t2b + tV = tN + kp707106781 * tU ; t1v = tN - kp707106781 * tU + t1o = t1m - t1n ; t1y = t1n + t1m + t1t = t15 - kp414213562 * t10 ; t16 = t10 + kp414213562 * t15 + t1h = t1b - kp414213562 * t1g ; t1s = t1g + kp414213562 * t1b + t1r = t1l + kp707106781 * t1o ; t1p = t1l - kp707106781 * t1o + t1q = t16 + t1h ; t1i = t16 - t1h ; t1w = t1t + t1s ; t1u = t1s - t1t + t1z = t1x + kp707106781 * t1y ; t1L = t1x - kp707106781 * t1y + t1M = t1B - kp414213562 * t1A ; t1C = t1A + kp414213562 * t1B + t1F = t1D - kp414213562 * t1E ; t1N = t1E + kp414213562 * t1D + t1P = t1H + kp707106781 * t1I ; t1J = t1H - kp707106781 * t1I + t1K = t1F - t1C ; t1G = t1C + t1F ; t1O = t1M - t1N ; t1Q = t1M + t1N + r15 = (t1z + kp923879532 * t1G) :+ (t1P + kp923879532 * t1Q) + r14 = (t1T + kp707106781 * t24) :+ (t2d + kp707106781 * t2e) + r13 = (tV + kp923879532 * t1i) :+ (t1r + kp923879532 * t1u) + r12 = (t2f + t2i) :+ (tv + tK) + r11 = (t1L + kp923879532 * t1O) :+ (t1J + kp923879532 * t1K) + r10 = (t29 + kp707106781 * t2c) :+ (t27 + kp707106781 * t28) + r9 = (t1v - kp923879532 * t1w) :+ (t1p - kp923879532 * t1q) + r8 = (tf - tu) :+ (t2j - t2k) + r7 = (t1z - kp923879532 * t1G) :+ (t1P - kp923879532 * t1Q) + r6 = (t1T - kp707106781 * t24) :+ (t2d - kp707106781 * t2e) + r5 = (tV - kp923879532 * t1i) :+ (t1r - kp923879532 * t1u) + r4 = (t2f - t2i) :+ (tK - tv) + r3 = (t1L - kp923879532 * t1O) :+ (t1J - kp923879532 * t1K) + r2 = (t29 - kp707106781 * t2c) :+ (t27 - kp707106781 * t28) + r1 = (t1v + kp923879532 * t1w) :+ (t1p + kp923879532 * t1q) + MV.unsafeWrite xsout 0 $ (tf + tu) :+ (t2j + t2k) + MV.unsafeWrite xsout 1 $ if sign == 1 then r1 else r15 + MV.unsafeWrite xsout 2 $ if sign == 1 then r2 else r14 + MV.unsafeWrite xsout 3 $ if sign == 1 then r3 else r13 + MV.unsafeWrite xsout 4 $ if sign == 1 then r4 else r12 + MV.unsafeWrite xsout 5 $ if sign == 1 then r5 else r11 + MV.unsafeWrite xsout 6 $ if sign == 1 then r6 else r10 + MV.unsafeWrite xsout 7 $ if sign == 1 then r7 else r9 + MV.unsafeWrite xsout 8 $ if sign == 1 then r8 else r8 + MV.unsafeWrite xsout 9 $ if sign == 1 then r9 else r7 + MV.unsafeWrite xsout 10 $ if sign == 1 then r10 else r6 + MV.unsafeWrite xsout 11 $ if sign == 1 then r11 else r5 + MV.unsafeWrite xsout 12 $ if sign == 1 then r12 else r4 + MV.unsafeWrite xsout 13 $ if sign == 1 then r13 else r3 + MV.unsafeWrite xsout 14 $ if sign == 1 then r14 else r2 + MV.unsafeWrite xsout 15 $ if sign == 1 then r15 else r1 + +-- | Length 32 hard-coded FFT. +kp980785280, kp198912367, kp831469612, kp668178637 :: Double +--kp923879532, kp707106781, kp414213562 :: Double +kp980785280 = 0.980785280403230449126182236134239036973933731 +kp198912367 = 0.198912367379658006911597622644676228597850501 +kp831469612 = 0.831469612302545237078788377617905756738560812 +kp668178637 = 0.668178637919298919997757686523080761552472251 +--kp923879532 = 0.923879532511286756128183189396788286822416626 +--kp707106781 = 0.707106781186547524400844362104849039284835938 +--kp414213562 = 0.414213562373095048801688724209698078569671875 +special32 :: Int -> MVCD s -> MVCD s -> ST s () +special32 sign xsin xsout = do + xr0 :+ xi0 <- MV.unsafeRead xsin 0 ; xr1 :+ xi1 <- MV.unsafeRead xsin 1 + xr2 :+ xi2 <- MV.unsafeRead xsin 2 ; xr3 :+ xi3 <- MV.unsafeRead xsin 3 + xr4 :+ xi4 <- MV.unsafeRead xsin 4 ; xr5 :+ xi5 <- MV.unsafeRead xsin 5 + xr6 :+ xi6 <- MV.unsafeRead xsin 6 ; xr7 :+ xi7 <- MV.unsafeRead xsin 7 + xr8 :+ xi8 <- MV.unsafeRead xsin 8 ; xr9 :+ xi9 <- MV.unsafeRead xsin 9 + xr10 :+ xi10 <- MV.unsafeRead xsin 10 ; xr11 :+ xi11 <- MV.unsafeRead xsin 11 + xr12 :+ xi12 <- MV.unsafeRead xsin 12 ; xr13 :+ xi13 <- MV.unsafeRead xsin 13 + xr14 :+ xi14 <- MV.unsafeRead xsin 14 ; xr15 :+ xi15 <- MV.unsafeRead xsin 15 + xr16 :+ xi16 <- MV.unsafeRead xsin 16 ; xr17 :+ xi17 <- MV.unsafeRead xsin 17 + xr18 :+ xi18 <- MV.unsafeRead xsin 18 ; xr19 :+ xi19 <- MV.unsafeRead xsin 19 + xr20 :+ xi20 <- MV.unsafeRead xsin 20 ; xr21 :+ xi21 <- MV.unsafeRead xsin 21 + xr22 :+ xi22 <- MV.unsafeRead xsin 22 ; xr23 :+ xi23 <- MV.unsafeRead xsin 23 + xr24 :+ xi24 <- MV.unsafeRead xsin 24 ; xr25 :+ xi25 <- MV.unsafeRead xsin 25 + xr26 :+ xi26 <- MV.unsafeRead xsin 26 ; xr27 :+ xi27 <- MV.unsafeRead xsin 27 + xr28 :+ xi28 <- MV.unsafeRead xsin 28 ; xr29 :+ xi29 <- MV.unsafeRead xsin 29 + xr30 :+ xi30 <- MV.unsafeRead xsin 30 ; xr31 :+ xi31 <- MV.unsafeRead xsin 31 + let t1x = xr0-xr16 ; t3 = xr0+xr16 ; t2R = xi0-xi16 ; t14 = xi0+xi16 + t2S = xr8-xr24 ; t6 = xr8+xr24 ; t1y = xi8-xi24 ; t17 = xi8+xi24 + t2T = t2R-t2S ; t3T = t2S+t2R ; t4r = t3-t6 ; t7 = t3+t6 + t3t = t1x-t1y ; t1z = t1x+t1y ; t18 = t14+t17 ; t4Z = t14-t17 + t1A = xr4-xr20 ; ta = xr4+xr20 ; t1B = xi4-xi20 ; t1b = xi4+xi20 + t1C = t1A+t1B ; t2U = t1B-t1A ; t1D = xr28-xr12 ; td = xr28+xr12 + t1E = xi28-xi12 ; t1e = xi28+xi12 ; te = ta+td ; t50 = td-ta + t1F = t1D-t1E ; t2V = t1D+t1E ; t4s = t1b-t1e ; t1f = t1b+t1e + t2W = t2U+t2V ; t3u = t2U-t2V ; t3U = t1F-t1C ; t1G = t1C+t1F + t1L = xr2-xr18 ; ti = xr2+xr18 ; t1I = xi2-xi18 ; t1j = xi2+xi18 + t1J = xr10-xr26 ; tl = xr10+xr26 ; t1M = xi10-xi26 ; t1m = xi10+xi26 + t3w = t1J+t1I ; t1K = t1I-t1J ; t4v = ti-tl ; tm = ti+tl + t3x = t1L-t1M ; t1N = t1L+t1M ; t4u = t1j-t1m ; t1n = t1j+t1m + t3X = t3x - kp414213562 * t3w ; t3y = t3w + kp414213562 * t3x + t2Z = t1N + kp414213562 * t1K ; t1O = t1K - kp414213562 * t1N + t53 = t4v+t4u ; t4w = t4u-t4v ; t1S = xr30-xr14 ; tp = xr30+xr14 + t1P = xi30-xi14 ; t1q = xi30+xi14 ; t1Q = xr6-xr22 ; ts = xr6+xr22 + t1T = xi6-xi22 ; t1t = xi6+xi22 ; t3z = t1Q+t1P ; t1R = t1P-t1Q + t4x = tp-ts ; tt = tp+ts ; t3A = t1S-t1T ; t1U = t1S+t1T + t4y = t1q-t1t ; t1u = t1q+t1t + t3W = t3A + kp414213562 * t3z ; t3B = t3z - kp414213562 * t3A + t2Y = t1U - kp414213562 * t1R ; t1V = t1R + kp414213562 * t1U + t52 = t4x-t4y ; t4z = t4x+t4y ; t2G = xr31-xr15 ; tN = xr31+xr15 + t4N = xi31+xi15 ; t2r = xi31-xi15 ; t2s = xr7-xr23 ; tQ = xr7+xr23 + t4O = xi7+xi23 ; t2J = xi7-xi23 ; t2x = xr3-xr19 ; tU = xr3+xr19 + t4T = xi3+xi19 ; t2w = xi3-xi19 ; t3O = t2s+t2r ; t2t = t2r-t2s + t2z = xr27-xr11 ; tX = xr27+xr11 ; t4U = xi27+xi11 ; t2C = xi27-xi11 + t3L = t2G-t2J ; t2K = t2G+t2J ; t4S = tN-tQ ; tR = tN+tQ + tY = tU+tX ; t4Q = tX-tU ; t4P = t4N-t4O ; t5G = t4N+t4O + t5H = t4T+t4U ; t4V = t4T-t4U ; t5F = tR-tY ; tZ = tR+tY + t5I = t5G-t5H ; t5X = t5G+t5H ; t2L = t2x+t2w ; t2y = t2w-t2x + t2D = t2z+t2C ; t2M = t2z-t2C ; t4R = t4P-t4Q ; t5k = t4Q+t4P + t3M = t2D-t2y ; t2E = t2y+t2D ; t5j = t4S+t4V ; t4W = t4S-t4V + t3P = t2L-t2M ; t2N = t2L+t2M ; t2f = xr1-xr17 ; ty = xr1+xr17 + t4C = xi1+xi17 ; t20 = xi1-xi17 ; t21 = xr9-xr25 ; tB = xr9+xr25 + t4D = xi9+xi25 ; t2i = xi9-xi25 ; t26 = xr5-xr21 ; tF = xr5+xr21 + t4I = xi5+xi21 ; t25 = xi5-xi21 ; t3H = t21+t20 ; t22 = t20-t21 + t28 = xr29-xr13 ; tI = xr29+xr13 ; t4J = xi29+xi13 ; t2b = xi29-xi13 + t3E = t2f-t2i ; t2j = t2f+t2i ; t4H = ty-tB ; tC = ty+tB + tJ = tF+tI ; t4F = tI-tF ; t4E = t4C-t4D ; t5B = t4C+t4D + t5C = t4I+t4J ; t4K = t4I-t4J ; t5A = tC-tJ ; tK = tC+tJ + t5D = t5B-t5C ; t5W = t5B+t5C ; t2k = t26+t25 ; t27 = t25-t26 + t2c = t28+t2b ; t2l = t28-t2b ; t4G = t4E-t4F ; t5h = t4F+t4E + t3F = t2c-t27 ; t2d = t27+t2c ; t5d = t4r+t4s ; t4t = t4r-t4s + t5g = t4H+t4K ; t4L = t4H-t4K ; t3I = t2k-t2l ; t2m = t2k+t2l + t4A = t4w-t4z ; t5o = t4w+t4z + t4X = t4R - kp414213562 * t4W ; t58 = t4W + kp414213562 * t4R + t59 = t4L - kp414213562 * t4G ; t4M = t4G + kp414213562 * t4L + t5b = t4t - kp707106781 * t4A ; t4B = t4t + kp707106781 * t4A + t5c = t59+t58 ; t5a = t58-t59 ; t5n = t50+t4Z ; t51 = t4Z-t50 + t54 = t52-t53 ; t5e = t53+t52 ; t56 = t4M+t4X ; t4Y = t4M-t4X + t57 = t51 + kp707106781 * t54 ; t55 = t51 - kp707106781 * t54 + t5i = t5g + kp414213562 * t5h ; t5s = t5h - kp414213562 * t5g + t5t = t5k + kp414213562 * t5j ; t5l = t5j - kp414213562 * t5k + t5r = t5d - kp707106781 * t5e ; t5f = t5d + kp707106781 * t5e + t5w = t5s+t5t ; t5u = t5s-t5t ; t5q = t5l-t5i ; t5m = t5i+t5l + t5v = t5n + kp707106781 * t5o ; t5p = t5n - kp707106781 * t5o + tf = t7+te ; t5x = t7-te ; t5y = t1n-t1u ; t1v = t1n+t1u + t5E = t5A+t5D ; t5Q = t5D-t5A ; t5R = t5F+t5I ; t5J = t5F-t5I + t5P = t5x-t5y ; t5z = t5x+t5y ; t5U = t5Q+t5R ; t5S = t5Q-t5R + t1g = t18+t1f ; t5L = t18-t1f ; t5M = tt-tm ; tu = tm+tt + t5O = t5J-t5E ; t5K = t5E+t5J ; t5T = t5M+t5L ; t5N = t5L-t5M + t5V = tf-tu ; tv = tf+tu ; t60 = t5W+t5X ; t5Y = t5W-t5X + t11 = tZ-tK ; t10 = tK+tZ ; t5Z = t1g+t1v ; t1w = t1g-t1v + t39 = t1z + kp707106781 * t1G ; t1H = t1z - kp707106781 * t1G + t1W = t1O-t1V ; t3k = t1O+t1V + t3j = t2T + kp707106781 * t2W ; t2X = t2T - kp707106781 * t2W + t30 = t2Y-t2Z ; t3a = t2Z+t2Y + t3d = t22 + kp707106781 * t2d ; t2e = t22 - kp707106781 * t2d + t37 = t1H - kp923879532 * t1W ; t1X = t1H + kp923879532 * t1W + t33 = t2X + kp923879532 * t30 ; t31 = t2X - kp923879532 * t30 + t2n = t2j - kp707106781 * t2m ; t3c = t2j + kp707106781 * t2m + t3g = t2t + kp707106781 * t2E ; t2F = t2t - kp707106781 * t2E + t2O = t2K - kp707106781 * t2N ; t3f = t2K + kp707106781 * t2N + t47 = t3t - kp707106781 * t3u ; t3v = t3t + kp707106781 * t3u + t35 = t2n - kp668178637 * t2e ; t2o = t2e + kp668178637 * t2n + t34 = t2O + kp668178637 * t2F ; t2P = t2F - kp668178637 * t2O + t3C = t3y-t3B ; t4i = t3y+t3B + t4h = t3T - kp707106781 * t3U ; t3V = t3T + kp707106781 * t3U + t38 = t35+t34 ; t36 = t34-t35 ; t32 = t2o+t2P ; t2Q = t2o-t2P + t41 = t3v - kp923879532 * t3C ; t3D = t3v + kp923879532 * t3C + t3Y = t3W-t3X ; t48 = t3X+t3W + t4b = t3E + kp707106781 * t3F ; t3G = t3E - kp707106781 * t3F + t3J = t3H - kp707106781 * t3I ; t4a = t3H + kp707106781 * t3I + t4e = t3L + kp707106781 * t3M ; t3N = t3L - kp707106781 * t3M + t45 = t3V + kp923879532 * t3Y ; t3Z = t3V - kp923879532 * t3Y + t42 = t3J - kp668178637 * t3G ; t3K = t3G + kp668178637 * t3J + t3Q = t3O - kp707106781 * t3P ; t4d = t3O + kp707106781 * t3P + t43 = t3Q + kp668178637 * t3N ; t3R = t3N - kp668178637 * t3Q + t4p = t47 + kp923879532 * t48 ; t49 = t47 - kp923879532 * t48 + t44 = t42-t43 ; t46 = t42+t43 ; t40 = t3R-t3K ; t3S = t3K+t3R + t4l = t4h - kp923879532 * t4i ; t4j = t4h + kp923879532 * t4i + t4n = t4b - kp198912367 * t4a ; t4c = t4a + kp198912367 * t4b + t4m = t4e + kp198912367 * t4d ; t4f = t4d - kp198912367 * t4e + t3n = t39 - kp923879532 * t3a ; t3b = t39 + kp923879532 * t3a + t4q = t4n+t4m ; t4o = t4m-t4n ; t4k = t4c+t4f ; t4g = t4c-t4f + t3r = t3j + kp923879532 * t3k ; t3l = t3j - kp923879532 * t3k + t3o = t3d - kp198912367 * t3c ; t3e = t3c + kp198912367 * t3d + t3h = t3f - kp198912367 * t3g ; t3p = t3g + kp198912367 * t3f + t3s = t3o+t3p ; t3q = t3o-t3p ; t3i = t3e+t3h ; t3m = t3h-t3e + r31 = (t3b + kp980785280 * t3i) :+ (t3r + kp980785280 * t3s) + r30 = (t5f + kp923879532 * t5m) :+ (t5v + kp923879532 * t5w) + r29 = (t3D + kp831469612 * t3S) :+ (t45 + kp831469612 * t46) + r28 = (t5z + kp707106781 * t5K) :+ (t5T + kp707106781 * t5U) + r27 = (t1X + kp831469612 * t2Q) :+ (t33 + kp831469612 * t36) + r26 = (t4B + kp923879532 * t4Y) :+ (t57 + kp923879532 * t5a) + r25 = (t49 + kp980785280 * t4g) :+ (t4l + kp980785280 * t4o) + r24 = (t5V+t5Y) :+ (t11+t1w) + r23 = (t3n + kp980785280 * t3q) :+ (t3l + kp980785280 * t3m) + r22 = (t5r + kp923879532 * t5u) :+ (t5p + kp923879532 * t5q) + r21 = (t41 + kp831469612 * t44) :+ (t3Z + kp831469612 * t40) + r20 = (t5P + kp707106781 * t5S) :+ (t5N + kp707106781 * t5O) + r19 = (t37 - kp831469612 * t38) :+ (t31 - kp831469612 * t32) + r18 = (t5b - kp923879532 * t5c) :+ (t55 - kp923879532 * t56) + r17 = (t4p - kp980785280 * t4q) :+ (t4j - kp980785280 * t4k) + r16 = (tv-t10) :+ (t5Z-t60) + r15 = (t3b - kp980785280 * t3i) :+ (t3r - kp980785280 * t3s) + r14 = (t5f - kp923879532 * t5m) :+ (t5v - kp923879532 * t5w) + r13 = (t3D - kp831469612 * t3S) :+ (t45 - kp831469612 * t46) + r12 = (t5z - kp707106781 * t5K) :+ (t5T - kp707106781 * t5U) + r11 = (t1X - kp831469612 * t2Q) :+ (t33 - kp831469612 * t36) + r10 = (t4B - kp923879532 * t4Y) :+ (t57 - kp923879532 * t5a) + r9 = (t49 - kp980785280 * t4g) :+ (t4l - kp980785280 * t4o) + r8 = (t5V-t5Y) :+ (t1w-t11) + r7 = (t3n - kp980785280 * t3q) :+ (t3l - kp980785280 * t3m) + r6 = (t5r - kp923879532 * t5u) :+ (t5p - kp923879532 * t5q) + r5 = (t41 - kp831469612 * t44) :+ (t3Z - kp831469612 * t40) + r4 = (t5P - kp707106781 * t5S) :+ (t5N - kp707106781 * t5O) + r3 = (t37 + kp831469612 * t38) :+ (t31 + kp831469612 * t32) + r2 = (t5b + kp923879532 * t5c) :+ (t55 + kp923879532 * t56) + r1 = (t4p + kp980785280 * t4q) :+ (t4j + kp980785280 * t4k) + MV.unsafeWrite xsout 0 $ (tv+t10) :+ (t5Z+t60) + MV.unsafeWrite xsout 1 $ if sign == 1 then r1 else r31 + MV.unsafeWrite xsout 2 $ if sign == 1 then r2 else r30 + MV.unsafeWrite xsout 3 $ if sign == 1 then r3 else r29 + MV.unsafeWrite xsout 4 $ if sign == 1 then r4 else r28 + MV.unsafeWrite xsout 5 $ if sign == 1 then r5 else r27 + MV.unsafeWrite xsout 6 $ if sign == 1 then r6 else r26 + MV.unsafeWrite xsout 7 $ if sign == 1 then r7 else r25 + MV.unsafeWrite xsout 8 $ if sign == 1 then r8 else r24 + MV.unsafeWrite xsout 9 $ if sign == 1 then r9 else r23 + MV.unsafeWrite xsout 10 $ if sign == 1 then r10 else r22 + MV.unsafeWrite xsout 11 $ if sign == 1 then r11 else r21 + MV.unsafeWrite xsout 12 $ if sign == 1 then r12 else r20 + MV.unsafeWrite xsout 13 $ if sign == 1 then r13 else r19 + MV.unsafeWrite xsout 14 $ if sign == 1 then r14 else r18 + MV.unsafeWrite xsout 15 $ if sign == 1 then r15 else r17 + MV.unsafeWrite xsout 16 $ if sign == 1 then r16 else r16 + MV.unsafeWrite xsout 17 $ if sign == 1 then r17 else r15 + MV.unsafeWrite xsout 18 $ if sign == 1 then r18 else r14 + MV.unsafeWrite xsout 19 $ if sign == 1 then r19 else r13 + MV.unsafeWrite xsout 20 $ if sign == 1 then r20 else r12 + MV.unsafeWrite xsout 21 $ if sign == 1 then r21 else r11 + MV.unsafeWrite xsout 22 $ if sign == 1 then r22 else r10 + MV.unsafeWrite xsout 23 $ if sign == 1 then r23 else r9 + MV.unsafeWrite xsout 24 $ if sign == 1 then r24 else r8 + MV.unsafeWrite xsout 25 $ if sign == 1 then r25 else r7 + MV.unsafeWrite xsout 26 $ if sign == 1 then r26 else r6 + MV.unsafeWrite xsout 27 $ if sign == 1 then r27 else r5 + MV.unsafeWrite xsout 28 $ if sign == 1 then r28 else r4 + MV.unsafeWrite xsout 29 $ if sign == 1 then r29 else r3 + MV.unsafeWrite xsout 30 $ if sign == 1 then r30 else r2 + MV.unsafeWrite xsout 31 $ if sign == 1 then r31 else r1 + +-- | Length 64 hard-coded FFT. +kp956940335, kp881921264, kp534511135, kp303346683 :: Double +kp995184726, kp773010453, kp820678790, kp098491403 :: Double +--kp980785280, kp831469612, kp668178637, kp198912367 :: Double +--kp923879532, kp707106781, kp414213562 :: Double +kp956940335 = 0.956940335732208864935797886980269969482849206 +kp881921264 = 0.881921264348355029712756863660388349508442621 +kp534511135 = 0.534511135950791641089685961295362908582039528 +kp303346683 = 0.303346683607342391675883946941299872384187453 +kp995184726 = 0.995184726672196886244836953109479921575474869 +kp773010453 = 0.773010453362736960810906609758469800971041293 +kp820678790 = 0.820678790828660330972281985331011598767386482 +kp098491403 = 0.098491403357164253077197521291327432293052451 +--kp980785280 = 0.980785280403230449126182236134239036973933731 +--kp831469612 = 0.831469612302545237078788377617905756738560812 +--kp668178637 = 0.668178637919298919997757686523080761552472251 +--kp198912367 = 0.198912367379658006911597622644676228597850501 +--kp923879532 = 0.923879532511286756128183189396788286822416626 +--kp707106781 = 0.707106781186547524400844362104849039284835938 +--kp414213562 = 0.414213562373095048801688724209698078569671875 +special64 :: Int -> MVCD s -> MVCD s -> ST s () +special64 sign xsin xsout = do + xr0 :+ xi0 <- MV.unsafeRead xsin 0 ; xr1 :+ xi1 <- MV.unsafeRead xsin 1 + xr2 :+ xi2 <- MV.unsafeRead xsin 2 ; xr3 :+ xi3 <- MV.unsafeRead xsin 3 + xr4 :+ xi4 <- MV.unsafeRead xsin 4 ; xr5 :+ xi5 <- MV.unsafeRead xsin 5 + xr6 :+ xi6 <- MV.unsafeRead xsin 6 ; xr7 :+ xi7 <- MV.unsafeRead xsin 7 + xr8 :+ xi8 <- MV.unsafeRead xsin 8 ; xr9 :+ xi9 <- MV.unsafeRead xsin 9 + xr10 :+ xi10 <- MV.unsafeRead xsin 10 ; xr11 :+ xi11 <- MV.unsafeRead xsin 11 + xr12 :+ xi12 <- MV.unsafeRead xsin 12 ; xr13 :+ xi13 <- MV.unsafeRead xsin 13 + xr14 :+ xi14 <- MV.unsafeRead xsin 14 ; xr15 :+ xi15 <- MV.unsafeRead xsin 15 + xr16 :+ xi16 <- MV.unsafeRead xsin 16 ; xr17 :+ xi17 <- MV.unsafeRead xsin 17 + xr18 :+ xi18 <- MV.unsafeRead xsin 18 ; xr19 :+ xi19 <- MV.unsafeRead xsin 19 + xr20 :+ xi20 <- MV.unsafeRead xsin 20 ; xr21 :+ xi21 <- MV.unsafeRead xsin 21 + xr22 :+ xi22 <- MV.unsafeRead xsin 22 ; xr23 :+ xi23 <- MV.unsafeRead xsin 23 + xr24 :+ xi24 <- MV.unsafeRead xsin 24 ; xr25 :+ xi25 <- MV.unsafeRead xsin 25 + xr26 :+ xi26 <- MV.unsafeRead xsin 26 ; xr27 :+ xi27 <- MV.unsafeRead xsin 27 + xr28 :+ xi28 <- MV.unsafeRead xsin 28 ; xr29 :+ xi29 <- MV.unsafeRead xsin 29 + xr30 :+ xi30 <- MV.unsafeRead xsin 30 ; xr31 :+ xi31 <- MV.unsafeRead xsin 31 + xr32 :+ xi32 <- MV.unsafeRead xsin 32 ; xr33 :+ xi33 <- MV.unsafeRead xsin 33 + xr34 :+ xi34 <- MV.unsafeRead xsin 34 ; xr35 :+ xi35 <- MV.unsafeRead xsin 35 + xr36 :+ xi36 <- MV.unsafeRead xsin 36 ; xr37 :+ xi37 <- MV.unsafeRead xsin 37 + xr38 :+ xi38 <- MV.unsafeRead xsin 38 ; xr39 :+ xi39 <- MV.unsafeRead xsin 39 + xr40 :+ xi40 <- MV.unsafeRead xsin 40 ; xr41 :+ xi41 <- MV.unsafeRead xsin 41 + xr42 :+ xi42 <- MV.unsafeRead xsin 42 ; xr43 :+ xi43 <- MV.unsafeRead xsin 43 + xr44 :+ xi44 <- MV.unsafeRead xsin 44 ; xr45 :+ xi45 <- MV.unsafeRead xsin 45 + xr46 :+ xi46 <- MV.unsafeRead xsin 46 ; xr47 :+ xi47 <- MV.unsafeRead xsin 47 + xr48 :+ xi48 <- MV.unsafeRead xsin 48 ; xr49 :+ xi49 <- MV.unsafeRead xsin 49 + xr50 :+ xi50 <- MV.unsafeRead xsin 50 ; xr51 :+ xi51 <- MV.unsafeRead xsin 51 + xr52 :+ xi52 <- MV.unsafeRead xsin 52 ; xr53 :+ xi53 <- MV.unsafeRead xsin 53 + xr54 :+ xi54 <- MV.unsafeRead xsin 54 ; xr55 :+ xi55 <- MV.unsafeRead xsin 55 + xr56 :+ xi56 <- MV.unsafeRead xsin 56 ; xr57 :+ xi57 <- MV.unsafeRead xsin 57 + xr58 :+ xi58 <- MV.unsafeRead xsin 58 ; xr59 :+ xi59 <- MV.unsafeRead xsin 59 + xr60 :+ xi60 <- MV.unsafeRead xsin 60 ; xr61 :+ xi61 <- MV.unsafeRead xsin 61 + xr62 :+ xi62 <- MV.unsafeRead xsin 62 ; xr63 :+ xi63 <- MV.unsafeRead xsin 63 + let t35 = xr0-xr32 ; t3 = xr0+xr32 ; t5Y = xi0-xi32 ; t26 = xi0+xi32 + t5X = xr16-xr48 ; t6 = xr16+xr48 ; t36 = xi16-xi48 ; t29 = xi16+xi48 + t39 = xr8-xr40 ; ta = xr8+xr40 ; t38 = xi8-xi40 ; t2d = xi8+xi40 + t7B = t35+t36 ; t37 = t35-t36 ; t3b = xr56-xr24 ; td = xr56+xr24 + t3c = xi56-xi24 ; t2g = xi56+xi24 ; t5Z = t5X+t5Y ; t8F = t5Y-t5X + taf = t3-t6 ; t7 = t3+t6 ; te = ta+td ; tbz = td-ta + tbA = t26-t29 ; t2a = t26+t29 ; t3d = t3b+t3c ; t60 = t3b-t3c + td9 = t7-te ; tf = t7+te ; tcB = tbA-tbz ; tbB = tbz+tbA + t61 = t39+t38 ; t3a = t38-t39 ; t2h = t2d+t2g ; tag = t2d-t2g + t7C = t61+t60 ; t62 = t60-t61 ; tdH = t2a-t2h ; t2i = t2a+t2h + tcb = taf-tag ; tah = taf+tag ; t8G = t3a+t3d ; t3e = t3a-t3d + t3j = xr4-xr36 ; ti = xr4+xr36 ; t3h = xi4-xi36 ; t2l = xi4+xi36 + t3g = xr20-xr52 ; tl = xr20+xr52 ; t3k = xi20-xi52 ; t2o = xi20+xi52 + t3q = xr60-xr28 ; tp = xr60+xr28 ; t3o = xi60-xi28 ; t2s = xi60+xi28 + tai = ti-tl ; tm = ti+tl ; t3n = xr12-xr44 ; ts = xr12+xr44 + t3r = xi12-xi44 ; t2v = xi12+xi44 ; tal = tp-ts ; tt = tp+ts + taj = t2l-t2o ; t2p = t2l+t2o ; tam = t2s-t2v ; t2w = t2s+t2v + tu = tm+tt ; tdI = tt-tm ; tak = tai+taj ; tbC = taj-tai + tbD = tal+tam ; tan = tal-tam ; t7F = t3h-t3g ; t3i = t3g+t3h + t3l = t3j-t3k ; t7E = t3j+t3k ; tda = t2p-t2w ; t2x = t2p+t2w + t65 = t3l-kp414213562*t3i ; t3m = t3i+kp414213562*t3l + t3s = t3q-t3r ; t7H = t3q+t3r ; t7I = t3o-t3n ; t3p = t3n+t3o + t8I = t7F-kp414213562*t7E ; t7G = t7E+kp414213562*t7F + t8J = t7I+kp414213562*t7H ; t7J = t7H-kp414213562*t7I + t64 = t3s+kp414213562*t3p ; t3t = t3p-kp414213562*t3s + t3H = xr2-xr34 ; ty = xr2+xr34 ; t3x = xi2-xi34 ; t2B = xi2+xi34 + t3w = xr18-xr50 ; tB = xr18+xr50 ; t3I = xi18-xi50 ; t2E = xi18+xi50 + t3C = xr58-xr26 ; tI = xr58+xr26 ; t3D = xi58-xi26 ; t2L = xi58+xi26 + t3z = xr10-xr42 ; tF = xr10+xr42 ; t3E = t3C-t3D ; t3K = t3C+t3D + t2I = xi10+xi42 ; t3A = xi10-xi42 ; tat = ty-tB ; tC = ty+tB + tJ = tF+tI ; taq = tI-tF ; t3L = t3A-t3z ; t3B = t3z+t3A + tdd = tC-tJ ; tK = tC+tJ ; tar = t2B-t2E ; t2F = t2B+t2E + tau = t2I-t2L ; t2M = t2I+t2L ; tce = tar-taq ; tas = taq+tar + tcf = tat-tau ; tav = tat+tau ; t7M = t3x-t3w ; t3y = t3w+t3x + t3F = t3B-t3E ; t7Q = t3B+t3E ; tdc = t2F-t2M ; t2N = t2F+t2M + t6G = t3y+kp707106781*t3F ; t3G = t3y-kp707106781*t3F + t7N = t3L+t3K ; t3M = t3K-t3L ; t3J = t3H-t3I ; t7P = t3H+t3I + t9k = t7M-kp707106781*t7N ; t7O = t7M+kp707106781*t7N + t9l = t7P-kp707106781*t7Q ; t7R = t7P+kp707106781*t7Q + t6H = t3J+kp707106781*t3M ; t3N = t3J-kp707106781*t3M + t5I = xr63-xr31 ; t1z = xr63+xr31 ; tb8 = xi63+xi31 ; t56 = xi63-xi31 + t53 = xr15-xr47 ; t1C = xr15+xr47 ; tb9 = xi15+xi47 ; t5L = xi15-xi47 + t5d = xr55-xr23 ; t1J = xr55+xr23 ; t5g = xi55-xi23 ; tbq = xi55+xi23 + t58 = xr7-xr39 ; t1G = xr7+xr39 ; t5N = t5d+t5g ; t5h = t5d-t5g + tbp = xi7+xi39 ; t5b = xi7-xi39 ; tbo = t1z-t1C ; t1D = t1z+t1C + t1K = t1G+t1J ; tb7 = t1J-t1G ; t5c = t58+t5b ; t5O = t5b-t58 + tdA = t1D-t1K ; t1L = t1D+t1K ; tbr = tbp-tbq ; tdw = tbp+tbq + tba = tb8-tb9 ; tdv = tb8+tb9 ; t8l = t56-t53 ; t57 = t53+t56 + tct = tbo-tbr ; tbs = tbo+tbr ; teo = tdv+tdw ; tdx = tdv-tdw + t5i = t5c-t5h ; t8x = t5c+t5h ; t8w = t5I+t5L ; t5M = t5I-t5L + t5P = t5N-t5O ; t8m = t5O+t5N + t6Y = t57+kp707106781*t5i ; t5j = t57-kp707106781*t5i + t6V = t5M+kp707106781*t5P ; t5Q = t5M-kp707106781*t5P + t9z = t8w-kp707106781*t8x ; t8y = t8w+kp707106781*t8x + tcw = tba-tb7 ; tbb = tb7+tba + t9C = t8l-kp707106781*t8m ; t8n = t8l+kp707106781*t8m + t40 = xr62-xr30 ; tN = xr62+xr30 ; t3Q = xi62-xi30 ; t2Q = xi62+xi30 + t3P = xr14-xr46 ; tQ = xr14+xr46 ; t41 = xi14-xi46 ; t2T = xi14+xi46 + t3V = xr54-xr22 ; tX = xr54+xr22 ; t3W = xi54-xi22 ; t30 = xi54+xi22 + t3S = xr6-xr38 ; tU = xr6+xr38 ; t3X = t3V-t3W ; t43 = t3V+t3W + t2X = xi6+xi38 ; t3T = xi6-xi38 ; taA = tN-tQ ; tR = tN+tQ + tY = tU+tX ; tax = tX-tU ; t44 = t3T-t3S ; t3U = t3S+t3T + tdf = tR-tY ; tZ = tR+tY ; tay = t2Q-t2T ; t2U = t2Q+t2T + taB = t2X-t30 ; t31 = t2X+t30 ; tch = tay-tax ; taz = tax+tay + tci = taA-taB ; taC = taA+taB ; t7T = t3Q-t3P ; t3R = t3P+t3Q + t3Y = t3U-t3X ; t7X = t3U+t3X ; tdg = t2U-t31 ; t32 = t2U+t31 + t6J = t3R+kp707106781*t3Y ; t3Z = t3R-kp707106781*t3Y + t7U = t44+t43 ; t45 = t43-t44 ; t42 = t40-t41 ; t7W = t40+t41 + t9n = t7T-kp707106781*t7U ; t7V = t7T+kp707106781*t7U + t9o = t7W-kp707106781*t7X ; t7Y = t7W+kp707106781*t7X + t6K = t42+kp707106781*t45 ; t46 = t42-kp707106781*t45 + t4P = xr1-xr33 ; t14 = xr1+xr33 ; taH = xi1+xi33 ; t4d = xi1-xi33 + t4a = xr17-xr49 ; t17 = xr17+xr49 ; taI = xi17+xi49 ; t4S = xi17-xi49 + t4k = xr57-xr25 ; t1e = xr57+xr25 ; t4n = xi57-xi25 ; taZ = xi57+xi25 + t4f = xr9-xr41 ; t1b = xr9+xr41 ; t4U = t4k+t4n ; t4o = t4k-t4n + taY = xi9+xi41 ; t4i = xi9-xi41 ; taX = t14-t17 ; t18 = t14+t17 + t1f = t1b+t1e ; taG = t1e-t1b ; t4j = t4f+t4i ; t4V = t4i-t4f + tdp = t18-t1f ; t1g = t18+t1f ; tb0 = taY-taZ ; tdl = taY+taZ + taJ = taH-taI ; tdk = taH+taI ; t82 = t4d-t4a ; t4e = t4a+t4d + tcm = taX-tb0 ; tb1 = taX+tb0 ; tej = tdk+tdl ; tdm = tdk-tdl + t4p = t4j-t4o ; t8e = t4j+t4o ; t8d = t4P+t4S ; t4T = t4P-t4S + t4W = t4U-t4V ; t83 = t4V+t4U + t6R = t4e+kp707106781*t4p ; t4q = t4e-kp707106781*t4p + t6O = t4T+kp707106781*t4W ; t4X = t4T-kp707106781*t4W + t9s = t8d-kp707106781*t8e ; t8f = t8d+kp707106781*t8e + tcp = taJ-taG ; taK = taG+taJ + t9v = t82-kp707106781*t83 ; t84 = t82+kp707106781*t83 + t4C = xr5-xr37 ; t1j = xr5+xr37 ; taL = xi5+xi37 ; t4K = xi5-xi37 + t4H = xr21-xr53 ; t1m = xr21+xr53 ; t85 = t4K-t4H ; t4L = t4H+t4K + taO = t1j-t1m ; t1n = t1j+t1m ; t4F = xi21-xi53 ; taM = xi21+xi53 + tdq = taL+taM ; taN = taL-taM ; t86 = t4C+t4F ; t4G = t4C-t4F + t4r = xr61-xr29 ; t1q = xr61+xr29 ; taR = xi61+xi29 ; t4z = xi61-xi29 + t4w = xr13-xr45 ; t1t = xr13+xr45 ; t88 = t4z-t4w ; t4A = t4w+t4z + taQ = t1q-t1t ; t1u = t1q+t1t ; t4u = xi13-xi45 ; taS = xi13+xi45 + tb2 = taO+taN ; taP = taN-taO ; tdr = taR+taS ; taT = taR-taS + t89 = t4r+t4u ; t4v = t4r-t4u ; tdn = t1u-t1n ; t1v = t1n+t1u + tb3 = taQ-taT ; taU = taQ+taT + t4Z = t4A-kp414213562*t4v ; t4B = t4v+kp414213562*t4A + tcq = tb2-tb3 ; tb4 = tb2+tb3 ; tek = tdq+tdr ; tds = tdq-tdr + t4M = t4G-kp414213562*t4L ; t4Y = t4L+kp414213562*t4G + t87 = t85-kp414213562*t86 ; t8g = t86+kp414213562*t85 + t6P = t4M+t4B ; t4N = t4B-t4M ; t6S = t4Y+t4Z ; t50 = t4Y-t4Z + t8h = t89-kp414213562*t88 ; t8a = t88+kp414213562*t89 + t9w = t8g-t8h ; t8i = t8g+t8h ; tcn = taU-taP ; taV = taP+taU + t9t = t8a-t87 ; t8b = t87+t8a ; t5v = xr3-xr35 ; t1O = xr3+xr35 + tbc = xi3+xi35 ; t5D = xi3-xi35 ; t5A = xr19-xr51 ; t1R = xr19+xr51 + t8o = t5D-t5A ; t5E = t5A+t5D ; tbf = t1O-t1R ; t1S = t1O+t1R + t5y = xi19-xi51 ; tbd = xi19+xi51 ; tdB = tbc+tbd ; tbe = tbc-tbd + t8p = t5v+t5y ; t5z = t5v-t5y ; t5k = xr59-xr27 ; t1V = xr59+xr27 + tbi = xi59+xi27 ; t5s = xi59-xi27 ; t5p = xr11-xr43 ; t1Y = xr11+xr43 + t8r = t5s-t5p ; t5t = t5p+t5s ; tbh = t1V-t1Y ; t1Z = t1V+t1Y + t5n = xi11-xi43 ; tbj = xi11+xi43 ; tbt = tbf+tbe ; tbg = tbe-tbf + tdC = tbi+tbj ; tbk = tbi-tbj ; t8s = t5k+t5n ; t5o = t5k-t5n + tdy = t1Z-t1S ; t20 = t1S+t1Z ; tbu = tbh-tbk ; tbl = tbh+tbk + t5S = t5t-kp414213562*t5o ; t5u = t5o+kp414213562*t5t + tcx = tbt-tbu ; tbv = tbt+tbu ; tep = tdB+tdC ; tdD = tdB-tdC + t5F = t5z-kp414213562*t5E ; t5R = t5E+kp414213562*t5z + t8q = t8o-kp414213562*t8p ; t8z = t8p+kp414213562*t8o + t6W = t5F+t5u ; t5G = t5u-t5F ; t6Z = t5R+t5S ; t5T = t5R-t5S + t8A = t8s-kp414213562*t8r ; t8t = t8r+kp414213562*t8s + t9D = t8z-t8A ; t8B = t8z+t8A ; tcu = tbl-tbg ; tbm = tbg+tbl + tef = tf-tu ; tv = tf+tu ; t10 = tK+tZ ; teu = tZ-tK + tel = tej-tek ; teE = tej+tek ; t9A = t8t-t8q ; t8u = t8q+t8t + teD = tv-t10 ; t11 = tv+t10 ; teF = teo+tep ; teq = teo-tep + tei = t1g-t1v ; t1w = t1g+t1v ; t21 = t1L+t20 ; ten = t1L-t20 + tet = t2i-t2x ; t2y = t2i+t2x ; teI = teE+teF ; teG = teE-teF + t23 = t21-t1w ; t22 = t1w+t21 ; t33 = t2N+t32 ; teg = t2N-t32 + t34 = t2y-t33 ; teH = t2y+t33 ; tex = tef-teg ; teh = tef+teg + teB = teu+tet ; tev = tet-teu ; tey = tel-tei ; tem = tei+tel + tdV = td9+tda ; tdb = td9-tda ; tdJ = tdH-tdI ; te5 = tdI+tdH + tez = ten+teq ; ter = ten-teq ; tdL = tdd+tdc ; tde = tdc-tdd + teA = tey-tez ; teC = tey+tez ; tew = ter-tem ; tes = tem+ter + tdh = tdf+tdg ; tdK = tdf-tdg ; tdE = tdA-tdD ; te1 = tdA+tdD + te2 = tdy+tdx ; tdz = tdx-tdy ; te6 = tde+tdh ; tdi = tde-tdh + teb = te2+kp414213562*te1 ; te3 = te1-kp414213562*te2 + tdZ = tdn+tdm ; tdo = tdm-tdn ; tdt = tdp-tds ; tdY = tdp+tds + tdW = tdL+tdK ; tdM = tdK-tdL + tdR = tdt-kp414213562*tdo ; tdu = tdo+kp414213562*tdt + tdT = tdb-kp707106781*tdi ; tdj = tdb+kp707106781*tdi + tea = tdZ-kp414213562*tdY ; te0 = tdY+kp414213562*tdZ + tdQ = tdE+kp414213562*tdz ; tdF = tdz-kp414213562*tdE + tdP = tdJ+kp707106781*tdM ; tdN = tdJ-kp707106781*tdM + tdS = tdQ-tdR ; tdU = tdR+tdQ ; tdO = tdu+tdF ; tdG = tdu-tdF + te9 = tdV-kp707106781*tdW ; tdX = tdV+kp707106781*tdW + te4 = te0+te3 ; te8 = te3-te0 + te7 = te5-kp707106781*te6 ; ted = te5+kp707106781*te6 + tee = tea+teb ; tec = tea-teb ; tbE = tbC+tbD ; tcc = tbC-tbD + tcC = tan-tak ; tao = tak+tan + tcF = tcf-kp414213562*tce ; tcg = tce+kp414213562*tcf + tcP = tcb-kp707106781*tcc ; tcd = tcb+kp707106781*tcc + tcZ = tcB-kp707106781*tcC ; tcD = tcB+kp707106781*tcC + tcj = tch-kp414213562*tci ; tcE = tci+kp414213562*tch + tcy = tcw-kp707106781*tcx ; tcV = tcw+kp707106781*tcx + tcW = tct+kp707106781*tcu ; tcv = tct-kp707106781*tcu + tcT = tcm+kp707106781*tcn ; tco = tcm-kp707106781*tcn + td0 = tcg+tcj ; tck = tcg-tcj + td4 = tcW+kp198912367*tcV ; tcX = tcV-kp198912367*tcW + tcr = tcp-kp707106781*tcq ; tcS = tcp+kp707106781*tcq + tcK = tcr-kp668178637*tco ; tcs = tco+kp668178637*tcr + tcQ = tcF+tcE ; tcG = tcE-tcF + tcJ = tcd-kp923879532*tck ; tcl = tcd+kp923879532*tck + td5 = tcT-kp198912367*tcS ; tcU = tcS+kp198912367*tcT + tcL = tcy+kp668178637*tcv ; tcz = tcv-kp668178637*tcy + tcN = tcD+kp923879532*tcG ; tcH = tcD-kp923879532*tcG + tcO = tcK+tcL ; tcM = tcK-tcL ; tcI = tcz-tcs ; tcA = tcs+tcz + td7 = tcP+kp923879532*tcQ ; tcR = tcP-kp923879532*tcQ + tcY = tcU-tcX ; td2 = tcU+tcX + td1 = tcZ+kp923879532*td0 ; td3 = tcZ-kp923879532*td0 + td6 = td4-td5 ; td8 = td5+td4 + tbH = tav+kp414213562*tas ; taw = tas-kp414213562*tav + tbR = tah+kp707106781*tao ; tap = tah-kp707106781*tao + tc1 = tbB+kp707106781*tbE ; tbF = tbB-kp707106781*tbE + taD = taz+kp414213562*taC ; tbG = taC-kp414213562*taz + tbw = tbs-kp707106781*tbv ; tbX = tbs+kp707106781*tbv + tbY = tbb+kp707106781*tbm ; tbn = tbb-kp707106781*tbm + tbV = taK+kp707106781*taV ; taW = taK-kp707106781*taV + tc2 = taw+taD ; taE = taw-taD + tc7 = tbY+kp198912367*tbX ; tbZ = tbX-kp198912367*tbY + tb5 = tb1-kp707106781*tb4 ; tbU = tb1+kp707106781*tb4 + tbN = tb5-kp668178637*taW ; tb6 = taW+kp668178637*tb5 + tbS = tbH+tbG ; tbI = tbG-tbH + tbP = tap-kp923879532*taE ; taF = tap+kp923879532*taE + tc6 = tbV-kp198912367*tbU ; tbW = tbU+kp198912367*tbV + tbM = tbw+kp668178637*tbn ; tbx = tbn-kp668178637*tbw + tbL = tbF+kp923879532*tbI ; tbJ = tbF-kp923879532*tbI + tbO = tbM-tbN ; tbQ = tbN+tbM ; tbK = tb6+tbx ; tby = tb6-tbx + tc5 = tbR-kp923879532*tbS ; tbT = tbR+kp923879532*tbS + tc0 = tbW+tbZ ; tc4 = tbZ-tbW + tc3 = tc1-kp923879532*tc2 ; tc9 = tc1+kp923879532*tc2 + tca = tc6+tc7 ; tc8 = tc6-tc7 + t3f = t37+kp707106781*t3e ; t6D = t37-kp707106781*t3e + t6E = t65+t64 ; t66 = t64-t65 + t6T = t6R-kp923879532*t6S ; t7k = t6R+kp923879532*t6S + t7h = t6D+kp923879532*t6E ; t6F = t6D-kp923879532*t6E + t7l = t6O+kp923879532*t6P ; t6Q = t6O-kp923879532*t6P + t70 = t6Y-kp923879532*t6Z ; t7n = t6Y+kp923879532*t6Z + t7o = t6V+kp923879532*t6W ; t6X = t6V-kp923879532*t6W + t77 = t6H-kp198912367*t6G ; t6I = t6G+kp198912367*t6H + t7x = t7l-kp098491403*t7k ; t7m = t7k+kp098491403*t7l + t7w = t7o+kp098491403*t7n ; t7p = t7n-kp098491403*t7o + t6L = t6J-kp198912367*t6K ; t76 = t6K+kp198912367*t6J + t63 = t5Z+kp707106781*t62 ; t73 = t5Z-kp707106781*t62 + t7s = t6I+t6L ; t6M = t6I-t6L + t7c = t6T-kp820678790*t6Q ; t6U = t6Q+kp820678790*t6T + t74 = t3m+t3t ; t3u = t3m-t3t + t7r = t73+kp923879532*t74 ; t75 = t73-kp923879532*t74 + t7i = t77+t76 ; t78 = t76-t77 + t7b = t6F-kp980785280*t6M ; t6N = t6F+kp980785280*t6M + t7f = t75+kp980785280*t78 ; t79 = t75-kp980785280*t78 + t71 = t6X-kp820678790*t70 ; t7d = t70+kp820678790*t6X + t7z = t7h+kp980785280*t7i ; t7j = t7h-kp980785280*t7i + t7q = t7m-t7p ; t7u = t7m+t7p ; t7g = t7c+t7d ; t7e = t7c-t7d + t72 = t6U+t71 ; t7a = t71-t6U + t7t = t7r+kp980785280*t7s ; t7v = t7r-kp980785280*t7s + t7y = t7w-t7x ; t7A = t7x+t7w + t7D = t7B+kp707106781*t7C ; t9h = t7B-kp707106781*t7C + t9i = t8I-t8J ; t8K = t8I+t8J + t9x = t9v-kp923879532*t9w ; t9Y = t9v+kp923879532*t9w + t9V = t9h-kp923879532*t9i ; t9j = t9h+kp923879532*t9i + t9Z = t9s+kp923879532*t9t ; t9u = t9s-kp923879532*t9t + t9E = t9C-kp923879532*t9D ; ta1 = t9C+kp923879532*t9D + ta2 = t9z+kp923879532*t9A ; t9B = t9z-kp923879532*t9A + t9L = t9l-kp668178637*t9k ; t9m = t9k+kp668178637*t9l + tab = t9Z-kp303346683*t9Y ; ta0 = t9Y+kp303346683*t9Z + taa = ta2+kp303346683*ta1 ; ta3 = ta1-kp303346683*ta2 + t9p = t9n-kp668178637*t9o ; t9K = t9o+kp668178637*t9n + t8H = t8F+kp707106781*t8G ; t9H = t8F-kp707106781*t8G + ta6 = t9m+t9p ; t9q = t9m-t9p + t9Q = t9x-kp534511135*t9u ; t9y = t9u+kp534511135*t9x + t9I = t7J-t7G ; t7K = t7G+t7J + ta5 = t9H-kp923879532*t9I ; t9J = t9H+kp923879532*t9I + t9W = t9L+t9K ; t9M = t9K-t9L + t9P = t9j-kp831469612*t9q ; t9r = t9j+kp831469612*t9q + t9T = t9J+kp831469612*t9M ; t9N = t9J-kp831469612*t9M + t9F = t9B-kp534511135*t9E ; t9R = t9E+kp534511135*t9B + tad = t9V+kp831469612*t9W ; t9X = t9V-kp831469612*t9W + ta4 = ta0-ta3 ; ta8 = ta0+ta3 ; t9U = t9Q+t9R ; t9S = t9Q-t9R + t9G = t9y+t9F ; t9O = t9F-t9y + ta7 = ta5+kp831469612*ta6 ; ta9 = ta5-kp831469612*ta6 + tac = taa-tab ; tae = tab+taa + t51 = t4X-kp923879532*t50 ; t6m = t4X+kp923879532*t50 + t6j = t3f+kp923879532*t3u ; t3v = t3f-kp923879532*t3u + t6n = t4q+kp923879532*t4N ; t4O = t4q-kp923879532*t4N + t5U = t5Q-kp923879532*t5T ; t6p = t5Q+kp923879532*t5T + t6q = t5j+kp923879532*t5G ; t5H = t5j-kp923879532*t5G + t69 = t3N+kp668178637*t3G ; t3O = t3G-kp668178637*t3N + t6y = t6n-kp303346683*t6m ; t6o = t6m+kp303346683*t6n + t6z = t6q+kp303346683*t6p ; t6r = t6p-kp303346683*t6q + t47 = t3Z+kp668178637*t46 ; t68 = t46-kp668178637*t3Z + t6u = t3O+t47 ; t48 = t3O-t47 + t6f = t51-kp534511135*t4O ; t52 = t4O+kp534511135*t51 + t6t = t63+kp923879532*t66 ; t67 = t63-kp923879532*t66 + t6k = t69+t68 ; t6a = t68-t69 + t6h = t3v-kp831469612*t48 ; t49 = t3v+kp831469612*t48 + t6d = t67+kp831469612*t6a ; t6b = t67-kp831469612*t6a + t5V = t5H-kp534511135*t5U ; t6e = t5U+kp534511135*t5H + t6x = t6j-kp831469612*t6k ; t6l = t6j+kp831469612*t6k + t6s = t6o+t6r ; t6w = t6r-t6o ; t6g = t6e-t6f ; t6i = t6f+t6e + t5W = t52-t5V ; t6c = t52+t5V + t6v = t6t-kp831469612*t6u ; t6B = t6t+kp831469612*t6u + t6C = t6y+t6z ; t6A = t6y-t6z + t8j = t8f-kp923879532*t8i ; t90 = t8f+kp923879532*t8i + t8X = t7D+kp923879532*t7K ; t7L = t7D-kp923879532*t7K + t91 = t84+kp923879532*t8b ; t8c = t84-kp923879532*t8b + t8C = t8y-kp923879532*t8B ; t93 = t8y+kp923879532*t8B + t94 = t8n+kp923879532*t8u ; t8v = t8n-kp923879532*t8u + t8N = t7R+kp198912367*t7O ; t7S = t7O-kp198912367*t7R + t9c = t91-kp098491403*t90 ; t92 = t90+kp098491403*t91 + t9d = t94+kp098491403*t93 ; t95 = t93-kp098491403*t94 + t7Z = t7V+kp198912367*t7Y ; t8M = t7Y-kp198912367*t7V + t98 = t7S+t7Z ; t80 = t7S-t7Z + t8T = t8j-kp820678790*t8c ; t8k = t8c+kp820678790*t8j + t97 = t8H+kp923879532*t8K ; t8L = t8H-kp923879532*t8K + t8Y = t8N+t8M ; t8O = t8M-t8N + t8V = t7L-kp980785280*t80 ; t81 = t7L+kp980785280*t80 + t8R = t8L+kp980785280*t8O ; t8P = t8L-kp980785280*t8O + t8D = t8v-kp820678790*t8C ; t8S = t8C+kp820678790*t8v + t9b = t8X-kp980785280*t8Y ; t8Z = t8X+kp980785280*t8Y + t96 = t92+t95 ; t9a = t95-t92 ; t8U = t8S-t8T ; t8W = t8T+t8S + t8E = t8k-t8D ; t8Q = t8k+t8D + t99 = t97-kp980785280*t98 ; t9f = t97+kp980785280*t98 + t9g = t9c+t9d ; t9e = t9c-t9d + r63 = (t8Z+kp995184726*t96) :+ (t9f+kp995184726*t9g) + r62 = (tbT+kp980785280*tc0) :+ (tc9+kp980785280*tca) + r61 = (t6l+kp956940335*t6s) :+ (t6B+kp956940335*t6C) + r60 = (tdX+kp923879532*te4) :+ (ted+kp923879532*tee) + r59 = (t9r+kp881921264*t9G) :+ (t9T+kp881921264*t9U) + r58 = (tcl+kp831469612*tcA) :+ (tcN+kp831469612*tcO) + r57 = (t6N+kp773010453*t72) :+ (t7f+kp773010453*t7g) + r56 = (teh+kp707106781*tes) :+ (teB+kp707106781*teC) + r55 = (t81+kp773010453*t8E) :+ (t8R+kp773010453*t8U) + r54 = (taF+kp831469612*tby) :+ (tbL+kp831469612*tbO) + r53 = (t49+kp881921264*t5W) :+ (t6d+kp881921264*t6g) + r52 = (tdj+kp923879532*tdG) :+ (tdP+kp923879532*tdS) + r51 = (t9X+kp956940335*ta4) :+ (ta9+kp956940335*tac) + r50 = (tcR+kp980785280*tcY) :+ (td3+kp980785280*td6) + r49 = (t7j+kp995184726*t7q) :+ (t7v+kp995184726*t7y) + r48 = (teD+teG) :+ (t23+t34) + r47 = (t9b+kp995184726*t9e) :+ (t99+kp995184726*t9a) + r46 = (tc5+kp980785280*tc8) :+ (tc3+kp980785280*tc4) + r45 = (t6x+kp956940335*t6A) :+ (t6v+kp956940335*t6w) + r44 = (te9+kp923879532*tec) :+ (te7+kp923879532*te8) + r43 = (t9P+kp881921264*t9S) :+ (t9N+kp881921264*t9O) + r42 = (tcJ+kp831469612*tcM) :+ (tcH+kp831469612*tcI) + r41 = (t7b+kp773010453*t7e) :+ (t79+kp773010453*t7a) + r40 = (tex+kp707106781*teA) :+ (tev+kp707106781*tew) + r39 = (t8V-kp773010453*t8W) :+ (t8P-kp773010453*t8Q) + r38 = (tbP-kp831469612*tbQ) :+ (tbJ-kp831469612*tbK) + r37 = (t6h-kp881921264*t6i) :+ (t6b-kp881921264*t6c) + r36 = (tdT-kp923879532*tdU) :+ (tdN-kp923879532*tdO) + r35 = (tad-kp956940335*tae) :+ (ta7-kp956940335*ta8) + r34 = (td7-kp980785280*td8) :+ (td1-kp980785280*td2) + r33 = (t7z-kp995184726*t7A) :+ (t7t-kp995184726*t7u) + r32 = (t11-t22) :+ (teH-teI) + r31 = (t8Z-kp995184726*t96) :+ (t9f-kp995184726*t9g) + r30 = (tbT-kp980785280*tc0) :+ (tc9-kp980785280*tca) + r29 = (t6l-kp956940335*t6s) :+ (t6B-kp956940335*t6C) + r28 = (tdX-kp923879532*te4) :+ (ted-kp923879532*tee) + r27 = (t9r-kp881921264*t9G) :+ (t9T-kp881921264*t9U) + r26 = (tcl-kp831469612*tcA) :+ (tcN-kp831469612*tcO) + r25 = (t6N-kp773010453*t72) :+ (t7f-kp773010453*t7g) + r24 = (teh-kp707106781*tes) :+ (teB-kp707106781*teC) + r23 = (t81-kp773010453*t8E) :+ (t8R-kp773010453*t8U) + r22 = (taF-kp831469612*tby) :+ (tbL-kp831469612*tbO) + r21 = (t49-kp881921264*t5W) :+ (t6d-kp881921264*t6g) + r20 = (tdj-kp923879532*tdG) :+ (tdP-kp923879532*tdS) + r19 = (t9X-kp956940335*ta4) :+ (ta9-kp956940335*tac) + r18 = (tcR-kp980785280*tcY) :+ (td3-kp980785280*td6) + r17 = (t7j-kp995184726*t7q) :+ (t7v-kp995184726*t7y) + r16 = (teD-teG) :+ (t34-t23) + r15 = (t9b-kp995184726*t9e) :+ (t99-kp995184726*t9a) + r14 = (tc5-kp980785280*tc8) :+ (tc3-kp980785280*tc4) + r13 = (t6x-kp956940335*t6A) :+ (t6v-kp956940335*t6w) + r12 = (te9-kp923879532*tec) :+ (te7-kp923879532*te8) + r11 = (t9P-kp881921264*t9S) :+ (t9N-kp881921264*t9O) + r10 = (tcJ-kp831469612*tcM) :+ (tcH-kp831469612*tcI) + r9 = (t7b-kp773010453*t7e) :+ (t79-kp773010453*t7a) + r8 = (tex-kp707106781*teA) :+ (tev-kp707106781*tew) + r7 = (t8V+kp773010453*t8W) :+ (t8P+kp773010453*t8Q) + r6 = (tbP+kp831469612*tbQ) :+ (tbJ+kp831469612*tbK) + r5 = (t6h+kp881921264*t6i) :+ (t6b+kp881921264*t6c) + r4 = (tdT+kp923879532*tdU) :+ (tdN+kp923879532*tdO) + r3 = (tad+kp956940335*tae) :+ (ta7+kp956940335*ta8) + r2 = (td7+kp980785280*td8) :+ (td1+kp980785280*td2) + r1 = (t7z+kp995184726*t7A) :+ (t7t+kp995184726*t7u) + MV.unsafeWrite xsout 0 $ (t11+t22) :+ (teH+teI) + MV.unsafeWrite xsout 1 $ if sign == 1 then r1 else r63 + MV.unsafeWrite xsout 2 $ if sign == 1 then r2 else r62 + MV.unsafeWrite xsout 3 $ if sign == 1 then r3 else r61 + MV.unsafeWrite xsout 4 $ if sign == 1 then r4 else r60 + MV.unsafeWrite xsout 5 $ if sign == 1 then r5 else r59 + MV.unsafeWrite xsout 6 $ if sign == 1 then r6 else r58 + MV.unsafeWrite xsout 7 $ if sign == 1 then r7 else r57 + MV.unsafeWrite xsout 8 $ if sign == 1 then r8 else r56 + MV.unsafeWrite xsout 9 $ if sign == 1 then r9 else r55 + MV.unsafeWrite xsout 10 $ if sign == 1 then r10 else r54 + MV.unsafeWrite xsout 11 $ if sign == 1 then r11 else r53 + MV.unsafeWrite xsout 12 $ if sign == 1 then r12 else r52 + MV.unsafeWrite xsout 13 $ if sign == 1 then r13 else r51 + MV.unsafeWrite xsout 14 $ if sign == 1 then r14 else r50 + MV.unsafeWrite xsout 15 $ if sign == 1 then r15 else r49 + MV.unsafeWrite xsout 16 $ if sign == 1 then r16 else r48 + MV.unsafeWrite xsout 17 $ if sign == 1 then r17 else r47 + MV.unsafeWrite xsout 18 $ if sign == 1 then r18 else r46 + MV.unsafeWrite xsout 19 $ if sign == 1 then r19 else r45 + MV.unsafeWrite xsout 20 $ if sign == 1 then r20 else r44 + MV.unsafeWrite xsout 21 $ if sign == 1 then r21 else r43 + MV.unsafeWrite xsout 22 $ if sign == 1 then r22 else r42 + MV.unsafeWrite xsout 23 $ if sign == 1 then r23 else r41 + MV.unsafeWrite xsout 24 $ if sign == 1 then r24 else r40 + MV.unsafeWrite xsout 25 $ if sign == 1 then r25 else r39 + MV.unsafeWrite xsout 26 $ if sign == 1 then r26 else r38 + MV.unsafeWrite xsout 27 $ if sign == 1 then r27 else r37 + MV.unsafeWrite xsout 28 $ if sign == 1 then r28 else r36 + MV.unsafeWrite xsout 29 $ if sign == 1 then r29 else r35 + MV.unsafeWrite xsout 30 $ if sign == 1 then r30 else r34 + MV.unsafeWrite xsout 31 $ if sign == 1 then r31 else r33 + MV.unsafeWrite xsout 32 $ if sign == 1 then r32 else r32 + MV.unsafeWrite xsout 33 $ if sign == 1 then r33 else r31 + MV.unsafeWrite xsout 34 $ if sign == 1 then r34 else r30 + MV.unsafeWrite xsout 35 $ if sign == 1 then r35 else r29 + MV.unsafeWrite xsout 36 $ if sign == 1 then r36 else r28 + MV.unsafeWrite xsout 37 $ if sign == 1 then r37 else r27 + MV.unsafeWrite xsout 38 $ if sign == 1 then r38 else r26 + MV.unsafeWrite xsout 39 $ if sign == 1 then r39 else r25 + MV.unsafeWrite xsout 40 $ if sign == 1 then r40 else r24 + MV.unsafeWrite xsout 41 $ if sign == 1 then r41 else r23 + MV.unsafeWrite xsout 42 $ if sign == 1 then r42 else r22 + MV.unsafeWrite xsout 43 $ if sign == 1 then r43 else r21 + MV.unsafeWrite xsout 44 $ if sign == 1 then r44 else r20 + MV.unsafeWrite xsout 45 $ if sign == 1 then r45 else r19 + MV.unsafeWrite xsout 46 $ if sign == 1 then r46 else r18 + MV.unsafeWrite xsout 47 $ if sign == 1 then r47 else r17 + MV.unsafeWrite xsout 48 $ if sign == 1 then r48 else r16 + MV.unsafeWrite xsout 49 $ if sign == 1 then r49 else r15 + MV.unsafeWrite xsout 50 $ if sign == 1 then r50 else r14 + MV.unsafeWrite xsout 51 $ if sign == 1 then r51 else r13 + MV.unsafeWrite xsout 52 $ if sign == 1 then r52 else r12 + MV.unsafeWrite xsout 53 $ if sign == 1 then r53 else r11 + MV.unsafeWrite xsout 54 $ if sign == 1 then r54 else r10 + MV.unsafeWrite xsout 55 $ if sign == 1 then r55 else r9 + MV.unsafeWrite xsout 56 $ if sign == 1 then r56 else r8 + MV.unsafeWrite xsout 57 $ if sign == 1 then r57 else r7 + MV.unsafeWrite xsout 58 $ if sign == 1 then r58 else r6 + MV.unsafeWrite xsout 59 $ if sign == 1 then r59 else r5 + MV.unsafeWrite xsout 60 $ if sign == 1 then r60 else r4 + MV.unsafeWrite xsout 61 $ if sign == 1 then r61 else r3 + MV.unsafeWrite xsout 62 $ if sign == 1 then r62 else r2 + MV.unsafeWrite xsout 63 $ if sign == 1 then r63 else r1 diff --git a/Numeric/FFT/Special/Primes.hs b/Numeric/FFT/Special/Primes.hs new file mode 100644 --- /dev/null +++ b/Numeric/FFT/Special/Primes.hs @@ -0,0 +1,339 @@ +module Numeric.FFT.Special.Primes + ( special3, special5, special7, special11, special13 ) where + +import Control.Monad.ST +import Data.IntMap.Strict (IntMap) +import qualified Data.IntMap.Strict as IM +import Data.Complex +import Data.Vector.Unboxed +import qualified Data.Vector.Unboxed.Mutable as MV + +import Numeric.FFT.Types +import Numeric.FFT.Utils + + +-- | Length 3 hard-coded FFT. +kp500000000, kp866025403 :: Double +kp866025403 = 0.866025403784438646763723170752936183471402627 +kp500000000 = 0.500000000000000000000000000000000000000000000 +special3 :: Int -> MVCD s -> MVCD s -> ST s () +special3 sign xsin xsout = do + xr0 :+ xi0 <- MV.unsafeRead xsin 0 ; xr1 :+ xi1 <- MV.unsafeRead xsin 1 + xr2 :+ xi2 <- MV.unsafeRead xsin 2 + let rp = xr1 + xr2 ; rm = xr1 - xr2 + ip = xi1 + xi2 ; im = xi1 - xi2 + tr = xr0 - kp500000000 * rp + ti = xi0 - kp500000000 * ip + r1 = (tr - kp866025403 * im) :+ (ti + kp866025403 * rm) + r2 = (tr + kp866025403 * im) :+ (ti - kp866025403 * rm) + MV.unsafeWrite xsout 0 $ (xr0 + rp) :+ (xi0 + ip) + MV.unsafeWrite xsout 1 $ if sign == 1 then r1 else r2 + MV.unsafeWrite xsout 2 $ if sign == 1 then r2 else r1 + +-- | Length 5 hard-coded FFT. +kp951056516, kp559016994, kp250000000, kp618033988 :: Double +kp951056516 = 0.951056516295153572116439333379382143405698634 +kp559016994 = 0.559016994374947424102293417182819058860154590 +kp250000000 = 0.250000000000000000000000000000000000000000000 +kp618033988 = 0.618033988749894848204586834365638117720309180 +special5 :: Int -> MVCD s -> MVCD s -> ST s () +special5 sign xsin xsout = do + xr0 :+ xi0 <- MV.unsafeRead xsin 0 ; xr1 :+ xi1 <- MV.unsafeRead xsin 1 + xr2 :+ xi2 <- MV.unsafeRead xsin 2 ; xr3 :+ xi3 <- MV.unsafeRead xsin 3 + xr4 :+ xi4 <- MV.unsafeRead xsin 4 + let ts = xr1 - xr4 ; t4 = xr1 + xr4 ; tt = xr2 - xr3 ; t7 = xr2 + xr3 + t8 = t4 + t7 ; ta = t4 - t7 ; te = xi1 - xi4 ; tm = xi1 + xi4 + tn = xi2 + xi3 ; th = xi2 - xi3 ; to = tm + tn ; tq = tm - tn + ti = te + kp618033988 * th ; tk = th - kp618033988 * te + t9 = xr0 - kp250000000 * t8 ; tu = ts + kp618033988 * tt + tw = tt - kp618033988 * ts ; tp = xi0 - kp250000000 * to + tb = t9 + kp559016994 * ta ; tj = t9 - kp559016994 * ta + tr = tp + kp559016994 * tq ; tv = tp - kp559016994 * tq + r4 = (tb + kp951056516 * ti) :+ (tr - kp951056516 * tu) + r3 = (tj - kp951056516 * tk) :+ (tv + kp951056516 * tw) + r2 = (tj + kp951056516 * tk) :+ (tv - kp951056516 * tw) + r1 = (tb - kp951056516 * ti) :+ (tr + kp951056516 * tu) + MV.unsafeWrite xsout 0 $ (xr0 + t8) :+ (xi0 + to) + MV.unsafeWrite xsout 1 $ if sign == 1 then r1 else r4 + MV.unsafeWrite xsout 2 $ if sign == 1 then r2 else r3 + MV.unsafeWrite xsout 3 $ if sign == 1 then r3 else r2 + MV.unsafeWrite xsout 4 $ if sign == 1 then r4 else r1 + +-- | Length 7 hard-coded FFT. +kp974927912, kp900968867, kp801937735 :: Double +kp692021471, kp356895867, kp554958132 :: Double +kp974927912 = 0.974927912181823607018131682993931217232785801 +kp900968867 = 0.900968867902419126236102319507445051165919162 +kp801937735 = 0.801937735804838252472204639014890102331838324 +kp692021471 = 0.692021471630095869627814897002069140197260599 +kp356895867 = 0.356895867892209443894399510021300583399127187 +kp554958132 = 0.554958132087371191422194871006410481067288862 +special7 :: Int -> MVCD s -> MVCD s -> ST s () +special7 sign xsin xsout = do + xr0 :+ xi0 <- MV.unsafeRead xsin 0 ; xr1 :+ xi1 <- MV.unsafeRead xsin 1 + xr2 :+ xi2 <- MV.unsafeRead xsin 2 ; xr3 :+ xi3 <- MV.unsafeRead xsin 3 + xr4 :+ xi4 <- MV.unsafeRead xsin 4 ; xr5 :+ xi5 <- MV.unsafeRead xsin 5 + xr6 :+ xi6 <- MV.unsafeRead xsin 6 + let tI = xr6 - xr1 ; t4 = xr1 + xr6 ; tG = xr4 - xr3 ; ta = xr3 + xr4 + tT = tI + kp554958132 * tG ; tp = ta - kp356895867 * t4 + tH = xr5 - xr2 ; t7 = xr2 + xr5 + tJ = tH - kp554958132 * tI ; tO = tG + kp554958132 * tH + tu = t7 - kp356895867 * ta ; tb = t4 - kp356895867 * t7 + tB = xi2 + xi5 ; tg = xi2 - xi5 + tC = xi3 + xi4 ; tm = xi3 - xi4 ; tA = xi1 + xi6 ; tj = xi1 - xi6 + tD = tB - kp356895867 * tC ; ts = tm + kp554958132 * tg + tL = tC - kp356895867 * tA ; tQ = tA - kp356895867 * tB + tx = tg - kp554958132 * tj ; tn = tj + kp554958132 * tm + tc = ta - kp692021471 * tb ; tU = tH + kp801937735 * tT + to = tg + kp801937735 * tn ; tR = tC - kp692021471 * tQ + td = xr0 - kp900968867 * tc ; tt = tj - kp801937735 * ts + tq = t7 - kp692021471 * tp ; tS = xi0 - kp900968867 * tR + tr = xr0 - kp900968867 * tq ; tP = tI - kp801937735 * tO + tM = tB - kp692021471 * tL ; ty = tm - kp801937735 * tx + tv = t4 - kp692021471 * tu ; tK = tG - kp801937735 * tJ + tN = xi0 - kp900968867 * tM ; tE = tA - kp692021471 * tD + tw = xr0 - kp900968867 * tv ; tF = xi0 - kp900968867 * tE + r6 = (td + kp974927912 * to) :+ (tS + kp974927912 * tU) + r5 = (tr + kp974927912 * tt) :+ (tN + kp974927912 * tP) + r4 = (tw + kp974927912 * ty) :+ (tF + kp974927912 * tK) + r3 = (tw - kp974927912 * ty) :+ (tF - kp974927912 * tK) + r2 = (tr - kp974927912 * tt) :+ (tN - kp974927912 * tP) + r1 = (td - kp974927912 * to) :+ (tS - kp974927912 * tU) + MV.unsafeWrite xsout 0 $ (xr0 + t4 + t7 + ta) :+ (xi0 + tA + tB + tC) + MV.unsafeWrite xsout 1 $ if sign == 1 then r1 else r6 + MV.unsafeWrite xsout 2 $ if sign == 1 then r2 else r5 + MV.unsafeWrite xsout 3 $ if sign == 1 then r3 else r4 + MV.unsafeWrite xsout 4 $ if sign == 1 then r4 else r3 + MV.unsafeWrite xsout 5 $ if sign == 1 then r5 else r2 + MV.unsafeWrite xsout 6 $ if sign == 1 then r6 else r1 + +-- | Length 11 hard-coded FFT. +kp989821441, kp959492973, kp918985947, kp876768831, kp830830026 :: Double +kp778434453, kp715370323, kp634356270, kp342584725, kp521108558 :: Double +kp989821441 = 0.989821441880932732376092037776718787376519372 +kp959492973 = 0.959492973614497389890368057066327699062454848 +kp918985947 = 0.918985947228994779780736114132655398124909697 +kp876768831 = 0.876768831002589333891339807079336796764054852 +kp830830026 = 0.830830026003772851058548298459246407048009821 +kp778434453 = 0.778434453334651800608337670740821884709317477 +kp715370323 = 0.715370323453429719112414662767260662417897278 +kp634356270 = 0.634356270682424498893150776899916060542806975 +kp342584725 = 0.342584725681637509502641509861112333758894680 +kp521108558 = 0.521108558113202722944698153526659300680427422 +special11 :: Int -> MVCD s -> MVCD s -> ST s () +special11 sign xsin xsout = do + xr0 :+ xi0 <- MV.unsafeRead xsin 0 ; xr1 :+ xi1 <- MV.unsafeRead xsin 1 + xr2 :+ xi2 <- MV.unsafeRead xsin 2 ; xr3 :+ xi3 <- MV.unsafeRead xsin 3 + xr4 :+ xi4 <- MV.unsafeRead xsin 4 ; xr5 :+ xi5 <- MV.unsafeRead xsin 5 + xr6 :+ xi6 <- MV.unsafeRead xsin 6 ; xr7 :+ xi7 <- MV.unsafeRead xsin 7 + xr8 :+ xi8 <- MV.unsafeRead xsin 8 ; xr9 :+ xi9 <- MV.unsafeRead xsin 9 + xr10 :+ xi10 <- MV.unsafeRead xsin 10 + let t1u = xr10 - xr1 ; t4 = xr1 + xr10 ; t1q = xr6 - xr5 ; tg = xr5 + xr6; + t1t = xr9 - xr2 ; t7 = xr2 + xr9 ; t1s = xr8 - xr3 ; ta = xr3 + xr8; + t25 = t1u + kp521108558 * t1q ; t1W = t1q + kp521108558 * t1s + tO = ta - kp342584725 * t4 ; th = t7 - kp342584725 * ta + td = xr4 + xr7 ; t1r = xr7 - xr4 + tP = tg - kp634356270 * tO ; t1X = t1t - kp715370323 * t1W + t26 = t1r + kp715370323 * t25 ; tF = t4 - kp342584725 * td + ti = td - kp634356270 * th ; t1N = t1r - kp521108558 * t1t + t1v = t1t - kp521108558 * t1u ; tG = t7 - kp634356270 * tF + tX = tg - kp342584725 * t7 ; t1O = t1q + kp715370323 * t1N + t1w = t1s - kp715370323 * t1v ; t1E = t1s + kp521108558 * t1r + tY = t4 - kp634356270 * tX ; t16 = td - kp342584725 * tg + t1F = t1u + kp715370323 * t1E ; t17 = ta - kp634356270 * t16 + to = xi3 - xi8 ; t1i = xi3 + xi8 ; t1k = xi5 + xi6 ; tA = xi5 - xi6 + t1h = xi2 + xi9 ; tr = xi2 - xi9 ; t1j = xi4 + xi7 ; tu = xi4 - xi7 + t20 = t1h - kp342584725 * t1i ; tK = tA + kp521108558 * to + tT = tu - kp521108558 * tr ; t1g = xi1 + xi10 ; tx = xi1 - xi10 + t21 = t1j - kp634356270 * t20 ; tU = tA + kp715370323 * tT + tL = tr - kp715370323 * tK ; tB = tx + kp521108558 * tA + t1R = t1g - kp342584725 * t1j ; t1I = t1i - kp342584725 * t1g + t1l = t1j - kp342584725 * t1k ; tC = tu + kp715370323 * tB + t1S = t1h - kp634356270 * t1R ; t1J = t1k - kp634356270 * t1I + t1m = t1i - kp634356270 * t1l ; t12 = to + kp521108558 * tu + t1z = t1k - kp342584725 * t1h ; t1b = tr - kp521108558 * tx + t13 = tx + kp715370323 * t12 ; t1A = t1g - kp634356270 * t1z + t1c = to - kp715370323 * t1b ; tj = t4 - kp778434453 * ti + tD = tr + kp830830026 * tC ; t22 = t1g - kp778434453 * t21 + t27 = t1t + kp830830026 * t26 ; tk = tg - kp876768831 * tj + tE = to + kp918985947 * tD ; t23 = t1k - kp876768831 * t22 + t28 = t1s + kp918985947 * t27 ; tl = xr0 - kp959492973 * tk + t1T = t1k - kp778434453 * t1S ; t24 = xi0 - kp959492973 * t23 + t1Y = t1u + kp830830026 * t1X ; t1U = t1i - kp876768831 * t1T + t1Z = t1r - kp918985947 * t1Y ; t1V = xi0 - kp959492973 * t1U + tH = tg - kp778434453 * tG ; tM = tx + kp830830026 * tL + tQ = td - kp778434453 * tP ; tI = ta - kp876768831 * tH + tN = tu - kp918985947 * tM ; tR = t7 - kp876768831 * tQ + tV = to - kp830830026 * tU ; tJ = xr0 - kp959492973 * tI + t1K = t1j - kp778434453 * t1J ; tS = xr0 - kp959492973 * tR + tW = tx - kp918985947 * tV ; t1L = t1h - kp876768831 * t1K + t1P = t1s - kp830830026 * t1O ; t1M = xi0 - kp959492973 * t1L + tZ = ta - kp778434453 * tY ; t14 = tA - kp830830026 * t13 + t1Q = t1u - kp918985947 * t1P ; t1B = t1i - kp778434453 * t1A + t10 = td - kp876768831 * tZ ; t15 = tr + kp918985947 * t14 + t11 = xr0 - kp959492973 * t10 ; t1C = t1j - kp876768831 * t1B + t1G = t1q - kp830830026 * t1F ; t1n = t1h - kp778434453 * t1m + t1D = xi0 - kp959492973 * t1C ; t1H = t1t + kp918985947 * t1G + t1o = t1g - kp876768831 * t1n ; t1x = t1r - kp830830026 * t1w + t18 = t7 - kp778434453 * t17 ; t1p = xi0 - kp959492973 * t1o + t1y = t1q - kp918985947 * t1x ; t19 = t4 - kp876768831 * t18 + t1d = tu - kp830830026 * t1c ; t1a = xr0 - kp959492973 * t19 + t1e = tA - kp918985947 * t1d + r10 = (tl + kp989821441 * tE) :+ (t24 + kp989821441 * t28) + r9 = (tJ - kp989821441 * tN) :+ (t1V - kp989821441 * t1Z) + r8 = (tS + kp989821441 * tW) :+ (t1M + kp989821441 * t1Q) + r7 = (t11 - kp989821441 * t15) :+ (t1D - kp989821441 * t1H) + r6 = (t1a + kp989821441 * t1e) :+ (t1p + kp989821441 * t1y) + r5 = (t1a - kp989821441 * t1e) :+ (t1p - kp989821441 * t1y) + r4 = (t11 + kp989821441 * t15) :+ (t1D + kp989821441 * t1H) + r3 = (tS - kp989821441 * tW) :+ (t1M - kp989821441 * t1Q) + r2 = (tJ + kp989821441 * tN) :+ (t1V + kp989821441 * t1Z) + r1 = (tl - kp989821441 * tE) :+ (t24 - kp989821441 * t28) + MV.unsafeWrite xsout 0 $ (xr0+t4+t7+ta+td+tg) :+ (xi0+t1g+t1h+t1i+t1j+t1k) + MV.unsafeWrite xsout 1 $ if sign == 1 then r1 else r10 + MV.unsafeWrite xsout 2 $ if sign == 1 then r2 else r9 + MV.unsafeWrite xsout 3 $ if sign == 1 then r3 else r8 + MV.unsafeWrite xsout 4 $ if sign == 1 then r4 else r7 + MV.unsafeWrite xsout 5 $ if sign == 1 then r5 else r6 + MV.unsafeWrite xsout 6 $ if sign == 1 then r6 else r5 + MV.unsafeWrite xsout 7 $ if sign == 1 then r7 else r4 + MV.unsafeWrite xsout 8 $ if sign == 1 then r8 else r3 + MV.unsafeWrite xsout 9 $ if sign == 1 then r9 else r2 + MV.unsafeWrite xsout 10 $ if sign == 1 then r10 else r1 + +-- | Length 13 hard-coded FFT. +kp875502302, kp520028571, kp575140729 :: Double +kp600477271, kp300462606, kp516520780 :: Double +kp968287244, kp503537032, kp251768516 :: Double +kp581704778, kp859542535, kp083333333 :: Double +kp957805992, kp522026385, kp853480001 :: Double +kp769338817, kp612264650, kp038632954 :: Double +kp302775637, kp514918778, kp686558370 :: Double +kp226109445, kp301479260 :: Double +--kp866025403, kp500000000 :: Double +kp875502302 = 0.875502302409147941146295545768755143177842006 +kp520028571 = 0.520028571888864619117130500499232802493238139 +kp575140729 = 0.575140729474003121368385547455453388461001608 +kp600477271 = 0.600477271932665282925769253334763009352012849 +kp300462606 = 0.300462606288665774426601772289207995520941381 +kp516520780 = 0.516520780623489722840901288569017135705033622 +kp968287244 = 0.968287244361984016049539446938120421179794516 +kp503537032 = 0.503537032863766627246873853868466977093348562 +kp251768516 = 0.251768516431883313623436926934233488546674281 +kp581704778 = 0.581704778510515730456870384989698884939833902 +kp859542535 = 0.859542535098774820163672132761689612766401925 +kp083333333 = 0.083333333333333333333333333333333333333333333 +kp957805992 = 0.957805992594665126462521754605754580515587217 +kp522026385 = 0.522026385161275033714027226654165028300441940 +kp853480001 = 0.853480001859823990758994934970528322872359049 +kp769338817 = 0.769338817572980603471413688209101117038278899 +kp612264650 = 0.612264650376756543746494474777125408779395514 +kp038632954 = 0.038632954644348171955506895830342264440241080 +kp302775637 = 0.302775637731994646559610633735247973125648287 +kp514918778 = 0.514918778086315755491789696138117261566051239 +kp686558370 = 0.686558370781754340655719594850823015421401653 +kp226109445 = 0.226109445035782405468510155372505010481906348 +kp301479260 = 0.301479260047709873958013540496673347309208464 +--kp866025403 = 0.866025403784438646763723170752936183471402627 +--kp500000000 = 0.500000000000000000000000000000000000000000000 +special13 :: Int -> MVCD s -> MVCD s -> ST s () +special13 sign xsin xsout = do + xr0 :+ xi0 <- MV.unsafeRead xsin 0 ; xr1 :+ xi1 <- MV.unsafeRead xsin 1 + xr2 :+ xi2 <- MV.unsafeRead xsin 2 ; xr3 :+ xi3 <- MV.unsafeRead xsin 3 + xr4 :+ xi4 <- MV.unsafeRead xsin 4 ; xr5 :+ xi5 <- MV.unsafeRead xsin 5 + xr6 :+ xi6 <- MV.unsafeRead xsin 6 ; xr7 :+ xi7 <- MV.unsafeRead xsin 7 + xr8 :+ xi8 <- MV.unsafeRead xsin 8 ; xr9 :+ xi9 <- MV.unsafeRead xsin 9 + xr10 :+ xi10 <- MV.unsafeRead xsin 10 ; xr11 :+ xi11 <- MV.unsafeRead xsin 11 + xr12 :+ xi12 <- MV.unsafeRead xsin 12 + let t2d = xr8 - xr5 ; tf = xr8 + xr5 ; ta = xr10 + xr4 ; tq = xr10 - xr4 + ty = kp500000000 * ta - xr12 + tb = xr12 + ta ; tr = xr9 - xr3 ; t5 = xr3 + xr9 ; t6 = xr1 + t5 + tx = xr1 - kp500000000 * t5 + ti = xr11 + xr6 ; tt = xr11 - xr6 ; tu = xr7 - xr2 ; tl = xr7 + xr2 + tc = t6 + tb ; t2n = t6 - tb ; t2b = ti - tl ; tm = ti + tl + t2e = tt + tu ; tv = tt - tu ; ts = tq - tr ; t2g = tr + tq + tz = tx - ty ; t2a = tx + ty + tA = tf - kp500000000 * tm + tn = tf + tm + t2f = t2d - kp500000000 * t2e + t2o = t2d + t2e ; to = tc + tn ; tH = tc - tn + t2h = t2f + kp866025403 * t2g ; t2k = t2f - kp866025403 * t2g + tE = tz - tA ; tB = tz + tA ; tF = ts - tv ; tw = ts + tv + t2j = t2a - kp866025403 * t2b ; t2c = t2a + kp866025403 * t2b + t1R = xi8 + xi5 ; tM = xi8 - xi5 ; t17 = xi10 + xi4 ; t10 = xi10 - xi4 + t18 = kp500000000 * t17 - xi12 + t1l = xi12 + t17 ; tX = xi9 - xi3 ; t14 = xi3 + xi9 ; t1k = xi1 + t14 + t15 = xi1 - kp500000000 * t14 + tP = xi11 - xi6 ; t1a = xi11 + xi6 ; t1b = xi7 + xi2 ; tS = xi7 - xi2 + t1Q = t1k + t1l ; t1m = t1k - t1l ; t11 = tX + t10 ; t1W = t10 - tX + t1X = tP - tS ; tT = tP + tS ; t1S = t1a + t1b ; t1c = t1a - t1b + t19 = t15 + t18 ; t1Z = t15 - t18 ; t1j = tM + tT + tU = tM - kp500000000 * tT + t1T = t1R + t1S + t20 = t1R - kp500000000 * t1S ; t12 = tU + kp866025403 * t11 + t1f = tU - kp866025403 * t11 + t21 = t1Z + t20 ; t24 = t1Z - t20 ; t27 = t1Q - t1T ; t1U = t1Q + t1T + t1g = t19 - kp866025403 * t1c ; t1d = t19 + kp866025403 * t1c + t25 = t1W - t1X ; t1Y = t1W + t1X + tC = tw + kp301479260 * tB ; t1x = tB - kp226109445 * tw + t1y = tF + kp686558370 * tE ; tG = tE - kp514918778 * tF + t1n = t1j - kp302775637 * t1m ; t1G = t1m + kp302775637 * t1j + t1u = t1d - kp038632954 * t12 ; t1e = t12 + kp038632954 * t1d + t1h = t1f + kp612264650 * t1g ; t1v = t1g - kp612264650 * t1f + t1J = t1x + kp769338817 * t1y ; t1z = t1x - kp769338817 * t1y + t1H = t1u - kp853480001 * t1v ; t1w = t1u + kp853480001 * t1v + t1I = t1G - kp522026385 * t1H ; t1O = t1H + kp957805992 * t1G + tp = xr0 - kp083333333 * to ; t1E = t1e + kp853480001 * t1h + t1i = t1e - kp853480001 * t1h ; t1q = tH - kp859542535 * tG + tI = tG + kp581704778 * tH ; t1o = t1i + kp957805992 * t1n + t1s = t1n - kp522026385 * t1i ; t1p = tp - kp251768516 * tC + tD = tp + kp503537032 * tC ; t1C = t1w - kp968287244 * t1z + t1A = t1w + kp968287244 * t1z ; tJ = tD + kp516520780 * tI + t1N = tD - kp516520780 * tI ; t1D = t1p - kp300462606 * t1q + t1r = t1p + kp300462606 * t1q ; t1t = t1r - kp575140729 * t1s + t1B = t1r + kp575140729 * t1s ; t1L = t1D - kp520028571 * t1E + t1F = t1D + kp520028571 * t1E ; t1K = t1I + kp875502302 * t1J + t1M = t1I - kp875502302 * t1J ; t2D = t21 - kp226109445 * t1Y + t22 = t1Y + kp301479260 * t21 ; t26 = t24 - kp514918778 * t25 + t2E = t25 + kp686558370 * t24 ; t2v = t2o - kp302775637 * t2n + t2p = t2n + kp302775637 * t2o ; t2i = t2c - kp038632954 * t2h + t2s = t2h + kp038632954 * t2c ; t2t = t2k + kp612264650 * t2j + t2l = t2j - kp612264650 * t2k ; t2F = t2D - kp769338817 * t2E + t2N = t2D + kp769338817 * t2E ; t2K = t2s + kp853480001 * t2t + t2u = t2s - kp853480001 * t2t ; t2w = t2u + kp957805992 * t2v + t2A = t2v - kp522026385 * t2u ; t1V = xi0 - kp083333333 * t1U + t2m = t2i - kp853480001 * t2l ; t2C = t2i + kp853480001 * t2l + t28 = t26 + kp581704778 * t27 ; t2y = t27 - kp859542535 * t26 + t2M = t2p - kp522026385 * t2m ; t2q = t2m + kp957805992 * t2p + t23 = t1V + kp503537032 * t22 ; t2x = t1V - kp251768516 * t22 + t2O = t2M - kp875502302 * t2N ; t2Q = t2M + kp875502302 * t2N + t2r = t23 + kp516520780 * t28 ; t29 = t23 - kp516520780 * t28 + t2z = t2x + kp300462606 * t2y ; t2J = t2x - kp300462606 * t2y + t2P = t2J + kp520028571 * t2K ; t2L = t2J - kp520028571 * t2K + t2B = t2z + kp575140729 * t2A ; t2H = t2z - kp575140729 * t2A + t2I = t2C + kp968287244 * t2F ; t2G = t2C - kp968287244 * t2F + r12 = (tJ - kp600477271 * t1o) :+ (t2r + kp600477271 * t2w) + r11 = (t1F + kp575140729 * t1K) :+ (t2L - kp575140729 * t2O) + r10 = (t1t + kp520028571 * t1A) :+ (t2B - kp520028571 * t2G) + r9 = (t1B + kp520028571 * t1C) :+ (t2H - kp520028571 * t2I) + r8 = (t1N + kp600477271 * t1O) :+ (t29 - kp600477271 * t2q) + r7 = (t1L + kp575140729 * t1M) :+ (t2P - kp575140729 * t2Q) + r6 = (t1F - kp575140729 * t1K) :+ (t2L + kp575140729 * t2O) + r5 = (t1N - kp600477271 * t1O) :+ (t29 + kp600477271 * t2q) + r4 = (t1t - kp520028571 * t1A) :+ (t2B + kp520028571 * t2G) + r3 = (t1B - kp520028571 * t1C) :+ (t2H + kp520028571 * t2I) + r2 = (t1L - kp575140729 * t1M) :+ (t2P + kp575140729 * t2Q) + r1 = (tJ + kp600477271 * t1o) :+ (t2r - kp600477271 * t2w) + MV.unsafeWrite xsout 0 $ (xr0 + to) :+ (xi0 + t1U) + MV.unsafeWrite xsout 1 $ if sign == 1 then r1 else r12 + MV.unsafeWrite xsout 2 $ if sign == 1 then r2 else r11 + MV.unsafeWrite xsout 3 $ if sign == 1 then r3 else r10 + MV.unsafeWrite xsout 4 $ if sign == 1 then r4 else r9 + MV.unsafeWrite xsout 5 $ if sign == 1 then r5 else r8 + MV.unsafeWrite xsout 6 $ if sign == 1 then r6 else r7 + MV.unsafeWrite xsout 7 $ if sign == 1 then r7 else r6 + MV.unsafeWrite xsout 8 $ if sign == 1 then r8 else r5 + MV.unsafeWrite xsout 9 $ if sign == 1 then r9 else r4 + MV.unsafeWrite xsout 10 $ if sign == 1 then r10 else r3 + MV.unsafeWrite xsout 11 $ if sign == 1 then r11 else r2 + MV.unsafeWrite xsout 12 $ if sign == 1 then r12 else r1 diff --git a/Numeric/FFT/Types.hs b/Numeric/FFT/Types.hs new file mode 100644 --- /dev/null +++ b/Numeric/FFT/Types.hs @@ -0,0 +1,70 @@ +module Numeric.FFT.Types + ( VCD, MVCD, VVCD, VMVCD, VVVCD, VI + , Direction (..), Plan (..), BaseTransform (..) + ) where + +import Data.IntMap.Strict (IntMap) +import Data.Vector.Unboxed +import qualified Data.Vector as V +import qualified Data.Vector.Unboxed.Mutable as MV +import Data.Complex + + +-- | Some useful type synonyms. +type VCD = Vector (Complex Double) +type MVCD s = MV.MVector s (Complex Double) +type VVCD = V.Vector VCD +type VVVCD = V.Vector VVCD +type VMVCD a = V.Vector (MVCD a) +type VI = Vector Int + +-- | Transform direction: 'Forward' is the normal FFT, 'Inverse' is +-- inverse FFT. +data Direction = Forward | Inverse deriving (Eq, Show) + +-- | A FFT plan. This depends only on the problem size and can be +-- pre-computed and reused to transform (and inverse transform) any +-- number of vectors of the given size. +data Plan = Plan { plDLInfo :: V.Vector (Int, Int, VVVCD, VVVCD) + -- ^ Size information and diagonal matrix entries + -- for Danielson-Lanczos recursive decomposition of + -- problem size. + , plPermute :: Maybe VI + -- ^ Input vector permutation to use before base + -- transformation and recursive Danielson-Lanczos + -- composition. + , plBase :: BaseTransform + -- ^ Base transformation used for each sub-vector + -- before performing recursive Danielson-Lanczos + -- steps to form the full FFT result. + } deriving (Eq, Show) + +-- | A "base transform" used at the "bottom" of the recursive +-- Cooley-Tukey decomposition of the input problem size: either a +-- simple DFT, a special hard-coded small problem size case, or a +-- Rader prime-length FFT invocation. +data BaseTransform = SpecialBase { baseSize :: Int } + -- ^ Hard-coded small-size base transform. + | DFTBase { baseSize :: Int + , dftWsFwd :: VCD + , dftWsInv :: VCD } + -- ^ Simple DFT base transform, giving problem + -- size and powers of roots of unity needed for + -- transform. + | RaderBase { baseSize :: Int + , raderOutPerm :: VI + , raderBFwd :: VCD + , raderBInv :: VCD + , raderConvSize :: Int + , raderConvPlan :: Plan } + -- ^ Prime-length Rader FFT base transform, + -- giving problem size, output index permutation + -- (the input index permutation is folded into + -- the main input permutation of the full + -- transform), pre-transformed Rader b sequence + -- for forward and inverse problems, the (padded + -- or not) problem size for Rader sequence + -- convolution and a sub-plan (either of size + -- baseSize-1 or the next larger power of two) + -- for computing the Rader convolution. + deriving (Eq, Show) diff --git a/Numeric/FFT/Utils.hs b/Numeric/FFT/Utils.hs new file mode 100644 --- /dev/null +++ b/Numeric/FFT/Utils.hs @@ -0,0 +1,199 @@ +module Numeric.FFT.Utils + ( omega, slicevecs, slicemvecs, primes, isPrime + , allFactors, factors + , primitiveRoot, invModN, log2, isPow2, dupperm, (%.%) + , compositions, makeComp, multisetPerms + , backpermuteM + ) where + +import Prelude hiding (all, concatMap, dropWhile, enumFromTo, + filter, head, length, map, maximum, null, reverse) +import qualified Prelude as P +import qualified Control.Monad as CM +import Control.Monad.ST +import Data.Bits +import Data.Complex +import Data.Vector.Unboxed +import qualified Data.Vector as V +import qualified Data.Vector.Unboxed.Mutable as MV +import Data.List (nub) +import qualified Data.List as L + +import Numeric.FFT.Types + + +-- | Roots of unity. +omega :: Int -> Complex Double +omega n = cis (2 * pi / fromIntegral n) + +-- | Slice a vector @v@ into equally sized parts, each of length @m@. +slicevecs :: Int -> VCD -> VVCD +slicevecs m v = V.map (\i -> slice (i * m) m v) $ + V.enumFromN 0 (length v `div` m) + +-- | Slice a mutable vector @v@ into equally sized parts, each of +-- length @m@. +slicemvecs :: Int -> MVCD a -> VMVCD a +slicemvecs m v = V.map (\i -> MV.slice (i * m) m v) $ + V.enumFromN 0 (MV.length v `div` m) + +-- | Determine primitive roots modulo n. +-- +-- From Wikipedia (https://en.wikipedia.org/wiki/Primitive_root_modulo_n): +-- +-- No simple general formula to compute primitive roots modulo n is +-- known. There are however methods to locate a primitive root that +-- are faster than simply trying out all candidates. If the +-- multiplicative order of a number m modulo n is equal to phi(n) (the +-- order of Z_n^x), then it is a primitive root. In fact the converse +-- is true: if m is a primitive root modulo n, then the multiplicative +-- order of m is phi(n). We can use this to test for primitive roots. +-- [Here, phi(n) is Euler's totient function, and Z_n^x is the +-- multiplicative group of integers modulo n.] +-- +-- First, compute phi(n). Then determine the different prime factors +-- of phi(n), say p1, ..., pk. Now, for every element m of Z_n^x, +-- compute +-- +-- m^(phi(n) / pi) mod n for i = 1, ..., k +-- +-- using a fast algorithm for modular exponentiation such as +-- exponentiation by squaring. A number m for which these k results +-- are all different from 1 is a primitive root. +-- +-- [In our case, n is restricted to being prime, and phi(p) = p - 1 +-- for prime p.] +-- +primitiveRoot :: Int -> Int +primitiveRoot p + | isPrime p = + let tot = p - 1 + -- ^ Euler's totient function for prime values. + totpows = map (tot `div`) $ fromList $ nub $ toList $ allFactors tot + -- ^ Powers to check. + check n = all (/=1) $ map (expt p n) totpows + -- ^ All powers are different from 1 => primitive root. + in fromIntegral $ head $ dropWhile (not . check) $ fromList [1..p-1] + | otherwise = error "Attempt to take primitive root of non-prime value" + +-- | Fast exponentation modulo n by squaring. +expt :: Int -> Int -> Int -> Int +expt n b pow = fromIntegral $ go pow + where bb = fromIntegral b + nb = fromIntegral n + go :: Int -> Integer + go p + | p == 0 = 1 + | p `mod` 2 == 1 = (bb * go (p - 1)) `mod` nb + | otherwise = let h = go (p `div` 2) in (h * h) `mod` nb + +-- | Find inverse element in multiplicative integer group modulo n. +invModN :: Int -> Int -> Int +invModN n g = head $ filter (\iv -> (g * iv) `mod` n == 1) $ enumFromTo 1 (n-1) + +-- | Prime sieve from Haskell wiki. +primes :: Integral a => [a] +primes = 2 : primes' + where primes' = sieve [3, 5 ..] 9 primes' + sieve (x:xs) q ps@ ~(p:t) + | x < q = x : sieve xs q ps + | True = sieve [n | n <- xs, rem n p /= 0] (P.head t^2) t + +-- | Naive primality testing. +isPrime :: Integral a => a -> Bool +isPrime n = n `P.elem` P.takeWhile (<= n) primes + +-- | Simple prime factorisation. +allFactors :: (Integral a, Unbox a) => a -> Vector a +allFactors n = fromList $ go n primes + where go cur pss@(p:ps) + | cur == p = [p] + | cur `mod` p == 0 = p : go (cur `div` p) pss + | otherwise = go cur ps + +-- | Simple prime factorisation: small factors only; largest/last +-- factor picked out as "special". +factors :: (Integral a, Unbox a) => a -> (a, Vector a) +factors n = let (lst, rest) = go n primes in (lst, fromList rest) + where go cur pss@(p:ps) + | cur == p = (p, []) + | cur `mod` p == 0 = let (lst, rest) = go (cur `div` p) pss + in (lst, p : rest) + | otherwise = go cur ps + +-- | Base-2 logarithm. +log2 :: Int -> Int +log2 1 = 0 +log2 n = 1 + log2 (n `div` 2) + +-- | Check for powers of two. +isPow2 :: Int -> Bool +isPow2 1 = True +isPow2 n + | n `mod` 2 == 0 = isPow2 $ n `div` 2 + | otherwise = False + +-- | Duplicate a sub-permutation to fill a given vector length. +dupperm :: Int -> VI -> VI +dupperm n p = + let sublen = length p + shift di = map (+(sublen * di)) p + in concatMap shift $ enumFromN 0 (n `div` sublen) + +-- | Composition of permutations. +(%.%) :: VI -> VI -> VI +p1 %.% p2 = backpermute p2 p1 + +-- | Generate all compositions of a given integer. +compositions :: Int -> V.Vector (Vector Int) +compositions 0 = V.empty +compositions n = let fs = allFactors n + in V.reverse $ V.map (makeComp fs) $ V.enumFromN 0 (2^(n-1)) + +-- | Generate a single composition of a given integer. +makeComp :: Vector Int -> Int -> Vector Int +makeComp fs i = fromList $ foldOps (toList fs) $ makeOps (length fs) i + where foldOps :: [Int] -> [Bool] -> [Int] + foldOps (f:fs) ops = go f fs ops + where go acc [] [] = [acc] + go acc (f:fs) (op:ops) = if op then go (acc * f) fs ops + else acc : go f fs ops + makeOps :: Int -> Int -> [Bool] + makeOps n i = P.replicate (n - 1 - P.length bs) False P.++ bs + where bs = P.dropWhile not $ P.reverse $ + P.map (testBit i) [0..bitSize i-1] + +-- | Generate all distinct permutations of a multiset in lexicographic +-- order. +multisetPerms :: Vector Int -> [Vector Int] +multisetPerms idp = sidp : L.unfoldr step sidp + where sidp = fromList $ L.sort $ toList idp + step v = case permStep v of + Nothing -> Nothing + Just p -> Just (p, p) + permStep :: Vector Int -> Maybe (Vector Int) + permStep v = + if null ks + then Nothing + else let k = maximum ks + ls = filter (\i -> v ! k < v ! i) $ enumFromN 0 n + l = maximum ls + in Just $ revEnd k (swap k l) + where n = length v + ks = filter (\i -> v ! i < v ! (i+1)) $ enumFromN 0 (n-1) + swap a b = generate n $ \i -> + if i == a then v ! b + else if i == b then v ! a + else v ! i + revEnd f vv = generate n $ \i -> + if i <= f + then vv ! i + else vv ! (n - i + f) + +-- | ST monad version of vector permutation. +backpermuteM :: Int -> VI -> MVCD s -> MVCD s -> ST s () +backpermuteM n perm vin vout = do + CM.forM_ [0..n-1] $ \i -> do + idx <- indexM perm i + x <- MV.unsafeRead vin idx + MV.unsafeWrite vout i x diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,4 @@ +arb-fft +======= + +Pure Haskell arbitrary length FFT library 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/arb-fft.cabal b/arb-fft.cabal new file mode 100644 --- /dev/null +++ b/arb-fft.cabal @@ -0,0 +1,76 @@ +name: arb-fft +version: 0.1.0.0 +synopsis: Pure Haskell arbitrary length FFT library +homepage: https://github.com/ian-ross/arb-fft +license: GPL-3 +license-file: LICENSE +author: Ian Ross +maintainer: ian@skybluetrades.net +copyright: Copyright (2013) Ian Ross +category: Math +build-type: Simple +extra-source-files: README.md +cabal-version: >=1.10 + +description: + This library provides a pure Haskell complex-to-complex Fast Fourier + Transform implementation for arbitrary length input vectors, using a + mixed-radix decimation-in-time algorithm with specialised + straight-line code for a range of base transform sizes, Rader's + algorithm for prime length base transforms, and an empirical + optimisation scheme to select a good problem decomposition. + . + + This package is probably primarily of pedagogical interest (FFTW is + about five times faster for most input sizes). There is a long + series of blog articles describing the development of the package, + indexed at . + +source-repository head + type: git + location: https://github.com/ian-ross/arb-fft + +Library + exposed-modules: Numeric.FFT + other-modules: Numeric.FFT.Plan + Numeric.FFT.Execute + Numeric.FFT.Special + Numeric.FFT.Special.PowersOfTwo + Numeric.FFT.Special.Primes + Numeric.FFT.Special.Miscellaneous + Numeric.FFT.Types + Numeric.FFT.Utils + ghc-options: -O2 -fllvm + ghc-prof-options: -auto-all -caf-all + build-depends: base >= 4.6 && < 5, + containers >= 0.5.0.0 && < 0.6, + criterion >= 0.8.0.0 && < 0.9, + directory >= 1.2.0.1 && < 1.3, + filepath >= 1.3.0.1 && < 1.4, + primitive >= 0.5.1.0 && < 0.6, + transformers >= 0.3.0.0 && < 0.4, + vector >= 0.10.9.1 && < 0.11 + default-language: Haskell2010 + +Test-Suite basic-test + type: exitcode-stdio-1.0 + hs-source-dirs: test + main-is: basic-test.hs + build-depends: arb-fft, + base >= 4.6 && < 5, + containers >= 0.5.0.0 && < 0.6, + vector >= 0.10.9.1 && < 0.11, + QuickCheck >= 2.6 && < 2.7, + tasty >= 0.3, + tasty-quickcheck >= 0.3 + default-language: Haskell2010 + +Executable profile-256 + hs-source-dirs: test + main-is: profile-256.hs + build-depends: arb-fft, + base >= 4.6 && < 5, + containers >= 0.5.0.0 && < 0.6, + vector >= 0.10.9.1 && < 0.11, + criterion + default-language: Haskell2010 diff --git a/test/basic-test.hs b/test/basic-test.hs new file mode 100644 --- /dev/null +++ b/test/basic-test.hs @@ -0,0 +1,70 @@ +{-# LANGUAGE ScopedTypeVariables, TypeSynonymInstances, FlexibleInstances #-} +module Main where + +import Prelude hiding ((++), length, map, maximum, sum, zipWith) +import qualified Prelude as P +import Test.Tasty +import Test.Tasty.QuickCheck +import Test.QuickCheck +import Test.QuickCheck.Monadic +import Control.Applicative ((<$>)) +import Data.Complex +import Data.Vector +import System.IO.Unsafe (unsafePerformIO) + +import Debug.Trace + +import Numeric.FFT + +main :: IO () +main = defaultMain tests + +tests :: TestTree +tests = testGroup "Tests" + [ testProperty "FFT vs. DFT" prop_dft_vs_fft + , testProperty "FFT/IFFT round-trip" prop_ifft + ] + +-- Clean up number display. +defuzz :: Vector (Complex Double) -> Vector (Complex Double) +defuzz = map (\(r :+ i) -> df r :+ df i) + where df x = if abs x < 1.0E-6 then 0 else x + +-- Check FFT against DFT. +check :: Vector (Complex Double) -> IO (Double, Vector (Complex Double)) +check v = do + tst <- fft v + let diff = defuzz $ zipWith (-) tst (basicDFT v) + return (maximum $ map magnitude diff, diff) + +-- QuickCheck property for FFT vs. DFT testing. +prop_dft_vs_fft :: Property +prop_dft_vs_fft = monadicIO $ do + v <- pick arbitrary + chk <- run $ check v + assert $ fst chk < 1.0E-6 + +-- QuickCheck property for inverse FFT round-trip testing. +prop_ifft :: Property +prop_ifft = monadicIO $ do + v <- pick arbitrary + fwd <- run $ fft v + bwd <- run $ ifft fwd + let diff = zipWith (-) v bwd + assert $ maximum (map magnitude diff) < 1.0E-6 + +-- Non-zero length arbitrary vectors. +instance Arbitrary (Vector (Complex Double)) where + arbitrary = fromList <$> listOf1 arbitrary + + +-- Roots of unity. +omega :: Int -> Complex Double +omega n = cis (2 * pi / fromIntegral n) + +-- Naïve DFT. +basicDFT :: Vector (Complex Double) -> Vector (Complex Double) +basicDFT h = generate n doone + where n = length h + w = omega n + doone i = sum $ zipWith (*) h $ generate n (\k -> w^^(i*k)) diff --git a/test/profile-256.hs b/test/profile-256.hs new file mode 100644 --- /dev/null +++ b/test/profile-256.hs @@ -0,0 +1,15 @@ +module Main where + +import Criterion.Main +import Data.Complex +import Data.Vector +import qualified Numeric.FFT as FFT + +tstvec :: Int -> Vector (Complex Double) +tstvec sz = generate sz (\i -> let ii = fromIntegral i + in sin (2*pi*ii/1024) + sin (2*pi*ii/511)) + +main :: IO () +main = do + p <- FFT.plan 256 + run (nf (FFT.fftWith p) (tstvec 256)) 1000