arb-fft (empty) → 0.1.0.0
raw patch · 15 files changed
+3671/−0 lines, 15 filesdep +QuickCheckdep +arb-fftdep +basesetup-changed
Dependencies added: QuickCheck, arb-fft, base, containers, criterion, directory, filepath, primitive, tasty, tasty-quickcheck, transformers, vector
Files
- LICENSE +675/−0
- Numeric/FFT.hs +40/−0
- Numeric/FFT/Execute.hs +228/−0
- Numeric/FFT/Plan.hs +312/−0
- Numeric/FFT/Special.hs +54/−0
- Numeric/FFT/Special/Miscellaneous.hs +779/−0
- Numeric/FFT/Special/PowersOfTwo.hs +808/−0
- Numeric/FFT/Special/Primes.hs +339/−0
- Numeric/FFT/Types.hs +70/−0
- Numeric/FFT/Utils.hs +199/−0
- README.md +4/−0
- Setup.hs +2/−0
- arb-fft.cabal +76/−0
- test/basic-test.hs +70/−0
- test/profile-256.hs +15/−0
+ LICENSE view
@@ -0,0 +1,675 @@+ GNU GENERAL PUBLIC LICENSE+ Version 3, 29 June 2007++ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>+ Everyone is permitted to copy and distribute verbatim copies+ of this license document, but changing it is not allowed.++ Preamble++ The GNU General Public License is a free, copyleft license for+software and other kinds of works.++ The licenses for most software and other practical works are designed+to take away your freedom to share and change the works. By contrast,+the GNU General Public License is intended to guarantee your freedom to+share and change all versions of a program--to make sure it remains free+software for all its users. We, the Free Software Foundation, use the+GNU General Public License for most of our software; it applies also to+any other work released this way by its authors. You can apply it to+your programs, too.++ When we speak of free software, we are referring to freedom, not+price. Our General Public Licenses are designed to make sure that you+have the freedom to distribute copies of free software (and charge for+them if you wish), that you receive source code or can get it if you+want it, that you can change the software or use pieces of it in new+free programs, and that you know you can do these things.++ To protect your rights, we need to prevent others from denying you+these rights or asking you to surrender the rights. Therefore, you have+certain responsibilities if you distribute copies of the software, or if+you modify it: responsibilities to respect the freedom of others.++ For example, if you distribute copies of such a program, whether+gratis or for a fee, you must pass on to the recipients the same+freedoms that you received. You must make sure that they, too, receive+or can get the source code. And you must show them these terms so they+know their rights.++ Developers that use the GNU GPL protect your rights with two steps:+(1) assert copyright on the software, and (2) offer you this License+giving you legal permission to copy, distribute and/or modify it.++ For the developers' and authors' protection, the GPL clearly explains+that there is no warranty for this free software. For both users' and+authors' sake, the GPL requires that modified versions be marked as+changed, so that their problems will not be attributed erroneously to+authors of previous versions.++ Some devices are designed to deny users access to install or run+modified versions of the software inside them, although the manufacturer+can do so. This is fundamentally incompatible with the aim of+protecting users' freedom to change the software. The systematic+pattern of such abuse occurs in the area of products for individuals to+use, which is precisely where it is most unacceptable. Therefore, we+have designed this version of the GPL to prohibit the practice for those+products. If such problems arise substantially in other domains, we+stand ready to extend this provision to those domains in future versions+of the GPL, as needed to protect the freedom of users.++ Finally, every program is threatened constantly by software patents.+States should not allow patents to restrict development and use of+software on general-purpose computers, but in those that do, we wish to+avoid the special danger that patents applied to a free program could+make it effectively proprietary. To prevent this, the GPL assures that+patents cannot be used to render the program non-free.++ The precise terms and conditions for copying, distribution and+modification follow.++ TERMS AND CONDITIONS++ 0. Definitions.++ "This License" refers to version 3 of the GNU General Public License.++ "Copyright" also means copyright-like laws that apply to other kinds of+works, such as semiconductor masks.+ + "The Program" refers to any copyrightable work licensed under this+License. Each licensee is addressed as "you". "Licensees" and+"recipients" may be individuals or organizations.++ To "modify" a work means to copy from or adapt all or part of the work+in a fashion requiring copyright permission, other than the making of an+exact copy. The resulting work is called a "modified version" of the+earlier work or a work "based on" the earlier work.++ A "covered work" means either the unmodified Program or a work based+on the Program.++ To "propagate" a work means to do anything with it that, without+permission, would make you directly or secondarily liable for+infringement under applicable copyright law, except executing it on a+computer or modifying a private copy. Propagation includes copying,+distribution (with or without modification), making available to the+public, and in some countries other activities as well.++ To "convey" a work means any kind of propagation that enables other+parties to make or receive copies. Mere interaction with a user through+a computer network, with no transfer of a copy, is not conveying.++ An interactive user interface displays "Appropriate Legal Notices"+to the extent that it includes a convenient and prominently visible+feature that (1) displays an appropriate copyright notice, and (2)+tells the user that there is no warranty for the work (except to the+extent that warranties are provided), that licensees may convey the+work under this License, and how to view a copy of this License. If+the interface presents a list of user commands or options, such as a+menu, a prominent item in the list meets this criterion.++ 1. Source Code.++ The "source code" for a work means the preferred form of the work+for making modifications to it. "Object code" means any non-source+form of a work.++ A "Standard Interface" means an interface that either is an official+standard defined by a recognized standards body, or, in the case of+interfaces specified for a particular programming language, one that+is widely used among developers working in that language.++ The "System Libraries" of an executable work include anything, other+than the work as a whole, that (a) is included in the normal form of+packaging a Major Component, but which is not part of that Major+Component, and (b) serves only to enable use of the work with that+Major Component, or to implement a Standard Interface for which an+implementation is available to the public in source code form. A+"Major Component", in this context, means a major essential component+(kernel, window system, and so on) of the specific operating system+(if any) on which the executable work runs, or a compiler used to+produce the work, or an object code interpreter used to run it.++ The "Corresponding Source" for a work in object code form means all+the source code needed to generate, install, and (for an executable+work) run the object code and to modify the work, including scripts to+control those activities. However, it does not include the work's+System Libraries, or general-purpose tools or generally available free+programs which are used unmodified in performing those activities but+which are not part of the work. For example, Corresponding Source+includes interface definition files associated with source files for+the work, and the source code for shared libraries and dynamically+linked subprograms that the work is specifically designed to require,+such as by intimate data communication or control flow between those+subprograms and other parts of the work.++ The Corresponding Source need not include anything that users+can regenerate automatically from other parts of the Corresponding+Source.++ The Corresponding Source for a work in source code form is that+same work.++ 2. Basic Permissions.++ All rights granted under this License are granted for the term of+copyright on the Program, and are irrevocable provided the stated+conditions are met. This License explicitly affirms your unlimited+permission to run the unmodified Program. The output from running a+covered work is covered by this License only if the output, given its+content, constitutes a covered work. This License acknowledges your+rights of fair use or other equivalent, as provided by copyright law.++ You may make, run and propagate covered works that you do not+convey, without conditions so long as your license otherwise remains+in force. You may convey covered works to others for the sole purpose+of having them make modifications exclusively for you, or provide you+with facilities for running those works, provided that you comply with+the terms of this License in conveying all material for which you do+not control copyright. Those thus making or running the covered works+for you must do so exclusively on your behalf, under your direction+and control, on terms that prohibit them from making any copies of+your copyrighted material outside their relationship with you.++ Conveying under any other circumstances is permitted solely under+the conditions stated below. Sublicensing is not allowed; section 10+makes it unnecessary.++ 3. Protecting Users' Legal Rights From Anti-Circumvention Law.++ No covered work shall be deemed part of an effective technological+measure under any applicable law fulfilling obligations under article+11 of the WIPO copyright treaty adopted on 20 December 1996, or+similar laws prohibiting or restricting circumvention of such+measures.++ When you convey a covered work, you waive any legal power to forbid+circumvention of technological measures to the extent such circumvention+is effected by exercising rights under this License with respect to+the covered work, and you disclaim any intention to limit operation or+modification of the work as a means of enforcing, against the work's+users, your or third parties' legal rights to forbid circumvention of+technological measures.++ 4. Conveying Verbatim Copies.++ You may convey verbatim copies of the Program's source code as you+receive it, in any medium, provided that you conspicuously and+appropriately publish on each copy an appropriate copyright notice;+keep intact all notices stating that this License and any+non-permissive terms added in accord with section 7 apply to the code;+keep intact all notices of the absence of any warranty; and give all+recipients a copy of this License along with the Program.++ You may charge any price or no price for each copy that you convey,+and you may offer support or warranty protection for a fee.++ 5. Conveying Modified Source Versions.++ You may convey a work based on the Program, or the modifications to+produce it from the Program, in the form of source code under the+terms of section 4, provided that you also meet all of these conditions:++ a) The work must carry prominent notices stating that you modified+ it, and giving a relevant date.++ b) The work must carry prominent notices stating that it is+ released under this License and any conditions added under section+ 7. This requirement modifies the requirement in section 4 to+ "keep intact all notices".++ c) You must license the entire work, as a whole, under this+ License to anyone who comes into possession of a copy. This+ License will therefore apply, along with any applicable section 7+ additional terms, to the whole of the work, and all its parts,+ regardless of how they are packaged. This License gives no+ permission to license the work in any other way, but it does not+ invalidate such permission if you have separately received it.++ d) If the work has interactive user interfaces, each must display+ Appropriate Legal Notices; however, if the Program has interactive+ interfaces that do not display Appropriate Legal Notices, your+ work need not make them do so.++ A compilation of a covered work with other separate and independent+works, which are not by their nature extensions of the covered work,+and which are not combined with it such as to form a larger program,+in or on a volume of a storage or distribution medium, is called an+"aggregate" if the compilation and its resulting copyright are not+used to limit the access or legal rights of the compilation's users+beyond what the individual works permit. Inclusion of a covered work+in an aggregate does not cause this License to apply to the other+parts of the aggregate.++ 6. Conveying Non-Source Forms.++ You may convey a covered work in object code form under the terms+of sections 4 and 5, provided that you also convey the+machine-readable Corresponding Source under the terms of this License,+in one of these ways:++ a) Convey the object code in, or embodied in, a physical product+ (including a physical distribution medium), accompanied by the+ Corresponding Source fixed on a durable physical medium+ customarily used for software interchange.++ b) Convey the object code in, or embodied in, a physical product+ (including a physical distribution medium), accompanied by a+ written offer, valid for at least three years and valid for as+ long as you offer spare parts or customer support for that product+ model, to give anyone who possesses the object code either (1) a+ copy of the Corresponding Source for all the software in the+ product that is covered by this License, on a durable physical+ medium customarily used for software interchange, for a price no+ more than your reasonable cost of physically performing this+ conveying of source, or (2) access to copy the+ Corresponding Source from a network server at no charge.++ c) Convey individual copies of the object code with a copy of the+ written offer to provide the Corresponding Source. This+ alternative is allowed only occasionally and noncommercially, and+ only if you received the object code with such an offer, in accord+ with subsection 6b.++ d) Convey the object code by offering access from a designated+ place (gratis or for a charge), and offer equivalent access to the+ Corresponding Source in the same way through the same place at no+ further charge. You need not require recipients to copy the+ Corresponding Source along with the object code. If the place to+ copy the object code is a network server, the Corresponding Source+ may be on a different server (operated by you or a third party)+ that supports equivalent copying facilities, provided you maintain+ clear directions next to the object code saying where to find the+ Corresponding Source. Regardless of what server hosts the+ Corresponding Source, you remain obligated to ensure that it is+ available for as long as needed to satisfy these requirements.++ e) Convey the object code using peer-to-peer transmission, provided+ you inform other peers where the object code and Corresponding+ Source of the work are being offered to the general public at no+ charge under subsection 6d.++ A separable portion of the object code, whose source code is excluded+from the Corresponding Source as a System Library, need not be+included in conveying the object code work.++ A "User Product" is either (1) a "consumer product", which means any+tangible personal property which is normally used for personal, family,+or household purposes, or (2) anything designed or sold for incorporation+into a dwelling. In determining whether a product is a consumer product,+doubtful cases shall be resolved in favor of coverage. For a particular+product received by a particular user, "normally used" refers to a+typical or common use of that class of product, regardless of the status+of the particular user or of the way in which the particular user+actually uses, or expects or is expected to use, the product. A product+is a consumer product regardless of whether the product has substantial+commercial, industrial or non-consumer uses, unless such uses represent+the only significant mode of use of the product.++ "Installation Information" for a User Product means any methods,+procedures, authorization keys, or other information required to install+and execute modified versions of a covered work in that User Product from+a modified version of its Corresponding Source. The information must+suffice to ensure that the continued functioning of the modified object+code is in no case prevented or interfered with solely because+modification has been made.++ If you convey an object code work under this section in, or with, or+specifically for use in, a User Product, and the conveying occurs as+part of a transaction in which the right of possession and use of the+User Product is transferred to the recipient in perpetuity or for a+fixed term (regardless of how the transaction is characterized), the+Corresponding Source conveyed under this section must be accompanied+by the Installation Information. But this requirement does not apply+if neither you nor any third party retains the ability to install+modified object code on the User Product (for example, the work has+been installed in ROM).++ The requirement to provide Installation Information does not include a+requirement to continue to provide support service, warranty, or updates+for a work that has been modified or installed by the recipient, or for+the User Product in which it has been modified or installed. Access to a+network may be denied when the modification itself materially and+adversely affects the operation of the network or violates the rules and+protocols for communication across the network.++ Corresponding Source conveyed, and Installation Information provided,+in accord with this section must be in a format that is publicly+documented (and with an implementation available to the public in+source code form), and must require no special password or key for+unpacking, reading or copying.++ 7. Additional Terms.++ "Additional permissions" are terms that supplement the terms of this+License by making exceptions from one or more of its conditions.+Additional permissions that are applicable to the entire Program shall+be treated as though they were included in this License, to the extent+that they are valid under applicable law. If additional permissions+apply only to part of the Program, that part may be used separately+under those permissions, but the entire Program remains governed by+this License without regard to the additional permissions.++ When you convey a copy of a covered work, you may at your option+remove any additional permissions from that copy, or from any part of+it. (Additional permissions may be written to require their own+removal in certain cases when you modify the work.) You may place+additional permissions on material, added by you to a covered work,+for which you have or can give appropriate copyright permission.++ Notwithstanding any other provision of this License, for material you+add to a covered work, you may (if authorized by the copyright holders of+that material) supplement the terms of this License with terms:++ a) Disclaiming warranty or limiting liability differently from the+ terms of sections 15 and 16 of this License; or++ b) Requiring preservation of specified reasonable legal notices or+ author attributions in that material or in the Appropriate Legal+ Notices displayed by works containing it; or++ c) Prohibiting misrepresentation of the origin of that material, or+ requiring that modified versions of such material be marked in+ reasonable ways as different from the original version; or++ d) Limiting the use for publicity purposes of names of licensors or+ authors of the material; or++ e) Declining to grant rights under trademark law for use of some+ trade names, trademarks, or service marks; or++ f) Requiring indemnification of licensors and authors of that+ material by anyone who conveys the material (or modified versions of+ it) with contractual assumptions of liability to the recipient, for+ any liability that these contractual assumptions directly impose on+ those licensors and authors.++ All other non-permissive additional terms are considered "further+restrictions" within the meaning of section 10. If the Program as you+received it, or any part of it, contains a notice stating that it is+governed by this License along with a term that is a further+restriction, you may remove that term. If a license document contains+a further restriction but permits relicensing or conveying under this+License, you may add to a covered work material governed by the terms+of that license document, provided that the further restriction does+not survive such relicensing or conveying.++ If you add terms to a covered work in accord with this section, you+must place, in the relevant source files, a statement of the+additional terms that apply to those files, or a notice indicating+where to find the applicable terms.++ Additional terms, permissive or non-permissive, may be stated in the+form of a separately written license, or stated as exceptions;+the above requirements apply either way.++ 8. Termination.++ You may not propagate or modify a covered work except as expressly+provided under this License. Any attempt otherwise to propagate or+modify it is void, and will automatically terminate your rights under+this License (including any patent licenses granted under the third+paragraph of section 11).++ However, if you cease all violation of this License, then your+license from a particular copyright holder is reinstated (a)+provisionally, unless and until the copyright holder explicitly and+finally terminates your license, and (b) permanently, if the copyright+holder fails to notify you of the violation by some reasonable means+prior to 60 days after the cessation.++ Moreover, your license from a particular copyright holder is+reinstated permanently if the copyright holder notifies you of the+violation by some reasonable means, this is the first time you have+received notice of violation of this License (for any work) from that+copyright holder, and you cure the violation prior to 30 days after+your receipt of the notice.++ Termination of your rights under this section does not terminate the+licenses of parties who have received copies or rights from you under+this License. If your rights have been terminated and not permanently+reinstated, you do not qualify to receive new licenses for the same+material under section 10.++ 9. Acceptance Not Required for Having Copies.++ You are not required to accept this License in order to receive or+run a copy of the Program. Ancillary propagation of a covered work+occurring solely as a consequence of using peer-to-peer transmission+to receive a copy likewise does not require acceptance. However,+nothing other than this License grants you permission to propagate or+modify any covered work. These actions infringe copyright if you do+not accept this License. Therefore, by modifying or propagating a+covered work, you indicate your acceptance of this License to do so.++ 10. Automatic Licensing of Downstream Recipients.++ Each time you convey a covered work, the recipient automatically+receives a license from the original licensors, to run, modify and+propagate that work, subject to this License. You are not responsible+for enforcing compliance by third parties with this License.++ An "entity transaction" is a transaction transferring control of an+organization, or substantially all assets of one, or subdividing an+organization, or merging organizations. If propagation of a covered+work results from an entity transaction, each party to that+transaction who receives a copy of the work also receives whatever+licenses to the work the party's predecessor in interest had or could+give under the previous paragraph, plus a right to possession of the+Corresponding Source of the work from the predecessor in interest, if+the predecessor has it or can get it with reasonable efforts.++ You may not impose any further restrictions on the exercise of the+rights granted or affirmed under this License. For example, you may+not impose a license fee, royalty, or other charge for exercise of+rights granted under this License, and you may not initiate litigation+(including a cross-claim or counterclaim in a lawsuit) alleging that+any patent claim is infringed by making, using, selling, offering for+sale, or importing the Program or any portion of it.++ 11. Patents.++ A "contributor" is a copyright holder who authorizes use under this+License of the Program or a work on which the Program is based. The+work thus licensed is called the contributor's "contributor version".++ A contributor's "essential patent claims" are all patent claims+owned or controlled by the contributor, whether already acquired or+hereafter acquired, that would be infringed by some manner, permitted+by this License, of making, using, or selling its contributor version,+but do not include claims that would be infringed only as a+consequence of further modification of the contributor version. For+purposes of this definition, "control" includes the right to grant+patent sublicenses in a manner consistent with the requirements of+this License.++ Each contributor grants you a non-exclusive, worldwide, royalty-free+patent license under the contributor's essential patent claims, to+make, use, sell, offer for sale, import and otherwise run, modify and+propagate the contents of its contributor version.++ In the following three paragraphs, a "patent license" is any express+agreement or commitment, however denominated, not to enforce a patent+(such as an express permission to practice a patent or covenant not to+sue for patent infringement). To "grant" such a patent license to a+party means to make such an agreement or commitment not to enforce a+patent against the party.++ If you convey a covered work, knowingly relying on a patent license,+and the Corresponding Source of the work is not available for anyone+to copy, free of charge and under the terms of this License, through a+publicly available network server or other readily accessible means,+then you must either (1) cause the Corresponding Source to be so+available, or (2) arrange to deprive yourself of the benefit of the+patent license for this particular work, or (3) arrange, in a manner+consistent with the requirements of this License, to extend the patent+license to downstream recipients. "Knowingly relying" means you have+actual knowledge that, but for the patent license, your conveying the+covered work in a country, or your recipient's use of the covered work+in a country, would infringe one or more identifiable patents in that+country that you have reason to believe are valid.+ + If, pursuant to or in connection with a single transaction or+arrangement, you convey, or propagate by procuring conveyance of, a+covered work, and grant a patent license to some of the parties+receiving the covered work authorizing them to use, propagate, modify+or convey a specific copy of the covered work, then the patent license+you grant is automatically extended to all recipients of the covered+work and works based on it.++ A patent license is "discriminatory" if it does not include within+the scope of its coverage, prohibits the exercise of, or is+conditioned on the non-exercise of one or more of the rights that are+specifically granted under this License. You may not convey a covered+work if you are a party to an arrangement with a third party that is+in the business of distributing software, under which you make payment+to the third party based on the extent of your activity of conveying+the work, and under which the third party grants, to any of the+parties who would receive the covered work from you, a discriminatory+patent license (a) in connection with copies of the covered work+conveyed by you (or copies made from those copies), or (b) primarily+for and in connection with specific products or compilations that+contain the covered work, unless you entered into that arrangement,+or that patent license was granted, prior to 28 March 2007.++ Nothing in this License shall be construed as excluding or limiting+any implied license or other defenses to infringement that may+otherwise be available to you under applicable patent law.++ 12. No Surrender of Others' Freedom.++ If conditions are imposed on you (whether by court order, agreement or+otherwise) that contradict the conditions of this License, they do not+excuse you from the conditions of this License. If you cannot convey a+covered work so as to satisfy simultaneously your obligations under this+License and any other pertinent obligations, then as a consequence you may+not convey it at all. For example, if you agree to terms that obligate you+to collect a royalty for further conveying from those to whom you convey+the Program, the only way you could satisfy both those terms and this+License would be to refrain entirely from conveying the Program.++ 13. Use with the GNU Affero General Public License.++ Notwithstanding any other provision of this License, you have+permission to link or combine any covered work with a work licensed+under version 3 of the GNU Affero General Public License into a single+combined work, and to convey the resulting work. The terms of this+License will continue to apply to the part which is the covered work,+but the special requirements of the GNU Affero General Public License,+section 13, concerning interaction through a network will apply to the+combination as such.++ 14. Revised Versions of this License.++ The Free Software Foundation may publish revised and/or new versions of+the GNU General Public License from time to time. Such new versions will+be similar in spirit to the present version, but may differ in detail to+address new problems or concerns.++ Each version is given a distinguishing version number. If the+Program specifies that a certain numbered version of the GNU General+Public License "or any later version" applies to it, you have the+option of following the terms and conditions either of that numbered+version or of any later version published by the Free Software+Foundation. If the Program does not specify a version number of the+GNU General Public License, you may choose any version ever published+by the Free Software Foundation.++ If the Program specifies that a proxy can decide which future+versions of the GNU General Public License can be used, that proxy's+public statement of acceptance of a version permanently authorizes you+to choose that version for the Program.++ Later license versions may give you additional or different+permissions. However, no additional obligations are imposed on any+author or copyright holder as a result of your choosing to follow a+later version.++ 15. Disclaimer of Warranty.++ THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY+APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR+PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM+IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.++ 16. Limitation of Liability.++ IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF+SUCH DAMAGES.++ 17. Interpretation of Sections 15 and 16.++ If the disclaimer of warranty and limitation of liability provided+above cannot be given local legal effect according to their terms,+reviewing courts shall apply local law that most closely approximates+an absolute waiver of all civil liability in connection with the+Program, unless a warranty or assumption of liability accompanies a+copy of the Program in return for a fee.++ END OF TERMS AND CONDITIONS++ How to Apply These Terms to Your New Programs++ If you develop a new program, and you want it to be of the greatest+possible use to the public, the best way to achieve this is to make it+free software which everyone can redistribute and change under these terms.++ To do so, attach the following notices to the program. It is safest+to attach them to the start of each source file to most effectively+state the exclusion of warranty; and each file should have at least+the "copyright" line and a pointer to where the full notice is found.++ <one line to give the program's name and a brief idea of what it does.>+ Copyright (C) <year> <name of author>++ This program is free software: you can redistribute it and/or modify+ it under the terms of the GNU General Public License as published by+ the Free Software Foundation, either version 3 of the License, or+ (at your option) any later version.++ This program is distributed in the hope that it will be useful,+ but WITHOUT ANY WARRANTY; without even the implied warranty of+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the+ GNU General Public License for more details.++ You should have received a copy of the GNU General Public License+ along with this program. If not, see <http://www.gnu.org/licenses/>.++Also add information on how to contact you by electronic and paper mail.++ If the program does terminal interaction, make it output a short+notice like this when it starts in an interactive mode:++ <program> Copyright (C) <year> <name of author>+ This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.+ This is free software, and you are welcome to redistribute it+ under certain conditions; type `show c' for details.++The hypothetical commands `show w' and `show c' should show the appropriate+parts of the General Public License. Of course, your program's commands+might be different; for a GUI interface, you would use an "about box".++ You should also get your employer (if you work as a programmer) or school,+if any, to sign a "copyright disclaimer" for the program, if necessary.+For more information on this, and how to apply and follow the GNU GPL, see+<http://www.gnu.org/licenses/>.++ The GNU General Public License does not permit incorporating your program+into proprietary programs. If your program is a subroutine library, you+may consider it more useful to permit linking proprietary applications with+the library. If this is what you want to do, use the GNU Lesser General+Public License instead of this License. But first, please read+<http://www.gnu.org/philosophy/why-not-lgpl.html>.+
+ Numeric/FFT.hs view
@@ -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
+ Numeric/FFT/Execute.hs view
@@ -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
+ Numeric/FFT/Plan.hs view
@@ -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
+ Numeric/FFT/Special.hs view
@@ -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
+ Numeric/FFT/Special/Miscellaneous.hs view
@@ -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+
+ Numeric/FFT/Special/PowersOfTwo.hs view
@@ -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
+ Numeric/FFT/Special/Primes.hs view
@@ -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
+ Numeric/FFT/Types.hs view
@@ -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)
+ Numeric/FFT/Utils.hs view
@@ -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
+ README.md view
@@ -0,0 +1,4 @@+arb-fft+=======++Pure Haskell arbitrary length FFT library
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ arb-fft.cabal view
@@ -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 <http://www.skybluetrades.net/haskell-fft-index.html>.++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
+ test/basic-test.hs view
@@ -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))
+ test/profile-256.hs view
@@ -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