packages feed

symtegration (empty) → 0.6.1

raw patch · 65 files changed

+6877/−0 lines, 65 filesdep +QuickCheckdep +addep +basesetup-changed

Dependencies added: QuickCheck, ad, base, containers, doctest-parallel, hspec, symtegration, text, text-show

Files

+ CHANGELOG.md view
@@ -0,0 +1,92 @@+# Changelog for `symtegration`++All notable changes to this project will be documented in this file.++The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),+and this project adheres to the+[Haskell Package Versioning Policy](https://pvp.haskell.org/).++## Unreleased++### 0.6.1 - 2025-01-30++*   Do not set `-Werror` by default in preparation for upload to Hackage.++### 0.6.0 - 2025-01-29++*   For rational function integration, use complex logarithms if we are not+    able to derive real function integrals.++*   Add function to map polynomial coefficients monadically.++### 0.5.0 - 2025-01-20++*   Implement integration by parts.++*   Attempt integration by substitution after factoring out constant factors.++*   Prefer positive integers in fraction denominators.++*   Cancel out common integer fractions in $\frac{1}{x} \times y$ as well.++## 0.4.0 - 2025-01-14++*   Integrate more rational functions.++    *   Find all real roots for integration involving solution of cubic equations.++    *   For integration involving solution of quartic equations,+        find real roots for more special cases.++*   Cancel out common integer factors in fractions.++*   Fewer parentheses in Haskell code output.++*   Fewer parentheses in LaTeX output.++*   Test with GHC 9.12.1.++## 0.3.0 - 2025-01-05++*   Implementation of Rioboo's algorithm.++    *   Supports integration of more rational functions.++    *   Integration of rational functions with rational number coefficients now+        only limited by finding solutions for polynomials.  As of yet, only+        rational functions which require solutions for polymials up to degree 2+        can be integrated.++*   Output `pi` as `\pi` in LaTeX.++## 0.2.0 - 2025-01-02++*   Integration of rational functions.++    *   Hermite reduction.++    *   Lazard-Rioboo-Trager integration.++*   Improvements to LaTeX output.++*   Remove simplification based on recursive heuristics,+    which were much more ad hoc.++*   Make `foldTerms` order consistent with simplification order,+    from lower to higher terms.++## 0.1.0 - 2024-12-24++*   Symbolic representation.++*   Simplification.++*   Basic integration support.++    *   Integration of polynomials.++    *   Integration of trigonometric functions.++    *   Integration of exponential and logarithmic functions.++    *   Integration by substitution.
+ LICENSE view
@@ -0,0 +1,201 @@+                                 Apache License+                           Version 2.0, January 2004+                        http://www.apache.org/licenses/++   TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION++   1. Definitions.++      "License" shall mean the terms and conditions for use, reproduction,+      and distribution as defined by Sections 1 through 9 of this document.++      "Licensor" shall mean the copyright owner or entity authorized by+      the copyright owner that is granting the License.++      "Legal Entity" shall mean the union of the acting entity and all+      other entities that control, are controlled by, or are under common+      control with that entity. For the purposes of this definition,+      "control" means (i) the power, direct or indirect, to cause the+      direction or management of such entity, whether by contract or+      otherwise, or (ii) ownership of fifty percent (50%) or more of the+      outstanding shares, or (iii) beneficial ownership of such entity.++      "You" (or "Your") shall mean an individual or Legal Entity+      exercising permissions granted by this License.++      "Source" form shall mean the preferred form for making modifications,+      including but not limited to software source code, documentation+      source, and configuration files.++      "Object" form shall mean any form resulting from mechanical+      transformation or translation of a Source form, including but+      not limited to compiled object code, generated documentation,+      and conversions to other media types.++      "Work" shall mean the work of authorship, whether in Source or+      Object form, made available under the License, as indicated by a+      copyright notice that is included in or attached to the work+      (an example is provided in the Appendix below).++      "Derivative Works" shall mean any work, whether in Source or Object+      form, that is based on (or derived from) the Work and for which the+      editorial revisions, annotations, elaborations, or other modifications+      represent, as a whole, an original work of authorship. For the purposes+      of this License, Derivative Works shall not include works that remain+      separable from, or merely link (or bind by name) to the interfaces of,+      the Work and Derivative Works thereof.++      "Contribution" shall mean any work of authorship, including+      the original version of the Work and any modifications or additions+      to that Work or Derivative Works thereof, that is intentionally+      submitted to Licensor for inclusion in the Work by the copyright owner+      or by an individual or Legal Entity authorized to submit on behalf of+      the copyright owner. For the purposes of this definition, "submitted"+      means any form of electronic, verbal, or written communication sent+      to the Licensor or its representatives, including but not limited to+      communication on electronic mailing lists, source code control systems,+      and issue tracking systems that are managed by, or on behalf of, the+      Licensor for the purpose of discussing and improving the Work, but+      excluding communication that is conspicuously marked or otherwise+      designated in writing by the copyright owner as "Not a Contribution."++      "Contributor" shall mean Licensor and any individual or Legal Entity+      on behalf of whom a Contribution has been received by Licensor and+      subsequently incorporated within the Work.++   2. Grant of Copyright License. Subject to the terms and conditions of+      this License, each Contributor hereby grants to You a perpetual,+      worldwide, non-exclusive, no-charge, royalty-free, irrevocable+      copyright license to reproduce, prepare Derivative Works of,+      publicly display, publicly perform, sublicense, and distribute the+      Work and such Derivative Works in Source or Object form.++   3. Grant of Patent License. Subject to the terms and conditions of+      this License, each Contributor hereby grants to You a perpetual,+      worldwide, non-exclusive, no-charge, royalty-free, irrevocable+      (except as stated in this section) patent license to make, have made,+      use, offer to sell, sell, import, and otherwise transfer the Work,+      where such license applies only to those patent claims licensable+      by such Contributor that are necessarily infringed by their+      Contribution(s) alone or by combination of their Contribution(s)+      with the Work to which such Contribution(s) was submitted. If You+      institute patent litigation against any entity (including a+      cross-claim or counterclaim in a lawsuit) alleging that the Work+      or a Contribution incorporated within the Work constitutes direct+      or contributory patent infringement, then any patent licenses+      granted to You under this License for that Work shall terminate+      as of the date such litigation is filed.++   4. Redistribution. You may reproduce and distribute copies of the+      Work or Derivative Works thereof in any medium, with or without+      modifications, and in Source or Object form, provided that You+      meet the following conditions:++      (a) You must give any other recipients of the Work or+          Derivative Works a copy of this License; and++      (b) You must cause any modified files to carry prominent notices+          stating that You changed the files; and++      (c) You must retain, in the Source form of any Derivative Works+          that You distribute, all copyright, patent, trademark, and+          attribution notices from the Source form of the Work,+          excluding those notices that do not pertain to any part of+          the Derivative Works; and++      (d) If the Work includes a "NOTICE" text file as part of its+          distribution, then any Derivative Works that You distribute must+          include a readable copy of the attribution notices contained+          within such NOTICE file, excluding those notices that do not+          pertain to any part of the Derivative Works, in at least one+          of the following places: within a NOTICE text file distributed+          as part of the Derivative Works; within the Source form or+          documentation, if provided along with the Derivative Works; or,+          within a display generated by the Derivative Works, if and+          wherever such third-party notices normally appear. The contents+          of the NOTICE file are for informational purposes only and+          do not modify the License. You may add Your own attribution+          notices within Derivative Works that You distribute, alongside+          or as an addendum to the NOTICE text from the Work, provided+          that such additional attribution notices cannot be construed+          as modifying the License.++      You may add Your own copyright statement to Your modifications and+      may provide additional or different license terms and conditions+      for use, reproduction, or distribution of Your modifications, or+      for any such Derivative Works as a whole, provided Your use,+      reproduction, and distribution of the Work otherwise complies with+      the conditions stated in this License.++   5. Submission of Contributions. Unless You explicitly state otherwise,+      any Contribution intentionally submitted for inclusion in the Work+      by You to the Licensor shall be under the terms and conditions of+      this License, without any additional terms or conditions.+      Notwithstanding the above, nothing herein shall supersede or modify+      the terms of any separate license agreement you may have executed+      with Licensor regarding such Contributions.++   6. Trademarks. This License does not grant permission to use the trade+      names, trademarks, service marks, or product names of the Licensor,+      except as required for reasonable and customary use in describing the+      origin of the Work and reproducing the content of the NOTICE file.++   7. Disclaimer of Warranty. Unless required by applicable law or+      agreed to in writing, Licensor provides the Work (and each+      Contributor provides its Contributions) on an "AS IS" BASIS,+      WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or+      implied, including, without limitation, any warranties or conditions+      of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A+      PARTICULAR PURPOSE. You are solely responsible for determining the+      appropriateness of using or redistributing the Work and assume any+      risks associated with Your exercise of permissions under this License.++   8. Limitation of Liability. In no event and under no legal theory,+      whether in tort (including negligence), contract, or otherwise,+      unless required by applicable law (such as deliberate and grossly+      negligent acts) or agreed to in writing, shall any Contributor be+      liable to You for damages, including any direct, indirect, special,+      incidental, or consequential damages of any character arising as a+      result of this License or out of the use or inability to use the+      Work (including but not limited to damages for loss of goodwill,+      work stoppage, computer failure or malfunction, or any and all+      other commercial damages or losses), even if such Contributor+      has been advised of the possibility of such damages.++   9. Accepting Warranty or Additional Liability. While redistributing+      the Work or Derivative Works thereof, You may choose to offer,+      and charge a fee for, acceptance of support, warranty, indemnity,+      or other liability obligations and/or rights consistent with this+      License. However, in accepting such obligations, You may act only+      on Your own behalf and on Your sole responsibility, not on behalf+      of any other Contributor, and only if You agree to indemnify,+      defend, and hold each Contributor harmless for any liability+      incurred by, or claims asserted against, such Contributor by reason+      of your accepting any such warranty or additional liability.++   END OF TERMS AND CONDITIONS++   APPENDIX: How to apply the Apache License to your work.++      To apply the Apache License to your work, attach the following+      boilerplate notice, with the fields enclosed by brackets "[]"+      replaced with your own identifying information. (Don't include+      the brackets!)  The text should be enclosed in the appropriate+      comment syntax for the file format. We also recommend that a+      file or class name and description of purpose be included on the+      same "printed page" as the copyright notice for easier+      identification within third-party archives.++   Copyright [yyyy] [name of copyright owner]++   Licensed under the Apache License, Version 2.0 (the "License");+   you may not use this file except in compliance with the License.+   You may obtain a copy of the License at++       http://www.apache.org/licenses/LICENSE-2.0++   Unless required by applicable law or agreed to in writing, software+   distributed under the License is distributed on an "AS IS" BASIS,+   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.+   See the License for the specific language governing permissions and+   limitations under the License.
+ README.md view
@@ -0,0 +1,122 @@+# Symtegration++This is a Haskell library intended to support symbolic integration of mathematical expressions.++It offers the following:++*   Symbolic integration of mathematical expressions.++    *   Integration of polynomials.++    *   Integration of trigonometric functions.++    *   Integration of exponential and logarithmic functions.++    *   Integration of ratios of two polynomials.++    *   Integration by substitution.++    *   Integration by parts.++*   Symbolic representation of mathematical expressions.++*   Utility functions to make it easier to read the mathematical expressions.+    For example, deriving equivalent Haskell code for a mathematical expression,+    and some support for simplifying symbolic representations.++[![Build](https://github.com/symtegration/symtegration/actions/workflows/build.yaml/badge.svg)](https://github.com/symtegration/symtegration/actions/workflows/build.yaml)+[![OpenSSF Best Practices](https://www.bestpractices.dev/projects/9864/badge)](https://www.bestpractices.dev/projects/9864)+[![OpenSSF Scorecard](https://api.scorecard.dev/projects/github.com/symtegration/symtegration/badge)](https://scorecard.dev/viewer/?uri=github.com/symtegration/symtegration)+[![codecov](https://codecov.io/gh/symtegration/symtegration/graph/badge.svg?token=CNBUMA1CKD)](https://codecov.io/gh/symtegration/symtegration)++## Integration++Mathematical expressions with either numeric coefficients or symbolic coefficients+can be integrated.  For example:++```haskell+>>> import Symtegration+>>> toHaskell <$> integrate "x" (4 * "x" ** 3 + 1)+Just "x + x ** 4"+>>> toHaskell <$> integrate "z" ("x" * "z" + "y")+Just "y * z + 1 / 2 * x * z ** 2"+```++Concrete numbers can also be computed from these integrals.  For example:++```haskell+>>> import Symtegration+>>> let (Just p) = integrate "x" (4 * "x" ** 3 + 1)+>>> fractionalEvaluate p (\case "x" -> Just (3 / 7 :: Rational))+Just (1110 % 2401)+```++### Symbolic integration in GHCi++With Symtegration, symbolic integration can be done within [GHCi].+When executing GHCi within the Symtegration project, it is best+to load only the `Symtegration` module to avoid name collisions,+so start GHCi without loading any modules.++```shell+$ stack ghci --no-load+```++Within GHCi, explicitly load the `Symtegration` module.+You can then proceed to symbolically integrate mathematical expressions+and compute approximate or exact values from these integrals.++```haskell+>>> :load Symtegration+>>> toHaskell <$> integrate "x" ("a" * "x" ** 4 + "x" + "b")+Just "b * x + 1 / 2 * x ** 2 + a * (x ** 5) / 5"+>>>+>>> let (Just p) = integrate "x" ("x" ** 2)+>>> evaluate p (\case "x" -> Just 1)+Just 0.3333333333333333+>>>+>>> fractionalEvaluate p (\case "x" -> Just (1 :: Rational))+Just (1 % 3)+```++[GHCi]: https://downloads.haskell.org/ghc/latest/docs/users_guide/ghci.html++### Symbolic integration in IHaskell++Symtegration can also be used in [IHaskell] to do symbolic integration.+Its use can be seen in an [example IHaskell notebook],+which you can try out by [running on mybinder.org].++[IHaskell]: https://github.com/IHaskell/IHaskell+[example IHaskell notebook]: https://github.com/chungyc/haskell-notebooks/blob/main/Symtegration.ipynb+[running on mybinder.org]: https://mybinder.org/v2/gh/chungyc/haskell-notebooks/HEAD?labpath=Symtegration.ipynb++## Changes++See [`CHANGELOG.md`] for what has changed.++[`CHANGELOG.md`]: CHANGELOG.md++## Code of conduct++Be nice; see [`CODE_OF_CONDUCT.md`] for details.++[`CODE_OF_CONDUCT.md`]: docs/CODE_OF_CONDUCT.md++## Security policy++See [`SECURITY.md`] for details.++[`SECURITY.md`]: docs/SECURITY.md++## Contributing++See [`CONTRIBUTING.md`] for details.++[`CONTRIBUTING.md`]: docs/CONTRIBUTING.md++## License++Apache 2.0; see [`LICENSE`] for details.++[`LICENSE`]: LICENSE
+ Setup.hs view
@@ -0,0 +1,3 @@+import Distribution.Simple++main = defaultMain
+ docs/CODE_OF_CONDUCT.md view
@@ -0,0 +1,128 @@+# Contributor Covenant Code of Conduct++## Our Pledge++We as members, contributors, and leaders pledge to make participation in our+community a harassment-free experience for everyone, regardless of age, body+size, visible or invisible disability, ethnicity, sex characteristics, gender+identity and expression, level of experience, education, socio-economic status,+nationality, personal appearance, race, religion, or sexual identity+and orientation.++We pledge to act and interact in ways that contribute to an open, welcoming,+diverse, inclusive, and healthy community.++## Our Standards++Examples of behavior that contributes to a positive environment for our+community include:++* Demonstrating empathy and kindness toward other people+* Being respectful of differing opinions, viewpoints, and experiences+* Giving and gracefully accepting constructive feedback+* Accepting responsibility and apologizing to those affected by our mistakes,+  and learning from the experience+* Focusing on what is best not just for us as individuals, but for the+  overall community++Examples of unacceptable behavior include:++* The use of sexualized language or imagery, and sexual attention or+  advances of any kind+* Trolling, insulting or derogatory comments, and personal or political attacks+* Public or private harassment+* Publishing others' private information, such as a physical or email+  address, without their explicit permission+* Other conduct which could reasonably be considered inappropriate in a+  professional setting++## Enforcement Responsibilities++Community leaders are responsible for clarifying and enforcing our standards of+acceptable behavior and will take appropriate and fair corrective action in+response to any behavior that they deem inappropriate, threatening, offensive,+or harmful.++Community leaders have the right and responsibility to remove, edit, or reject+comments, commits, code, wiki edits, issues, and other contributions that are+not aligned to this Code of Conduct, and will communicate reasons for moderation+decisions when appropriate.++## Scope++This Code of Conduct applies within all community spaces, and also applies when+an individual is officially representing the community in public spaces.+Examples of representing our community include using an official e-mail address,+posting via an official social media account, or acting as an appointed+representative at an online or offline event.++## Enforcement++Instances of abusive, harassing, or otherwise unacceptable behavior may be+reported to the community leaders responsible for enforcement at+dev@chungyc.org.+All complaints will be reviewed and investigated promptly and fairly.++All community leaders are obligated to respect the privacy and security of the+reporter of any incident.++## Enforcement Guidelines++Community leaders will follow these Community Impact Guidelines in determining+the consequences for any action they deem in violation of this Code of Conduct:++### 1. Correction++**Community Impact**: Use of inappropriate language or other behavior deemed+unprofessional or unwelcome in the community.++**Consequence**: A private, written warning from community leaders, providing+clarity around the nature of the violation and an explanation of why the+behavior was inappropriate. A public apology may be requested.++### 2. Warning++**Community Impact**: A violation through a single incident or series+of actions.++**Consequence**: A warning with consequences for continued behavior. No+interaction with the people involved, including unsolicited interaction with+those enforcing the Code of Conduct, for a specified period of time. This+includes avoiding interactions in community spaces as well as external channels+like social media. Violating these terms may lead to a temporary or+permanent ban.++### 3. Temporary Ban++**Community Impact**: A serious violation of community standards, including+sustained inappropriate behavior.++**Consequence**: A temporary ban from any sort of interaction or public+communication with the community for a specified period of time. No public or+private interaction with the people involved, including unsolicited interaction+with those enforcing the Code of Conduct, is allowed during this period.+Violating these terms may lead to a permanent ban.++### 4. Permanent Ban++**Community Impact**: Demonstrating a pattern of violation of community+standards, including sustained inappropriate behavior,  harassment of an+individual, or aggression toward or disparagement of classes of individuals.++**Consequence**: A permanent ban from any sort of public interaction within+the community.++## Attribution++This Code of Conduct is adapted from the [Contributor Covenant][homepage],+version 2.0, available at+https://www.contributor-covenant.org/version/2/0/code_of_conduct.html.++Community Impact Guidelines were inspired by [Mozilla's code of conduct+enforcement ladder](https://github.com/mozilla/diversity).++[homepage]: https://www.contributor-covenant.org++For answers to common questions about this code of conduct, see the FAQ at+https://www.contributor-covenant.org/faq. Translations are available at+https://www.contributor-covenant.org/translations.
+ docs/CONTRIBUTING.md view
@@ -0,0 +1,94 @@+# How to contribute++## Before you begin++### Review community guidelines++This project follows the [Contributor Covenant Code of Conduct].++[Contributor Covenant Code of Conduct]: CODE_OF_CONDUCT.md++### Review license++Any contributions are to be licensed under the [Apache-2.0 license].+Review the license to determine whether you are willing to license+any contributions under the same license.++[Apache-2.0 license]: ../LICENSE++### Background material++_[Symbolic Integration I: Transcendental Functions]_ by Manuel Bronstein+is a primary reference for the algorithms used by this project.++[Symbolic Integration I: Transcendental Functions]: https://doi.org/10.1007/b138171++## Contribution process++### Code reviews++All external contributions require review.+[GitHub pull requests] are used for this purpose.++[GitHub pull requests]: https://docs.github.com/en/pull-requests++### Coding standards++Code should be pure to the extent possible, and partial functions should be avoided.+User-visible entities should be documented with [Haddock], including examples if feasible.+[HLint] should report no issues, and formatting should be according to [Ormolu].++All changes should be accompanied by corresponding tests.+Code should be tested with property-based tests to the extent possible.+This project uses [Hspec] and [QuickCheck] for testing.+Examples in the Haddock documentation are tested using [`doctest-parallel`].++[Haddock]: https://haskell-haddock.readthedocs.io/+[HLint]: https://github.com/ndmitchell/hlint+[Ormolu]: https://github.com/tweag/ormolu+[Hspec]: https://hspec.github.io/+[QuickCheck]: https://hackage.haskell.org/package/QuickCheck+[`doctest-parallel`]: https://github.com/martijnbastiaan/doctest-parallel++All warnings are enabled for builds.+If a certain warning is unavoidable, it should only be disabled on a per file basis.+While the warnings are not errors by default, code with compiler warnings will not+be merged, and the continuous build upgrades these to errors.+To upgrade compiler warnings to errors locally, use the `--pedantic` flag.++```bash+$ stack build --pedantic+$ stack test --pedantic+```++### Dependencies++This project aims to avoid using too many heavy dependencies.+Care should be taken not to add dependencies casually.+If the same thing can be done with some additional code in the project,+then adding a dependency should be avoided.+This is more important the more heavy a dependency is or the less maintained it is.++### Releases++When releasing, these files should be updated:++*   [`CHANGELOG.md`] with user-visible changes.++*   [`package.yaml`] with the new version.  There should be at least one+    subsequent `stack build` to update [`symtegration.cabal`] as well.++Versioning is based on [semantic versioning] and the [Haskell package versioning policy].+When there are differences between the two policies, the latter takes precedence.++Lower version bounds for dependencies should be verified by setting the versions+to the lowest minor versions in the Cabal configuration and checking that+builds and tests are still successful.  These changes to the Cabal configuration+are only for confirming that the lower bounds are still valid, and should not+be submitted to the repository.++[`CHANGELOG.md`]: ../CHANGELOG.md+[`package.yaml`]: ../package.yaml+[`symtegration.cabal`]: ../symtegration.cabal+[semantic versioning]: https://semver.org/+[Haskell package versioning policy]: https://pvp.haskell.org/
+ docs/SECURITY.md view
@@ -0,0 +1,10 @@+# Security Policy++## Reporting a vulnerability++To report a security bug, submit details such as what the vulnerability is,+what risks it may entail, and how to reproduce it [via GitHub].+I will try to respond and deal with the issue within 2 weeks.+Disclosure policy is to disclose a vulnerability within 90 days of it being reported.++[via GitHub]: https://github.com/chungyc/symtegration/security/advisories/new
+ src/Symtegration.hs view
@@ -0,0 +1,112 @@+-- |+-- Module: Symtegration+-- Description: Library for symbolic integration of mathematical expressions.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+--+-- Symtegration is a library for symbolic integration of mathematical expressions.+-- For normal use, this is the only module which needs to be loaded.+-- Other modules are used for finer control over what happens,+-- or for supporting the work that yet other modules do.+--+-- For example, with \(\int (4x^3 + 1) \, dx = x^4 + x\):+--+-- >>> import Symtegration+-- >>> toHaskell <$> integrate "x" (4 * "x" ** 3 + 1)+-- Just "x + x ** 4"+--+-- For another example, with \(\int (xz+y) \, dz = \frac{xz^2}{2} + yz\):+--+-- >>> import Symtegration+-- >>> toHaskell <$> integrate "z" ("x" * "z" + "y")+-- Just "y * z + 1 / 2 * x * z ** 2"+module Symtegration+  ( -- * Symbolic representation+    Expression,++    -- * Integration+    integrate,++    -- * Differentiation+    differentiate,++    -- * Computation+    evaluate,+    fractionalEvaluate,+    toFunction,++    -- * Conversion+    toHaskell,+    toLaTeX,++    -- * Simplification++    -- | When using only this module, explicitly simplifying mathematical expressions+    -- should usually not be necessary, since the exported functions automatically+    -- simplify results as appropriate.  One may want to explicitly simplify+    -- mathematical expressions when used with other packages, however,+    -- such as when using [Numeric.AD](https://hackage.haskell.org/package/ad)+    -- directly for differentiation.+    simplify,+    tidy,+  )+where++import Data.Text (Text)+import Symtegration.Differentiation (differentiate)+import Symtegration.Integration qualified as Integration+import Symtegration.Symbolic (Expression, evaluate, fractionalEvaluate, toFunction)+import Symtegration.Symbolic.Haskell (toHaskell)+import Symtegration.Symbolic.LaTeX (toLaTeX)+import Symtegration.Symbolic.Simplify (simplify, simplifyForVariable)+import Symtegration.Symbolic.Simplify.Tidy (tidy)++-- |+-- Returns the indefinite integral of a mathematical expression given+-- its symbolic representation.  It will return 'Nothing' if it is+-- unable to derive an integral.  The indefinite integral will be+-- simplified to a certain extent.+--+-- For example, with \(\int (4x^3 + 1) \, dx = x^4 + x\)+-- where all the coefficients are numbers:+--+-- >>> toHaskell <$> integrate "x" (4 * "x" ** 3 + 1)+-- Just "x + x ** 4"+--+-- It can also return indefinite integrals when the coefficients+-- are symbolic, as with \(\int (xz+y) \, dz = \frac{xz^2}{2} + yz\):+--+-- >>> toHaskell <$> integrate "z" ("x" * "z" + "y")+-- Just "y * z + 1 / 2 * x * z ** 2"+--+-- === __Definite integrals__+--+-- If the indefinite integral \(F = \int f(x) \, dx\) is continuous,+-- then the definite integral is+--+-- \[ \int_a^b f(x) \, dx = F(b) - F(a) \]+--+-- This is /not/ true in general if \(F\) is not continuous in the integral interval.+-- Care must be taken when computing a definite integral from an indefinite integral+-- which is not continuous.  For example, an indefinite integral such as the following+--+-- \[ \int f \, dx = F = \sum_{\alpha \mid 4 \alpha^2 + 1 = 0} \alpha \log (x^3 +2 \alpha x^2 - 3x - 4\alpha) \]+--+-- uses complex logarithms, where \(\alpha = \pm \frac{i}{2}\) and there are discontinuities at \(x=-\sqrt{2}\) and \(x=\sqrt{2}\).+--+-- Definite integrals for such cases can be handled by integrating over continuous intervals separately.+-- For example,+--+-- \[ \int_1^2 f \, dx = \left( F(2) - \lim_{x \rightarrow \sqrt{2}^+} F(x) \right) + \left( \lim_{x \rightarrow \sqrt{2}^-} F(x) - F(1) \right) \]+--+-- Symtegration will return real function integrals if it can,+-- but may return complex function integrals instead if it is unable to.+integrate ::+  -- | The symbol representing the variable being integrated over.+  Text ->+  -- | The mathematical expression being integrated.+  Expression ->+  -- | The indefinite integral, if derived.+  Maybe Expression+integrate var expr = tidy . simplifyForVariable var <$> Integration.integrate var expr
+ src/Symtegration/Differentiation.hs view
@@ -0,0 +1,42 @@+-- |+-- Module: Symtegration.Differentiation+-- Description: Differentiate mathematical expressions.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+--+-- Differentiate symbolic representations of mathematical expressions.+-- This module does not actually implement differentiation,+-- but is rather a thin wrapper over "Numeric.AD" providing+-- derivatives for 'Expression' with some simplification applied.+module Symtegration.Differentiation (differentiate) where++import Data.Text (Text)+import Numeric.AD.Rank1.Forward+import Symtegration.Symbolic+import Symtegration.Symbolic.Simplify++-- $setup+-- >>> import Symtegration.Symbolic.Haskell++-- | Differentiates a mathematical expression.+--+-- >>> toHaskell $ differentiate "x" $ "x" ** 2+-- "2 * x"+-- >>> toHaskell $ differentiate "x" $ "a" * sin "x"+-- "a * cos x"+--+-- This uses [Numeric.AD](https://hackage.haskell.org/package/ad).+differentiate ::+  -- | Symbol representing the variable.+  Text ->+  -- | Symbolic representation of the mathematical expression to differentiate.+  Expression ->+  -- | The derivative.+  Expression+differentiate v e = tidy $ simplifyForVariable v $ diff f $ Symbol v+  where+    f = toFunction e assign+    assign x+      | v == x = id+      | otherwise = const $ auto $ Symbol x
+ src/Symtegration/Integration.hs view
@@ -0,0 +1,47 @@+-- |+-- Module: Symtegration.Integration+-- Description: Symbolically integrates mathematical expressions.+-- Copyright: Copyright 2025 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Integration (integrate) where++import Data.Foldable (asum)+import Data.Text (Text)+import Symtegration.Integration.Exponential qualified as Exponential+import Symtegration.Integration.Parts qualified as Parts+import Symtegration.Integration.Powers qualified as Powers+import Symtegration.Integration.Rational qualified as Rational+import Symtegration.Integration.Substitution qualified as Substitution+import Symtegration.Integration.Sum qualified as Sum+import Symtegration.Integration.Term qualified as Term+import Symtegration.Integration.Trigonometric qualified as Trigonometric+import Symtegration.Symbolic+import Symtegration.Symbolic.Simplify++-- |+-- Return the indefinite integral of a mathematical expression given+-- its symbolic representation.  It will return 'Nothing' if it is+-- unable to derive an integral.  This will not apply any simplification.+integrate :: Text -> Expression -> Maybe Expression+integrate v e = asum $ map (\f -> f v e') withTermSum+  where+    e' = simplifyForVariable v e++-- | Functions which directly integrate.+base :: [Text -> Expression -> Maybe Expression]+base = [Powers.integrate, Exponential.integrate, Trigonometric.integrate, Rational.integrate]++-- | Includes integration of a term using other integration functions.+withTerm :: [Text -> Expression -> Maybe Expression]+withTerm =+  base+    ++ [ Term.integrate base,+         Substitution.integrate base,+         Parts.integrate [Term.integrate base],+         Term.integrate [Substitution.integrate base, Parts.integrate [Term.integrate base]]+       ]++-- | Includes integration of a sum of terms.+withTermSum :: [Text -> Expression -> Maybe Expression]+withTermSum = withTerm ++ [Sum.integrate withTerm]
+ src/Symtegration/Integration/Exponential.hs view
@@ -0,0 +1,41 @@+-- |+-- Module: Symtegration.Integration.Exponential+-- Description: Basic integration of exponential and logarithmic functions.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+--+-- Supports basic integration of exponential and logarithmic functions.+-- This does not support the integration of anything else,+-- even if it is trivial like integrating a constant.+module Symtegration.Integration.Exponential (integrate) where++import Data.Text (Text)+import Symtegration.Symbolic++-- $setup+-- >>> import Symtegration.Symbolic.Haskell++-- | Integrates exponential and logarithmic functions required by the 'Floating' type class.+--+-- >>> toHaskell <$> integrate "x" (exp "x")+-- Just "exp x"+-- >>> toHaskell <$> integrate "x" (log "x")+-- Just "x * log x - x"+integrate :: Text -> Expression -> Maybe Expression+integrate _ (Number _) = Nothing+integrate _ (Symbol _) = Nothing+integrate v e@(Exp' (Symbol s))+  | v == s = Just e+  | otherwise = Nothing+integrate v (Log' e@(Symbol s))+  | v == s = Just $ (e :*: Log' e) :-: e+  | otherwise = Nothing+integrate v e@(Number n :**: Symbol s)+  | v == s = Just $ c :*: e+  | otherwise = Nothing+  where+    c = Number 1 :/: Log' (Number n)+integrate v (LogBase' (Number n) (Symbol s))+  | v == s = fmap (\x -> x :/: Log' (Number n)) $ integrate v $ Log' (Symbol s)+integrate _ _ = Nothing
+ src/Symtegration/Integration/Factor.hs view
@@ -0,0 +1,69 @@+-- |+-- Module: Symtegration.Integration.Factor+-- Description: Factor a term into constant and non-constant parts.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Integration.Factor (factor, isConstant) where++import Data.Text (Text)+import Symtegration.Symbolic+import Symtegration.Symbolic.Simplify++-- $setup+-- >>> import Symtegration++-- | Factor a multiplicative term into a constant portion and the variable-dependent portion.+-- E.g., \(2a x \sin x\) into \(2a\) and \(x \sin x\) when the variable is \(x\).+--+-- >>> let s (x, y) = (toHaskell $ simplify x, toHaskell $ simplify y)+-- >>> s $ factor "x" $ 2 * ("a" * sin "x")+-- ("2 * a","sin x")+-- >>> s $ factor "x" $ "a" / "x"+-- ("a","1 / x")+--+-- Assumes algebraic ring ordering has been applied to the term.+factor ::+  -- | Symbol for the variable.+  Text ->+  -- | Term to separate into constant and non-constant portions.+  Expression ->+  (Expression, Expression)+factor _ e@(Number _) = (e, Number 1)+factor v e@(Symbol s) | v == s = (Number 1, e) | otherwise = (e, Number 1)+factor v e@(UnaryApply _ x) | isConstant v x = (e, Number 1) | otherwise = (Number 1, e)+factor v e@(x :*: (y :*: z))+  | isConstant v x, isConstant v y, isConstant v z = (e, Number 1)+  | isConstant v x, isConstant v y = (simplifyForVariable v $ x :*: (y :*: c), z')+  | isConstant v x = (simplifyForVariable v $ x :*: d, y')+  | otherwise = (Number 1, e)+  where+    (c, z') = factor v z+    (d, y') = factor v $ y :*: z+factor v e@(x :*: y)+  | isConstant v x, isConstant v y = (e, Number 1)+  | isConstant v x = (x, y)+  | otherwise = (Number 1, e)+factor v (x :/: y) = (simplify $ constX :/: constY, simplify $ varX :/: varY)+  where+    (constX, varX) = factor v x+    (constY, varY) = factor v y+factor v e | isConstant v e = (e, Number 1) | otherwise = (Number 1, e)++-- | Returns whether an expression contains the variable.+--+-- >>> isConstant "x" $ 1 + "x"+-- False+-- >>> isConstant "x" $ 1 + "a"+-- True+isConstant ::+  -- | Symbol for the variable.+  Text ->+  -- | Expression to check.+  Expression ->+  -- | Whether the expression is a constant.+  Bool+isConstant _ (Number _) = True+isConstant v (Symbol s) = s /= v+isConstant v (UnaryApply _ x) = isConstant v x+isConstant v (BinaryApply _ x y) = isConstant v x && isConstant v y
+ src/Symtegration/Integration/Parts.hs view
@@ -0,0 +1,64 @@+-- |+-- Module: Symtegration.Integration.Parts+-- Description: Integration by parts.+-- Copyright: Copyright 2025 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Integration.Parts (integrate) where++import Control.Applicative (asum, (<|>))+import Data.Text (Text)+import Symtegration.Differentiation+import Symtegration.Symbolic+import Symtegration.Symbolic.Simplify++-- $setup+-- >>> import Symtegration.Integration.Powers qualified as Powers+-- >>> import Symtegration.Integration.Trigonometric qualified as Trigonometric+-- >>> import Symtegration.Symbolic.Haskell+-- >>> import Symtegration.Symbolic.Simplify++-- | Integrates by parts.+--+-- Specifically, if for+--+-- \[ \int f g \, dx \]+--+-- it is the case that we can find \(F = \int f \, dx\) and \(\int F \frac{dg}{dx} \, dx\),+-- then we can derive the integral as+--+-- \[ \int f g \, dx = F g - \int F \frac{dg}{dx} \, dx \]+--+-- >>> let directMethods = [Powers.integrate, Trigonometric.integrate]+-- >>> toHaskell . simplify <$> integrate directMethods "x" ("x" * cos "x")+-- Just "x * sin x + cos x"+integrate ::+  -- | Integration algorithms to try on the parts.+  [Text -> Expression -> Maybe Expression] ->+  -- | Symbol for the variable.+  Text ->+  -- | Expression to integrate.+  Expression ->+  -- | Integral, if derived.+  Maybe Expression+integrate fs v (x :*: y) = integrate' fs v x y <|> integrate' fs v y x+integrate _ _ _ = Nothing++-- | The actual work of integrating by parts, except it tries the parts in only one order.+integrate' ::+  -- | Integration algorithms to try on the parts.+  [Text -> Expression -> Maybe Expression] ->+  -- | Symbol for the variable.+  Text ->+  -- | The part to be integrated.+  Expression ->+  -- | The part to be differentiated.+  Expression ->+  -- | Integral, if derived.+  Maybe Expression+integrate' fs v x y = do+  ix <- integrate'' x+  iixdy <- integrate'' $ simplifyForVariable v $ ix * differentiate v y+  return $ ix * y - iixdy+  where+    integrate'' z = asum $ map (\f -> f v z) fs
+ src/Symtegration/Integration/Powers.hs view
@@ -0,0 +1,45 @@+-- |+-- Module: Symtegration.Integration.Powers+-- Description: Integration of arbitrary powers of a variable.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Integration.Powers (integrate) where++import Data.Text (Text)+import Symtegration.Symbolic++-- $setup+-- >>> import Symtegration.Symbolic.Haskell+-- >>> import Symtegration.Symbolic.Simplify++-- | Integrates powers of a variable.+-- In other words, expressions of the form \(x^c\),+-- where \(c\) is a constant.+--+-- >>> toHaskell . simplify <$> integrate "x" "x"+-- Just "1 / 2 * x ** 2"+-- >>> toHaskell . simplify <$> integrate "x" ("x" ** (1/2))+-- Just "(2 * x ** (3 / 2)) / 3"+-- >>> toHaskell . simplify <$> integrate "x" ("x" ** (-1))+-- Just "log x"+integrate :: Text -> Expression -> Maybe Expression+integrate v (1 :/: Symbol s) =+  integrate v $ Symbol s :**: Number (-1)+integrate v (x :**: (Negate' (Number n :/: Number m))) =+  integrate v $ x :**: (Number (-n) :/: Number m)+integrate v (x :**: (Negate' (Number n))) =+  integrate v $ x :**: Number (-n)+integrate v e@(Number _) = Just $ e :*: Symbol v+integrate v e@(Symbol v')+  | v == v' = Just $ (Number 1 :/: Number 2) :*: (e :**: 2)+  | otherwise = Just $ e :*: Symbol v+integrate v (x@(Symbol s) :**: Number n)+  | s == v, -1 <- n = Just $ Log' x+  | s == v = Just $ (x :**: Number (n + 1)) :/: Number (n + 1)+  | otherwise = Nothing+integrate _ (_ :**: (_ :/: Number 0)) = Nothing+integrate v (x@(Symbol s) :**: y@(Number _ :/: Number _))+  | s == v = Just $ (x :**: (y :+: 1)) :/: (y :+: 1)+  | otherwise = Nothing+integrate _ _ = Nothing
+ src/Symtegration/Integration/Rational.hs view
@@ -0,0 +1,607 @@+-- |+-- Module: Symtegration.Integration.Rational+-- Description: Integration of rational functions.+-- Copyright: Copyright 2025 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+--+-- Integrates rational functions.+-- Rational functions are ratios of two polynomials, not functions of rational numbers.+-- Only rational number coefficients are supported.+module Symtegration.Integration.Rational+  ( -- * Integration+    integrate,++    -- * Algorithms++    -- | Algorithms used for integrating rational functions.+    hermiteReduce,+    rationalIntegralLogTerms,+    complexLogTermToAtan,+    complexLogTermToRealTerm,++    -- * Support++    -- | Functions and types useful when integrating rational functions.+    toRationalFunction,+    RationalFunction (..),+  )+where++import Data.Foldable (asum)+import Data.List (find, intersect)+import Data.Monoid (Sum (..))+import Data.Text (Text)+import Symtegration.Polynomial hiding (integrate)+import Symtegration.Polynomial qualified as Polynomial+import Symtegration.Polynomial.Indexed+import Symtegration.Polynomial.Solve+import Symtegration.Polynomial.Symbolic+import Symtegration.Symbolic+import Symtegration.Symbolic.Simplify++-- $setup+-- >>> :set -w+-- >>> import Symtegration.Polynomial hiding (integrate)+-- >>> import Symtegration.Polynomial.Indexed+-- >>> import Symtegration.Symbolic.Haskell+-- >>> import Symtegration.Symbolic.Simplify++-- | Integrate a ratio of two polynomials with rational number coefficients.+--+-- For example,+--+-- >>> let p = "x" ** 7 - 24 * "x" ** 4 - 4 * "x" ** 2 + 8 * "x" - 8+-- >>> let q = "x" ** 8 + 6 * "x" ** 6 + 12 * "x" ** 4 + 8 * "x" ** 2+-- >>> toHaskell . simplify <$> integrate "x" (p / q)+-- Just "3 / (2 + x ** 2) + (4 + 8 * x ** 2) / (4 * x + 4 * x ** 3 + x ** 5) + log x"+--+-- so that+--+-- \[\int \frac{x^7-24x^4-4x^2+8x-8}{x^8+6x^6+12x^4+8x^2} \, dx = \frac{3}{x^2+2} + \frac{8x^2+4}{x^5+4x^3+4x} + \log x\]+--+-- For another example,+--+-- >>> let f = 36 / ("x" ** 5 - 2 * "x" ** 4 - 2 * "x" ** 3 + 4 * "x" ** 2 + "x" - 2)+-- >>> toHaskell . simplify <$> integrate "x" f+-- Just "(-4) * log (8 + 8 * x) + 4 * log (16 + (-8) * x) + (6 + 12 * x) / ((-1) + x ** 2)"+--+-- so that+--+-- \[\int \frac{36}{x^5-2x^4-2x^3+4x^2+x-2} \, dx = \frac{12x+6}{x^2-1} + 4 \log \left( x - 2 \right) - 4 \log \left( x + 1 \right)\]+--+-- This function will attempt to find a real function integral if it can,+-- but if it cannot, it will try to find an integral which includes complex logarithms.+integrate :: Text -> Expression -> Maybe Expression+integrate v e+  | (x :/: y) <- e',+    (Just n) <- fromExpression (forVariable v) x,+    (Just d) <- fromExpression (forVariable v) y,+    d /= 0 =+      integrate' n d+  | otherwise = Nothing+  where+    e' = simplifyForVariable v e+    integrate' n d = (+) reduced . (+) poly <$> logs+      where+        -- Integrals directly from Hermite reduction.+        (g, h) = hermiteReduce $ toRationalFunction n d+        reduced = sum $ map fromRationalFunction g++        -- Integrate polynomials left over from the Hermite reduction.+        RationalFunction numer denom = h+        (q, r) = numer `divide` denom+        poly = toExpression v toRationalCoefficient $ Polynomial.integrate q++        -- Derive the log terms in the integral.+        h' = toRationalFunction r denom+        logTerms = rationalIntegralLogTerms h'+        logs = asum [realLogs, complexLogs] :: Maybe Expression++        -- Try to integrate into real functions first.+        realLogs+          | (Just terms) <- logTerms = sum <$> toMaybeList (map (complexLogTermToRealExpression v) terms)+          | otherwise = Nothing++        -- If it cannot be integrated into real functions, allow complex logarithms.+        complexLogs+          | (Just terms) <- logTerms = sum <$> toMaybeList (map (complexLogTermToComplexExpression v) terms)+          | otherwise = Nothing++        fromRationalFunction (RationalFunction u w) = u' / w'+          where+            u' = toExpression v toRationalCoefficient u+            w' = toExpression v toRationalCoefficient w++-- | Represents the ratio of two polynomials with rational number coefficients.+data RationalFunction = RationalFunction IndexedPolynomial IndexedPolynomial+  deriving (Eq)++instance Show RationalFunction where+  show (RationalFunction n d) = "(" <> show n <> ") / (" <> show d <> ")"++-- | The numerator and denominator in the results+-- for '(+)', '(-)', '(*)', and 'negate' will be coprime.+instance Num RationalFunction where+  (RationalFunction x y) + (RationalFunction u v) =+    toRationalFunction (x * v + u * y) (y * v)++  (RationalFunction x y) - (RationalFunction u v) =+    toRationalFunction (x * v - u * y) (y * v)++  (RationalFunction x y) * (RationalFunction u v) =+    toRationalFunction (x * u) (y * v)++  abs = id++  signum 0 = 0+  signum _ = 1++  fromInteger n = RationalFunction (fromInteger n) 1++instance Fractional RationalFunction where+  fromRational q = RationalFunction (scale q 1) 1+  recip (RationalFunction p q) = RationalFunction q p++-- | Form a rational function from two polynomials.+-- The polynomials will be reduced so that the numerator and denominator are coprime.+toRationalFunction ::+  -- | Numerator.+  IndexedPolynomial ->+  -- | Denominator.+  IndexedPolynomial ->+  RationalFunction+toRationalFunction x 0 = RationalFunction x 0+toRationalFunction x y = RationalFunction x' y'+  where+    g = monic $ greatestCommonDivisor x y+    (x', _) = x `divide` g+    (y', _) = y `divide` g++-- | Applies Hermite reduction to a rational function.+-- Returns a list of rational functions whose sums add up to the integral+-- and a rational function which remains to be integrated.+-- Only rational functions with rational number coefficients and+-- where the numerator and denominator are coprime are supported.+--+-- Specifically, for rational function \(f = \frac{A}{D}\),+-- where \(A\) and \(D\) are coprime polynomials, then for return value @(gs, h)@,+-- the sum of @gs@ is equal to \(g\) and @h@ is equal to \(h\) in the following:+--+-- \[ \frac{A}{D} = \frac{dg}{dx} + h \]+--+-- This is equivalent to the following:+--+-- \[ \int \frac{A}{D} \, dx = g + \int h \, dx \]+--+-- If preconditions are satisfied, i.e., \(D \neq 0\) and \(A\) and \(D\) are coprime,+-- then \(h\) will have a squarefree denominator.+--+-- For example,+--+-- >>> let p = power 7 - 24 * power 4 - 4 * power 2 + 8 * power 1 - 8 :: IndexedPolynomial+-- >>> let q = power 8 + 6 * power 6 + 12 * power 4 + 8 * power 2 :: IndexedPolynomial+-- >>> hermiteReduce $ toRationalFunction p q+-- ([(3) / (x^2 + 2),(8x^2 + 4) / (x^5 + 4x^3 + 4x)],(1) / (x))+--+-- so that+--+-- \[\int \frac{x^7-24x^4-4x^2+8x-8}{x^8+6x^6+12x^4+8x^2} \, dx = \frac{3}{x^2+2}+\frac{8x^2+4}{x^5+4x^3+4x}+\int \frac{1}{x} \, dx\]+--+-- \(g\) is returned as a list of rational functions which sum to \(g\)+-- instead of a single rational function, because the former could sometimes+-- be simpler to read.+hermiteReduce :: RationalFunction -> ([RationalFunction], RationalFunction)+hermiteReduce h@(RationalFunction _ 0) = ([], h)+hermiteReduce h@(RationalFunction x y)+  | (Just z) <- reduce x [] common = z+  | otherwise = ([], h) -- Should never happen, but a fallback if it does.+  where+    common = monic $ greatestCommonDivisor y $ differentiate y+    (divisor, _) = y `divide` common+    reduce a g d+      | degree d > 0 = do+          let d' = monic $ greatestCommonDivisor d $ differentiate d+          let (d'', _) = d `divide` d'+          let (d''', _) = (divisor * differentiate d) `divide` d+          (b, c) <- diophantineEuclidean (-d''') d'' a+          let (b', _) = (differentiate b * divisor) `divide` d''+          let a' = c - b'+          let g' = toRationalFunction b d : g+          reduce a' g' d'+      | otherwise = Just (g, toRationalFunction a divisor)++-- | For rational function \(\frac{A}{D}\), where \(\deg(A) < \deg(D)\),+-- and \(D\) is non-zero, squarefree, and coprime with \(A\),+-- returns the components which form the logarithmic terms of \(\int \frac{A}{D} \, dx\).+-- Specifically, when a list of \((Q_i(t), S_i(t, x))\) is returned,+-- where \(Q_i(t)\) are polynomials of \(t\) and \(S_i(t, x)\) are polynomials of \(x\)+-- with coefficients formed from polynomials of \(t\), then+--+-- \[+-- \int \frac{A}{D} \, dx = \sum_{i=1}^n \sum_{a \in \{t \mid Q_i(t) = 0\}} a \log \left(S_i(a,x)\right)+-- \]+--+-- For example,+--+-- >>> let p = power 4 - 3 * power 2 + 6 :: IndexedPolynomial+-- >>> let q = power 6 - 5 * power 4 + 5 * power 2 + 4 :: IndexedPolynomial+-- >>> let f = toRationalFunction p q+-- >>> let gs = rationalIntegralLogTerms f+-- >>> length <$> gs+-- Just 1+-- >>> fst . head <$> gs+-- Just x^2 + (1 % 4)+-- >>> foldTerms (\e c -> show (e, c) <> " ") . snd . head <$> gs+-- Just "(0,792x^2 + (-16)) (1,(-2440)x^3 + 32x) (2,(-400)x^2 + 7) (3,800x^3 + (-14)x) "+--+-- so it is the case that+--+-- \[+-- \int \frac{x^4-3x^2+6}{x^6-5x^4+5x^2+4} \, dx+-- = \sum_{a \mid a^2+\frac{1}{4} = 0} a \log \left( (800a^3-14a)x^3+(-400a^2+7)x^2+(-2440a^3+32a)x + 792a^2-16 \right)+-- \]+--+-- It may return 'Nothing' if \(\frac{A}{D}\) is not in the expected form.+rationalIntegralLogTerms ::+  RationalFunction ->+  Maybe [(IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial)]+rationalIntegralLogTerms (RationalFunction a d) = do+  -- For A/D, get the resultant and subresultant polynomial remainder sequence+  -- for D and (A - t * D').+  let sa = mapCoefficients fromRational a+  let sd = mapCoefficients fromRational d+  let t = RationalFunction (power 1) 1+  let (resultant, prs) = subresultant sd $ sa - scale t (differentiate sd)++  -- Turn rational functions into polynomials if possible.+  -- When the preconditions are satisfied, these should all be polynomials.+  sd' <- mapCoefficientsM toPoly sd+  resultant' <- toPoly resultant+  prs' <- toMaybeList $ map (mapCoefficientsM toPoly) prs :: Maybe [IndexedPolynomialWith IndexedPolynomial]++  -- Derive what make up the log terms in the integral.+  let qs = squarefree resultant' :: [IndexedPolynomial]+  let terms = zipWith (toTerm sd' prs') [1 ..] qs++  -- Ignore log terms which end up being multiples of 0 = log 1.+  return $ filter ((/=) 1 . snd) terms+  where+    toTerm ::+      IndexedPolynomialWith IndexedPolynomial ->+      [IndexedPolynomialWith IndexedPolynomial] ->+      Int ->+      IndexedPolynomial ->+      (IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial)+    toTerm sd prs i q+      | degree q == 0 = (q, 1)+      | i == degree d = (q, sd)+      | (Just r) <- find ((==) i . degree) prs = derive q r+      | otherwise = (q, 1)++    derive ::+      IndexedPolynomial ->+      IndexedPolynomialWith IndexedPolynomial ->+      (IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial)+    derive q s = (q, s')+      where+        as = squarefree $ leadingCoefficient s+        s' = foldl scalePoly s (zip ([1 ..] :: [Int]) as)+          where+            scalePoly x (j, u) =+              getSum $ foldTerms (reduceTerm (monic $ greatestCommonDivisor u q ^ j)) x+            reduceTerm v e c = Sum $ scale (exactDivide c v) $ power e+            exactDivide u v = r+              where+                (r, _) = u `divide` v++-- | Given polynomials \(A\) and \(B\),+-- return a sum \(f\) of inverse tangents such that the following is true.+--+-- \[+-- \frac{df}{dx} = \frac{d}{dx} i \log \left( \frac{A + iB}{A - iB} \right)+-- \]+--+-- This allows integrals to be evaluated with only real-valued functions.+-- It also avoids the discontinuities in real-valued indefinite integrals which may result+-- when the integral uses logarithms with complex arguments.+--+-- For example,+--+-- >>> toHaskell $ simplify $ complexLogTermToAtan "x" (power 3 - 3 * power 1) (power 2 - 2)+-- "2 * atan x + 2 * atan ((x + (-3) * x ** 3 + x ** 5) / 2) + 2 * atan (x ** 3)"+--+-- so it is the case that+--+-- \[ \frac{d}{dx} \left( i \log \left( \frac{(x^3-3x) + i(x^2-2)}{(x^3-3x) - i(x^2-2)} \right) \right) =+-- \frac{d}{dx} \left( 2 \tan^{-1} \left(\frac{x^5-3x^3+x}{2}\right) + 2 \tan^{-1} \left(x^3\right) + 2 \tan^{-1} x \right) \]+complexLogTermToAtan ::+  -- | Symbol for the variable.+  Text ->+  -- | Polynomial \(A\).+  IndexedPolynomial ->+  -- | Polynomial \(B\).+  IndexedPolynomial ->+  -- | Sum \(f\) of inverse tangents.+  Expression+complexLogTermToAtan v a b+  | r == 0 = 2 * atan (a' / b')+  | degree a < degree b = complexLogTermToAtan v (-b) a+  | otherwise = 2 * atan (s' / g') + complexLogTermToAtan v d c+  where+    (_, r) = a `divide` b+    (d, c, g) = extendedEuclidean b (-a)+    a' = toExpression v toRationalCoefficient a+    b' = toExpression v toRationalCoefficient b+    g' = toExpression v toRationalCoefficient g+    s' = toExpression v toRationalCoefficient $ a * d + b * c++-- | For the ingredients of a complex logarithm, return the ingredients of an equivalent real function in terms of an indefinite integral.+--+-- Specifically, for polynomials \(\left(R(t), S(t,x)\right)\) such that+--+-- \[+-- \frac{df}{dx} = \frac{d}{dx} \sum_{\alpha \in \{ t \mid R(t) = 0 \}} \left( \alpha \log \left( S(\alpha,x) \right) \right)+-- \]+--+-- then with return value \(\left( \left(P(u,v), Q(u,v)\right), \left(A(u,v,x), B(u,v,x)\right) \right)\),+-- and a return value \(g_{uv}\) from 'complexLogTermToAtan' for \(A(u,v)\) and \(B(u,v)\), the real function is+--+-- \[+-- \frac{df}{dx} = \frac{d}{dx} \left(+-- \sum_{(a,b) \in \{(u,v) \in (\mathbb{R}, \mathbb{R}) \mid P(u,v)=Q(u,v)=0, b > 0\}}+--   \left( a \log \left( A(a,b,x)^2 + B(a,b,x)^2 \right) + b g_{ab}(x) \right)+-- + \sum_{a \in \{t \in \mathbb{R} \mid R(t)=0 \}} \left( a \log (S(a,x)) \right)+-- \right)+-- \]+--+-- The return value are polynomials \(\left( (P,Q), (A,B) \right)\), where+--+-- * \(P\) is a \(u\)-polynomial, i.e., a polynomial with variable \(u\), with coefficients which are \(v\)-polynomials.+--+-- * \(Q\) is a \(u\)-polynomial, with coefficients which are \(v\)-polynomials.+--+-- * \(A\) is an \(x\)-polynomial, with coefficients which are \(u\)-polynomials, which in turn have coefficients with \(v\)-polynomials.+--+-- * \(B\) is an \(x\)-polynomial, with coefficients which are \(u\)-polynomials, which in turn have coefficients with \(v\)-polynomials.+--+-- For example,+--+-- >>> let r = 4 * power 2 + 1 :: IndexedPolynomial+-- >>> let s = power 3 + scale (2 * power 1) (power 2) - 3 * power 1 - scale (4 * power 1) 1 :: IndexedPolynomialWith IndexedPolynomial+-- >>> complexLogTermToRealTerm (r, s)+-- (([(0,(-4)x^2 + 1),(2,4)],[(1,8x)]),([(0,[(1,(-4))]),(1,[(0,(-3))]),(2,[(1,2)]),(3,[(0,1)])],[(0,[(0,(-4)x)]),(2,[(0,2x)])]))+--+-- While the return value may be hard to parse, this means:+--+-- \[+-- \begin{align*}+-- P & = 4u^2 - 4v^2 + 1 \\+-- Q & = 8uv \\+-- A & = x^3 + 2ux^2 - 3x - 4u \\+-- B & = 2vx^2 - 4v+-- \end{align*}+-- \]+complexLogTermToRealTerm ::+  (IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial) ->+  ( (IndexedPolynomialWith IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial),+    (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial), IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))+  )+complexLogTermToRealTerm (q, s) = ((qp, qq), (sp, sq))+  where+    -- For all of the following, i is the imaginary number.+    -- We use an i polynomial instead of Complex to represent complex numbers+    -- because the Complex a is not an instance of the Num class unless a is+    -- an instance of the RealFloat class.++    -- We use polynomial coefficients to introduce a separate variable.+    -- An alternative would have been to use Expression coefficients,+    -- but this would require a guarantee that we can rewrite an Expression+    -- down to the degree where we can tease apart the real and imaginary parts+    -- in a complex number.++    -- Compute q(u+iv) as an i polynomial with coefficients+    -- of u polynomials with coefficients+    -- of v polynomials with rational coefficients.+    q' = getSum $ foldTerms reduceImaginary $ getSum $ foldTerms fromTerm q+      where+        fromTerm :: Int -> Rational -> Sum (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial))+        fromTerm e c = Sum $ c' * (u + i * v) ^ e+          where+            c' = scale (scale (scale c 1) 1) 1+        i = power 1+        u = scale (power 1) 1+        v = scale (scale (power 1) 1) 1+    -- q' == qp + i * qq+    (qp, qq) = (coefficient q' 0, coefficient q' 1)++    -- Compute s(u+iv,x) as an i polynomial with coefficients+    -- of x polynomials with coefficients+    -- of u polynomials with coefficients+    -- of v polynomials with rational coefficients.+    s' = getSum $ foldTerms reduceImaginary $ getSum $ foldTerms fromTerm s+      where+        fromTerm :: Int -> IndexedPolynomial -> Sum (IndexedPolynomialWith (IndexedPolynomialWith (IndexedPolynomialWith IndexedPolynomial)))+        fromTerm e c = Sum $ c' * x ^ e+          where+            c' = getSum $ foldTerms fromCoefficient c+            fromCoefficient e' c'' = Sum $ c''' * (u + i * v) ^ e'+              where+                c''' = scale (scale (scale (scale c'' 1) 1) 1) 1+        i = power 1+        x = scale (power 1) 1+        u = scale (scale (power 1) 1) 1+        v = scale (scale (scale (power 1) 1) 1) 1+    -- s' = sp + i * sq+    (sp, sq) = (coefficient s' 0, coefficient s' 1)++    -- For terms in polynomials of i, reduce them to the form x or i*x.+    reduceImaginary :: (Eq a, Num a) => Int -> a -> Sum (IndexedPolynomialWith a)+    reduceImaginary e c = Sum $ case e `mod` 4 of+      0 -> c'+      1 -> c' * i+      2 -> c' * (-1)+      3 -> c' * (-i)+      _ -> 0 -- Not possible.+      where+        i = power 1+        c' = scale c 1++-- | For the ingredients of a complex logarithm, return an equivalent real function in terms of an indefinite integral.+--+-- Specifically, for polynomials \(\left(R(t), S(t,x)\right)\) such that+--+-- \[+-- \frac{df}{dx} = \frac{d}{dx} \sum_{\alpha \in \{ t \mid R(t) = 0 \}} \left( \alpha \log \left( S(\alpha,x) \right) \right)+-- \]+--+-- a symbolic representation for \(f\) will be returned.  See 'complexLogTermToRealTerm' for specifics as to how \(f\) is derived.+complexLogTermToRealExpression ::+  -- | Symbol for the variable.+  Text ->+  -- | Polynomials \(R(t)\) and \(S(t,x)\).+  (IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial) ->+  -- | Expression for the real function \(f\).+  Maybe Expression+complexLogTermToRealExpression v (r, s)+  | (Just xys) <- solveBivariatePolynomials p q,+    (Just h) <- f xys,+    (Just zs) <- toRationalList (solve r) =+      Just $ sum h + g zs+  | otherwise = Nothing+  where+    ((p, q), (a, b)) = complexLogTermToRealTerm (r, s)++    f :: [(Rational, Rational)] -> Maybe [Expression]+    f xys = toMaybeList $ do+      (x, y) <- filter ((> 0) . snd) xys+      let flatten'' = mapCoefficients (toExpr (fromRational y) fromRational) -- v-polynomials into Expressions.+      let flatten' = mapCoefficients (toExpr (fromRational x) id . flatten'') -- u-polynomials into Expressions.+      let flatten = toExpr (Symbol v) id . flatten' -- x-polynomials into Expressions.+      -- a and b flattened into Expressions.+      let a' = flatten a+      let b' = flatten b+      -- a and b flattened into x-polynomials with rational number coefficients.+      return $ do+        a'' <- convertCoefficients $ flatten' a+        b'' <- convertCoefficients $ flatten' b+        return $ fromRational x * log (a' * a' + b' * b') + fromRational y * complexLogTermToAtan v a'' b''++    g zs = sum $ do+      z <- zs+      let s' = mapCoefficients (toExpr (fromRational z) fromRational) s+      return $ fromRational z * Log' (toExpression v toSymbolicCoefficient s')++    toRationalList :: Maybe [Expression] -> Maybe [Rational]+    toRationalList Nothing = Nothing+    toRationalList (Just []) = Just []+    toRationalList (Just (x : xs))+      | (Just x'') <- convert (simplify x'), (Just xs'') <- xs' = Just $ x'' : xs''+      | otherwise = Nothing+      where+        x' = simplify x+        xs' = toRationalList $ Just xs++    -- Convert a simplified Expression into a rational number.+    convert (Number n) = Just $ fromIntegral n+    convert (Number n :/: Number m) = Just $ fromIntegral n / fromIntegral m+    convert _ = Nothing++    -- Convert polynomial with Expression coefficients into a polynomial with rational number coefficients.+    convertCoefficients :: IndexedPolynomialWith Expression -> Maybe IndexedPolynomial+    convertCoefficients x = sum . map (\(e, c) -> scale c (power e)) <$> toMaybeList (foldTerms (\e c -> [(e,) <$> convert (simplify c)]) x)++    -- Turns a polynomial into an Expression.+    -- Function h is used to turn the coefficient into an Expression.+    toExpr x h u = getSum $ foldTerms (\e'' c -> Sum $ h c * (x ** Number (fromIntegral e''))) u++-- | From the ingredients of a complex logarithm, return the expression for the complex algorithm.+-- Specifically, for polynomials \(\left(Q(t), S(t,x)\right)\),+-- a symbolic representation for the following will be returned.+--+-- \[+-- \sum_{\alpha \in \{ t \mid Q(t) = 0 \}} \left( \alpha \log \left( S(\alpha,x) \right) \right)+-- \]+complexLogTermToComplexExpression ::+  -- | Symbol for the variable.+  Text ->+  -- | Polynomials \(Q(t)\) and \(S(t,x)\).+  (IndexedPolynomial, IndexedPolynomialWith IndexedPolynomial) ->+  -- | Expression for the logarithm.+  Maybe Expression+complexLogTermToComplexExpression v (q, s) = do+  as <- complexSolve q+  let terms = do+        a <- as+        let s' = mapCoefficients (collapse a) s+        let s'' = toExpression v toSymbolicCoefficient s'+        return $ a * log s''+  return $ sum terms+  where+    -- Collapse a polynomial coefficient of a polynomial into an expression with the variable substituted.+    -- E.g., turn (t+2)x+1 into (3+2)x+1 for t=3.+    collapse a c' = getSum $ foldTerms (\e c -> Sum $ fromRational c * a ** fromIntegral e) c'++-- | Returns the roots for two variables in two polynomials.+--+-- Only supports rational roots.  If not all real roots are rational, then it will return 'Nothing'.+-- Returning all real roots would be preferable, but this is not supported at this time.+--+-- If the function cannot derive the roots otherwise, either, 'Nothing' will be returned as well.+solveBivariatePolynomials ::+  IndexedPolynomialWith IndexedPolynomial ->+  IndexedPolynomialWith IndexedPolynomial ->+  Maybe [(Rational, Rational)]+solveBivariatePolynomials p q = do+  let p' = toRationalFunctionCoefficients p+  let q' = toRationalFunctionCoefficients q+  resultant <- toPoly $ fst $ subresultant p' q'+  vs' <- solve resultant+  vs <- toMaybeList $ map (convert . simplify) vs'+  concat <$> toMaybeList (map solveForU vs)+  where+    toRationalFunctionCoefficients = mapCoefficients (`toRationalFunction` 1)++    -- For each v, returns list of (u,v) such that P(u,v)=Q(u,v)=0.+    solveForU :: Rational -> Maybe [(Rational, Rational)]+    solveForU v+      | 0 <- p' = do+          -- Any u will make p'=0 true, so we only need to solve p'.+          u <- map (convert . simplify) <$> solve q'+          map (,v) <$> toMaybeList u+      | 0 <- q' = do+          -- Any u will make q'=0 true, so we only need to solve p'.+          u <- map (convert . simplify) <$> solve p'+          map (,v) <$> toMaybeList u+      | otherwise = do+          up <- map (convert . simplify) <$> solve p'+          uq <- map (convert . simplify) <$> solve q'+          up' <- toMaybeList up+          uq' <- toMaybeList uq+          return $ map (,v) $ up' `intersect` uq'+      where+        p' = mapCoefficients (getSum . foldTerms (\e c -> Sum $ c * v ^ e)) p+        q' = mapCoefficients (getSum . foldTerms (\e c -> Sum $ c * v ^ e)) q++    -- Turn a simplified Expression into a rational number if possible.+    convert :: Expression -> Maybe Rational+    convert (Number n) = Just $ fromIntegral n+    convert (Number n :/: Number m) = Just $ fromIntegral n / fromIntegral m+    convert _ = Nothing++-- | Turn the rational function into a polynomial if possible.+toPoly :: RationalFunction -> Maybe IndexedPolynomial+toPoly (RationalFunction p q)+  | degree q == 0, q /= 0 = Just p'+  | otherwise = Nothing+  where+    p' = scale (1 / leadingCoefficient q) p++-- | If there are any nothings, then turn the list into nothing.+-- Otherwise, turn it into the list of just the elements.+toMaybeList :: [Maybe a] -> Maybe [a]+toMaybeList [] = Just []+toMaybeList (Nothing : _) = Nothing+toMaybeList (Just x : xs)+  | (Just xs') <- toMaybeList xs = Just (x : xs')+  | otherwise = Nothing
+ src/Symtegration/Integration/Substitution.hs view
@@ -0,0 +1,76 @@+-- |+-- Module: Symtegration.Integration.Substitution+-- Description: Integration by substitution.+-- Copyright: Copyright 2025 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Integration.Substitution (integrate) where++import Data.Foldable (asum)+import Data.Text (Text)+import Symtegration.Differentiation+import Symtegration.Integration.Factor+import Symtegration.Symbolic++-- $setup+-- >>> import Symtegration.Symbolic.Haskell+-- >>> import Symtegration.Symbolic.Simplify++-- | Integrates by substitution.+--+-- Specifically, if for+--+-- \[ \int f(g(x)) h(x) \, dx\]+--+-- it is the case that \(\frac{dg(x)}{dx} = h(x)\), then compute \(\int f(v) \, dv\) and substitute with \(v=g(x)\).+--+-- >>> import Symtegration.Integration.Trigonometric qualified as Trigonometric+-- >>> toHaskell <$> simplify <$> integrate [Trigonometric.integrate] "x" (sin ("a" * "x" + 1))+-- Just "(-1) * 1 / a * cos (1 + a * x)"+integrate ::+  -- | Integration algorithms to try after substitution.+  [Text -> Expression -> Maybe Expression] ->+  -- | Symbol for the variable.+  Text ->+  -- | Expression to integrate.+  Expression ->+  -- | Integral, if derived.+  Maybe Expression+integrate fs v (x :*: UnaryApply func y)+  | Number 0 <- d = Nothing -- Argument is constant.+  | x' == y',+    -- Re-use v as the variable, as it is the one symbol guaranteed not to appear outside the argument.+    Just e <- integrateSubstitution fs v (UnaryApply func (Symbol v)) =+      Just $ (c :/: d) :*: substitute e (\s -> if s == v then Just y else Nothing)+  | otherwise = Nothing+  where+    (c, x') = factor v x+    (d, y') = factor v $ differentiate v y+integrate fs v (e@(UnaryApply _ _) :*: x) = integrate fs v $ x :*: e+integrate fs v e@(UnaryApply _ _) = integrate fs v $ Number 1 :*: e+integrate fs v (x :*: BinaryApply func y z)+  -- Re-use v as the variable, as it is the one symbol guaranteed not to appear outside the argument.+  | c /= Number 0,+    x' == y',+    isConstant v z,+    Just e <- integrateSubstitution fs v (BinaryApply func (Symbol v) z) =+      Just $ (b :/: c) :*: substitute e (\s -> if s == v then Just y else Nothing)+  | d /= Number 0,+    x' == z',+    isConstant v y,+    Just e <- integrateSubstitution fs v (BinaryApply func y (Symbol v)) =+      Just $ (b :/: d) :*: substitute e (\s -> if s == v then Just z else Nothing)+  | otherwise = Nothing+  where+    (b, x') = factor v x+    (c, y') = factor v $ differentiate v y+    (d, z') = factor v $ differentiate v z+integrate fs v (e@(BinaryApply _ _ _) :*: x) = integrate fs v $ x :*: e+integrate fs v e@(BinaryApply func _ _)+  | func /= Multiply = integrate fs v $ Number 1 :*: e+  | otherwise = Nothing+integrate _ _ _ = Nothing++-- | Use the given functions to integrate the given expression.+integrateSubstitution :: [Text -> Expression -> Maybe Expression] -> Text -> Expression -> Maybe Expression+integrateSubstitution fs v e = asum $ map (\f -> f v e) fs
+ src/Symtegration/Integration/Sum.hs view
@@ -0,0 +1,50 @@+-- |+-- Module: Symtegration.Integration.Sum+-- Description: Integrates the sum of multiple terms in an expression.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Integration.Sum (integrate) where++import Data.Foldable (asum)+import Data.Text (Text)+import Symtegration.Symbolic++-- $setup+-- >>> import Symtegration.Symbolic.Haskell+-- >>> import Symtegration.Symbolic.Simplify++-- | Integrate term by term and returns the sum, using direct methods on each term.+--+-- >>> import Symtegration.Integration.Powers qualified as P+-- >>> import Symtegration.Integration.Trigonometric qualified as T+-- >>> let f = "x" + sin "x"+-- >>> P.integrate "x" f+-- Nothing+-- >>> T.integrate "x" f+-- Nothing+-- >>> let g = integrate [P.integrate, T.integrate] "x" f+-- >>> toHaskell . simplify <$> g+-- Just "(-1) * cos x + 1 / 2 * x ** 2"+integrate ::+  -- | Functions for directly integrating each term.+  [Text -> Expression -> Maybe Expression] ->+  -- | The variable being integrated over.+  Text ->+  -- | The expression being integrated.+  Expression ->+  -- | The integral, if successful.+  Maybe Expression+integrate fs v (Negate' x) =+  UnaryApply Negate <$> integrate fs v x+integrate fs v (x :-: y) =+  integrate fs v (x :+: Negate' y)+integrate fs v (x@(_ :+: _) :+: y@(_ :+: _)) =+  BinaryApply Add <$> integrate fs v x <*> integrate fs v y+integrate fs v (x@(_ :+: _) :+: y) =+  BinaryApply Add <$> integrate fs v x <*> asum [f v y | f <- fs]+integrate fs v (x :+: y@(_ :+: _)) =+  BinaryApply Add <$> asum [f v x | f <- fs] <*> integrate fs v y+integrate fs v (x :+: y) =+  BinaryApply Add <$> asum [f v x | f <- fs] <*> asum [f v y | f <- fs]+integrate _ _ _ = Nothing
+ src/Symtegration/Integration/Term.hs view
@@ -0,0 +1,43 @@+-- |+-- Module: Symtegration.Integration.Term+-- Description: Integrates a single term.+-- Copyright: Copyright 2025 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Integration.Term (integrate) where++import Data.Foldable (asum)+import Data.Text (Text)+import Symtegration.Integration.Factor+import Symtegration.Symbolic+import Symtegration.Symbolic.Simplify++-- $setup+-- >>> import Symtegration.Symbolic.Haskell+-- >>> import Symtegration.Symbolic.Simplify++-- | Integrate a single term, separating out the constant factor and+-- applying direct integration methods to the non-constant factor.+--+-- >>> import Symtegration.Integration.Trigonometric qualified as T+-- >>> let f = "a" * sin "x"+-- >>> T.integrate "x" f+-- Nothing+-- >>> let g = integrate [T.integrate] "x" f+-- >>> toHaskell . simplify <$> g+-- Just "(-1) * a * cos x"+--+-- Assumes the expression has had algebraic ring ordering applied.+integrate ::+  -- | Functions for directly integrating the non-constant factor.+  [Text -> Expression -> Maybe Expression] ->+  -- | The variable being integrated over.+  Text ->+  -- | The expression being integrated.+  Expression ->+  -- | The integral, if successful.+  Maybe Expression+integrate fs v e = asum $ map (\f -> (:*:) c <$> f v u) fs+  where+    e' = simplifyForVariable v e+    (c, u) = factor v e'
+ src/Symtegration/Integration/Trigonometric.hs view
@@ -0,0 +1,64 @@+-- |+-- Module: Symtegration.Integration.Trigonometric+-- Description: Basic integration of trigonometric functions.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+--+-- Supports basic integration of trigonometric functions.+-- This does not support the integration of anything else,+-- even if it is trivial like integrating a constant.+module Symtegration.Integration.Trigonometric (integrate) where++import Data.Text (Text)+import Symtegration.Symbolic++-- $setup+-- >>> import Symtegration.Symbolic.Haskell++-- | Integrates trigonometric functions required by the 'Floating' type class.+--+-- >>> toHaskell <$> integrate "x" (sin "x")+-- Just "negate (cos x)"+-- >>> toHaskell <$> integrate "x" (cos "x")+-- Just "sin x"+integrate :: Text -> Expression -> Maybe Expression+integrate _ (Number _) = Nothing+integrate _ (Symbol _) = Nothing+integrate v (Sin' x@(Symbol s))+  | s == v = Just $ Negate' $ Cos' x+  | otherwise = Nothing+integrate v (Cos' x@(Symbol s))+  | s == v = Just $ Sin' x+  | otherwise = Nothing+integrate v (Tan' x@(Symbol s))+  | s == v = Just $ Negate' $ Log' $ Abs' $ Cos' x+  | otherwise = Nothing+integrate v (Asin' x@(Symbol s))+  | s == v = Just $ (x :*: Asin' x) :+: Sqrt' (1 :-: (x :**: 2))+  | otherwise = Nothing+integrate v (Acos' x@(Symbol s))+  | s == v = Just $ (x :*: Acos' x) :-: Sqrt' (1 :-: (x :**: 2))+  | otherwise = Nothing+integrate v (Atan' x@(Symbol s))+  | s == v = Just $ (x :*: Atan' x) :-: (Log' ((x :**: 2) :+: 1) :/: 2)+  | otherwise = Nothing+integrate v (Sinh' x@(Symbol s))+  | s == v = Just $ Cosh' x+  | otherwise = Nothing+integrate v (Cosh' x@(Symbol s))+  | s == v = Just $ Sinh' x+  | otherwise = Nothing+integrate v (Tanh' x@(Symbol s))+  | s == v = Just $ Log' $ Cosh' x+  | otherwise = Nothing+integrate v (Asinh' x@(Symbol s))+  | s == v = Just $ (x :*: Asinh' x) :-: Sqrt' ((x :**: 2) + 1)+  | otherwise = Nothing+integrate v (Acosh' x@(Symbol s))+  | s == v = Just $ (x :*: Acosh' x) :-: (Sqrt' (x :+: 1) :*: Sqrt' (x :-: 1))+  | otherwise = Nothing+integrate v (Atanh' x@(Symbol s))+  | s == v = Just $ (x :*: Atanh' x) :+: (Log' (1 :-: (x :**: 2)) :/: 2)+  | otherwise = Nothing+integrate _ _ = Nothing
+ src/Symtegration/Numeric.hs view
@@ -0,0 +1,45 @@+-- |+-- Module: Symtegration.Numeric+-- Description: Numerical algorithms that are useful for implementing symbolic integration.+-- Copyright: Copyright 2025 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+--+-- This module contains numerical algorithms that are useful to more than one module,+-- ultimately for the purpose of symbolic integration of mathematical expressions.+-- By numerical algorithms here, we means algorithm that work on pure numbers and not symbols.+-- The algorithms should still return exact results.+module Symtegration.Numeric (root) where++-- | Compute the integer root to the given power.+-- I.e., find \(m\) such that \(m^k = n\).+--+-- >>> root 27 3+-- Just 3+-- >>> root (-27) 3+-- Just (-3)+-- >>> root 2 2+-- Nothing+root ::+  -- | Number \(n\) whose root we want.+  Integer ->+  -- | The power \(k\) of the root.+  Integer ->+  -- | The root \(m\).+  Maybe Integer+root 0 _ = Just 0+root 1 _ = Just 1+root n k+  | k < 0 = Nothing+  | GT <- compare n 0 = search n 1 n+  | LT <- compare n 0, odd k = (* (-1)) <$> search (-n) 1 (-n)+  | otherwise = Nothing+  where+    search m low hi+      | low >= hi, c /= EQ = Nothing+      | EQ <- c = Just mid+      | GT <- c = search m low (mid - 1)+      | LT <- c = search m (mid + 1) hi+      where+        mid = (low + hi) `div` 2+        c = compare (mid ^ k) m
+ src/Symtegration/Polynomial.hs view
@@ -0,0 +1,455 @@+-- |+-- Module: Symtegration.Polynomial+-- Description: Polynomials for Symtegration.+-- Copyright: Copyright 2025 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+--+-- This modules defines a type class that concrete types representing polynomials+-- should be an instance of.  It includes important algorithms operating on+-- polynomials.  In particular, algorithms for polynomial division and+-- the extended Euclidean algorithm are included.+module Symtegration.Polynomial+  ( -- * Polynomials+    Polynomial (..),+    monic,+    mapCoefficients,+    mapCoefficientsM,++    -- * Algorithms+    divide,+    pseudoDivide,+    extendedEuclidean,+    diophantineEuclidean,+    greatestCommonDivisor,+    subresultant,+    differentiate,+    integrate,+    squarefree,+  )+where++import Data.Monoid (Sum (..))++-- $setup+-- >>> import Data.Ratio ((%), denominator, numerator)+-- >>> import Symtegration.Symbolic+-- >>> import Symtegration.Symbolic.Simplify+-- >>> import Symtegration.Polynomial.Indexed++-- | Polynomials must support the operations specified in this type class.+-- All powers must be non-negative.+class (Integral e, Num c) => Polynomial p e c where+  -- | Returns the degree of a given polynomial.+  --+  -- The following returns 9 for the highest term in \(3x^9 + 2x^4 + x\):+  --+  -- >>> degree (3 * power 9 + 2 * power 4 + power 1 :: IndexedPolynomial)+  -- 9+  degree :: p e c -> e++  -- | Returns the coefficient for the term with the given power.+  --+  -- The following returns 4 from the \(4x^3\) term in \(x^4 + 4x^3 + 3\):+  --+  -- >>> coefficient (power 4 + 4 * power 3 + 3 :: IndexedPolynomial) 3+  -- 4 % 1+  coefficient :: p e c -> e -> c++  -- | Returns the leading coefficient.+  --+  -- The following returns 6 from the \(6x^3\) term in \(6x^3 + 2x^2\):+  --+  -- >>> leadingCoefficient (6 * power 3 + 2 * power 2 :: IndexedPolynomial)+  -- 6 % 1+  --+  -- The leading coefficient is never zero unless the polynomial itself is zero.+  leadingCoefficient :: p e c -> c++  -- | Returns the polynomial without the leading term.+  --+  -- >>> deleteLeadingTerm (2 * power 3 + power 1 + 2 :: IndexedPolynomial)+  -- x + 2+  deleteLeadingTerm :: p e c -> p e c++  -- | Fold the terms, i.e., the powers and coefficients, using the given monoid.+  -- Only terms with non-zero coefficients will be folded.+  -- Folding is ordered from lower to higher terms.+  --+  -- For example with \(3x^5 - 2x + 7\),+  --+  -- >>> foldTerms (\e c -> show (e, c)) (3 * power 5 - 2 * power 1 + 7 :: IndexedPolynomial)+  -- "(0,7 % 1)(1,(-2) % 1)(5,3 % 1)"+  foldTerms :: (Monoid m) => (e -> c -> m) -> p e c -> m++  -- | Multiplies a polynomial by a scalar.+  --+  -- The following divides \(6x + 2\) by 2:+  --+  -- >>> scale (1 % 2) (6 * power 1 + 2 :: IndexedPolynomial)+  -- 3x + 1+  scale :: c -> p e c -> p e c++  -- | Returns a single term with the variable raised to the given power.+  --+  -- The following is equivalent to \(x^5\):+  --+  -- >>> power 5 :: IndexedPolynomial+  -- x^5+  power :: e -> p e c++-- | Scale the polynomial so that its leading coefficient is one.+--+-- >>> monic $ 4 * power 2 + 4 * power 1 + 4 :: IndexedPolynomial+-- x^2 + x + 1+--+-- The exception is when the polynomial is zero.+--+-- >>> monic 0 :: IndexedPolynomial+-- 0+monic :: (Polynomial p e c, Eq c, Fractional c) => p e c -> p e c+monic p+  | leadingCoefficient p == 0 = p+  | otherwise = scale (1 / leadingCoefficient p) p++-- | Maps the coefficients in a polynomial to form another polynomial.+--+-- For example, it can be used to convert a polynomial with 'Rational' coefficients+-- into a polynomial with 'Expression' coefficients.+--+-- >>> let p = 2 * power 1 + 1 :: IndexedPolynomial+-- >>> let q = mapCoefficients fromRational p :: IndexedSymbolicPolynomial+-- >>> simplify $ coefficient q 1+-- Number 2+--+-- Note that only non-zero coefficients are mapped.+mapCoefficients ::+  (Polynomial p e c, Polynomial p e c', Num (p e c), Num (p e c')) =>+  (c -> c') ->+  p e c ->+  p e c'+mapCoefficients f p = getSum $ foldTerms convertTerm p+  where+    convertTerm e c = Sum $ scale (f c) (power e)++-- | Maps the coefficients in a polynomial to form another polynomial, but in a monad.+-- Specifically, it maps each coefficient in a monadic action,+-- and collects the products of each result and power.+--+-- For example, with the 'Maybe' monad:+--+-- >>> let f q | denominator q == 1 = Just q | otherwise = Nothing+-- >>> let p = scale 2 (power 2) + scale 3 (power 1) :: IndexedPolynomial+-- >>> mapCoefficientsM f p+-- Just 2x^2 + 3x+-- >>> let q = scale (1/2) (power 2) + scale 3 (power 1) :: IndexedPolynomial+-- >>> mapCoefficientsM f q+-- Nothing+--+-- As an another example, with the 'Either' monad:+--+-- >>> let f q | denominator q == 1 = Right q | otherwise = Left "not integer"+-- >>> let p = scale 2 (power 2) + scale 3 (power 1) :: IndexedPolynomial+-- >>> mapCoefficientsM f p+-- Right 2x^2 + 3x+-- >>> let q = scale (1/2) (power 2) + scale 3 (power 1) :: IndexedPolynomial+-- >>> mapCoefficientsM f q+-- Left "not integer"+--+-- Note that only non-zero coefficients are mapped.+mapCoefficientsM ::+  (Polynomial p e c, Polynomial p e c', Num (p e c), Num (p e c'), Monad m) =>+  (c -> m c') ->+  p e c ->+  m (p e c')+mapCoefficientsM f p = sum <$> mapM f' terms+  where+    terms = foldTerms (\e c -> [(e, c)]) p+    f' (e, c) = do+      c' <- f c+      return $ scale c' $ power e++-- | Polynomial division.  It returns the quotient polynomial and the remainder polynomial.+--+-- For example, dividing \(p = x^3-12x^2-42\) by \(q = x^2 - 2x + 1\)+-- returns \(x-10\) as the quotient and \(-21x-32\) as the remainder,+-- since \(p = (x-10)q -21x - 32\):+--+-- >>> let p = power 3 - 12 * power 2 - 42 :: IndexedPolynomial+-- >>> let q = power 2 - 2 * power 1 + 1 :: IndexedPolynomial+-- >>> divide p q+-- (x + (-10),(-21)x + (-32))+divide ::+  (Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>+  -- | Dividend polynomial being divided.+  p e c ->+  -- | Divisor polynomial dividing the dividend.+  p e c ->+  -- | Quotient and remainder.+  (p e c, p e c)+divide p q = go 0 p+  where+    go quotient remainder+      | remainder /= 0, delta >= 0 = go (quotient + t) (remainder' - qt')+      | otherwise = (quotient, remainder)+      where+        delta = degree remainder - degree q+        t = scale (leadingCoefficient remainder / leadingCoefficient q) $ power delta+        -- remainder and q * t will have the same leading coefficients.+        -- Subtract them without the leading terms.+        -- Not necessary for purely numeric coefficients,+        -- but guarantees the cancellation of the leading terms when coefficients are symbolic.+        remainder' = deleteLeadingTerm remainder+        qt' = deleteLeadingTerm $ q * t++-- | Polynomial pseudo-division.  It returns the pseudo-quotient and pseudo-remainder polynomials.+--+-- Equivalent to \(b^{\delta+1} p\) divided by \(q\),+-- where \(p\) and \(q\) are polynomials with integer coefficients,+-- \(b\) is the leading coefficient of \(q\) and \(\delta=\max(-1, \deg(p) - \deg(q))\).+-- This guarantees the pseudo-quotient and pseudo-remainder exist,+-- even when the quotient and remainder do not when only integer coefficients are allowed.+--+-- For example, with \(p = 3x^3 + x^2 + x + 5\) and \(q = 5x^2 - 3x + 1\),+-- it is the case that \(5^2p = (15x + 14)q + (52x + 111)\):+--+-- >>> let p = 3 * power 3 + power 2 + power 1 + 5 :: IndexedPolynomial+-- >>> let q = 5 * power 2 - 3 * power 1 + 1 :: IndexedPolynomial+-- >>> pseudoDivide p q+-- (15x + 14,52x + 111)+pseudoDivide ::+  (Polynomial p e c, Eq (p e c), Num (p e c), Num c) =>+  -- | Dividend polynomial being pseudo-divided.+  p e c ->+  -- | Divisor polynomial pseudo-dividing the dividend.+  p e c ->+  -- | Pseudo-quotient and pseudo-remainder.+  (p e c, p e c)+pseudoDivide p q+  | degree p < degree q = (0, p)+  | otherwise = go (1 + degree p - degree q) 0 p+  where+    b = leadingCoefficient q+    go n quotient remainder+      | remainder /= 0, delta >= 0 = go (n - 1) quotient' remainder'+      | otherwise = (scale (b ^ n) quotient, scale (b ^ n) remainder)+      where+        delta = degree remainder - degree q+        t = scale (leadingCoefficient remainder) (power delta)+        quotient' = scale b quotient + t+        -- Subtract with the leading terms deleted.+        -- The leading terms cancel out numerically,+        -- but guarantee cancellation when the coefficients are symbolic.+        remainder' = deleteLeadingTerm (scale b remainder) - deleteLeadingTerm (t * q)++-- | The extended Euclidean algorithm.  For polynomials \(p\) and \(q\),+-- it returns the greatest common divisor between \(p\) and \(q\).+-- It also returns \(s\) and \(t\) such that \(sp+tq = \gcd(p,q)\).+--+-- For example, for \(p=2x^5-2x\) and \(q=x^4-2x^2+1\), it is the case+-- that \(\gcd(p,q)=-x^2+1\) and \((-\frac{1}{4}x) p + (\frac{1}{2}x^2 + 1) q = -x^2+1\):+--+-- >>> let p = 2 * power 5 - 2 * power 1 :: IndexedPolynomial+-- >>> let q = power 4 - 2 * power 2 + 1 :: IndexedPolynomial+-- >>> extendedEuclidean p q+-- (((-1) % 4)x,(1 % 2)x^2 + 1,(-1)x^2 + 1)+extendedEuclidean ::+  (Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>+  -- | Polynomial \(p\).+  p e c ->+  -- | Polynomial \(q\).+  p e c ->+  -- | \(s\), \(t\), and \(\gcd(p,q)\).+  (p e c, p e c, p e c)+extendedEuclidean u v = descend u v 1 0 0 1+  where+    descend g 0 s t _ _ = (s, t, g)+    descend a b a1 a2 b1 b2 = descend b r b1 b2 r1 r2+      where+        (q, r) = divide a b+        r1 = a1 - q * b1+        r2 = a2 - q * b2++-- | Solves \(sa + tb = c\) for given polynomials \(a\), \(b\), and \(c\).+-- It will be the case that either \(s=0\) or+-- the degree of \(s\) will be less than the degree of \(b\).+--+-- >>> let a = power 4 - 2 * power 3 - 6 * power 2 + 12 * power 1 + 15 :: IndexedPolynomial+-- >>> let b = power 3 + power 2 - 4 * power 1 - 4 :: IndexedPolynomial+-- >>> let c = power 2 - 1 :: IndexedPolynomial+-- >>> diophantineEuclidean a b c+-- Just (((-1) % 5)x^2 + (4 % 5)x + ((-3) % 5),(1 % 5)x^3 + ((-7) % 5)x^2 + (16 % 5)x + (-2))+--+-- If there is no such \((s,t)\), then 'Nothing' is returned.+diophantineEuclidean ::+  (Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>+  -- | Polynomial \(a\).+  p e c ->+  -- | Polynomial \(b\).+  p e c ->+  -- | Polynomial \(c\).+  p e c ->+  -- | \((s,t)\) such that \(sa + tb = c\).+  Maybe (p e c, p e c)+diophantineEuclidean a b c+  | r /= 0 = Nothing+  | s' /= 0, degree s' >= degree b = Just (r', t' + q' * a)+  | otherwise = Just (s', t')+  where+    (s, t, g) = extendedEuclidean a b+    (q, r) = divide c g+    s' = q * s+    t' = q * t+    (q', r') = divide s' b++-- | Returns the greatest common divisor btween two polynomials.+--+-- Convenient wrapper over 'extendedEuclidean' which only returns the greatest common divisor.+greatestCommonDivisor ::+  (Polynomial p e c, Eq (p e c), Num (p e c), Fractional c) =>+  -- | Polynomial \(p\).+  p e c ->+  -- | Polynomial \(q\).+  p e c ->+  -- | \(\gcd(p,q)\).+  p e c+greatestCommonDivisor p q = g+  where+    (_, _, g) = extendedEuclidean p q++-- | Returns the resultant and the subresultant polynomial remainder sequence for the given polynomials.+--+-- >>> subresultant (power 2 + 1) (power 2 - 1 :: IndexedPolynomial)+-- (4 % 1,[x^2 + 1,x^2 + (-1),(-2),0])+-- >>> subresultant (2 * power 2 - 3 * power 1 + 1) (5 * power 2 + power 1 - 6 :: IndexedPolynomial)+-- (0 % 1,[2x^2 + (-3)x + 1,5x^2 + x + (-6),17x + (-17),0])+-- >>> subresultant (power 3 + 2 * power 2 + 3 * power 1 + 4) (5 * power 2 + 6 * power 1 + 7 :: IndexedPolynomial)+-- (832 % 1,[x^3 + 2x^2 + 3x + 4,5x^2 + 6x + 7,16x + 72,832,0])+--+-- === __Reference__+--+-- See sections 1.4 and 1.5 in+-- [/Symbolic Integration I: Transcendental Functions/](https://doi.org/10.1007/b138171)+-- by Manuel Bronstein for the definition of resultants, subresultants,+-- polynomial remainder sequences, and subresultant polynomial remainder sequences.+subresultant ::+  (Polynomial p e c, Eq (p e c), Num (p e c), Num e, Fractional c) =>+  -- | First element in the remainder sequence.+  p e c ->+  -- | Second element in the remainder sequence.+  p e c ->+  -- | The resultant and the subresultant polynomial remainder sequence.+  (c, [p e c])+subresultant p q+  | degree p >= degree q = (resultantFromSequence rs betas, rs)+  | otherwise = ((-1) ^ (degree q * degree p) * resultant, prs)+  where+    (rs, betas) = subresultantRemainderSequence (p, q) gamma beta+    gamma = -1+    beta = (-1) ^ (1 + delta)+    delta = degree p - degree q++    (resultant, prs) = subresultant q p++-- | Derives the subresultant polynomial remainder sequence for 'subresultant'.+-- Constructs \(\gamma_i\), \(\beta_i\), and the remainder sequence as it goes along.+-- Returns the remainder sequence and the sequence of \(\beta_i\).+subresultantRemainderSequence ::+  (Polynomial p e c, Eq (p e c), Num (p e c), Num e, Fractional c) =>+  -- | The previous and current remainders in the sequence.+  (p e c, p e c) ->+  -- | \(\gamma_i\) as defined for the subresultant PRS.+  c ->+  -- | \(\beta_i\) as defined for the subresultant PRS.+  c ->+  -- | Polynomial remainder sequence and sequence of \(\beta_i\).+  ([p e c], [c])+subresultantRemainderSequence (rprev, rcurr) gamma beta+  | rcurr /= 0 = (rprev : rs, beta : betas)+  | otherwise = ([rprev, rcurr], [beta])+  where+    (rs, betas) = subresultantRemainderSequence (rcurr, rnext) gamma' beta'+    (_, r) = pseudoDivide rprev rcurr+    rnext = scale (1 / beta) r+    lc = leadingCoefficient rcurr+    delta = degree rprev - degree rcurr+    delta' = degree rcurr - degree rnext+    gamma' = ((-lc) ^ delta) * (gamma ^^ (1 - delta))+    beta' = (-lc) * (gamma' ^ delta')++-- | Constructs the resultant based on the subresultant polynomial remainder sequence+-- and the sequence of \(\beta_i\) used to construct the subresultant PRS.+resultantFromSequence ::+  (Polynomial p e c, Eq (p e c), Num (p e c), Num e, Fractional c) =>+  -- | Subresultant polynomial remainder sequence.+  [p e c] ->+  -- | Sequence of \(\beta_i\) used for deriving the subresultant PRS.+  [c] ->+  -- | Resultant.+  c+resultantFromSequence rs betas = go rs betas 1 1+  where+    go (r : r' : r'' : rs') (beta : betas') c s+      | [] <- rs', degree r' > 0 = 0+      | [] <- rs', degree r == 1 = leadingCoefficient r'+      | [] <- rs' = s * c * leadingCoefficient r' ^ degree r+      | otherwise = go (r' : r'' : rs') betas' c' s'+      where+        s' | odd (degree r), odd (degree r') = -s | otherwise = s+        c' = c * ((beta / (lc ^ (1 + delta))) ^ degree r') * (lc ^ (degree r - degree r''))+        lc = leadingCoefficient r'+        delta = degree r - degree r'+    go _ _ _ _ = 0++-- | Returns the derivative of the given polynomial.+--+-- >>> differentiate (power 2 + power 1 :: IndexedPolynomial)+-- 2x + 1+differentiate :: (Polynomial p e c, Num (p e c), Num c) => p e c -> p e c+differentiate p = getSum $ foldTerms diffTerm p+  where+    diffTerm 0 _ = Sum 0+    diffTerm e c = Sum $ scale (fromIntegral e * c) $ power (e - 1)++-- | Returns the integral of the given polynomial.+--+-- >>> integrate (power 2 + power 1 :: IndexedPolynomial)+-- (1 % 3)x^3 + (1 % 2)x^2+integrate :: (Polynomial p e c, Num (p e c), Fractional c) => p e c -> p e c+integrate p = getSum $ foldTerms integrateTerm p+  where+    integrateTerm e c = Sum $ scale (c / (1 + fromIntegral e)) $ power (e + 1)++-- | Returns the squarefree factorization of the given polynomial.+--+-- Specifically, for a polynomial \(p\), find \([p_1, p_2, \ldots, p_n]\) such that+--+-- \[ p = \sum_{k=1}^n p_k^k \]+--+-- where all \(p_k\) are squarefree, i.e., there is no polynomial \(q\) such that \(q^2 = p_k\).+--+-- For example, the squarefree factorization of \(x^8 + 6x^6 + 12x^4 + 8x^2\)+-- is \(x^2 (x^2 + 2)^3\):+--+-- >>> squarefree (power 8 + 6 * power 6 + 12 * power 4 + 8 * power 2 :: IndexedPolynomial)+-- [1,x,x^2 + 2]+squarefree :: (Polynomial p e c, Eq (p e c), Num (p e c), Eq c, Fractional c) => p e c -> [p e c]+squarefree 0 = [0]+squarefree p+  | (x : xs) <- factor u v = scale c x : xs+  | otherwise = [scale c 1]+  where+    c = leadingCoefficient p+    q = scale (1 / c) p+    q' = differentiate q+    g = monic $ greatestCommonDivisor q q'+    (u, _) = q `divide` g+    (v, _) = q' `divide` g+    factor s y+      | z == 0 = [s]+      | otherwise = f : factor s' y'+      where+        z = y - differentiate s+        f = monic $ greatestCommonDivisor s z+        (s', _) = s `divide` f+        (y', _) = z `divide` f
+ src/Symtegration/Polynomial/Indexed.hs view
@@ -0,0 +1,111 @@+-- |+-- Module: Symtegration.Polynomial.Indexed+-- Description: A polynomial representation mapping the power of each term to its coefficient.+-- Copyright: Copyright 2025 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Polynomial.Indexed+  ( IndexedPolynomial,+    IndexedSymbolicPolynomial,+    IndexedPolynomialWith,+  )+where++import Data.IntMap (IntMap)+import Data.IntMap qualified as IntMap+import Data.List (intersperse)+import Data.Maybe (fromMaybe)+import Data.Ratio (denominator, numerator)+import Data.Text (unpack)+import Symtegration.Polynomial+import Symtegration.Symbolic+import TextShow++-- | Polynomial representation which maps the power of each term to its coefficient.+-- Exponents are represented with 'Int', while coefficients are represented with 'Rational'.+-- It is an instance of the 'Polynomial' type class.+type IndexedPolynomial = IndexedPolynomialWith Rational++-- | Polynomial representation which maps the power of each term to its coefficient.+-- Exponents are represented with 'Int', while coefficients are represented with 'Expression'.+-- It is an instance of the 'Polynomial' type class.+type IndexedSymbolicPolynomial = IndexedPolynomialWith Expression++-- | Polynomial representation which maps the power of each term to its coefficient.+-- Exponents are represented with 'Int'.  Coefficients have a type as specified by the type parameter.+-- These types are an instance of the 'Polynomial' type class.+type IndexedPolynomialWith a = P Int a++-- | Type with two type parameters so that it can be an instance of 'Polynomial'.+-- The first type parameter is not involved in the data constructor;+-- it is used to set the exponent type for 'Polynomial'.+newtype P a b = P (IntMap b) deriving (Eq)++instance Show (P Int Rational) where+  show = unpack . showt++instance TextShow (P Int Rational) where+  showb (P m)+    | IntMap.null m = "0"+    | otherwise =+        mconcat $+          intersperse " + " $+            map showTerm $+              IntMap.toDescList m+    where+      showTerm (0, 1) = "1"+      showTerm (0, c) = showCoefficient c+      showTerm (1, c) = showCoefficient c <> "x"+      showTerm (e, 1) = "x^" <> showb e+      showTerm (e, c) = showCoefficient c <> "x^" <> showb e+      showCoefficient r+        | 1 <- r = mempty+        | 1 <- denominator r, r > 0 = showb $ numerator r+        | 1 <- denominator r, r < 0 = showbParen True $ showb $ numerator r+        | otherwise = showbParen True $ showb r++instance (Polynomial p e c, TextShow (p e c)) => Show (IndexedPolynomialWith (p e c)) where+  show = unpack . showt++instance (Polynomial p e c, TextShow (p e c)) => TextShow (IndexedPolynomialWith (p e c)) where+  showb (P m)+    | IntMap.null m = "0"+    | otherwise = showb $ IntMap.toList m++instance Show (P Int Expression) where+  show = unpack . showt++instance TextShow (P Int Expression) where+  showb (P m)+    | IntMap.null m = "0"+    | otherwise = showb $ IntMap.toList m++instance (Eq a, Num a) => Num (P Int a) where+  (P p) + (P q) = P $ filterNonzero $ IntMap.unionWith (+) p q++  (P p) * (P q) = P $ filterNonzero $ IntMap.foldlWithKey' accumulate IntMap.empty p+    where+      accumulate m e c = IntMap.unionWith (+) m $ multiplyTerm e c+      multiplyTerm e c = IntMap.mapKeysMonotonic (+ e) $ IntMap.map (* c) q++  abs = id+  signum 0 = 0+  signum _ = 1+  fromInteger 0 = P IntMap.empty+  fromInteger n = P $ IntMap.singleton 0 $ fromInteger n+  negate (P m) = P $ IntMap.map negate m++-- | Get rid of zero coefficients to ensure that zero coefficients do not cause+-- two polynomials represented by an 'IntMap' are not considered different.+filterNonzero :: (Eq a, Num a) => IntMap a -> IntMap a+filterNonzero = IntMap.filter (/= 0)++instance (Eq a, Num a) => Polynomial P Int a where+  degree (P m) = maybe 0 fst $ IntMap.lookupMax m+  coefficient (P m) k = fromMaybe 0 $ IntMap.lookup k m+  leadingCoefficient (P m) = maybe 0 snd $ IntMap.lookupMax m+  deleteLeadingTerm (P m) = P $ IntMap.deleteMax m+  foldTerms f (P m) = IntMap.foldMapWithKey f m+  scale 0 _ = P IntMap.empty+  scale x (P m) = P $ IntMap.map (* x) m+  power n = P $ IntMap.singleton (fromIntegral n) 1
+ src/Symtegration/Polynomial/Solve.hs view
@@ -0,0 +1,260 @@+-- |+-- Module: Symtegration.Polynomial.Solve+-- Description: Derive the roots of polynomials with rational coefficients.+-- Copyright: Copyright 2025 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+--+-- This module supports deriving exact solutions to polynomial equations.+-- It cannot derive solutions for all polynomials; it will only return those which it can.+module Symtegration.Polynomial.Solve (solve, complexSolve) where++import Data.List (nub)+import Symtegration.Polynomial+import Symtegration.Polynomial.Indexed+import Symtegration.Symbolic++-- $setup+-- >>> import Symtegration+-- >>> import Symtegration.Polynomial++-- | Derive the roots for the given polynomial.  Only real roots are returned.+--+-- >>> map (toHaskell . simplify) <$> solve (2 * power 1 - 6)+-- Just ["3"]+--+-- >>> map (toHaskell . simplify) <$> solve (power 2 - 4)+-- Just ["2","-2"]+--+-- Returns 'Nothing' if the function does not know how to derive the roots.+solve :: IndexedPolynomial -> Maybe [Expression]+solve p+  | degree p == 1 = solveLinear (c 1) (c 0)+  | degree p == 2 = solveQuadratic (c 2) (c 1) (c 0)+  | degree p == 3 = solveCubic (c 3) (c 2) (c 1) (c 0)+  | degree p == 4 = solveQuartic (c 4) (c 3) (c 2) (c 1) (c 0)+  | otherwise = Nothing+  where+    c = coefficient p++-- | Returns the real root for a polynomial of degree 1.+solveLinear :: Rational -> Rational -> Maybe [Expression]+solveLinear a b = Just [fromRational ((-b) / a)]++-- | Returns the real roots for a polynomial of degree 2.+solveQuadratic :: Rational -> Rational -> Rational -> Maybe [Expression]+solveQuadratic a b c+  | sq == 0 = Just [fromRational $ (-b) / (2 * a)]+  | sq > 0 =+      Just+        [ ((-b') + sq' ** (1 / 2)) / (2 * a'),+          ((-b') - sq' ** (1 / 2)) / (2 * a')+        ]+  | otherwise = Just []+  where+    sq = b * b - 4 * a * c+    sq' = fromRational sq+    a' = fromRational a+    b' = fromRational b++-- | Returns the real roots for a polynomial of degree 3.+solveCubic :: Rational -> Rational -> Rational -> Rational -> Maybe [Expression]+solveCubic a b c d = map restore <$> depressedRoots+  where+    restore x = x - fromRational b / (3 * fromRational a)+    depressedRoots = solveDepressedCubic p q+    p = (3 * a * c - b ^ two) / (3 * a ^ two)+    q = (2 * b ^ three - 9 * a * b * c + 27 * a ^ two * d) / (27 * a ^ three)+    two = 2 :: Int+    three = 3 :: Int++-- | Solve depressed cubic equations of the form \(x^3 + px + q = 0\).+-- Only returns real roots.+--+-- #### References+--+-- * [Wikipedia](https://en.wikipedia.org/wiki/Cubic_equation#Trigonometric_and_hyperbolic_solutions)+-- * [Wolfram MathWorld](https://mathworld.wolfram.com/CubicFormula.html)+solveDepressedCubic :: Rational -> Rational -> Maybe [Expression]+solveDepressedCubic 0 q+  | q < 0 = Just [fromRational (-q) ** (1 / 3)]+  | otherwise = Just [negate $ fromRational q ** (1 / 3)]+solveDepressedCubic p q+  | s < 0 =+      let c = 2 * sqrt (-(p' / 3))+          theta = acos (3 / 2 * q' / p' * sqrt (-(3 / p'))) / 3+       in Just [c * cos theta, c * cos (theta - 2 * pi / 3), c * cos (theta - 4 * pi / 3)]+  | p < 0,+    s > 0 =+      Just [(-2) * signum q' * sqrt (-(p' / 3)) * cosh (acosh ((-3) / 2 * abs q' / p' * sqrt (-(3 / p'))) / 3)]+  | s == 0 = Just [fromRational (3 * q / p), fromRational ((-3) / 2 * q / p)]+  | p > 0 =+      Just [(-2) * sqrt (p' / 3) * sinh (asinh (3 / 2 * q' / p' * sqrt (3 / p')) / 3)]+  | otherwise = Nothing+  where+    s = 4 * p ^ (3 :: Int) + 27 * q ^ (2 :: Int)+    p' = fromRational p+    q' = fromRational q++-- | Returns the real roots for a polynomial of degree 4.+solveQuartic :: Rational -> Rational -> Rational -> Rational -> Rational -> Maybe [Expression]+solveQuartic a b 0 0 0+  | b /= 0 = Just [0, fromRational $ -(b / a)]+  | otherwise = Just [0]+solveQuartic a b c 0 0+  | (Just xs) <- solveQuadratic a b c = Just $ nub $ 0 : xs+  | otherwise = Nothing+solveQuartic a b c d 0+  | (Just xs) <- solveCubic a b c d = Just $ nub $ 0 : xs+  | otherwise = Nothing+solveQuartic a 0 0 0 b+  | a > 0, b > 0 = Just []+  | a < 0, b < 0 = Just []+  | b == 0 = Just [0]+  | otherwise = Just [x, -x]+  where+    x = fromRational ((-b) / a) ** (1 / 4)+solveQuartic a 0 b 0 c+  | sq < 0 = Just []+  | sq == 0, st < 0 = Just []+  | sq == 0 = Just [sqrt st', -sqrt st']+  | a > 0, sq > 0, b > 0, sq > b * b = Just [sqrt x1, -sqrt x1]+  | a < 0, sq > 0, b < 0, sq > b * b = Just [sqrt x2, -sqrt x2]+  | a > 0, sq > 0, b < 0, sq < b * b = Just [sqrt x1, -sqrt x1, sqrt x2, -sqrt x2]+  | a < 0, sq > 0, b > 0, sq < b * b = Just [sqrt x1, -sqrt x1, sqrt x2, -sqrt x2]+  | otherwise = Nothing+  where+    sq = b * b - 4 * a * c+    st = (-b) / (2 * a)+    sq' = fromRational sq+    st' = fromRational st+    a' = fromRational a+    b' = fromRational b+    x1 = ((-b') + sqrt sq') / (2 * a')+    x2 = ((-b') - sqrt sq') / (2 * a')+solveQuartic _ _ _ _ _ = Nothing++-- | Derive the roots for the given polynomial.+-- All roots are returned, including complex roots.+--+-- >>> map (toHaskell . simplify) <$> complexSolve (2 * power 1 - 6)+-- Just ["3"]+--+-- >>> map (toHaskell . simplify) <$> complexSolve (power 2 + 1)+-- Just ["(-1) ** (1 / 2)","(-1) * (-1) ** (1 / 2)"]+--+-- >>> map (toHaskell . simplify) <$> complexSolve (power 3 + 1)+-- Just ["-1","(-1) * ((-1) + (-1) ** (1 / 2) * 3 ** (1 / 2)) / 2","(-1) * ((-1) + (-1) * (-1) ** (1 / 2) * 3 ** (1 / 2)) / 2"]+--+-- Returns 'Nothing' if the function does not know how to derive the roots.+complexSolve :: IndexedPolynomial -> Maybe [Expression]+complexSolve p+  | degree p == 1 = complexSolveLinear (c 1) (c 0)+  | degree p == 2 = complexSolveQuadratic (c 2) (c 1) (c 0)+  | degree p == 3 = complexSolveCubic (c 3) (c 2) (c 1) (c 0)+  | degree p == 4 = complexSolveQuartic (c 4) (c 3) (c 2) (c 1) (c 0)+  | otherwise = Nothing+  where+    c = coefficient p++-- | Returns the roots for a polynomial of degree 1.+complexSolveLinear :: Rational -> Rational -> Maybe [Expression]+complexSolveLinear a b = Just [fromRational $ (-b) / a]++-- | Returns the roots for a polynomial of degree 2.+complexSolveQuadratic :: Rational -> Rational -> Rational -> Maybe [Expression]+complexSolveQuadratic a b c+  | sq == 0 = Just [p]+  | otherwise = Just [p + q, p - q]+  where+    sq = b * b - 4 * a * c+    p = fromRational $ (-b) / (2 * a)+    q = sqrt (fromRational sq) / fromRational (2 * a)++-- | Returns the roots for a polynomial of degree 3.+complexSolveCubic :: Rational -> Rational -> Rational -> Rational -> Maybe [Expression]+complexSolveCubic _ 0 0 0 = Just [0]+complexSolveCubic a b 0 0 = Just [0, fromRational $ (-b) / a]+complexSolveCubic a b c 0+  | Just xs <- complexSolveQuadratic a b c = Just $ nub $ 0 : xs+  | otherwise = Just [0]+complexSolveCubic a b c d = map restore <$> complexSolveDepressedCubic p q+  where+    restore t = t - fromRational (b / (3 * a))+    p = (3 * a * c - b * b) / (3 * a * a)+    q = (2 * b * b * b - 9 * a * b * c + 27 * a * a * d) / (27 * a * a * a)++-- | Solve depressed cubic equations of the form \(x^3 + px + q = 0\).+--+-- #### References+--+-- * [Wikipedia](https://en.wikipedia.org/wiki/Cubic_equation)+complexSolveDepressedCubic :: Rational -> Rational -> Maybe [Expression]+complexSolveDepressedCubic p q+  | discriminant == 0, p == 0 = Just [0]+  | discriminant == 0 = Just $ map fromRational $ nub [3 * q / p, (-3) / 2 * q / p]+  | p == 0 = Just [x * e | let x = fromRational (-q) ** (1 / 3), e <- [1, e1, e2]]+  | otherwise =+      Just+        [ c - fromRational p / (3 * c),+          c * e1 - fromRational p / (3 * c * e1),+          c * e2 - fromRational p / (3 * c * e2)+        ]+  where+    discriminant = -(4 * p * p * p + 27 * q * q)+    c = (fromRational (-(q / 2)) + sqrt (fromRational (q * q / 4 + p * p * p / 27))) ** (1 / 3)+    e1 = (-1 + sqrt (-3)) / 2+    e2 = (-1 - sqrt (-3)) / 2++-- | Returns the roots for a polynomial of degree 4.+complexSolveQuartic :: Rational -> Rational -> Rational -> Rational -> Rational -> Maybe [Expression]+complexSolveQuartic _ 0 0 0 0 = Just [0]+complexSolveQuartic a b 0 0 0 = Just $ nub [0, fromRational $ -(b / a)]+complexSolveQuartic a b c 0 0+  | Just xs <- complexSolveQuadratic a b c = Just $ nub $ 0 : xs+  | otherwise = Just [0]+complexSolveQuartic a b c d 0+  | Just xs <- complexSolveCubic a b c d = Just $ nub $ 0 : xs+  | otherwise = Just [0]+complexSolveQuartic a 0 b 0 c = concatMap restore <$> complexSolveQuadratic a b c+  where+    restore 0 = [0]+    restore x = [sqrt x, -sqrt x]+complexSolveQuartic a b c d e = map restore <$> complexSolveDepressedQuartic p q r+  where+    restore x = x - fromRational (b / (4 * a))++    p = (-3) * b ^ two / (8 * a ^ two) + c / a+    q = b ^ three / (8 * a ^ three) - b * c / (2 * a ^ two) + d / a+    r = (-3) * b ^ four / (256 * a ^ four) + c * b ^ two / (16 * a ^ three) - b * d / (4 * a ^ two) + e / a++    two = 2 :: Int+    three = 3 :: Int+    four = 4 :: Int++-- | Returns the roots for a depressed quartic equation \(x^4+ax^2+bx+c=0\).+-- Complex numbers roots are included.+--+-- #### References+--+-- * [Wikipedia](https://en.wikipedia.org/wiki/Quartic_equation#The_general_case)+complexSolveDepressedQuartic :: Rational -> Rational -> Rational -> Maybe [Expression]+complexSolveDepressedQuartic a 0 c = concatMap restore <$> complexSolveQuadratic 1 a c+  where+    restore 0 = [0]+    restore x = [sqrt x, -sqrt x]+complexSolveDepressedQuartic a b c = do+  -- Get any cubic root of the following cubic equation.+  ys <- complexSolveCubic 2 (-a) (-(2 * c)) (a * c - b * b / 4)+  y <- case ys of x : _ -> Just x; [] -> Nothing++  -- Because b /= 0, it is the case that s /= 0.+  let s = sqrt $ 2 * y - fromRational a+  let t = (-2) * y - fromRational a++  return+    [ (1 / 2) * (-s + sqrt (t + 2 * fromRational b / s)),+      (1 / 2) * (-s - sqrt (t + 2 * fromRational b / s)),+      (1 / 2) * (s + sqrt (t - 2 * fromRational b / s)),+      (1 / 2) * (s - sqrt (t - 2 * fromRational b / s))+    ]
+ src/Symtegration/Polynomial/Symbolic.hs view
@@ -0,0 +1,158 @@+-- |+-- Module: Symtegration.Polynomial.Symbolic+-- Description: Conversion between data structures storing general mathematical expressions and those specialized for storing polynomials.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Polynomial.Symbolic+  ( -- * Converting expression to polynomial+    fromExpression,+    forVariable,+    withSymbolicCoefficients,++    -- * Converting polynomial to expression+    toExpression,+    toRationalCoefficient,+    toSymbolicCoefficient,+  )+where++import Data.Maybe (fromMaybe)+import Data.Monoid (Sum (..))+import Data.Ratio (denominator, numerator)+import Data.Text (Text)+import Symtegration.Polynomial+import Symtegration.Symbolic++-- $setup+-- >>> import Symtegration+-- >>> import Symtegration.Polynomial+-- >>> import Symtegration.Polynomial.Indexed++-- | Converts an 'Expression' into a 'Polynomial'.+-- 'Nothing' will be returned if the conversion is not possible.+--+-- Specify the symbol representing the variable for the polynomial with 'forVariable'.+-- For example,+--+-- >>> fromExpression (forVariable "x") (("x" + 4) ** 3) :: Maybe IndexedPolynomial+-- Just x^3 + 12x^2 + 48x + 64+--+-- By default, symbols other than the variable for the polynomial are not allowed.+-- To use symbols representing constants, use 'withSymbolicCoefficients' as well.+-- Note that the polynomial type the expression is being converted into+-- must be able to handle symbolic mathematical expressions for the coefficients.+-- For example,+--+-- >>> let expr = ("a" + "b") * "x" + "c" :: Expression+-- >>> let (Just p) = fromExpression (withSymbolicCoefficients (forVariable "x")) expr :: Maybe IndexedSymbolicPolynomial+-- >>> toHaskell $ simplify $ coefficient p 1+-- "a + b"+--+-- The expressions which can be converted must only use 'negate', '(+)', '(*)', '(-)',+-- '(/)' with only numbers, coefficients which do not contain the variable,+-- '(**)' with a non-negative integral exponent, and expressions formed thereof.+fromExpression ::+  (Polynomial p e c, Num (p e c), Fractional c) =>+  (Text -> Maybe (p e c), Expression -> Maybe c) ->+  Expression ->+  Maybe (p e c)+fromExpression _ (Number n) = Just $ fromInteger n+fromExpression (cf, _) (Symbol x) = cf x+fromExpression t (Negate' x) = negate <$> fromExpression t x+fromExpression t (x :+: y) = (+) <$> fromExpression t x <*> fromExpression t y+fromExpression t (x :*: y) = (*) <$> fromExpression t x <*> fromExpression t y+fromExpression t (x :-: y) = (-) <$> fromExpression t x <*> fromExpression t y+fromExpression t (x :**: (Number n))+  | n >= 0 = (^ n) <$> fromExpression t x+  | otherwise = Nothing+fromExpression _ (_ :**: _) = Nothing+fromExpression _ (_ :/: Number 0) = Nothing+fromExpression _ (Number n :/: Number m) = Just $ scale r 1+  where+    r = fromInteger n / fromInteger m+fromExpression (_, eval) e+  | Just e' <- eval e = Just $ scale e' 1+  | otherwise = Nothing++-- | Specifies the symbol representing the variable for 'fromExpression'.+forVariable ::+  (Polynomial p e c, Num (p e c), Fractional c) =>+  Text ->+  (Text -> Maybe (p e c), Expression -> Maybe c)+forVariable v = (fromSymbol, toCoefficient)+  where+    fromSymbol s+      | v == s = Just $ power 1+      | otherwise = Nothing++    toCoefficient (Symbol _) = Nothing+    toCoefficient (Number n) = Just $ fromInteger n+    toCoefficient (Negate' x) = negate <$> toCoefficient x+    toCoefficient (Abs' x) = abs <$> toCoefficient x+    toCoefficient (Signum' x) = signum <$> toCoefficient x+    toCoefficient (x :+: y) = (+) <$> toCoefficient x <*> toCoefficient y+    toCoefficient (x :*: y) = (*) <$> toCoefficient x <*> toCoefficient y+    toCoefficient (x :-: y) = (-) <$> toCoefficient x <*> toCoefficient y+    toCoefficient (_ :/: 0) = Nothing+    toCoefficient (x :/: y) = (/) <$> toCoefficient x <*> toCoefficient y+    toCoefficient (x :**: (Number n)) = (^^ n) <$> toCoefficient x+    toCoefficient _ = Nothing++-- | Specifies that non-variable symbols are allowed for 'fromExpression'.+-- The coefficients will be represented by 'Expression' values.+withSymbolicCoefficients ::+  (Polynomial p e Expression, Num (p e Expression), Integral e) =>+  (Text -> Maybe (p e Expression), Expression -> Maybe Expression) ->+  (Text -> Maybe (p e Expression), Expression -> Maybe Expression)+withSymbolicCoefficients (fromSymbol, _) = (fromSymbol', toCoefficient)+  where+    fromSymbol' s = Just $ fromMaybe (scale (Symbol s) 1) (fromSymbol s)++    toCoefficient e@(Symbol s)+      | Nothing <- fromSymbol s = Just e+      | otherwise = Nothing+    toCoefficient e@(Number _) = Just e+    toCoefficient (UnaryApply func x) = UnaryApply func <$> toCoefficient x+    toCoefficient (BinaryApply func x y) = BinaryApply func <$> x' <*> y'+      where+        x' = toCoefficient x+        y' = toCoefficient y++-- | Converts a 'Polynomial' into an 'Expression'.+-- The symbol which will represent the variable is the first argument.+--+-- How the coefficients are converted must also be specified.+-- To evaluate the coefficients to an exact rational number,+-- use 'toRationalCoefficient'.  For example,+--+-- >>> let (Just p) = fromExpression (forVariable "x") (3 * "x"**4 + 1) :: Maybe IndexedPolynomial+-- >>> toHaskell $ simplify $ toExpression "x" toRationalCoefficient p+-- "1 + 3 * x ** 4"+--+-- To evaluate the coefficients symbolically, use 'toSymbolicCoefficient'.+--+-- >>> let (Just p) = fromExpression (withSymbolicCoefficients (forVariable "x")) (("a"+"b") * "x"**4 + 1) :: Maybe IndexedSymbolicPolynomial+-- >>> toHaskell $ simplify $ toExpression "x" toSymbolicCoefficient p+-- "1 + x ** 4 * (a + b)"+toExpression :: (Polynomial p e c) => Text -> (c -> Expression) -> p e c -> Expression+toExpression x cf p = getSum $ foldTerms convert p+  where+    convert 0 c = Sum $ cf c+    convert e c = Sum $ cf c * xp+      where+        xp = Symbol x ** Number (fromIntegral e)++-- | Specifies that coefficients are numbers for 'toExpression'.+toRationalCoefficient :: (Real c) => c -> Expression+toRationalCoefficient c+  | d == 1 = Number n+  | otherwise = Number n :/: Number d+  where+    r = toRational c+    n = fromInteger $ numerator r+    d = fromInteger $ denominator r++-- | Specifies that coefficients are symbolic for 'toExpression'.+toSymbolicCoefficient :: Expression -> Expression+toSymbolicCoefficient = id
+ src/Symtegration/Symbolic.hs view
@@ -0,0 +1,423 @@+{-# LANGUAGE DerivingVia #-}+{-# LANGUAGE PatternSynonyms #-}++-- |+-- Module: Symtegration.Symbolic+-- Description: Library for symbolically representing mathematical expressions.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Symbolic+  ( -- * Representation+    Expression (..),+    UnaryFunction (..),+    BinaryFunction (..),++    -- * Manipulation+    substitute,++    -- * Computation+    evaluate,+    fractionalEvaluate,+    toFunction,+    getUnaryFunction,+    getBinaryFunction,++    -- * Pattern synonyms++    -- | Pattern synonyms are defined to make it more convenient to pattern match on 'Expression'.++    -- ** Constants+    pattern Pi',++    -- ** Unary functions+    pattern Negate',+    pattern Abs',+    pattern Signum',+    pattern Exp',+    pattern Log',+    pattern Sqrt',+    pattern Sin',+    pattern Cos',+    pattern Tan',+    pattern Asin',+    pattern Acos',+    pattern Atan',+    pattern Sinh',+    pattern Cosh',+    pattern Tanh',+    pattern Asinh',+    pattern Acosh',+    pattern Atanh',++    -- ** Binary functions+    pattern (:+:),+    pattern (:*:),+    pattern (:-:),+    pattern (:/:),+    pattern (:**:),+    pattern LogBase',+  )+where++import Data.Ratio+import Data.String (IsString, fromString)+import Data.Text+import GHC.Generics (Generic)+import TextShow (TextShow)+import TextShow.Generic (FromGeneric (..))++-- $setup+-- >>> import Symtegration++-- | Symbolic representation of a mathematical expression.+-- It is an instance of the 'Num', 'Fractional', and 'Floating' type classes,+-- so normal Haskell expressions can be used, although the expressions+-- are limited to using the functions defined by these type classses.+-- The type is also an instance of the 'IsString' type class,+-- so symbols can be expressed as Haskell string with the @OverloadedStrings@ extension.+-- The structure of these values is intended to be visible.+--+-- >>> 2 :: Expression+-- Number 2+-- >>> "x" :: Expression+-- Symbol "x"+-- >>> 2 + sin "x" :: Expression+-- BinaryApply Add (Number 2) (UnaryApply Sin (Symbol "x"))+--+-- A somewhat more concise representation can be obtained using 'Symtegration.toHaskell':+--+-- >>> toHaskell $ 2 * "y" + sin "x"+-- "2 * y + sin x"+data Expression+  = -- | Represents a concrete number.+    Number Integer+  | -- | Represents a symbol, which could either be a variable or a constant.+    Symbol Text+  | -- | Represents the application of an unary function.+    UnaryApply UnaryFunction Expression+  | -- | Represents the application of a binary function.+    BinaryApply BinaryFunction Expression Expression+  deriving+    ( -- | Structural equality, not semantic equality.+      -- E.g., @"a" - "a" /= 0@.+      Eq,+      Show,+      Read,+      Generic+    )+  deriving (TextShow) via FromGeneric Expression++pattern Pi' :: Expression+pattern Pi' = Symbol "pi"++-- | Symbolic representation for unary functions.+data UnaryFunction+  = -- | 'negate'+    Negate+  | -- | 'abs'+    Abs+  | -- | 'signum'+    Signum+  | -- | 'exp'+    Exp+  | -- | 'log'+    Log+  | -- | 'sqrt'+    Sqrt+  | -- | 'sin'+    Sin+  | -- | 'cos'+    Cos+  | -- | 'tan'+    Tan+  | -- | 'asin'+    Asin+  | -- | 'acos'+    Acos+  | -- | 'atan'+    Atan+  | -- | 'sinh'+    Sinh+  | -- | 'cosh'+    Cosh+  | -- | 'tanh'+    Tanh+  | -- | 'asinh'+    Asinh+  | -- | 'acosh'+    Acosh+  | -- | 'atanh'+    Atanh+  deriving (Eq, Enum, Bounded, Show, Read, Generic)+  deriving (TextShow) via FromGeneric UnaryFunction++pattern Negate', Abs', Signum', Exp', Log', Sqrt', Sin', Cos', Tan', Asin', Acos', Atan', Sinh', Cosh', Tanh', Asinh', Acosh', Atanh' :: Expression -> Expression+pattern Negate' x = UnaryApply Negate x+pattern Abs' x = UnaryApply Abs x+pattern Signum' x = UnaryApply Signum x+pattern Exp' x = UnaryApply Exp x+pattern Log' x = UnaryApply Log x+pattern Sqrt' x = UnaryApply Sqrt x+pattern Sin' x = UnaryApply Sin x+pattern Cos' x = UnaryApply Cos x+pattern Tan' x = UnaryApply Tan x+pattern Asin' x = UnaryApply Asin x+pattern Acos' x = UnaryApply Acos x+pattern Atan' x = UnaryApply Atan x+pattern Sinh' x = UnaryApply Sinh x+pattern Cosh' x = UnaryApply Cosh x+pattern Tanh' x = UnaryApply Tanh x+pattern Asinh' x = UnaryApply Asinh x+pattern Acosh' x = UnaryApply Acosh x+pattern Atanh' x = UnaryApply Atanh x++-- | Symbolic representation for binary functions.+data BinaryFunction+  = -- | '(+)'+    Add+  | -- | '(*)'+    Multiply+  | -- | '(-)'+    Subtract+  | -- | '(/)'+    Divide+  | -- | '(**)'+    Power+  | -- | 'logBase'+    LogBase+  deriving (Eq, Enum, Bounded, Show, Read, Generic)+  deriving (TextShow) via FromGeneric BinaryFunction++pattern (:+:), (:*:), (:-:), (:/:), (:**:), LogBase' :: Expression -> Expression -> Expression+pattern x :+: y = BinaryApply Add x y+pattern x :*: y = BinaryApply Multiply x y+pattern x :-: y = BinaryApply Subtract x y+pattern x :/: y = BinaryApply Divide x y+pattern x :**: y = BinaryApply Power x y+pattern LogBase' x y = BinaryApply LogBase x y++instance IsString Expression where+  fromString = Symbol . fromString++instance Num Expression where+  (+) = BinaryApply Add+  (-) = BinaryApply Subtract+  (*) = BinaryApply Multiply+  negate = UnaryApply Negate+  abs = UnaryApply Abs+  signum = UnaryApply Signum+  fromInteger = Number++instance Fractional Expression where+  (/) = BinaryApply Divide+  fromRational q | d == 1 = n | otherwise = BinaryApply Divide n d+    where+      n = Number $ numerator q+      d = Number $ denominator q++instance Floating Expression where+  pi = Symbol "pi"+  exp = UnaryApply Exp+  log = UnaryApply Log+  sqrt = UnaryApply Sqrt+  (**) = BinaryApply Power+  logBase = BinaryApply LogBase+  sin = UnaryApply Sin+  cos = UnaryApply Cos+  tan = UnaryApply Tan+  asin = UnaryApply Asin+  acos = UnaryApply Acos+  atan = UnaryApply Atan+  sinh = UnaryApply Sinh+  cosh = UnaryApply Cosh+  tanh = UnaryApply Tanh+  asinh = UnaryApply Asinh+  acosh = UnaryApply Acosh+  atanh = UnaryApply Atanh++-- | Returns a function corresponding to the symbolic representation of an unary function.+--+-- >>> (getUnaryFunction Cos) pi == (cos pi :: Double)+-- True+getUnaryFunction :: (Floating a) => UnaryFunction -> (a -> a)+getUnaryFunction Negate = negate+getUnaryFunction Abs = abs+getUnaryFunction Signum = signum+getUnaryFunction Exp = exp+getUnaryFunction Log = log+getUnaryFunction Sqrt = sqrt+getUnaryFunction Sin = sin+getUnaryFunction Cos = cos+getUnaryFunction Tan = tan+getUnaryFunction Asin = asin+getUnaryFunction Acos = acos+getUnaryFunction Atan = atan+getUnaryFunction Sinh = sinh+getUnaryFunction Cosh = cosh+getUnaryFunction Tanh = tanh+getUnaryFunction Asinh = asinh+getUnaryFunction Acosh = acosh+getUnaryFunction Atanh = atanh++-- | Returns a function corresponding to the symbolic representation of a binary function.+--+-- >>> (getBinaryFunction Add) 2 5 == (2 + 5 :: Double)+-- True+getBinaryFunction :: (Floating a) => BinaryFunction -> (a -> a -> a)+getBinaryFunction Add = (+)+getBinaryFunction Multiply = (*)+getBinaryFunction Subtract = (-)+getBinaryFunction Divide = (/)+getBinaryFunction Power = (**)+getBinaryFunction LogBase = logBase++-- | Substitute the symbols with the corresponding expressions they are mapped to.+-- The symbols will be replaced as is; there is no special treatment if the+-- expression they are replaced by also contains the same symbol.+--+-- >>> toHaskell $ substitute ("x" + "y") (\case "x" -> Just ("a" * "b"); "y" -> Just 4)+-- "a * b + 4"+substitute ::+  -- | Expression to apply substitution.+  Expression ->+  -- | Maps symbols to expressions they are to be substituted with.+  (Text -> Maybe Expression) ->+  -- | Expression with substitution applied.+  Expression+substitute e@(Number _) _ = e+substitute e@(Symbol s) f+  | (Just x) <- f s = x+  | otherwise = e+substitute (UnaryApply func x) f = UnaryApply func (substitute x f)+substitute (BinaryApply func x y) f = BinaryApply func (substitute x f) (substitute y f)++-- | Calculates the value for a mathematical expression for a given assignment of values to symbols.+--+-- For example, when \(x=5\), then \(2x+1=11\).+--+-- >>> evaluate (2 * "x" + 1) (\case "x" -> Just 5)+-- Just 11.0+--+-- All symbols except for @"pi"@ in a mathematical expression must be assigned a value.+-- Otherwise, a value cannot be computed.+--+-- >>> evaluate (2 * "x" + 1) (const Nothing)+-- Nothing+--+-- The symbol @"pi"@ is always used to represent \(\pi\),+-- and any assignment to @"pi"@ will be ignored.+-- For example, the following is \(\pi - \pi\), not \(100 - \pi\).+--+-- >>> evaluate ("pi" - pi) (\case "x" -> Just 100)+-- Just 0.0+evaluate ::+  (Floating a) =>+  -- | Mathematical expression to evaluate.+  Expression ->+  -- | Maps symbols to concrete values.+  (Text -> Maybe a) ->+  -- | Evaluation result.+  Maybe a+evaluate (Number n) _ = Just $ fromInteger n+evaluate (Symbol "pi") _ = Just pi+evaluate (Symbol x) m = m x+evaluate (UnaryApply fun expr) m = fmap f v+  where+    f = getUnaryFunction fun+    v = evaluate expr m+evaluate (BinaryApply fun expr1 expr2) m = f <$> v1 <*> v2+  where+    f = getBinaryFunction fun+    v1 = evaluate expr1 m+    v2 = evaluate expr2 m++-- |+-- Evaluates a mathematical expression with only operations available to 'Fractional' values.+-- In particular, this allows exact evaluations with 'Rational' values.+-- 'Nothing' will be returned if a function not supported by all 'Fractional' values+-- is used by the mathematical expression.+--+-- As an exception, the '(**)' operator is allowed with constant integer exponents,+-- even though '(**)' is not a function applicable to all 'Fractional' types.+--+-- For example,+--+-- >>> let p = 1 / (3 * "x"**5 - 2 * "x" + 1) :: Expression+-- >>> fractionalEvaluate p (\case "x" -> Just (2 / 7 :: Rational))+-- Just (16807 % 7299)+--+-- Compare against 'evaluate', which cannot even use 'Rational' computations+-- because 'Rational' is not an instance of the 'Floating' type class:+--+-- >>> evaluate p (\case "x" -> Just (2 / 7 :: Double))+-- Just 2.3026441978353196+fractionalEvaluate ::+  (Eq a, Fractional a) =>+  -- | Mathematical expression to evaluate.+  Expression ->+  -- | Maps symbols to concrete values.+  (Text -> Maybe a) ->+  -- | Evaluation result.+  Maybe a+fractionalEvaluate (Number n) _ = Just $ fromInteger n+fractionalEvaluate (Symbol x) m = m x+fractionalEvaluate (Negate' x) m = negate <$> fractionalEvaluate x m+fractionalEvaluate (Abs' x) m = abs <$> fractionalEvaluate x m+fractionalEvaluate (Signum' x) m = signum <$> fractionalEvaluate x m+fractionalEvaluate (x :+: y) m = (+) <$> fractionalEvaluate x m <*> fractionalEvaluate y m+fractionalEvaluate (x :-: y) m = (-) <$> fractionalEvaluate x m <*> fractionalEvaluate y m+fractionalEvaluate (x :*: y) m = (*) <$> fractionalEvaluate x m <*> fractionalEvaluate y m+fractionalEvaluate (x :/: y) m+  | Just 0 <- y' = Nothing+  | otherwise = (/) <$> x' <*> y'+  where+    x' = fractionalEvaluate x m+    y' = fractionalEvaluate y m+fractionalEvaluate (x :**: (Number n)) m = (^^ n) <$> fractionalEvaluate x m+fractionalEvaluate _ _ = Nothing++-- | Returns a function based on a given expression.  This requires+-- a specification of how a symbol maps the argument to a value+-- to be used in its place.+--+-- For example, the symbol "x" could use the argument as is as its value.+-- I.e., "x" can be mapped to a function which maps the argument to itself.+--+-- >>> let f = toFunction ("x" ** 2 + 1) (\case "x" -> id) :: Double -> Double+-- >>> f 3  -- 3 ** 2 + 1+-- 10.0+-- >>> f 10  -- 10 ** 2 + 1+-- 101.0+--+-- For another example, "x" could map the first element from a tuple argument,+-- and "y" could map the second element from the tuple argument.  I.e.,+-- for a tuple argument to the function, the first element will be used as "x"+-- and the second element will be used as "y".+--+-- >>> let m = \case "x" -> (\(x,_) -> x); "y" -> (\(_,y) -> y)+-- >>> let g = toFunction ("x" + 2 * "y") m :: (Double, Double) -> Double+-- >>> g (3,4)  -- 3 + 2 * 4+-- 11.0+-- >>> g (7,1)  -- 7 + 2 * 1+-- 9.0+toFunction ::+  (Floating b) =>+  -- | The expression to be converted into a function.+  Expression ->+  -- | Maps how the argument to the function should be mapped to a value for a symbol.+  -- E.g., "x" could map the first element in a tuple as the value to use in its place.+  (Text -> (a -> b)) ->+  -- | The function generated from the expression.+  (a -> b)+toFunction (Number n) _ = const $ fromInteger n+toFunction (Symbol s) m = m s+toFunction (UnaryApply func x) m = f . g+  where+    f = getUnaryFunction func+    g = toFunction x m+toFunction (BinaryApply func x y) m = \v -> f (g v) (h v)+  where+    f = getBinaryFunction func+    g = toFunction x m+    h = toFunction y m
+ src/Symtegration/Symbolic/Haskell.hs view
@@ -0,0 +1,118 @@+-- |+-- Module: Symtegration.Symbolic.Haskell+-- Description: Converts a symbolic representation of a mathematical expression into equivalent Haskell code.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+--+-- Support for converting symbolic representations of mathematical expressions+-- into equivalent Haskell code.+module Symtegration.Symbolic.Haskell+  ( toHaskell,++    -- * Support functions+    getUnaryFunctionText,+    getBinaryFunctionText,+  )+where++import Data.Text+import Symtegration.Symbolic+import TextShow (showt)++-- $setup+-- >>> import Symtegration.Symbolic++-- | Converts an 'Expression' into an equivalent Haskell expression.+--+-- >>> toHaskell $ BinaryApply Add (Number 1) (Number 3)+-- "1 + 3"+-- >>> toHaskell $ 1 + 3+-- "1 + 3"+--+-- Symbols are included without quotation.+--+-- >>> toHaskell $ ("x" + "y") * 4+-- "(x + y) * 4"+toHaskell :: Expression -> Text+toHaskell (Number n) = showt n+toHaskell (Symbol t) = t+toHaskell (UnaryApply fun x) = funcText <> " " <> asArg x+  where+    funcText = getUnaryFunctionText fun+toHaskell (LogBase' x y) = funcText <> " " <> asArg x <> " " <> asArg y+  where+    funcText = getBinaryFunctionText LogBase+toHaskell (x :+: y) = asAddArg x <> " + " <> asAddArg y+toHaskell (x :-: y@(_ :+: _)) = asAddArg x <> " - " <> asArg y+toHaskell (x :-: y@(_ :-: _)) = asAddArg x <> " - " <> asArg y+toHaskell (x :-: y) = asAddArg x <> " - " <> asAddArg y+toHaskell (x :*: y) = asMultiplyArg x <> " * " <> asMultiplyArg y+toHaskell (BinaryApply op x y) = asArg x <> " " <> opText <> " " <> asArg y+  where+    opText = getBinaryFunctionText op++-- | Converts an 'Expression' to Haskell code appropriate for use as an argument.+-- In other words, show numbers and symbols as is, while surrounding everything+-- else in parentheses.+asArg :: Expression -> Text+asArg x@(Number n)+  | n >= 0 = toHaskell x+  | otherwise = "(" <> toHaskell x <> ")"+asArg x@(Symbol _) = toHaskell x+asArg x = par $ toHaskell x++-- | Converts an 'Expression' to an argument appropriate for addition.+asAddArg :: Expression -> Text+asAddArg x@(Number _) = asArg x+asAddArg x@(Symbol _) = asArg x+-- No operation has lower precedence than addition,+-- and addition is commutative, so no parentheses are needed.+asAddArg x = toHaskell x++-- | Converts an 'Expression' to an argument appropriate for multiplication.+asMultiplyArg :: Expression -> Text+asMultiplyArg x@(Number _) = asArg x+asMultiplyArg x@(Symbol _) = asArg x+asMultiplyArg x@(_ :+: _) = par $ toHaskell x+asMultiplyArg x@(_ :-: _) = par $ toHaskell x+-- No other operation has lower precedence than multiplication,+-- and multiplication is commutative, so no parentheses are needed.+asMultiplyArg x = toHaskell x++-- | Surrounds text by parentheses.+par :: Text -> Text+par s = "(" <> s <> ")"++-- | Returns the corresponding Haskell function name.+getUnaryFunctionText :: UnaryFunction -> Text+getUnaryFunctionText Negate = "negate"+getUnaryFunctionText Abs = "abs"+getUnaryFunctionText Signum = "signum"+getUnaryFunctionText Exp = "exp"+getUnaryFunctionText Log = "log"+getUnaryFunctionText Sqrt = "sqrt"+getUnaryFunctionText Sin = "sin"+getUnaryFunctionText Cos = "cos"+getUnaryFunctionText Tan = "tan"+getUnaryFunctionText Asin = "asin"+getUnaryFunctionText Acos = "acos"+getUnaryFunctionText Atan = "atan"+getUnaryFunctionText Sinh = "sinh"+getUnaryFunctionText Cosh = "cosh"+getUnaryFunctionText Tanh = "tanh"+getUnaryFunctionText Asinh = "asinh"+getUnaryFunctionText Acosh = "acosh"+getUnaryFunctionText Atanh = "atanh"++-- | Returns the corresponding Haskell function name.+--+-- For binary operators, it will be the infix form.+-- In other words, @"+"@ will be returned for 'Add', not @"(+)"@.+getBinaryFunctionText :: BinaryFunction -> Text+getBinaryFunctionText Add = "+"+getBinaryFunctionText Multiply = "*"+getBinaryFunctionText Subtract = "-"+getBinaryFunctionText Divide = "/"+getBinaryFunctionText Power = "**"+getBinaryFunctionText LogBase = "logBase"
+ src/Symtegration/Symbolic/LaTeX.hs view
@@ -0,0 +1,137 @@+-- |+-- Module: Symtegration.Symbolic.LaTeX+-- Description: Converts a symbolic representation of a mathematical expression into equivalent LaTeX text.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+--+-- Support for converting symbolic representations of mathematical expressions+-- into equivalent LaTeX text.+module Symtegration.Symbolic.LaTeX (toLaTeX) where++import Data.Text (Text)+import Symtegration.Symbolic+import TextShow (showt)++-- | Converts an 'Expression' into an equivalent LaTeX expression.+--+-- >>> toLaTeX $ exp 5+-- "e^{5}"+--+-- Symbols are included without quotation.+--+-- >>> toLaTeX $ exp "x"+-- "e^{x}"+-- >>> toLaTeX $ "x" + 4 * sin "y"+-- "x + 4 \\sin y"+--+-- Since the text for symbols are included as is, we can also include LaTeX symbols:+--+-- >>> toLaTeX $ exp "\\delta_0"+-- "e^{\\delta_0}"+toLaTeX :: Expression -> Text+toLaTeX (Number n) = showt n+toLaTeX (Symbol "pi") = "\\pi"+toLaTeX (Symbol s) = s+toLaTeX (UnaryApply func x) = unary func x+toLaTeX (BinaryApply func x y) = binary func x y++-- | Converts unary functions into LaTeX.+unary :: UnaryFunction -> Expression -> Text+unary Negate x@(_ :+: _) = "-" <> asArg x+unary Negate x@(_ :-: _) = "-" <> asArg x+unary Negate x@(Negate' _) = "-" <> asArg x+unary Negate x = "-" <> toLaTeX x+unary Abs x = "\\left\\lvert " <> toLaTeX x <> " \\right\\rvert"+unary Signum x = "\\mathrm{signum}" <> par (toLaTeX x)+unary Exp x = "e^" <> brace (toLaTeX x)+unary Log x = "\\log " <> asNamedFunctionArg x+unary Sqrt x = "\\sqrt" <> brace (toLaTeX x)+unary Sin x = "\\sin " <> asNamedFunctionArg x+unary Cos x = "\\cos " <> asNamedFunctionArg x+unary Tan x = "\\tan " <> asNamedFunctionArg x+unary Asin x = "\\sin^{-1} " <> asNamedFunctionArg x+unary Acos x = "\\cos^{-1} " <> asNamedFunctionArg x+unary Atan x = "\\tan^{-1} " <> asNamedFunctionArg x+unary Sinh x = "\\sinh " <> asNamedFunctionArg x+unary Cosh x = "\\cosh " <> asNamedFunctionArg x+unary Tanh x = "\\tanh " <> asNamedFunctionArg x+unary Asinh x = "\\sinh^{-1} " <> asNamedFunctionArg x+unary Acosh x = "\\cosh^{-1} " <> asNamedFunctionArg x+unary Atanh x = "\\tanh^{-1} " <> asNamedFunctionArg x++-- | Converts binary functions into LaTeX.+binary :: BinaryFunction -> Expression -> Expression -> Text+binary Add x (Negate' y) = binary Subtract x y+binary Add x y = asAddInitialArg x <> " + " <> asAddTrailingArg y+binary Multiply x@(_ :*: Number _) y@(Number _ :*: _) = toLaTeX x <> " \\times " <> toLaTeX y+binary Multiply x@(Number _) y@(Number _ :*: _) = toLaTeX x <> " \\times " <> toLaTeX y+binary Multiply x@(_ :*: Number _) y@(Number _) = toLaTeX x <> " \\times " <> toLaTeX y+binary Multiply x y@(Number _) = asMultiplyArg x <> " \\times " <> asArg y+binary Multiply x@(Abs' _) y = toLaTeX x <> " " <> asMultiplyArg y+binary Multiply x@(Signum' _) y = toLaTeX x <> " " <> asMultiplyArg y+binary Multiply x@(Exp' _) y = toLaTeX x <> " " <> asMultiplyArg y+binary Multiply x@(Sqrt' _) y = toLaTeX x <> " " <> asMultiplyArg y+binary Multiply x@(UnaryApply _ _) y@(Symbol _) = par (toLaTeX x) <> " " <> asMultiplyArg y+binary Multiply x@(LogBase' _ _) y = par (toLaTeX x) <> " " <> asMultiplyArg y+binary Multiply x y = asMultiplyArg x <> " " <> asMultiplyArg y+binary Subtract x y@(Negate' _) = asAddInitialArg x <> " - " <> asArg y+binary Subtract x y@(_ :+: _) = asAddInitialArg x <> " - " <> asArg y+binary Subtract x y@(_ :-: _) = asAddInitialArg x <> " - " <> asArg y+binary Subtract x y = asAddInitialArg x <> " - " <> asAddTrailingArg y+binary Divide x y = "\\frac" <> brace (toLaTeX x) <> brace (toLaTeX y)+binary Power x y = asArg x <> "^" <> brace (toLaTeX y)+binary LogBase x y = "\\log_" <> brace (toLaTeX x) <> asNamedFunctionArg y++asArg :: Expression -> Text+asArg e@(Number n) | n >= 0 = toLaTeX e | otherwise = par $ toLaTeX e+asArg e@(Symbol _) = toLaTeX e+asArg e@(Negate' _) = par $ toLaTeX e+asArg e@(UnaryApply _ _) = toLaTeX e+asArg e@(_ :/: _) = toLaTeX e+asArg e@(Number _ :**: _) = par $ toLaTeX e+asArg e@(_ :**: _) = toLaTeX e+asArg e = par $ toLaTeX e++asAddInitialArg :: Expression -> Text+asAddInitialArg e@(Number _) = toLaTeX e+asAddInitialArg e@(Symbol _) = toLaTeX e+asAddInitialArg e@(Negate' _) = toLaTeX e+asAddInitialArg (x :+: y) = asAddInitialArg x <> " + " <> asAddTrailingArg y+asAddInitialArg (x :-: y@(Negate' _)) = asAddInitialArg x <> " - " <> asArg y+asAddInitialArg (x :-: y@(_ :+: _)) = asAddInitialArg x <> " - " <> asArg y+asAddInitialArg (x :-: y@(_ :-: _)) = asAddInitialArg x <> " - " <> asArg y+asAddInitialArg (x :-: y) = asAddInitialArg x <> " - " <> toLaTeX y+asAddInitialArg e = asAddTrailingArg e++asAddTrailingArg :: Expression -> Text+asAddTrailingArg e@(Number _) = asArg e+asAddTrailingArg e@(Symbol _) = asArg e+asAddTrailingArg e@(Negate' _) = asArg e+asAddTrailingArg e = toLaTeX e++asMultiplyArg :: Expression -> Text+asMultiplyArg e@(Number _) = asArg e+asMultiplyArg e@(Symbol _) = asArg e+asMultiplyArg e@(Negate' _) = asArg e+asMultiplyArg e@(UnaryApply _ _) = toLaTeX e+asMultiplyArg e@(_ :+: _) = asArg e+asMultiplyArg e@(_ :-: _) = asArg e+asMultiplyArg e@(BinaryApply _ _ _) = toLaTeX e++-- For arguments to named functions such as "sin" which do not always delimit their arguments.+-- E.g., it is preferred that "1 + sin x" be "1 + sin x" and not "1 + (sin x)",+-- but we want "cos (sin x)" to be "cos (sin x)" and not "cos sin x".+asNamedFunctionArg :: Expression -> Text+asNamedFunctionArg e@(Exp' _) = asArg e+asNamedFunctionArg e@(Abs' _) = asArg e+asNamedFunctionArg e@(Sqrt' _) = asArg e+asNamedFunctionArg e@(UnaryApply _ _) = par $ toLaTeX e+asNamedFunctionArg e@(LogBase' _ _) = par $ toLaTeX e+asNamedFunctionArg e = asArg e++par :: Text -> Text+par s = "\\left(" <> s <> "\\right)"++brace :: Text -> Text+brace s = "{" <> s <> "}"
+ src/Symtegration/Symbolic/Simplify.hs view
@@ -0,0 +1,61 @@+-- |+-- Module: Symtegration.Symbolic.Simplify+-- Description: Simplifes symbolic representations of mathematical expressions.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+--+-- Supports the simplification of the symbolic representation for a mathematical expression.+-- This is aimed towards making it easier to find common factors for the purpose of integration.+-- It requires the specification of which symbol represents the variable.+module Symtegration.Symbolic.Simplify (simplify, tidy, simplifyForVariable) where++import Data.Text (Text)+import Symtegration.Symbolic+import Symtegration.Symbolic.Simplify.AlgebraicRingOrder qualified as AlgebraicRingOrder+import Symtegration.Symbolic.Simplify.Fraction qualified as Fraction+import Symtegration.Symbolic.Simplify.NumericFolding qualified as NumericFolding+import Symtegration.Symbolic.Simplify.SymbolicFolding qualified as SymbolicFolding+import Symtegration.Symbolic.Simplify.Tidy++-- $setup+-- >>> import Symtegration.Symbolic.Haskell++-- | Simplifies symbolic representations of mathematical expressions.+--+-- All addition and multiplication will be associated to the left.+-- The simplification is done with an eye towards making it+-- easier to find common factors.+--+-- >>> toHaskell $ simplify $ 4 - "x" + "a" * "x" ** 3 + 2 * "x" - 3+-- "1 + x + a * x ** 3"+simplify :: Expression -> Expression+simplify = simplifyForVariable ""++-- | Simplifies symbolic representations of mathematical expressions+-- with special consideration for a particular variable.+--+-- All addition and multiplication will be associated to the left.+-- Terms with higher orders of the variable will appear later.+-- The simplification is done with an eye towards making it+-- easier to find common factors.+--+-- >>> toHaskell $ simplifyForVariable "x" $ 1 + "a" * "x" ** 3 + "x"+-- "1 + x + a * x ** 3"+-- >>> toHaskell $ simplifyForVariable "x" $ "a" ** 143 + "x" + "b" ** 2+-- "a ** 143 + b ** 2 + x"+-- >>> toHaskell $ simplifyForVariable "x" $ "a" * "x" + "x" + "b ** 2" + "x" ** 2+-- "b ** 2 + x + a * x + x ** 2"+simplifyForVariable ::+  -- | Symbol for the variable.+  Text ->+  -- | Expression to be simplified.+  Expression ->+  -- | Simplified expression.+  Expression+simplifyForVariable v e+  | e == e' = e+  | otherwise = simplifyForVariable v e' -- Another round.+  where+    e' = f e+    f = Fraction.simplify . NumericFolding.simplify . SymbolicFolding.simplify . AlgebraicRingOrder.order v
+ src/Symtegration/Symbolic/Simplify/AlgebraicRingOrder.hs view
@@ -0,0 +1,176 @@+-- |+-- Module: Symtegration.Symbolic.Simplify.AlgebraicRingOrder+-- Description: Order terms in a mathematical expression according to a deterministic order.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Symbolic.Simplify.AlgebraicRingOrder (order) where++import Data.List (sortBy)+import Data.Set qualified as Set+import Data.Text (Text)+import Symtegration.Symbolic++-- $setup+-- >>> import Symtegration.Symbolic.Haskell++-- | Order terms in an mathematical expression.+--+-- Terms will be ordered according to a deterministic set of rules.+-- The re-ordering aims to make it easier to identify common factors and terms.+-- Terms with higher integral powers of the variable are sorted later.+-- Addition and multiplication will be re-arranged to associate to the left.+--+-- >>> toHaskell $ order "x" $ "x" + 1+-- "1 + x"+-- >>> toHaskell $ order "x" $ 2 + 3 * "x"**2 + "x"+-- "2 + x + 3 * x ** 2"+order ::+  -- | Symbol representing the variable.+  Text ->+  -- | Expression to be ordered.+  Expression ->+  -- | Expression with terms re-ordered.+  Expression+order _ e@(Number _) = e+order _ e@(Symbol _) = e+order v (UnaryApply func x) = UnaryApply func $ order v x+order v (x :/: y) = order v x :/: order v y+order v (x :**: y) = order v x :**: order v y+order v (LogBase' x y) = LogBase' (order v x) (order v y)+order v e = fromAddList $ sortBy (compareExpressions v) orderedAddTerms+  where+    terms = toAddMultiplyList v e+    orderedAddTerms = map (fromMultiplyList . sortBy (compareExpressions v)) terms++-- | Gather additive terms formed out of multiplicative terms.+-- No particular ordering should be expected.+toAddMultiplyList :: Text -> Expression -> [[Expression]]+toAddMultiplyList v (x@(_ :+: _) :+: y@(_ :+: _)) = toAddMultiplyList v x ++ toAddMultiplyList v y+toAddMultiplyList v (x@(_ :+: _) :+: y) = toMultiplyList v y : toAddMultiplyList v x+toAddMultiplyList v (x :+: y@(_ :+: _)) = toMultiplyList v x : toAddMultiplyList v y+toAddMultiplyList v (x :+: y) = map (toMultiplyList v) [x, y]+toAddMultiplyList v (x :-: y) = toAddMultiplyList v (x :+: (Number (-1) :*: y))+toAddMultiplyList v x = [toMultiplyList v x]++-- | Gather multiplicative terms.+-- No particular ordering should be expected.+toMultiplyList :: Text -> Expression -> [Expression]+toMultiplyList v (x@(_ :*: _) :*: y@(_ :*: _)) = toMultiplyList v x ++ toMultiplyList v y+toMultiplyList v (x@(_ :*: _) :*: y) = y : toMultiplyList v x+toMultiplyList v (x :*: y@(_ :*: _)) = x : toMultiplyList v y+toMultiplyList v (x :*: y) = [order v x, order v y]+toMultiplyList _ x@(Number _) = [x]+toMultiplyList _ x@(Symbol _) = [x]+toMultiplyList v (Negate' x) = Number (-1) : toMultiplyList v x+toMultiplyList v (UnaryApply func x) = [UnaryApply func $ order v x]+toMultiplyList v (BinaryApply func x y) = [BinaryApply func (order v x) (order v y)]++-- | Convert a list of sub-expressions for a multiplicative term into a single expression.+fromMultiplyList :: [Expression] -> Expression+fromMultiplyList [] = Number 1+fromMultiplyList [x] = x+fromMultiplyList (x : xs) = x :*: fromMultiplyList xs++-- | Convert a list of sub-expressions for an additive term into a single expression.+fromAddList :: [Expression] -> Expression+fromAddList [] = Number 0+fromAddList [x] = x+fromAddList (x : xs) = x :+: fromAddList xs++-- | Defines a total order among expressions.+-- In particular, higher integral powers of the variable are ordered later.+compareExpressions :: Text -> Expression -> Expression -> Ordering+compareExpressions v x y+  | (Just LT) <- compareDegree = LT+  | (Just GT) <- compareDegree = GT+  | LT <- comparePseudoDegree = LT+  | GT <- comparePseudoDegree = GT+  | LT <- compareSymbolCount = LT+  | GT <- compareSymbolCount = GT+  | LT <- compareOp = LT+  | GT <- compareOp = GT+  | Number n <- x, Number m <- y = compare n m+  | Symbol s <- x, Symbol r <- y = compare s r+  | UnaryApply _ x' <- x, UnaryApply _ y' <- y = compareExpressions v x' y'+  | BinaryApply _ x' x'' <- x,+    BinaryApply _ y' y'' <- y =+      case compareExpressions v x' y' of+        EQ -> compareExpressions v x'' y''+        c -> c+  | otherwise = EQ+  where+    compareDegree = do+      xd <- degree v x+      yd <- degree v y+      case (xd, yd) of+        (0, 0) -> return EQ+        (0, _) -> return LT+        (_, 0) -> return GT+        _ -> return $ compare xd yd+    comparePseudoDegree = compare (pseudoDegree v x) (pseudoDegree v y)+    compareSymbolCount = compare (symbolCount x) (symbolCount y)+    compareOp = compare (expressionOrder v x) (expressionOrder v y)++-- | The integral power of the variable for a particular expression.+degree :: Text -> Expression -> Maybe Integer+degree _ (Number _) = Just 0+degree v (Symbol s) | v == s = Just 1 | otherwise = Just 0+degree v (Negate' x) = degree v x+degree v (x :+: y) = max <$> degree v x <*> degree v y+degree v (x :-: y) = max <$> degree v x <*> degree v y+degree v (x :*: y) = (+) <$> degree v x <*> degree v y+degree v (x :/: y) = (-) <$> degree v x <*> degree v y+degree v (x :**: (Number n)) = (n *) <$> degree v x+degree v (x :**: Negate' y) = degree v $ x :**: y+degree _ _ = Nothing++-- | The number of times the variable appears in an expression.+-- Used as part of a somewhat arbitrary ordering.+pseudoDegree :: Text -> Expression -> Integer+pseudoDegree _ (Number _) = 0+pseudoDegree v (Symbol s) | v == s = 1 | otherwise = 0+pseudoDegree v (Negate' x) = pseudoDegree v x+pseudoDegree v (UnaryApply _ x) = pseudoDegree v x+pseudoDegree v (BinaryApply _ x y) = pseudoDegree v x + pseudoDegree v y++symbolCount :: Expression -> Int+symbolCount x = Set.size $ collect x+  where+    collect (Number _) = Set.empty+    collect (Symbol s) = Set.singleton s+    collect (UnaryApply _ u) = collect u+    collect (BinaryApply _ u v) = Set.union (collect u) (collect v)++-- | A fixed order between functions and operators.+-- Ignores the actual values the functins and operators are given.+expressionOrder :: Text -> Expression -> Int+expressionOrder _ (Number _) = 0+-- constant symbol has expressionOrder 1+expressionOrder _ (UnaryApply Negate _) = 2+expressionOrder _ (UnaryApply Signum _) = 3+expressionOrder _ (UnaryApply Abs _) = 4+expressionOrder _ (BinaryApply Add _ _) = 5+expressionOrder _ (BinaryApply Subtract _ _) = 6+expressionOrder _ (BinaryApply Multiply _ _) = 7+expressionOrder _ (BinaryApply Divide _ _) = 8+expressionOrder _ (BinaryApply Power _ _) = 9+expressionOrder _ (UnaryApply Sqrt _) = 10+expressionOrder _ (UnaryApply Exp _) = 11+expressionOrder _ (UnaryApply Log _) = 12+expressionOrder _ (BinaryApply LogBase _ _) = 13+expressionOrder _ (UnaryApply Sin _) = 14+expressionOrder _ (UnaryApply Cos _) = 15+expressionOrder _ (UnaryApply Tan _) = 16+expressionOrder _ (UnaryApply Asin _) = 17+expressionOrder _ (UnaryApply Acos _) = 18+expressionOrder _ (UnaryApply Atan _) = 19+expressionOrder _ (UnaryApply Sinh _) = 20+expressionOrder _ (UnaryApply Cosh _) = 21+expressionOrder _ (UnaryApply Tanh _) = 22+expressionOrder _ (UnaryApply Asinh _) = 23+expressionOrder _ (UnaryApply Acosh _) = 24+expressionOrder _ (UnaryApply Atanh _) = 25+expressionOrder v (Symbol s)+  | v == s = 26+  | otherwise = 1
+ src/Symtegration/Symbolic/Simplify/Fraction.hs view
@@ -0,0 +1,70 @@+-- |+-- Module: Symtegration.Symbolic.Simplify.Fraction+-- Description: Cancel out common numeric factors in fractions.+-- Copyright: Copyright 2025 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Symbolic.Simplify.Fraction (simplify) where++import Symtegration.Symbolic++-- $setup+-- >>> import Symtegration.Symbolic+-- >>> import Symtegration.Symbolic.Haskell++-- | Cancel out common numeric factors in fractions.+--+-- >>> toHaskell $ simplify $ 10 / 20+-- "1 / 2"+--+-- >>> toHaskell $ simplify $ Number (-15) / Number (-10)+-- "3 / 2"+--+-- >>> toHaskell $ simplify $ (12 * "x") / (4 * "y")+-- "(3 * x) / (1 * y)"+--+-- >> toHaskell $ simplify $ (15 * "x" + 20 * "y") / (5 * "z" - 35 * "u")+-- "(3 * x + 4 * y) / (1 * z - 7 * u)"+--+-- Assumes numeric folding and algebraic ring ordering has been applied.+simplify :: Expression -> Expression+simplify e@(_ :/: Number 0) = e+simplify (Number n :/: Number m)+  | m > 0 = Number (n `div` g) :/: Number (m `div` g)+  | otherwise = Number ((-n) `div` g) :/: Number ((-m) `div` g)+  where+    g = gcd n m+simplify (x :/: y) = divideFactor g x' :/: divideFactor g y'+  where+    g+      | (Number n) <- y, n < 0 = negate $ gcd (commonFactor x') n+      | otherwise = gcd (commonFactor x') (commonFactor y')+    x' = simplify x+    y' = simplify y+simplify ((1 :/: x) :*: y) = (1 :/: divideFactor g x') :*: divideFactor g y'+  where+    g = gcd (commonFactor x') (commonFactor y')+    x' = simplify x+    y' = simplify y+simplify (UnaryApply func x) = UnaryApply func $ simplify x+simplify (BinaryApply func x y) = BinaryApply func (simplify x) (simplify y)+simplify e = e++-- | Finds a common factor which multiplies each term in an expression.+-- Ignores terms not in algebraic ring ordering or includes direct negations of numbers.+commonFactor :: Expression -> Integer+commonFactor (Number n) = n+commonFactor (x :+: y) = gcd (commonFactor x) (commonFactor y)+commonFactor (x :-: y) = gcd (commonFactor x) (commonFactor y)+commonFactor (Number n :*: _) = n+commonFactor _ = 1++-- | Divides each term in an expression by a common factor.+-- Specialized for dividing factors found by 'commonFactor.+divideFactor :: Integer -> Expression -> Expression+divideFactor 0 e = e+divideFactor 1 e = e+divideFactor g (Number n) = Number $ n `div` g+divideFactor g (x :+: y) = divideFactor g x :+: divideFactor g y+divideFactor g (Number n :*: x) = Number (n `div` g) :*: x+divideFactor g e = e :/: Number g
+ src/Symtegration/Symbolic/Simplify/NumericFolding.hs view
@@ -0,0 +1,183 @@+-- |+-- Module: Symtegration.Symbolic.Simplify.NumericFolding+-- Description: Constant folding of numeric constants.+-- Copyright: Copyright 2025 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+--+-- This merges numeric terms as much as it can to simplify expressions.+-- Simplifications are finitely equivalent; i.e., any calculation with+-- finite inputs should result in the equivalent finite input.+-- The changes will also be exact, and no numeric constant will be replaced+-- by an approximate floating-point number.+module Symtegration.Symbolic.Simplify.NumericFolding (simplify) where++import Symtegration.Numeric (root)+import Symtegration.Symbolic++-- $setup+-- >>> import Symtegration.Symbolic.Haskell++-- | Simplifies computations involving numeric constants.+-- Basically, it computes as much as it can as long as any change is exact.+--+-- >>> toHaskell $ simplify $ 1 + 4+-- "5"+-- >>> toHaskell $ simplify $ 8 ** (1/3)+-- "2"+-- >>> toHaskell $ simplify $ 7 ** (1/3)+-- "7 ** (1 / 3)"+-- >>> toHaskell $ simplify $ 5 * 10 * "x"+-- "50 * x"+--+-- It will replace subtraction by addition and square roots by powers of \(\frac{1}{2}\).+simplify :: Expression -> Expression+simplify e@(Number _) = e+simplify e@(Symbol _) = e+simplify (UnaryApply func x) = unary $ UnaryApply func $ simplify x+simplify (BinaryApply func x y) = binary $ BinaryApply func (simplify x) (simplify y)++-- | Simplify computations involving numeric constants in unary expressions.+-- The arguments should already have been simplified.+unary :: Expression -> Expression+unary (Negate' (Number n)) = Number (-n)+unary (Negate' (Number n :/: Number m))+  | m < 0 = simplify $ Number n :/: Number (-m)+  | otherwise = simplify $ Number (-n) :/: Number m+unary (Abs' (Number n)) = Number $ abs n+unary (Signum' (Number n)) = Number $ signum n+unary (Exp' x) = simplifyExp x+unary (Log' x) = simplifyLog x+unary (Sqrt' x) = simplify $ x :**: (Number 1 :/: Number 2)+unary (Sin' x) = simplifySin x+unary (Cos' x) = simplifyCos x+unary (Tan' x) = simplifyTan x+unary e = e++-- | Simplify computations involving numeric constants in binary expressions.+-- The arguments should already have been simplified.+binary :: Expression -> Expression+-- Fold addition.+binary (Number 0 :+: x) = x+binary (x :+: Number 0) = x+binary (Number n :+: Number m) = Number (n + m)+binary ((Number n :/: Number m) :+: Number k) = reduceRatio (n + m * k) m+binary (Number n :+: (Number m :/: Number k)) = reduceRatio (n * k + m) k+binary ((Number n :/: Number m) :+: (Number k :/: Number l)) = reduceRatio (n * l + k * m) (m * l)+binary ((x :+: Number n) :+: Number m) = Number (n + m) :+: x+binary ((Number n :+: x) :+: Number m) = Number (n + m) :+: x+binary (Number n :+: (x :+: Number m)) = Number (n + m) :+: x+binary (Number n :+: (Number m :+: x)) = Number (n + m) :+: x+-- Fold multiplication.+binary (Number 0 :*: _) = Number 0+binary (_ :*: Number 0) = Number 0+binary (Number 1 :*: x) = x+binary (x :*: Number 1) = x+binary (Number n :*: Number m) = Number (n * m)+binary (Number n :*: (Number m :/: Number k)) = reduceRatio (n * m) k+binary ((Number n :/: Number m) :*: Number k) = reduceRatio (n * k) m+binary ((Number n :/: Number m) :*: (Number k :/: Number l)) = reduceRatio (n * k) (m * l)+binary ((x :*: Number n) :*: Number m) = Number (n * m) :*: x+binary ((Number n :*: x) :*: Number m) = Number (n * m) :*: x+binary (Number n :*: (x :*: Number m)) = Number (n * m) :*: x+binary (Number n :*: (Number m :*: x)) = Number (n * m) :*: x+binary e@(Number n :*: (x :/: Number m)) | m /= 0, m == n = x | otherwise = e+binary e@((x :/: Number n) :*: Number m) | n /= 0, m == n = x | otherwise = e+binary (x@(Number _) :*: (y@(Number _ :/: Number _) :*: z)) = simplify (x :*: y) :*: z+binary (x@(Number _ :/: Number _) :*: (y@(Number _ :/: Number _) :*: z)) = simplify (x :*: y) :*: z+-- Subtractions are turned into addition.+binary (x :-: y) = simplify $ x :+: Negate' y+-- Fold division.+binary e@(_ :/: (_ :/: 0)) = e+binary (x :/: (y :/: z)) = simplify $ (x :*: z) :/: y+binary e@((_ :/: 0) :/: _) = e+binary e@((_ :/: _) :/: 0) = e+binary ((x :/: y) :/: z) = simplify $ x :/: (y :*: z)+binary (Number n :/: Number m) = reduceRatio n m+-- Fold exponentiation.+binary e@(Number 0 :**: Number 0) = e+binary (Number _ :**: Number 0) = Number 1+binary (Number 1 :**: _) = Number 1+binary (Number n :**: Number m)+  | m >= 0 = Number (n ^ m)+  | otherwise = Number 1 :/: Number (n ^ (-m))+binary ((Number n :/: Number m) :**: Number k)+  | k >= 0 = Number (n ^ k) :/: Number (m ^ k)+  | otherwise = Number (m ^ (-k)) :/: Number (n ^ (-k))+binary e@(Number n :**: c@(Number m :/: Number k))+  | (Just l) <- root n k, m >= 0 = Number (l ^ m)+  | (Just l) <- root n k, m < 0 = 1 :/: Number (l ^ (-m))+  | n < 0, n /= -1, even k = (-1) ** c * simplify (Number (-n) ** c)+  | otherwise = e+binary e@((Number n :/: Number m) :**: (Number k :/: Number l))+  | (Just n', Just m') <- (root n l, root m l) = (Number n' :/: Number m') :**: Number k+  | otherwise = e+-- Turn LogBase into Log.+binary (LogBase' b x) = simplify $ Log' x :/: Log' b+binary e = e++-- | Simplify integer ratios.  Basically turns them into integers if possible,+-- and if not, reduce the fractions so that the denominator and numerator+-- do not have a common factor.+reduceRatio :: Integer -> Integer -> Expression+reduceRatio n 0 = Number n :/: Number 0+reduceRatio n 1 = Number n+reduceRatio n m+  | m == d = Number (n `div` m)+  | m == -d = Number (n `div` m)+  | n < 0, m < 0 = Number (-(n `div` d)) :/: Number (-(m `div` d))+  | otherwise = Number (n `div` d) :/: Number (m `div` d)+  where+    d = gcd n m++-- | Simplify an exponential of Euler's number.  I.e., simplify \(e^X\).+-- Only the exponent is given as an argument, while the return value is+-- the full simplified expression.+simplifyExp :: Expression -> Expression+simplifyExp (Number 0) = Number 1+simplifyExp (Log' x) = x+simplifyExp e = Exp' e++-- | Simplify a logarithm.  I.e., simplify \(log X\).+-- Only the parameter \(X\) is given as an argument, while the return value is+-- the full simplified expression.+simplifyLog :: Expression -> Expression+simplifyLog (Number 1) = Number 0+simplifyLog (Exp' x) = x+simplifyLog e = Log' e++-- | Simplify a sine.  I.e., simplify \(\sin X\).+-- Only the parameter \(X\) is given as an argument, while the return value is+-- the full simplified expression.+simplifySin :: Expression -> Expression+simplifySin (Number 0) = 0+simplifySin (Number _ :*: Pi') = 0+simplifySin (Pi' :*: Number _) = 0+simplifySin ((Number n :/: 2) :*: Pi')+  | even n = 0+  | odd ((n - 1) `div` 2) = 1+  | otherwise = -1+simplifySin (Pi' :*: (Number n :/: 2))+  | even n = 0+  | odd ((n - 1) `div` 2) = 1+  | otherwise = -1+simplifySin e = Sin' e++-- | Simplify a cosine.  I.e., simplify \(\cos X\).+-- Only the parameter \(X\) is given as an argument, while the return value is+-- the full simplified expression.+simplifyCos :: Expression -> Expression+simplifyCos (Number 0) = 1+simplifyCos (Number n :*: Pi') | even n = 1 | odd n = -1+simplifyCos (Pi' :*: Number n) | even n = 1 | odd n = -1+-- Any 2k/2 would have been simplified to k already.+simplifyCos ((Number _ :/: 2) :*: Pi') = 0+simplifyCos (Pi' :*: (Number _ :/: 2)) = 0+simplifyCos e = Cos' e++-- | Simplify a tangent.  I.e., simplify \(\tan X\).+-- Only the parameter \(X\) is given as an argument, while the return value is+-- the full simplified expression.+simplifyTan :: Expression -> Expression+simplifyTan (Number 0) = 0+simplifyTan e = Tan' e
+ src/Symtegration/Symbolic/Simplify/SymbolicFolding.hs view
@@ -0,0 +1,112 @@+-- |+-- Module: Symtegration.Symbolic.Simplify.SymbolicFolding+-- Description: Folding of symbolic terms.+-- Copyright: Copyright 2025 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+--+-- This merges symbolic terms as much as it can to simplify expressions.+-- Simplifications are finitely equivalent; i.e., any calculation with+-- finite inputs should result in the equivalent finite input.+module Symtegration.Symbolic.Simplify.SymbolicFolding (simplify) where++import Symtegration.Symbolic++-- | Folds symbolic terms as much as it can to simplify expressions.+--+-- Assumes algebraic ring ordering has been applied.+simplify :: Expression -> Expression+simplify e@(Number _) = e+simplify e@(Symbol _) = e+simplify (UnaryApply func x) = unary $ UnaryApply func $ simplify x+simplify (BinaryApply func x y) = binary $ BinaryApply func (simplify x) (simplify y)++-- | Folds symbolic terms for unary expressions.+--+-- The arguments should already have been simplified.+unary :: Expression -> Expression+unary (Negate' (Negate' x)) = x+unary (Negate' x) = (-1) * x+unary e = e++-- | Folds symbolic terms for binary expressions.+--+-- The arguments should already have been simplified.+binary :: Expression -> Expression+-- Fold addition.+binary e@(x :+: Negate' y)+  | x == y = Number 0+  | otherwise = e+binary e@(Negate' x :+: y)+  | x == y = Number 0+  | otherwise = e+binary (Number 0 :+: x) = x+binary (x :+: Number 0) = x+binary e@((Number n :*: x) :+: ((Number m :*: y) :+: z))+  | x == y = (Number (m + n) :*: x) :+: z+  | otherwise = e+binary e@((Number n :*: x) :+: (y :+: z))+  | x == y = (Number (n + 1) :*: x) :+: z+  | otherwise = e+binary e@(x :+: ((Number n :*: y) :+: z))+  | x == y = (Number (n + 1) :*: x) :+: z+  | otherwise = e+binary e@((Number n :*: x) :+: (Number m :*: y))+  | x == y = Number (n + m) :*: x+  | otherwise = e+binary e@(x :+: (Number n :*: y))+  | x == y = Number (n + 1) :*: x+  | otherwise = e+binary e@((Number n :*: x) :+: y)+  | x == y = Number (n + 1) :*: x+  | otherwise = e+binary e@(x :+: (y :+: z))+  | x == y = Number 2 :*: x :+: z+  | otherwise = e+binary e@(x :+: y)+  | x == y = Number 2 :*: x+  | otherwise = e+-- Fold multiplication.+binary (Number 0 :*: _) = Number 0+binary (_ :*: Number 0) = Number 0+binary (x :*: Number 1) = x+binary (Number 1 :*: x) = x+binary e@(x :*: (y :**: Number n))+  | x == y = x :**: Number (n + 1)+  | otherwise = e+binary e@((x :**: y) :*: (x' :**: y'))+  | x == x' = x :**: (y :+: y')+  | otherwise = e+binary e@(x :*: ((y :**: Number n) :*: z))+  | x == y = (x :**: Number (n + 1)) :*: z+  | otherwise = e+binary e@((x :**: Number n) :*: (y :*: z))+  | x == y = (x :**: Number (n + 1)) :*: z+  | otherwise = e+binary e@(x :*: (y :*: z))+  | x == y = (x :**: Number 2) :*: z+  | otherwise = e+binary e@(x :*: y)+  | x == y = x :**: Number 2+  | otherwise = e+-- Fold division.+binary (x :/: (y :/: z)) = (x :*: z) :/: y+binary ((x :/: y) :/: z) = x :/: (y :*: z)+binary (x :/: Number 1) = x+binary (x :/: Number (-1)) = (-1) * x+-- Fold powers.+binary (_ :**: Number 0) = Number 1+binary (1 :**: _) = Number 1+binary (x :**: Number 1) = x+binary (Negate' x :**: Number n)+  | even n = x :**: Number n+  | otherwise = Negate' (x :**: Number n)+binary ((Number (-1) :*: x) :**: Number n)+  | even n = x :**: Number n+  | otherwise = Number (-1) :*: (x :**: Number n)+binary ((x :**: y) :**: z) = x :**: (y :*: z)+-- Fold subtraction.+binary e@(x :-: y)+  | x == y = Number 0+  | otherwise = e+binary e = e
+ src/Symtegration/Symbolic/Simplify/Tidy.hs view
@@ -0,0 +1,72 @@+-- |+-- Module: Symtegtarion.Symbolic.Simplify.Tidy+-- Description: Tidy up a simplified mathematical expression.+-- Copyright: Copyright 2025 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Symbolic.Simplify.Tidy (tidy) where++import Symtegration.Symbolic++-- $setup+-- >>> import Symtegration.Symbolic+-- >>> import Symtegration.Symbolic.Haskell++-- | Tidies up expressions for nicer output.+--+-- Assumes that other simplifications have been applied first.+-- In fact, it may undo changes that made other simplifications easier.+--+-- ==== __What is tidied up__+--+-- This section shows examples of what this function tidies up.+--+-- >>> toHaskell $ tidy $ "x" + negate "y"+-- "x - y"+--+-- >>> toHaskell $ tidy $ "x" + Number (-2) * "y"+-- "x - 2 * y"+--+-- >>> toHaskell $ tidy $ Number (-1) / Number 2+-- "negate (1 / 2)"+--+-- >>> toHaskell $ tidy $ Number (-1) * "x"+-- "negate x"+--+-- >>> toHaskell $ tidy $ (-"x") * "y"+-- "negate (x * y)"+--+-- >>> toHaskell $ tidy $ "x" * (-"y")+-- "negate (x * y)"+--+-- >>> toHaskell $ tidy $ (-"x") * (-"y")+-- "x * y"+--+-- >>> toHaskell $ tidy $ "x" + ((-"y") + "z")+-- "x - y + z"+--+-- >>> toHaskell $ tidy $ "x" ** (1/2)+-- "sqrt x"+tidy :: Expression -> Expression+tidy (UnaryApply func x) = unary $ UnaryApply func $ tidy x+tidy (BinaryApply func x y) = binary $ BinaryApply func (tidy x) (tidy y)+tidy e = e++unary :: Expression -> Expression+unary e = e++binary :: Expression -> Expression+binary (x :+: (Negate' y)) = x :-: y+binary (Number (-1) :*: x) = Negate' x+binary e@(Number n :*: x)+  | n < 0 = Negate' (Number (-n) :*: x)+  | otherwise = e+binary e@(Number n :/: x)+  | n < 0 = Negate' (Number (-n) :/: x)+  | otherwise = e+binary (Negate' x :*: Negate' y) = x :*: y+binary (Negate' x :*: y) = Negate' $ x :*: y+binary (x :*: Negate' y) = Negate' $ x :*: y+binary (x :+: (Negate' y :+: z)) = (x :-: y) :+: z+binary (x :**: (Number 1 :/: Number 2)) = sqrt x+binary e = e
+ symtegration.cabal view
@@ -0,0 +1,156 @@+cabal-version: 1.12++-- This file has been generated from package.yaml by hpack version 0.37.0.+--+-- see: https://github.com/sol/hpack++name:           symtegration+version:        0.6.1+synopsis:       Library for symbolic integration of mathematical expressions.+description:    Symtegration is a library providing symbolic integration of mathematical expressions.+                .+                For example,+                .+                >>> import Symtegration+                >>> toHaskell <$> integrate "x" (4 * "x" ** 3 + 1)+                Just "x + x ** 4"+                .+                See the "Symtegration" module for the main interface.+category:       Mathematics, Symbolic Computation+homepage:       https://symtegration.dev/+bug-reports:    https://github.com/symtegration/symtegration/issues+author:         Yoo Chung+maintainer:     dev@chungyc.org+copyright:      Copyright 2025 Yoo Chung+license:        Apache-2.0+license-file:   LICENSE+build-type:     Simple+tested-with:+    GHC == 9.12.1 || == 9.10.1 || == 9.8.4 || == 9.6.6+extra-source-files:+    CHANGELOG.md+    LICENSE+    README.md+    docs/CODE_OF_CONDUCT.md+    docs/CONTRIBUTING.md+    docs/SECURITY.md++source-repository head+  type: git+  location: https://github.com/symtegration/symtegration++library+  exposed-modules:+      Symtegration+      Symtegration.Differentiation+      Symtegration.Integration+      Symtegration.Integration.Exponential+      Symtegration.Integration.Factor+      Symtegration.Integration.Parts+      Symtegration.Integration.Powers+      Symtegration.Integration.Rational+      Symtegration.Integration.Substitution+      Symtegration.Integration.Sum+      Symtegration.Integration.Term+      Symtegration.Integration.Trigonometric+      Symtegration.Numeric+      Symtegration.Polynomial+      Symtegration.Polynomial.Indexed+      Symtegration.Polynomial.Solve+      Symtegration.Polynomial.Symbolic+      Symtegration.Symbolic+      Symtegration.Symbolic.Haskell+      Symtegration.Symbolic.LaTeX+      Symtegration.Symbolic.Simplify+      Symtegration.Symbolic.Simplify.AlgebraicRingOrder+      Symtegration.Symbolic.Simplify.Fraction+      Symtegration.Symbolic.Simplify.NumericFolding+      Symtegration.Symbolic.Simplify.SymbolicFolding+      Symtegration.Symbolic.Simplify.Tidy+  other-modules:+      Paths_symtegration+  hs-source-dirs:+      src+  default-extensions:+      LambdaCase+      OverloadedStrings+  ghc-options: -Wall+  build-depends:+      ad ==4.5.*+    , base >=4.18 && <4.22+    , containers >=0.6 && <0.8+    , text >=2.0 && <2.2+    , text-show >=3.10 && <3.12+  default-language: GHC2021++test-suite examples+  type: exitcode-stdio-1.0+  main-is: test/Examples.hs+  other-modules:+      Paths_symtegration+  default-extensions:+      LambdaCase+      OverloadedStrings+  ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N+  build-depends:+      ad ==4.5.*+    , base >=4.18 && <4.22+    , containers >=0.6 && <0.8+    , doctest-parallel ==0.3.*+    , symtegration+    , text >=2.0 && <2.2+    , text-show >=3.10 && <3.12+  default-language: GHC2021++test-suite spec+  type: exitcode-stdio-1.0+  main-is: Spec.hs+  other-modules:+      Symtegration.ErrorDouble+      Symtegration.FiniteDouble+      Symtegration.Integration.ExponentialSpec+      Symtegration.Integration.FactorSpec+      Symtegration.Integration.PartsSpec+      Symtegration.Integration.PowersSpec+      Symtegration.Integration.Properties+      Symtegration.Integration.RationalSpec+      Symtegration.Integration.SubstitutionSpec+      Symtegration.Integration.SumSpec+      Symtegration.Integration.TermSpec+      Symtegration.Integration.TrigonometricSpec+      Symtegration.IntegrationSpec+      Symtegration.NumericSpec+      Symtegration.Polynomial.Indexed.Arbitrary+      Symtegration.Polynomial.IndexedSpec+      Symtegration.Polynomial.SolveSpec+      Symtegration.Polynomial.SymbolicSpec+      Symtegration.PolynomialSpec+      Symtegration.Symbolic.Arbitrary+      Symtegration.Symbolic.HaskellSpec+      Symtegration.Symbolic.LaTeXSpec+      Symtegration.Symbolic.Simplify.AlgebraicRingOrderSpec+      Symtegration.Symbolic.Simplify.FractionSpec+      Symtegration.Symbolic.Simplify.NumericFoldingSpec+      Symtegration.Symbolic.Simplify.Properties+      Symtegration.Symbolic.Simplify.SymbolicFoldingSpec+      Symtegration.Symbolic.Simplify.TidySpec+      Symtegration.SymbolicSpec+      Paths_symtegration+  hs-source-dirs:+      test+  default-extensions:+      LambdaCase+      OverloadedStrings+  ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N+  build-tool-depends:+      hspec-discover:hspec-discover ==2.11.*+  build-depends:+      QuickCheck >=2.14 && <2.16+    , ad ==4.5.*+    , base >=4.18 && <4.22+    , containers >=0.6 && <0.8+    , hspec ==2.11.*+    , symtegration+    , text >=2.0 && <2.2+    , text-show >=3.10 && <3.12+  default-language: GHC2021
+ test/Examples.hs view
@@ -0,0 +1,7 @@+module Main (main) where++import System.Environment (getArgs)+import Test.DocTest (mainFromCabal)++main :: IO ()+main = mainFromCabal "symtegration" =<< getArgs
+ test/Spec.hs view
@@ -0,0 +1,1 @@+{-# OPTIONS_GHC -F -pgmF hspec-discover #-}
+ test/Symtegration/ErrorDouble.hs view
@@ -0,0 +1,119 @@+-- |+-- Description: Floating-point numbers with error ranges.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+--+-- Floating-point numbers with error bars.+-- Basically each number is a pair of 'Double' values denoting a range.+-- These are used to determine whether an expression is too sensitive+-- to small divergences in floating-point computations.  By avoiding such+-- expressions, one can avoid situations where a mathematically equivalent+-- reformulation of a mathematical expression can end up with vastly different results.+module Symtegration.ErrorDouble+  ( DoubleWithError,+    sensitiveFunction,+    sensitiveExpression,+  )+where++import Data.Foldable1 qualified as Foldable1+import Data.List.NonEmpty (NonEmpty (..))+import Data.Text (Text)+import Symtegration.Symbolic++-- | Floating-point number with a range of values simulating floating-point divergences.+newtype DoubleWithError = DE (Double, Double) deriving (Eq, Ord, Show)++-- | Relative size of error to introduce to an individual 'Double' value.+errorSize :: Double+errorSize = 1e-5++-- | The amount of relative error we will tolerate.+errorTolerance :: Double+errorTolerance = 1e-3++-- | Add some error to a 'Double' value.+includeError :: Double -> DoubleWithError+includeError 0 = DE (-errorSize, errorSize)+includeError x+  | u <= v = DE (u, v)+  | otherwise = DE (v, u)+  where+    u = x * (1 - errorSize)+    v = x * (1 + errorSize)++-- | The relative size of error present in a 'DoubleWithError' value.+relativeError :: DoubleWithError -> Double+relativeError (DE (u, 0)) = abs u+relativeError (DE (0, v)) = abs v+relativeError (DE (u, v)) = abs (u - v) / (abs u + abs v)++-- | Returns whether the given function is sensitive at the given value.+-- I.e., given an error in the value, whether the error will grow too great in the result.+sensitiveFunction :: (DoubleWithError -> DoubleWithError) -> Double -> Bool+sensitiveFunction f x = not isNotSensitive+  where+    y = f $ includeError x++    -- We want to say it is too sensitive if either is NaN, so don't use >= directly.+    isNotSensitive = relativeError y < errorTolerance++-- | Returns whether the given expression is sensitive at the given assignment of values.+-- I.e., given an error in the values, whether the error will grow too great in the result.+sensitiveExpression :: Expression -> (Text -> Maybe Double) -> Bool+sensitiveExpression e m = not isNotSensitive+  where+    y = evaluate e t+    t s+      | (Just x') <- m s = Just $ includeError x'+      | otherwise = Just $ includeError 0++    -- We want to say it is too sensitive if either is NaN, so don't use >= directly.+    isNotSensitive+      | (Just y') <- y = relativeError y' < errorTolerance+      | otherwise = False++binOp :: (Double -> Double -> Double) -> DoubleWithError -> DoubleWithError -> DoubleWithError+binOp f (DE (u, v)) (DE (u', v')) = DE (Foldable1.minimum bounds, Foldable1.maximum bounds)+  where+    bounds = f u u' :| [f u v', f v u', f v v']++unOp :: (Double -> Double) -> DoubleWithError -> DoubleWithError+unOp f (DE (u, v)) = DE (min u' v', max u' v')+  where+    u' = f u+    v' = f v++instance Num DoubleWithError where+  (+) = binOp (+)+  (-) = binOp (-)+  (*) = binOp (*)+  negate = unOp negate+  abs = unOp abs+  signum = unOp signum+  fromInteger = includeError . fromInteger++instance Fractional DoubleWithError where+  (/) = binOp (/)+  recip = unOp recip+  fromRational = includeError . fromRational++instance Floating DoubleWithError where+  pi = includeError pi+  exp = unOp exp+  log = unOp log+  (**) = binOp (**)+  logBase = binOp logBase+  sin = unOp sin+  cos = unOp cos+  tan = unOp tan+  asin = unOp asin+  acos = unOp acos+  atan = unOp atan+  sinh = unOp sinh+  cosh = unOp cosh+  tanh = unOp tanh+  asinh = unOp asinh+  acosh = unOp acosh+  atanh = unOp atanh
+ test/Symtegration/FiniteDouble.hs view
@@ -0,0 +1,154 @@+-- |+-- Description: A variant of 'Double' without infinities or multiple zeroes.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+--+-- 'FiniteDouble' is a variant of 'Double' which avoids sensitivities+-- which result in what would otherwise be equivalent mathematical functions+-- result in significantly different results.  Basically, it ensures that+-- any finite value resulting from a calculation on finite values does+-- not involve any infinities during the calculation.+--+-- A value of NaN compares equal to any other NaN, which makes it possible+-- to check whether two supposedly equivalent functions both return NaN.+--+-- The functions operating on 'FiniteDouble' are not allowed to return+-- infinities, which prevents seemingly equivalent functions from returning+-- completely different results.  For the same reason, only a positive zero+-- is allowed to be returned.+--+-- These are examples of seemingly equivalent functions which can return+-- significantly different finite results, which 'FiniteDouble' prevents:+--+-- *   @(atan (m / (0 * z)))@ and @(atan (m / 0))@+--+-- *   @d / (0 - (cosh (exp ((logBase f c) / (e * 0)))))@ and+--     @d / (0 - (cosh (exp ((logBase f c) / 0))))@+--+-- *   @tanh (s ** (f / (0 * k)))@ and @tanh (s ** (f / 0))@+module Symtegration.FiniteDouble (FiniteDouble (..), isFinite, Exact (..), Near (..)) where++import Test.QuickCheck++-- | A variant of 'Double' which only allows finite+newtype FiniteDouble = FiniteDouble Double++instance Show FiniteDouble where+  show (FiniteDouble x) = show x++instance Eq FiniteDouble where+  (FiniteDouble x) == (FiniteDouble y)+    | isNaN x && isNaN y = True+    | otherwise = x == y++instance Ord FiniteDouble where+  (FiniteDouble x) <= (FiniteDouble y)+    | isNaN x && isNaN y = True+    | otherwise = x <= y++instance Num FiniteDouble where+  (+) = binOp (+)+  (*) = binOp (*)+  (-) = binOp (-)+  abs = unaryOp abs+  signum = unaryOp signum+  fromInteger n = FiniteDouble $ fromInteger n++instance Fractional FiniteDouble where+  (/) = binOp (/)+  fromRational q = FiniteDouble $ fromRational q++instance Floating FiniteDouble where+  pi = FiniteDouble pi+  exp = unaryOp exp+  log = unaryOp log+  sin = unaryOp sin+  cos = unaryOp cos+  asin = unaryOp asin+  acos = unaryOp acos+  atan = unaryOp atan+  sinh = unaryOp sinh+  cosh = unaryOp cosh+  asinh = unaryOp asinh+  acosh = unaryOp acosh+  atanh = unaryOp atanh+  (**) = binOp (**)++instance Real FiniteDouble where+  toRational (FiniteDouble x) = toRational x++instance RealFrac FiniteDouble where+  properFraction (FiniteDouble x) = (n, FiniteDouble f)+    where+      (n, f) = properFraction x++instance RealFloat FiniteDouble where+  floatRadix (FiniteDouble x) = floatRadix x+  floatDigits (FiniteDouble x) = floatDigits x+  floatRange (FiniteDouble x) = floatRange x+  decodeFloat (FiniteDouble x) = decodeFloat x+  encodeFloat x y = FiniteDouble $ encodeFloat x y+  isNaN (FiniteDouble x) = isNaN x || isInfinite x+  isInfinite (FiniteDouble x) = isInfinite x+  isDenormalized (FiniteDouble x) = isDenormalized x+  isNegativeZero (FiniteDouble x) = isNegativeZero x+  isIEEE (FiniteDouble x) = isIEEE x++instance Arbitrary FiniteDouble where+  arbitrary = FiniteDouble <$> arbitrary+  shrink (FiniteDouble x) = FiniteDouble <$> shrink x++-- | Returns whether a number is finite and not a NaN.+isFinite :: FiniteDouble -> Bool+isFinite (FiniteDouble x)+  | isNaN x = False+  | isInfinite x = False+  | otherwise = True++binOp :: (Double -> Double -> Double) -> FiniteDouble -> FiniteDouble -> FiniteDouble+binOp op (FiniteDouble x) (FiniteDouble y)+  | isInfinite v = FiniteDouble nan+  | -0 <- v = FiniteDouble 0+  | otherwise = FiniteDouble v+  where+    v = x `op` y+    nan = 0 / 0++unaryOp :: (Double -> Double) -> FiniteDouble -> FiniteDouble+unaryOp op (FiniteDouble x)+  | isInfinite v = FiniteDouble nan+  | -0 <- v = FiniteDouble 0+  | otherwise = FiniteDouble v+  where+    v = op x+    nan = 0 / 0++-- | Wrapper type over 'FiniteDouble' so that the same return values are compared as equal.+-- In other words, a NaN compared to a NaN will be considered equal.+-- Used for comparing that two implementations apply the exact same operations.+newtype Exact = Exact FiniteDouble deriving (Show)++instance Eq Exact where+  (Exact x) == (Exact y)+    | isNaN x && isNaN y = True+    | otherwise = x == y++-- | Wrapper type for comparing whether 'FiniteDouble' values are close enough.+-- Intended for testing whether two supposedly equivalent functions return+-- values which are close enough.+newtype Near = Near FiniteDouble deriving (Show)++instance Eq Near where+  (==) (Near (FiniteDouble x)) (Near (FiniteDouble y))+    | isNaN x && isNaN y = True+    | isInfinite x || isInfinite y = x == y+    | x == 0 || y == 0 || x == (-0) || y == (-0) = x - y < threshold+    | otherwise = (x - y) / y < threshold+    where+      threshold = 1e-3++instance Ord Near where+  compare (Near x'@(FiniteDouble x)) (Near y'@(FiniteDouble y))+    | x' == y' = EQ+    | otherwise = compare x y
+ test/Symtegration/Integration/ExponentialSpec.hs view
@@ -0,0 +1,54 @@+-- |+-- Description: Tests basic integration of exponential and logarithmic functions.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Integration.ExponentialSpec (spec) where++import Data.Map qualified as Map+import Data.Text (Text)+import Symtegration.Integration.Exponential+import Symtegration.Integration.Properties+import Symtegration.Symbolic+import Symtegration.Symbolic.Arbitrary+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck++spec :: Spec+spec = parallel $ do+  modifyMaxSuccess (* 10) $+    prop "consistent with derivative of integral" $ \(F e) x ->+      antiderivativeProperty integrate (Map.singleton var x) e x++  describe "ignores constants" $ do+    prop "with exponential" $+      forAll (arbitrarySymbol `suchThat` (/= Symbol var)) $ \c ->+        integrate var (exp c) `shouldBe` Nothing++    prop "with logarithm" $+      forAll (arbitrarySymbol `suchThat` (/= Symbol var)) $ \c ->+        integrate var (log c) `shouldBe` Nothing++    prop "with power of number" $ \n ->+      forAll (arbitrarySymbol `suchThat` (/= Symbol var)) $ \c ->+        integrate var (Number n ** c) `shouldBe` Nothing++    prop "with logarithm with base" $ \n ->+      forAll (arbitrarySymbol `suchThat` (/= Symbol var)) $ \c ->+        integrate var (logBase (Number n) c) `shouldBe` Nothing++newtype F = F Expression deriving (Eq, Show)++instance Arbitrary F where+  arbitrary =+    F+      <$> oneof+        [ pure $ Exp' (Symbol var),+          pure $ Log' (Symbol var),+          (:**:) <$> fmap Number arbitrarySizedNatural <*> pure (Symbol var),+          LogBase' <$> fmap Number arbitrarySizedNatural <*> pure (Symbol var)+        ]++var :: Text+var = "x"
+ test/Symtegration/Integration/FactorSpec.hs view
@@ -0,0 +1,70 @@+-- |+-- Description: Tests for Symtegration.Integration.Factor+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Integration.FactorSpec (spec) where++import Data.Text (Text, unpack)+import Symtegration.Integration.Factor+import Symtegration.Symbolic+import Symtegration.Symbolic.Arbitrary+import Symtegration.Symbolic.Haskell+import Symtegration.Symbolic.Simplify+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck++spec :: Spec+spec = parallel $ do+  describe "isConstant" $ do+    prop "for constant expression" $+      forAll genConstant $ \e ->+        isConstant var e `shouldBe` True++    prop "not for non-constant expression" $+      forAll genVariable $ \e ->+        isConstant var e `shouldBe` False++  describe "factor" $ do+    prop "into non-constant and constant factors" $+      forAll genVariable $ \e ->+        counterexample ("e = " <> unpack (toHaskell $ simplifyForVariable var e)) $+          factor var (simplifyForVariable var e)+            `shouldSatisfy` (\(x, y) -> isConstant var x && (not (isConstant var y) || y == Number 1))++    prop "variable portion has no multiplicative constant" $+      forAll genVariable $ \e ->+        counterexample ("e = " <> unpack (toHaskell $ simplifyForVariable var e)) $+          factor var (simplifyForVariable var e)+            `shouldSatisfy` (\(_, x) -> notConstantFactors x || x == Number 1)++notConstantFactors :: Expression -> Bool+notConstantFactors (x :*: y) = notConstantFactors x && notConstantFactors y+notConstantFactors x = not (isConstant var x)++genConstant :: Gen Expression+genConstant = sized $ \case+  0 -> oneof [arbitraryNumber, arbitrarySymbol `suchThat` (/= Symbol var)]+  n ->+    frequency+      [ (1, arbitraryNumber),+        (1, arbitrarySymbol `suchThat` (/= Symbol var)),+        (10, resize (max 0 (n - 1)) $ UnaryApply <$> arbitrary <*> genConstant),+        (10, resize (n `div` 2) $ BinaryApply <$> arbitrary <*> genConstant <*> genConstant)+      ]++genVariable :: Gen Expression+genVariable = sized $ \case+  0 -> pure (Symbol var)+  n ->+    frequency+      [ (1, pure (Symbol var)),+        (10, resize (max 0 (n - 1)) $ UnaryApply <$> arbitrary <*> genVariable),+        (5, resize (n `div` 2) $ BinaryApply <$> arbitrary <*> genVariable <*> genConstant),+        (5, resize (n `div` 2) $ BinaryApply <$> arbitrary <*> genConstant <*> genVariable),+        (5, resize (n `div` 2) $ BinaryApply <$> arbitrary <*> genVariable <*> genVariable)+      ]++var :: Text+var = "x"
+ test/Symtegration/Integration/PartsSpec.hs view
@@ -0,0 +1,41 @@+-- |+-- Description: Tests for Symtegration.Integration.Parts+-- Copyright: Copyright 2025 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Integration.PartsSpec (spec) where++import Data.Map qualified as Map+import Data.Text (Text)+import Symtegration.Integration.Parts+import Symtegration.Integration.Powers qualified as Powers+import Symtegration.Integration.Properties+import Symtegration.Integration.Term qualified as Term+import Symtegration.Symbolic+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck++spec :: Spec+spec = parallel $ do+  describe "integrates by parts" $ do+    prop "for powers" $+      forAll genParts $ \e x ->+        antiderivativeProperty+          (integrate [Powers.integrate, Term.integrate [Powers.integrate]])+          (Map.singleton var x)+          e+          x++-- | Generate an expression which can be integrated by parts.+genParts :: Gen Expression+genParts = do+  n <- arbitrarySizedNatural+  m <- arbitrarySizedNatural+  -- A product of two powers has some non-negligible chance to be integrated by parts.+  return $ x ** Number n * x ** Number m+  where+    x = Symbol var++var :: Text+var = "x"
+ test/Symtegration/Integration/PowersSpec.hs view
@@ -0,0 +1,54 @@+-- |+-- Description: Tests of Symtegration.Integration.Powers+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Integration.PowersSpec (spec) where++import Data.Map qualified as Map+import Data.Ratio (denominator, numerator)+import Data.Text (Text)+import Symtegration.Integration.Powers+import Symtegration.Integration.Properties+import Symtegration.Symbolic+import Symtegration.Symbolic.Arbitrary+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck++spec :: Spec+spec = parallel $ do+  prop "consistent with derivative of integral" $ \(Pow e) x ->+    antiderivativeProperty integrate (Map.singleton var x) e x++  prop "integrates constant symbol" $+    forAll (arbitrarySymbol `suchThat` (/= Symbol var)) $ \c ->+      integrate var c `shouldSatisfy` flip elem (map Just [Symbol var * c, c * Symbol var])++  describe "ignores constants" $ do+    prop "with integer power" $ \n ->+      forAll (arbitrarySymbol `suchThat` (/= Symbol var)) $ \c ->+        integrate var (c :*: Number n) `shouldBe` Nothing++    prop "with fraction" $ \n m ->+      forAll (arbitrarySymbol `suchThat` (/= Symbol var)) $ \c ->+        integrate var (c :*: (Number n :/: Number m)) `shouldBe` Nothing++newtype Pow = Pow Expression deriving (Eq, Show)++instance Arbitrary Pow where+  arbitrary =+    Pow+      <$> frequency+        [ (1, pure $ Symbol var :**: Number (-1)),+          (2, (\n -> Symbol var :**: Number n) <$> genExponent),+          (10, (\(m, n) -> Symbol var :**: (Number m :/: Number n)) <$> genFractionalExponent)+        ]+    where+      genExponent = resize 4 arbitrarySizedIntegral+      genFractionalExponent = resize 4 $ do+        q <- arbitrarySizedFractional+        return (numerator q, denominator q)++var :: Text+var = "x"
+ test/Symtegration/Integration/Properties.hs view
@@ -0,0 +1,55 @@+-- |+-- Description: Provides general properties that can be used to testing various integration algorithms.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Integration.Properties (antiderivativeProperty) where++import Data.Map (Map)+import Data.Map qualified as Map+import Data.Text (Text)+import Data.Text qualified as Text+import Numeric.AD+import Symtegration.ErrorDouble+import Symtegration.FiniteDouble+import Symtegration.Symbolic+import Symtegration.Symbolic.Arbitrary+import Symtegration.Symbolic.Haskell+import Test.Hspec+import Test.QuickCheck++-- | Tests the property that a function should be consistent+-- with the derivative of its integral.+antiderivativeProperty ::+  (Text -> Expression -> Maybe Expression) ->+  Map Text Double ->+  Expression ->+  Double ->+  Property+antiderivativeProperty integrate m e x =+  not (Map.null m) ==> forAll (elements $ Map.keys m) $ \v -> check (integrate v e) v+  where+    check Nothing _ = label "integration fail" True+    check (Just integrated) v =+      not (sensitiveExpression e (assign m)) && not (sensitiveExpression integrated (assign m)) ==>+        isFinite (FiniteDouble $ f x) && isFinite (FiniteDouble $ f' x) ==>+          label "integration success" $+            counterexample ("derivative = " <> Text.unpack (toHaskell e)) $+              counterexample ("antiderivative = " <> Text.unpack (toHaskell integrated)) $+                Near (FiniteDouble (f' x)) `shouldBe` Near (FiniteDouble (f x))+      where+        -- The original function and the derivative of the integral should behave similarly.+        --+        -- These are (Double -> Double).  It seems Numeric.AD does not like FiniteDouble.+        f = toFunction e (replace v)+        f' = diff (toFunction integrated (replaceForDiff v))++    -- Map all but the variable symbol to concrete numbers.+    replace var s+      | s == var = id+      | (Just z) <- Map.lookup s m = const z+      | otherwise = const 0+    replaceForDiff var s+      | s == var = id+      | (Just z) <- Map.lookup s m = const $ auto z+      | otherwise = const $ auto 0
+ test/Symtegration/Integration/RationalSpec.hs view
@@ -0,0 +1,68 @@+-- |+-- Description: Tests of Symtegration.Integration.Rational.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Integration.RationalSpec (spec) where++import Data.Map qualified as Map+import Data.Text (Text)+import Symtegration.Integration.Properties+import Symtegration.Integration.Rational+import Symtegration.Polynomial hiding (integrate)+import Symtegration.Polynomial.Indexed+import Symtegration.Polynomial.Indexed.Arbitrary ()+import Symtegration.Polynomial.Symbolic+import Symtegration.Symbolic+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck++spec :: Spec+spec = parallel $ do+  describe "integrate" $ modifyMaxSuccess (* 10) $ do+    prop "consistent with derivative of integral" $ \(Rat e) x ->+      antiderivativeProperty integrate (Map.singleton var x) e x++  describe "toRationalFunction" $ do+    prop "has coprime numerator and denominator" $ \(NonZero p) (NonZero q) ->+      let coprime (RationalFunction p' q') =+            degree (greatestCommonDivisor p' q') == 0+       in toRationalFunction p q `shouldSatisfy` coprime++  describe "hermiteReduce" $ do+    prop "h has squarefree denominator" $ \(NonZero p) (NonZero q) ->+      let r@(_, h) = hermiteReduce $ toRationalFunction p q+          RationalFunction _ d = h+       in counterexample (show r) $+            greatestCommonDivisor d (differentiate d) `shouldSatisfy` ((==) 0 . degree)++    prop "adds back to original rational function" $ \(NonZero p) (NonZero q) ->+      let f = toRationalFunction p q+          r@(gs, h) = hermiteReduce $ toRationalFunction p q++          -- Manually derive derivative of g = sum gs.+          RationalFunction x y = sum gs+          x' = y * differentiate x - x * differentiate y+          y' = y * y+          g' = toRationalFunction x' y'++          -- With leading coefficients factored out and numerator and denominator coprime,+          -- the representation of a rational function should be unique.+          rep (RationalFunction u v) =+            (leadingCoefficient u / leadingCoefficient v, RationalFunction (monic u) (monic v))+       in counterexample (show r) $ rep (g' + h) `shouldBe` rep f++-- | For generating arbitrary rational functions with rational number coefficients.+newtype Rat = Rat Expression deriving (Eq, Show)++instance Arbitrary Rat where+  arbitrary = resize 6 $ do+    p <- arbitrary :: Gen IndexedPolynomial+    q <- arbitrary `suchThat` (/= 0) :: Gen IndexedPolynomial+    let p' = toExpression var toRationalCoefficient p+    let q' = toExpression var toRationalCoefficient q+    return $ Rat $ p' / q'++var :: Text+var = "x"
+ test/Symtegration/Integration/SubstitutionSpec.hs view
@@ -0,0 +1,52 @@+-- |+-- Description: Tests for Symtegration.Integration.Substitution+-- Copyright: Copyright 2025 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Integration.SubstitutionSpec (spec) where++import Data.Map qualified as Map+import Data.Text (Text)+import Symtegration.Integration.Powers qualified as Powers+import Symtegration.Integration.Properties+import Symtegration.Integration.Substitution+import Symtegration.Integration.Trigonometric qualified as Trigonometric+import Symtegration.Symbolic+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck++spec :: Spec+spec = parallel $ do+  describe "integrates by substitution" $ do+    prop "powers and trigonometric functions mixed" $+      forAll genExpression $ \e x ->+        antiderivativeProperty+          (integrate [Powers.integrate, Trigonometric.integrate])+          (Map.singleton var x)+          e+          x++-- | Generate an expression which combines polynomials and trigonometric functions.+genExpression :: Gen Expression+genExpression = sized $ \case+  0 -> oneof leaves+  n ->+    frequency $+      [(1, g) | g <- leaves]+        ++ [ (1, resize (max 0 (n - 1)) $ Negate' <$> genExpression),+             (10, resize (max 0 (n - 1)) $ Sin' <$> genExpression),+             (10, resize (max 0 (n - 1)) $ Cos' <$> genExpression),+             (10, resize (max 0 (n - 1)) $ Tan' <$> genExpression),+             (10, resize (n `div` 2) $ (:+:) <$> genExpression <*> genExpression),+             (10, resize (n `div` 2) $ (:-:) <$> genExpression <*> genExpression)+           ]+  where+    leaves =+      [ Number <$> arbitrary,+        pure $ Symbol var,+        (:+:) (Symbol var) . Number <$> choose (2, 6)+      ]++var :: Text+var = "x"
+ test/Symtegration/Integration/SumSpec.hs view
@@ -0,0 +1,52 @@+-- |+-- Description: Tests for Symtegration.Integration.Sum+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Integration.SumSpec (spec) where++import Data.Map qualified as Map+import Data.Text (Text)+import Symtegration.Integration.Powers qualified as Powers+import Symtegration.Integration.Properties+import Symtegration.Integration.Sum+import Symtegration.Integration.Trigonometric qualified as Trigonometric+import Symtegration.Symbolic+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck++spec :: Spec+spec = parallel $ do+  describe "integrates and adds by term" $ do+    prop "powers and trigonometric functions mixed" $+      forAll genExpression $ \e x ->+        antiderivativeProperty+          (integrate [Powers.integrate, Trigonometric.integrate])+          (Map.singleton var x)+          e+          x++-- | Generate an expression which adds polynomials and trigonometric functions together.+genExpression :: Gen Expression+genExpression = sized $ \case+  0 -> oneof leaves+  n ->+    frequency $+      [(1, g) | g <- leaves]+        ++ [ (1, resize (max 0 (n - 1)) $ Negate' <$> genExpression),+             (10, resize (n `div` 2) $ (:+:) <$> genExpression <*> genExpression),+             (10, resize (n `div` 2) $ (:-:) <$> genExpression <*> genExpression)+           ]+  where+    leaves =+      [ Number <$> arbitrary,+        pure $ Symbol var,+        pure $ Sin' $ Symbol var,+        pure $ Cos' $ Symbol var,+        pure $ Tan' $ Symbol var,+        (:+:) (Symbol var) . Number <$> choose (2, 6)+      ]++var :: Text+var = "x"
+ test/Symtegration/Integration/TermSpec.hs view
@@ -0,0 +1,59 @@+-- |+-- Description: Tests for Symtegration.Integration.Sum+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Integration.TermSpec (spec) where++import Data.Map qualified as Map+import Data.Text (Text, unpack)+import Symtegration.Integration.Powers qualified as Powers+import Symtegration.Integration.Properties+import Symtegration.Integration.Term+import Symtegration.Integration.Trigonometric qualified as Trigonometric+import Symtegration.Symbolic+import Symtegration.Symbolic.Arbitrary+import Symtegration.Symbolic.Haskell+import Symtegration.Symbolic.Simplify+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck++spec :: Spec+spec = parallel $ do+  describe "integrates term" $ do+    prop "for constant multiplied by simple term" $ \x ->+      forAll genConstant $ \c ->+        forAll genVariableTerm $ \e ->+          let e' = simplifyForVariable var $ c :*: e+              fs = [Powers.integrate, Trigonometric.integrate]+           in counterexample ("e' = " <> unpack (toHaskell e')) $+                antiderivativeProperty (integrate fs) (Map.singleton var x) e' x++-- | Expression with no variable.+genConstant :: Gen Expression+genConstant = sized $ \case+  0 -> Number <$> arbitrarySizedNatural+  n ->+    frequency+      [ (1, Number <$> arbitrarySizedNatural),+        (10, resize (max 0 (n - 1)) $ Exp' <$> genConstant),+        (10, resize (max 0 (n - 1)) $ Negate' <$> genConstant),+        (10, resize (n `div` 2) $ (:+:) <$> genConstant <*> genConstant),+        (10, resize (n `div` 2) $ (:*:) <$> genConstant <*> genConstant)+      ]++-- | Variable terms that the basic integration algorithms can integrate.+genVariableTerm :: Gen Expression+genVariableTerm =+  oneof+    [ pure $ Number 1,+      pure $ Symbol var,+      (:**:) (Symbol var) <$> arbitraryNumber,+      pure $ Sin' (Symbol var),+      pure $ Cos' (Symbol var),+      pure $ Tan' (Symbol var)+    ]++var :: Text+var = "x"
+ test/Symtegration/Integration/TrigonometricSpec.hs view
@@ -0,0 +1,49 @@+-- |+-- Description: Tests basic integration of trigonometric functions.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Integration.TrigonometricSpec (spec) where++import Data.Map qualified as Map+import Data.Text (Text)+import Symtegration.Integration.Properties+import Symtegration.Integration.Trigonometric+import Symtegration.Symbolic+import Symtegration.Symbolic.Arbitrary+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck++spec :: Spec+spec = parallel $ do+  modifyMaxSuccess (* 10) $+    prop "consistent with derivative of integral" $ \(Trig e) x ->+      antiderivativeProperty integrate (Map.singleton var x) e x++  prop "ignores constant symbols" $ \(Trig e) ->+    forAll (arbitrarySymbol `suchThat` (/= Symbol var)) $ \c ->+      integrate var (substitute e (\x -> if x == var then Just c else Nothing)) `shouldBe` Nothing++newtype Trig = Trig Expression deriving (Eq, Show)++instance Arbitrary Trig where+  arbitrary = Trig <$> elements [f (Symbol var) | f <- candidates]+    where+      candidates =+        [ Sin',+          Cos',+          Tan',+          Asin',+          Acos',+          Atan',+          Sinh',+          Cosh',+          Tanh',+          Asinh',+          Acosh',+          Atanh'+        ]++var :: Text+var = "x"
+ test/Symtegration/IntegrationSpec.hs view
@@ -0,0 +1,52 @@+-- |+-- Description: General testing of specific integration algorithms with numeric coefficients.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.IntegrationSpec (spec) where++import Data.Map qualified as Map+import Data.Text (Text)+import Symtegration.FiniteDouble+import Symtegration.Integration+import Symtegration.Integration.Exponential qualified as Exponential+import Symtegration.Integration.Powers qualified as Powers+import Symtegration.Integration.Properties qualified as Properties+import Symtegration.Integration.Rational qualified as Rational+import Symtegration.Integration.Trigonometric qualified as Trigonometric+import Symtegration.Symbolic+import Symtegration.Symbolic.Arbitrary+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck++spec :: Spec+spec = parallel $ do+  -- Each integration algorithm should have their own tests,+  -- where they focus the input expressions which are generated.+  -- These tests are for checking whether they could have problems+  -- with expressions they do not focus on.+  modifyMaxSuccess (* 10) $ context "for any expression" $ do+    describe "integral consistent with derivative" $ do+      prop "for integration of powers" $+        antiderivativeProperty Powers.integrate++      prop "for trigonometric integration" $+        antiderivativeProperty Trigonometric.integrate++      prop "for integration of exponential and logarithmic functions" $+        antiderivativeProperty Exponential.integrate++      prop "for rational functions" $+        antiderivativeProperty Rational.integrate++      prop "for general integration" $+        antiderivativeProperty integrate++antiderivativeProperty ::+  (Text -> Expression -> Maybe Expression) ->+  Complete ->+  Double ->+  Property+antiderivativeProperty f (Complete e m) =+  Properties.antiderivativeProperty f (Map.map (\(FiniteDouble z) -> z) m) e
+ test/Symtegration/NumericSpec.hs view
@@ -0,0 +1,33 @@+-- |+-- Description: Tests for Symtegration.Numeric.+-- Copyright: Copyright 2025 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.NumericSpec (spec) where++import Symtegration.Numeric+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck++spec :: Spec+spec = parallel $ do+  describe "root" $ do+    prop "finds root of zero" $ \(Positive e) -> root 0 e `shouldBe` Just 0++    prop "finds root of one" $ \(Positive e) -> root 1 e `shouldBe` Just 1++    prop "finds positive root" $ \(Positive x) (Positive e) ->+      root (x ^ e) e `shouldBe` Just x++    prop "finds negative root for odd power" $ \(Negative x) (Positive e) ->+      odd e ==>+        root (x ^ e) e `shouldBe` Just x++    prop "finds nothing for negative number with even power" $ \(Negative x) (Positive e) ->+      even e ==>+        root x e `shouldBe` Nothing++    prop "finds nothing when there is no integer root" $ \(Positive x) (Positive e) ->+      even e ==>+        root (x ^ e + 1) e `shouldBe` Nothing
+ test/Symtegration/Polynomial/Indexed/Arbitrary.hs view
@@ -0,0 +1,63 @@+{-# OPTIONS_GHC -fno-warn-orphans #-}++-- |+-- Description: Generate arbitrary instances of 'IndexedPolynomial'.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Polynomial.Indexed.Arbitrary where++import Symtegration.Polynomial+import Symtegration.Polynomial.Indexed+import Symtegration.Symbolic.Arbitrary ()+import Test.QuickCheck hiding (scale)++instance Arbitrary IndexedPolynomial where+  arbitrary = sized $ \case+    0 ->+      frequency+        [ (50, pure (power 1)),+          (10, scale <$> resize 3 arbitrary `suchThat` (/= 0) <*> pure 1),+          (1, pure 0)+        ]+    n ->+      frequency+        [ (1, resize 0 arbitrary),+          (10, resize (n `div` 2) $ (+) <$> arbitrary <*> arbitrary),+          (10, resize (n `div` 2) $ (*) <$> arbitrary `suchThat` (/= 0) <*> arbitrary `suchThat` (/= 0))+        ]++  shrink p+    | 0 <- degree p = []+    | otherwise = [p - scale c (power k) | k <- [0 .. degree p], let c = coefficient p k, c /= 0]++instance+  (Polynomial p e c, Arbitrary (p e c), Eq (p e c), Num (p e c), Eq c) =>+  Arbitrary (IndexedPolynomialWith (p e c))+  where+  arbitrary = sized $ \case+    0 -> frequency [(10, pure (power 1)), (1, scale <$> resize 4 arbitrary <*> pure 1)]+    n ->+      frequency+        [ (1, resize 0 arbitrary),+          (10, resize (n `div` 2) $ (+) <$> arbitrary <*> arbitrary),+          (10, resize (n `div` 2) $ (*) <$> arbitrary <*> arbitrary)+        ]++  shrink p+    | 0 <- degree p = []+    | otherwise = [p - scale c (power k) | k <- [0 .. degree p], let c = coefficient p k, c /= 0]++instance Arbitrary IndexedSymbolicPolynomial where+  arbitrary = sized $ \case+    0 -> frequency [(10, pure (power 1)), (1, scale <$> arbitrary <*> pure 1)]+    n ->+      frequency+        [ (1, frequency [(10, pure (power 1)), (1, scale <$> arbitrary <*> pure 1)]),+          (10, resize (n `div` 2) $ (+) <$> arbitrary <*> arbitrary),+          (10, resize (n `div` 2) $ (*) <$> arbitrary <*> arbitrary)+        ]++  shrink p+    | 0 <- degree p = []+    | otherwise = [p - scale c (power k) | k <- [0 .. degree p], let c = coefficient p k, c /= 0]
+ test/Symtegration/Polynomial/IndexedSpec.hs view
@@ -0,0 +1,102 @@+-- |+-- Description: Tests Symtegration.Poynomial.Indexed.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Polynomial.IndexedSpec (spec) where++import Data.List (dropWhileEnd)+import Symtegration.Polynomial+import Symtegration.Polynomial.Indexed+import Symtegration.Polynomial.Indexed.Arbitrary ()+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck hiding (scale)++spec :: Spec+spec = parallel $ describe "IndexedPolynomial" $ do+  prop "fromInteger" $ \n ->+    let p = fromInteger n :: IndexedPolynomial+     in conjoin+          [ degree p `shouldBe` 0,+            coefficient p 0 `shouldBe` fromInteger n,+            leadingCoefficient p `shouldBe` fromInteger n+          ]++  prop "from power" $ \n ->+    let p = power n :: IndexedPolynomial+     in conjoin+          [ degree p `shouldBe` n,+            coefficient p n `shouldBe` 1,+            leadingCoefficient p `shouldBe` 1,+            [coefficient p k | k <- [0 .. degree p - 1]] `shouldSatisfy` all (== 0)+          ]++  prop "from series of powers" $ \cs ->+    let cs' = dropWhileEnd (== 0) cs+        p = foldl accumulate 0 (zip [0 ..] cs') :: IndexedPolynomial+        accumulate p' (e, c) = p' + scale c (power e)+     in not (null cs') ==> getCoefficients p `shouldBe` cs'++  prop "leading coefficients match" $ \p ->+    leadingCoefficient p `shouldBe` coefficient (p :: IndexedPolynomial) (degree p)++  describe "addition" $ do+    prop "adds numbers" $ \m n ->+      fromInteger m + fromInteger n `shouldBe` (fromInteger (m + n) :: IndexedPolynomial)++    prop "adds number and polynomial" $ \m p c ->+      let leadingTerm = scale c (power $ 1 + degree p)+          p' = p + leadingTerm :: IndexedPolynomial+       in fromInteger m + p' `shouldBe` (fromInteger m + p) + leadingTerm++    prop "adds polynomials" $ \p q c ->+      let leadingTerm = scale c (power $ 1 + degree p)+          p' = p + leadingTerm :: IndexedPolynomial+       in p' + q `shouldBe` (p + q) + leadingTerm++  describe "multiplication" $ do+    prop "multiplies numbers" $ \m n ->+      fromInteger m * fromInteger n `shouldBe` (fromInteger (m * n) :: IndexedPolynomial)++    prop "multiplies number and polynomial" $ \m p c ->+      let leadingTerm = scale c (power $ 1 + degree p)+          p' = p + leadingTerm :: IndexedPolynomial+       in fromInteger m * p' `shouldBe` (fromInteger m * p) + fromInteger m * leadingTerm++    prop "multiplies polynomials" $ \p q c ->+      let leadingTerm = scale c (power $ 1 + degree p)+          p' = p + leadingTerm :: IndexedPolynomial+       in p' * q `shouldBe` (p * q) + (leadingTerm * q)++  describe "subtraction" $ do+    prop "is same as adding negation" $ \p q ->+      let q' = negate q :: IndexedPolynomial+       in p - q `shouldBe` p + q'++  describe "negate" $ do+    prop "negates coefficients" $ \p ->+      getCoefficients (negate p) `shouldBe` map negate (getCoefficients p)++  describe "signum" $ do+    it "is zero for zero" $ do+      signum (0 :: IndexedPolynomial) `shouldBe` 0++    prop "is either one or negative one" $ \(NonZero p) ->+      signum (p :: IndexedPolynomial) `shouldSatisfy` (\x -> x == 1 || x == -1)++    prop "is consistent with abs" $ \p ->+      abs p * signum (p :: IndexedPolynomial) `shouldBe` p++  describe "show" $ do+    prop "is total for IndexedPolynomial" $ \p -> total (show (p :: IndexedPolynomial))++    prop "is total for IndexedSymbolicPolynomial" $ \p ->+      total (show (p :: IndexedSymbolicPolynomial))++    prop "is total for IndexedPolynomialWith IndexedPolynomial" $ \p ->+      total (show (p :: IndexedPolynomialWith IndexedPolynomial))++-- | Returns the coefficients of the given polynomial, in ascending order of the power.+getCoefficients :: IndexedPolynomial -> [Rational]+getCoefficients p = [coefficient p k | k <- [0 .. degree p]]
+ test/Symtegration/Polynomial/SolveSpec.hs view
@@ -0,0 +1,243 @@+-- |+-- Description: Tests Symtegration.Poynomial.Solve.+-- Copyright: Copyright 2025 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Polynomial.SolveSpec (spec) where++import Data.Complex+import Data.List (nub, sort)+import Data.Monoid (Sum (..))+import Symtegration.FiniteDouble+import Symtegration.Polynomial+import Symtegration.Polynomial.Indexed+import Symtegration.Polynomial.Solve+import Symtegration.Polynomial.Symbolic+import Symtegration.Symbolic+import Symtegration.Symbolic.Haskell+import Symtegration.Symbolic.Simplify+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck hiding (scale)++spec :: Spec+spec = parallel $ do+  describe "solve" $ do+    describe "linear polynomials" $ do+      prop "found roots are roots" $ \(NonZero a) b ->+        let p = scale a (power 1) + scale b (power 0)+         in correctlySolves p++      prop "finds root" $ \(NonZero a) x ->+        let p = scale a 1 * (power 1 - scale x 1)+         in counterexample (show p) $+              solve p `shouldBe` Just [fromRational x]++    describe "quadratic polynomials" $ do+      prop "found roots are roots" $ \(NonZero a) b c ->+        let p = scale a (power 2) + scale b (power 1) + scale c (power 0)+         in correctlySolves p++      prop "finds all roots" $ \(NonZero a) x y ->+        let p = scale a 1 * (power 1 - scale x 1) * (power 1 - scale y 1)+         in counterexample (show p) $+              if x == y+                then toFiniteDoubleRoots (solve p) `shouldBe` Just (toFiniteDoubles [x])+                else toFiniteDoubleRoots (solve p) `shouldBe` Just (toFiniteDoubles [x, y])++      prop "does not find real roots" $ \(NonZero a) b c ->+        let p = scale a (power 2) + scale b (power 1) + scale c 1+            sq = b * b - 4 * a * c+         in sq < 0 ==> solve p `shouldBe` Just []++    describe "cubic polynomials" $ do+      modifyMaxSuccess (* 10) $+        prop "found roots are roots" $ \(NonZero a) b c d ->+          let p = scale a (power 3) + scale b (power 2) + scale c (power 1) + scale d 1+           in correctlySolves p++      prop "with zero lower order terms" $ \(NonZero a) ->+        let p = scale a (power 3)+         in correctlySolves p++      prop "with zero discriminant" $ \u ->+        let p = -(3 * u * u)+            q = 2 * u * u * u+            r = power 3 + scale p (power 1) + scale q 1+         in conjoin+              [ counterexample ("p = " <> show p <> ", q = " <> show q) $+                  4 * p * p * p + 27 * q * q === 0,+                correctlySolves r+              ]++      modifyMaxSuccess (* 100) $+        prop "finds roots" $ \(NonZero a) x y z ->+          let p = scale a 1 * (power 1 - scale x 1) * (power 1 - scale y 1) * (power 1 - scale z 1)+              roots = nub [x, y, z]+           in counterexample (show p) $+                toFiniteDoubleRoots (solve p) `shouldBe` Just (toFiniteDoubles roots)++    describe "quartic polynomials" $ do+      modifyMaxSuccess (* 10) $+        prop "found roots are roots" $ \(NonZero a) b c d e ->+          let p = scale a (power 4) + scale b (power 3) + scale c (power 2) + scale d (power 1) + scale e 1+           in correctlySolves p++      describe "special cases" $ do+        prop "ax^4 + bx^3 = 0" $ \(NonZero a) b ->+          let p = scale a (power 4) + scale b (power 3)+           in correctlySolves p++        prop "ax^4 + bx^3 + cx^2 = 0" $ \(NonZero a) b c ->+          let p = scale a (power 4) + scale b (power 3) + scale c (power 2)+           in correctlySolves p++        prop "ax^4 + bx^3 + cx^2 + dx = 0" $ \(NonZero a) b c d ->+          let p = scale a (power 4) + scale b (power 3) + scale c (power 2) + scale d (power 1)+           in correctlySolves p++        prop "ax^4 + b = 0" $ \(NonZero a) b ->+          let p = scale a (power 4) + scale b 1+           in correctlySolves p++        modifyMaxSuccess (* 10) $+          prop "ax^4 + bx^2 + c = 0" $ \(NonZero a) b c ->+            let p = scale a (power 4) + scale b (power 2) + scale c 1+             in correctlySolves p++      modifyMaxSuccess (* 1000) $+        prop "finds all real roots when any found" $ \(NonZero a) x y z w ->+          let p = scale a $ product [power 1 - scale v 1 | v <- [x, y, z, w]]+              roots = nub [x, y, z, w]+           in counterexample (show p) $+                case solve p of+                  Nothing -> label "not solved" True+                  xs@(Just _) ->+                    label "solved" $+                      toFiniteDoubleRoots xs `shouldBe` Just (toFiniteDoubles roots)++  describe "complexSolve" $ do+    describe "linear polynomials" $ do+      prop "finds root" $ \(NonZero a) x ->+        let p = scale a 1 * (power 1 - scale x 1)+         in counterexample (show p) $+              complexSolve p `shouldBe` Just [fromRational x]++    describe "quadratic polynomials" $ do+      prop "finds real solutions" $ \(NonZero a) b c ->+        let p = scale a (power 2) + scale b (power 1) + scale c 1+         in filter (/= Near (0 / 0)) <$> toFiniteDoubleRoots (complexSolve p)+              `shouldBe` toFiniteDoubleRoots (solve p)++    describe "cubic polynomials" $ do+      modifyMaxSuccess (* 100) $+        prop "finds roots" $ \(NonZero a) b c d ->+          let p = scale a (power 3) + scale b (power 2) + scale c (power 1) + scale d 1+           in consistentWithComplexRoots p (complexSolve p)++      describe "special cases" $ do+        prop "ax^3 = 0" $ \(NonZero a) ->+          let p = scale a (power 3)+           in consistentWithComplexRoots p (complexSolve p)++        prop "ax^3 + bx^2 = 0" $ \(NonZero a) (NonZero b) ->+          let p = scale a (power 3) + scale b (power 2)+           in consistentWithComplexRoots p (complexSolve p)++        prop "ax^3 + bx^2 + cx= 0" $ \(NonZero a) (NonZero b) (NonZero c) ->+          let p = scale a (power 3) + scale b (power 2) + scale c (power 1)+           in consistentWithComplexRoots p (complexSolve p)++    describe "quartic polynomials" $ do+      modifyMaxSuccess (* 100) $+        prop "finds roots" $ \(NonZero a) b c d e ->+          let p = scale a (power 4) + scale b (power 3) + scale c (power 2) + scale d (power 1) + scale e 1+           in consistentWithComplexRoots p (complexSolve p)++      describe "special cases" $ do+        prop "ax^4 = 0" $ \(NonZero a) ->+          let p = scale a (power 4)+           in consistentWithComplexRoots p (complexSolve p)++        prop "ax^4 + bx^3 = 0" $ \(NonZero a) (NonZero b) ->+          let p = scale a (power 4) + scale b (power 3)+           in consistentWithComplexRoots p (complexSolve p)++        prop "ax^4 + bx^3 + cx^2 = 0" $ \(NonZero a) (NonZero b) (NonZero c) ->+          let p = scale a (power 4) + scale b (power 3) + scale c (power 2)+           in consistentWithComplexRoots p (complexSolve p)++        prop "ax^4 + bx^3 + cx^2 + dx = 0" $ \(NonZero a) (NonZero b) (NonZero c) (NonZero d) ->+          let p = scale a (power 4) + scale b (power 3) + scale c (power 2) + scale d (power 1)+           in consistentWithComplexRoots p (complexSolve p)++        prop "ax^4 + bx^2 + c = 0" $ \(NonZero a) b c ->+          let p = scale a (power 4) + scale b (power 2) + scale c (power 0)+           in consistentWithComplexRoots p (complexSolve p)++-- | Passes if either all the roots found are indeed roots of the polynomial+-- or solutions could not be derived.+correctlySolves :: IndexedPolynomial -> Property+correctlySolves p =+  counterexample (show p) $+    counterexample (show $ map (toHaskell . simplify) <$> roots) $+      label (case roots of Nothing -> "not solved"; Just _ -> "solved") $+        roots `shouldSatisfy` areRoots p+  where+    roots = solve p++-- | Whether x is a root of p.+isRoot :: IndexedPolynomial -> Expression -> Bool+isRoot p x+  | (Just x') <- evaluate x (const Nothing) = Near (f x') == Near 0+  | otherwise = False+  where+    p' = toExpression "x" toRationalCoefficient p+    f = toFunction p' (\case "x" -> id; _ -> undefined)++-- | Whether the given roots are indeed roots of the given polynomial,+-- or if roots could not be found.+areRoots :: IndexedPolynomial -> Maybe [Expression] -> Bool+areRoots _ Nothing = True+areRoots p (Just xs) = all (isRoot p) xs++-- | Evaluate an expression into a floating-point value for comparisons.+eval :: Expression -> Near+eval e+  | (Just x) <- evaluate e (const Nothing) = Near x+  | otherwise = Near $ 0 / 0 -- not a number++-- | Convert a potential list of polynomial root solutions into floating-point values for comparisons.+toFiniteDoubleRoots :: Maybe [Expression] -> Maybe [Near]+toFiniteDoubleRoots = fmap (sort . map eval)++-- | Convert a list of rational numbers into floating-point values for comparisons.+toFiniteDoubles :: [Rational] -> [Near]+toFiniteDoubles = sort . map (Near . fromRational)++-- | Evaluate an expression to a concrete complex number.+complexEval :: Expression -> Complex Double+complexEval expr+  | (Just x) <- evaluate expr (const Nothing) = x+  | otherwise = 0 / 0 -- not a number++-- | Evaluate a polynomial with a complex number substituted in the variable.+complexPolyEval :: IndexedPolynomial -> Complex Double -> Complex Double+complexPolyEval p x = getSum $ foldTerms (\e c -> Sum $ fromRational c * x ** fromIntegral e) p++-- | Check whether the given polynomial is consistent with the given solutions,+-- which may include complex numbers.+consistentWithComplexRoots :: IndexedPolynomial -> Maybe [Expression] -> Property+consistentWithComplexRoots p roots =+  label (rootsLabel roots) $+    counterexample (show roots') $+      map (complexPolyEval p) <$> roots' `shouldSatisfy` closeEnough+  where+    roots' = map complexEval <$> roots++    rootsLabel Nothing = "did not solve"+    rootsLabel (Just xs) = "root count = " <> show (length xs)++    -- We do not check for sensitive functions, so use a generous error bound.+    closeEnough (Just xs) = all ((< 1) . magnitude) xs+    closeEnough Nothing = True
+ test/Symtegration/Polynomial/SymbolicSpec.hs view
@@ -0,0 +1,67 @@+-- |+-- Description: Tests for Symtegration.Polynomial.Symbolic+-- Copyright: Copyright 2025 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Polynomial.SymbolicSpec (spec) where++import Symtegration.Polynomial+import Symtegration.Polynomial.Indexed+import Symtegration.Polynomial.Indexed.Arbitrary ()+import Symtegration.Polynomial.Symbolic+import Symtegration.Symbolic+import Symtegration.Symbolic.Arbitrary+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck++spec :: Spec+spec = parallel $ do+  describe "fromExpression" $ do+    describe "with rational number coefficients" $ do+      prop "is inverse of toExpression" $ \p (SymbolText s) ->+        let e = toExpression s toRationalCoefficient (p :: IndexedPolynomial)+            p' = fromExpression (forVariable s) e+         in counterexample ("p = " <> show p) $+              counterexample ("p'" <> show p') $+                -- With exact rational coefficients, the polynomial representation of+                -- a particular polynomial is unique.+                p' `shouldBe` Just p++      prop "from number" $ \(SymbolText s) n ->+        fromExpression (forVariable s) (Number n)+          `shouldBe` Just (fromInteger n :: IndexedPolynomial)++      prop "from symbol" $ \(SymbolText s) ->+        fromExpression (forVariable s) (Symbol s)+          `shouldBe` Just (power 1 :: IndexedPolynomial)++      prop "from symbol with exponent" $ \(SymbolText s) (Positive n) ->+        n > 1 ==>+          fromExpression (forVariable s) (Symbol s :**: Number n)+            `shouldBe` Just (power (fromIntegral n) :: IndexedPolynomial)++    describe "with symbolic coefficients" $ do+      prop "from number" $ \(SymbolText s) n ->+        fromExpression (withSymbolicCoefficients (forVariable s)) (Number n)+          `shouldBe` Just (fromInteger n :: IndexedSymbolicPolynomial)++      prop "from symbol" $ \(SymbolText s) ->+        fromExpression (withSymbolicCoefficients (forVariable s)) (Symbol s)+          `shouldBe` Just (power 1 :: IndexedSymbolicPolynomial)++      prop "from symbol with exponent" $ \(SymbolText s) (Positive n) ->+        n > 1 ==>+          let e = Symbol s :**: Number n+              p = fromExpression (withSymbolicCoefficients $ forVariable s) e+           in reduceSymbolicCoefficients <$> p `shouldBe` Just (power $ fromIntegral n)++-- | Reduce symbolic coefficients into rational number coefficients.+-- The representation of polynomials with rational number coefficients is unique,+-- which make them easier to compare.+reduceSymbolicCoefficients :: IndexedSymbolicPolynomial -> IndexedPolynomial+reduceSymbolicCoefficients = mapCoefficients reduce+  where+    reduce e+      | (Just x) <- fractionalEvaluate e (const Nothing) = x+      | otherwise = 0
+ test/Symtegration/PolynomialSpec.hs view
@@ -0,0 +1,170 @@+-- |+-- Description: Tests Symtegration.Polynomial.+-- Copyright: Copyright 2025 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.PolynomialSpec (spec) where++import Symtegration.Polynomial+import Symtegration.Polynomial.Indexed+import Symtegration.Polynomial.Indexed.Arbitrary ()+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck hiding (scale)++spec :: Spec+spec = parallel $ do+  describe "monic" $ do+    prop "is zero for zero" $+      monic 0 `shouldBe` (0 :: IndexedPolynomial)++    prop "has leading coefficient of one" $ \p ->+      p /= 0 ==> leadingCoefficient (monic p :: IndexedPolynomial) `shouldBe` 1++    prop "is rational multiple of original polynomial" $ \p ->+      let p' = monic p :: IndexedPolynomial+          (q, r) = p `divide` p'+       in counterexample (show p') $+            conjoin [r `shouldBe` 0, degree q `shouldBe` 0]++  describe "mapCoefficients" $ do+    prop "scales" $ \p x ->+      p /= 0 && x /= 0 ==>+        let q = mapCoefficients (* x) p :: IndexedPolynomial+         in conjoin+              [ monic p === monic q,+                leadingCoefficient p * x === leadingCoefficient q+              ]++  describe "mapCoefficientsM" $ do+    prop "with Maybe" $ \p (Fun _ f) ->+      let q = mapCoefficientsM (f :: Rational -> Maybe Rational) (p :: IndexedPolynomial)+          p' = filter (\(_, c) -> c /= 0) <$> mapM (\(e, c) -> (e,) <$> f c) (toList p)+          toList = foldTerms (\e c -> [(e, c)])+       in toList <$> q `shouldBe` p'++  describe "polynomial algorithms" $ do+    describe "division" $ do+      prop "matches multiplication" $ \a b ->+        degree b /= 0 ==>+          let (q, r) = divide a b+           in b * q + r `shouldBe` (a :: IndexedPolynomial)++      prop "remainder has smaller degree than divisor" $ \a b ->+        degree b > 0 ==>+          let (_, r) = divide a b+           in degree r `shouldSatisfy` (< degree (b :: IndexedPolynomial))++    describe "pseudo-division" $ do+      prop "matches division for integer coefficients" $ \a b ->+        b /= 0 ==>+          let delta = max (-1) (degree a - degree b)+              x = leadingCoefficient b ^ (1 + delta)+           in pseudoDivide a b `shouldBe` divide (scale x a) (b :: IndexedPolynomial)++    describe "extended Euclidean algorithm" $ do+      prop "gets common divisor" $ \a b ->+        let (_, _, g :: IndexedPolynomial) = extendedEuclidean a b+         in conjoin (map (\p -> let (_, r) = divide p g in r `shouldBe` 0) [a, b])++      prop "coefficients generate greatest common divisor" $ \a b ->+        let (s, t, g :: IndexedPolynomial) = extendedEuclidean a b+         in s * a + t * b `shouldBe` g++      prop "any sa+tb must be multiple of gcd a b" $ \a b s t ->+        let (_, _, g :: IndexedPolynomial) = extendedEuclidean a b+         in snd (divide (s * a + t * b) g) `shouldBe` 0++    describe "diophantine extended Euclidean algorithm" $ do+      prop "solves for (s,t)" $ \a b c ->+        degree a > 0 && degree b > 0 && degree c > 0 ==>+          let p Nothing =+                label "no solution" $+                  snd (c `divide` greatestCommonDivisor a b) `shouldSatisfy` (/= 0)+              p (Just (s, t)) =+                label "solved" $+                  counterexample ("(s,t) = " <> show (s, t)) $+                    conjoin+                      [ s * a + t * b === (c :: IndexedPolynomial),+                        disjoin [s === 0, property $ degree s < degree b]+                      ]+           in p (diophantineEuclidean a b c)++    describe "greatest common divisor" $ do+      prop "is consistent with extended Euclidean algorithm" $ \a b ->+        let (_, _, g :: IndexedPolynomial) = extendedEuclidean a b+         in greatestCommonDivisor a b `shouldBe` g++    describe "subresultant polynomial remainder sequence" $ do+      prop "resultant is zero iff gcd has non-zero degree" $ \a b ->+        b /= 0 ==>+          let (resultant, _) = subresultant a (b :: IndexedPolynomial)+           in resultant == 0 `shouldBe` degree (greatestCommonDivisor a b) > 0++      modifyMaxSize (const 25) $+        prop "resultant has expected value" $+          forAll (arbitrarySizedFractional `suchThat` (/= 0)) $ \a ->+            forAll (arbitrarySizedFractional `suchThat` (/= 0)) $ \b ->+              forAll (listOf1 arbitrarySizedFractional) $ \as ->+                forAll (listOf1 arbitrarySizedFractional) $ \bs ->+                  let x = scale a $ product [power 1 - scale t 1 | t <- as] :: IndexedPolynomial+                      y = scale b $ product [power 1 - scale t 1 | t <- bs] :: IndexedPolynomial+                      r@(resultant, _) = subresultant x y+                      resultant' = a ^ length bs * b ^ length as * product [u - v | u <- as, v <- bs]+                   in counterexample (show r) $+                        resultant `shouldBe` resultant'++      prop "is polynomial remainder sequence" $ \a b ->+        let (_, prs) = subresultant a b+            -- Whether z is a numeric multiple of prem(x, y).+            fromPseudoRemainder (x, y, z)+              | y == 0 = z == 0+              | prem == 0 = z == 0+              | otherwise = degree q == 0 && r == 0+              where+                (_, prem) = pseudoDivide x y+                (q, r) = divide z (prem :: IndexedPolynomial)+            isPolynomialRemainderSequence xs =+              all fromPseudoRemainder $ zip3 xs (drop 1 xs) (drop 2 xs)+         in prs `shouldSatisfy` isPolynomialRemainderSequence++      prop "has zero as last element in sequence" $ \a b ->+        let (_, prs) = subresultant a (b :: IndexedPolynomial)+         in counterexample (show prs) $+              drop (length prs - 1) prs `shouldBe` [0]++    describe "differentiation" $ do+      prop "computes derivative of constant" $ \c ->+        differentiate (scale c (power 0) :: IndexedPolynomial) `shouldBe` 0++      prop "computes derivative of integral power" $ \(Positive e) c ->+        differentiate (scale c (power e) :: IndexedPolynomial)+          `shouldBe` scale (fromIntegral e * c) (power (e - 1))++      prop "computes derivative of compound polynomials" $ \a b ->+        differentiate (a + b :: IndexedPolynomial)+          `shouldBe` differentiate a + differentiate b++    describe "integration" $ do+      prop "computes integral of integral power" $ \(NonNegative e) c ->+        integrate (scale c (power e) :: IndexedPolynomial)+          `shouldBe` scale (c / (1 + fromIntegral e)) (power (e + 1))++      prop "computes integral of compound polynomials" $ \a b ->+        integrate (a + b :: IndexedPolynomial)+          `shouldBe` integrate a + integrate b++    describe "squarefree factorization" $ do+      prop "divides polynomial" $ \p ->+        let qs = squarefree p :: [IndexedPolynomial]+         in counterexample (show qs) $+              conjoin $+                map (\q -> counterexample (show q) $ snd (p `divide` q) === 0) qs++      prop "multiplies to polynomial" $ \p ->+        let qs = squarefree p :: [IndexedPolynomial]+            prod :: Int -> [IndexedPolynomial] -> IndexedPolynomial+            prod _ [] = 1+            prod k (x : xs) = x ^ k * prod (k + 1) xs+         in counterexample (show qs) $+              prod 1 qs `shouldBe` p
+ test/Symtegration/Symbolic/Arbitrary.hs view
@@ -0,0 +1,175 @@+{-# OPTIONS_GHC -fno-warn-orphans #-}++-- |+-- Description: QuickCheck Arbitrary instances for generating Symtegration.Symbolic values.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Symbolic.Arbitrary+  ( Simple (..),+    Compound (..),+    Complete (..),+    SymbolMap (..),+    SymbolText (..),+    arbitraryNumber,+    arbitrarySymbol,+    arbitraryUnaryFunction,+    arbitraryBinaryFunction,+    arbitrarySymbolText,+    shrinkSymbolText,+    arbitrarySymbolMap,+    shrinkSymbolMap,+    assign,+  )+where++import Data.Map (Map)+import Data.Map qualified as Map+import Data.Set (Set)+import Data.Set qualified as S+import Data.String (fromString)+import Data.Text (Text)+import Data.Text qualified as Text+import Symtegration.ErrorDouble+import Symtegration.FiniteDouble+import Symtegration.Symbolic+import Test.QuickCheck++instance Arbitrary Expression where+  arbitrary = sized $ \n -> case n of+    0 -> oneof [arbitraryNumber, arbitrarySymbol]+    _ ->+      frequency+        [ (1, arbitraryNumber),+          (1, arbitrarySymbol),+          (4, resize (max 0 (n - 1)) arbitraryUnaryFunction),+          (8, resize (n `div` 2) arbitraryBinaryFunction)+        ]++  shrink (Number n) = Number <$> shrink n+  shrink (Symbol s) = Symbol <$> shrinkSymbolText s+  shrink (UnaryApply func x) = x : (UnaryApply func <$> shrink x)+  shrink (BinaryApply func x y) =+    x : y : [BinaryApply func x' y' | (x', y') <- shrink (x, y)]++instance Arbitrary UnaryFunction where+  arbitrary = chooseEnum (minBound, maxBound)++instance Arbitrary BinaryFunction where+  arbitrary = chooseEnum (minBound, maxBound)++-- | Generates simple symbolic mathematical expressions.+-- Specically, those which represent a single symbol or a single number.+newtype Simple = Simple Expression deriving (Eq, Show)++instance Arbitrary Simple where+  arbitrary = Simple <$> oneof [arbitraryNumber, arbitrarySymbol]++-- | Generates a compound symbolic mathematical expression.+-- Specifically, either a unary function application or a binary function application.+newtype Compound = Compound Expression deriving (Eq, Show)++instance Arbitrary Compound where+  arbitrary = Compound <$> oneof [arbitraryUnaryFunction, arbitraryBinaryFunction]+  shrink (Compound e) = Compound <$> filter isCompound (shrink e)+    where+      isCompound (Number _) = False+      isCompound (Symbol _) = False+      isCompound _ = True++-- | Generates arbitrary expressions with a complete assignment of numbers to symbols.+-- The assignment of symbols to values will only contain symbols appearing in the expression.+-- Use the 'assign' function to turn the map into a function.+data Complete = Complete Expression (Map Text FiniteDouble) deriving (Eq, Show)++instance Arbitrary Complete where+  arbitrary = do+    expr <- arbitrary+    vals <- infiniteList+    let symbols = gatherSymbols expr+    let assignment = Map.fromList $ zip (S.toList symbols) vals+    if not (sensitiveExpression expr (assign assignment))+      -- Only use expressions where slight divergences do not result in huge errors.+      then return $ Complete expr (Map.map FiniteDouble assignment)+      -- If we do not have such an expression, try again.+      else arbitrary++  shrink (Complete e m) = [Complete e' (restrict m e') | e' <- shrink e]+    where+      -- Keep symbol assignments still relevant to a shrinked expression.+      restrict xs x = Map.restrictKeys xs $ gatherSymbols x++-- | Gather the symbols appearing in an expression.+gatherSymbols :: Expression -> Set Text+gatherSymbols (Number _) = S.empty+gatherSymbols (Symbol s) = S.singleton s+gatherSymbols (UnaryApply _ x) = gatherSymbols x+gatherSymbols (BinaryApply _ x y) = S.union (gatherSymbols x) (gatherSymbols y)++-- | Generates a random assignment from symbols to values.+-- Use the 'assign' function to turn it into a function.+newtype SymbolMap a = SymbolMap (Map Text a) deriving (Eq, Show)++instance (Arbitrary a) => Arbitrary (SymbolMap a) where+  arbitrary = SymbolMap <$> arbitrarySymbolMap+  shrink (SymbolMap m) = SymbolMap <$> shrinkSymbolMap m++-- | Generates random readable symbol.+newtype SymbolText = SymbolText Text deriving (Eq, Show)++instance Arbitrary SymbolText where+  arbitrary = SymbolText <$> arbitrarySymbolText+  shrink (SymbolText s) = SymbolText <$> shrinkSymbolText s++-- | Generate a random number.+arbitraryNumber :: Gen Expression+arbitraryNumber = Number <$> arbitrary++-- | Generate a random symbol with only letters.+arbitrarySymbol :: Gen Expression+arbitrarySymbol = Symbol <$> arbitrarySymbolText++-- | Generate a random expression with an unary function application.+arbitraryUnaryFunction :: Gen Expression+arbitraryUnaryFunction = UnaryApply <$> arbitrary <*> arbitrary++-- | Generate a random expression with a binary function application.+arbitraryBinaryFunction :: Gen Expression+arbitraryBinaryFunction = BinaryApply <$> arbitrary <*> arbitrary <*> arbitrary++-- | Generate a random map from readable symbols to values.+arbitrarySymbolMap :: (Arbitrary a) => Gen (Map Text a)+arbitrarySymbolMap = Map.fromList <$> listOf assocs+  where+    assocs = do+      s <- arbitrarySymbolText+      x <- arbitrary+      return (s, x)++-- | Shrinks a map from readable symbols to values.+shrinkSymbolMap :: (Arbitrary a) => Map Text a -> [Map Text a]+shrinkSymbolMap = shrinkMapBy Map.fromList Map.toList (shrinkList shrinkAssoc)+  where+    shrinkAssoc (s, x) = do+      s' <- shrinkSymbolText s+      x' <- shrink x+      return (s', x')++-- | Generate random text that is appropriate as a readable symbol.+-- They will be short, since what exactly are in the symbols is usually not important.+-- Does not generate the special symbol "pi".+arbitrarySymbolText :: Gen Text+arbitrarySymbolText = resize 3 $ fromString <$> listOf1 (choose ('a', 'z')) `suchThat` (/= "pi")++-- | Shrinks readable symbols.+shrinkSymbolText :: Text -> [Text]+shrinkSymbolText s =+  -- Exclude empty text and s itself.+  drop 1 $ reverse $ drop 1 $ Text.tails s++-- | For creating a function which assigns symbols to values+-- based on the given map, which are easier to generate with+-- specific properties and easier to show than a function itself.+-- Shorthand for writing @assign m@ instead of @flip Map.lookup m@.+assign :: Map Text a -> Text -> Maybe a+assign = flip Map.lookup
+ test/Symtegration/Symbolic/HaskellSpec.hs view
@@ -0,0 +1,135 @@+-- |+-- Description: Tests for Symtegration.Symbolic.Haskell+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Symbolic.HaskellSpec (spec) where++import Data.String (fromString)+import Data.Text (Text, toLower)+import Symtegration.Symbolic+import Symtegration.Symbolic.Arbitrary+import Symtegration.Symbolic.Haskell+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck+import TextShow (showt)++spec :: Spec+spec = parallel $ do+  describe "toHaskell" $ do+    prop "converts for number" $ \n ->+      toHaskell (Number n) `shouldBe` showt n++    prop "converts for symbol" $ \(PrintableString s) ->+      toHaskell (Symbol $ fromString s) `shouldBe` fromString s++    describe "converts for unary function" $ do+      prop "with non-negative number" $ \func (NonNegative n) ->+        toHaskell (UnaryApply func $ Number n)+          `shouldBe` getUnaryFunctionText func <> " " <> showt n++      prop "with negative number" $ \func (Negative n) ->+        toHaskell (UnaryApply func $ Number n)+          `shouldBe` getUnaryFunctionText func <> " " <> par (showt n)++      prop "with symbol" $ \func s ->+        toHaskell (UnaryApply func $ Symbol $ fromString s)+          `shouldBe` getUnaryFunctionText func <> " " <> fromString s++      prop "with compound argument" $ \func (Compound e) ->+        toHaskell (UnaryApply func e)+          `shouldBe` getUnaryFunctionText func <> " " <> par (toHaskell e)++    describe "converts for binary function" $ do+      prop "logBase with non-negative numbers" $ \(NonNegative m) (NonNegative n) ->+        toHaskell (BinaryApply LogBase (Number m) (Number n))+          `shouldBe` "logBase " <> showt m <> " " <> showt n++      prop "logBase with negative numbers" $ \(Negative m) (Negative n) ->+        toHaskell (BinaryApply LogBase (Number m) (Number n))+          `shouldBe` "logBase " <> par (showt m) <> " " <> par (showt n)++      prop "logBase with symbols" $ \s r ->+        toHaskell (BinaryApply LogBase (Symbol $ fromString s) (Symbol $ fromString r))+          `shouldBe` "logBase " <> fromString s <> " " <> fromString r++      prop "logBase with compound arguments" $ \(Compound e1) (Compound e2) ->+        toHaskell (BinaryApply LogBase e1 e2)+          `shouldBe` "logBase " <> par (toHaskell e1) <> " " <> par (toHaskell e2)++      prop "operators with non-negative numbers" $ \op (NonNegative m) (NonNegative n) ->+        op /= LogBase ==>+          toHaskell (BinaryApply op (Number m) (Number n))+            `shouldBe` showt m <> " " <> getBinaryFunctionText op <> " " <> showt n++      prop "operators with negative numbers" $ \op (Negative m) (Negative n) ->+        op /= LogBase ==>+          toHaskell (BinaryApply op (Number m) (Number n))+            `shouldBe` par (showt m) <> " " <> getBinaryFunctionText op <> " " <> par (showt n)++      prop "operators with symbols" $ \op s r ->+        op /= LogBase ==>+          toHaskell (BinaryApply op (Symbol $ fromString s) (Symbol $ fromString r))+            `shouldBe` fromString s <> " " <> getBinaryFunctionText op <> " " <> fromString r++      prop "addition with compound arguments" $ \(Compound e1) (Compound e2) ->+        let text1 = toHaskell e1+            text2 = toHaskell e2+            t = toHaskell $ e1 :+: e2+         in t `shouldBe` text1 <> " + " <> text2++      prop "multiplication with compound arguments" $ \(Compound e1) (Compound e2) ->+        let text1 = toHaskell e1+            text2 = toHaskell e2+            multiply x y = x <> " * " <> y+            t = toHaskell $ e1 :*: e2+         in t `shouldBe` case (e1, e2) of+              (_ :+: _, _ :+: _) -> par text1 `multiply` par text2+              (_ :+: _, _ :-: _) -> par text1 `multiply` par text2+              (_ :-: _, _ :+: _) -> par text1 `multiply` par text2+              (_ :-: _, _ :-: _) -> par text1 `multiply` par text2+              (_ :+: _, _) -> par text1 `multiply` text2+              (_ :-: _, _) -> par text1 `multiply` text2+              (_, _ :+: _) -> text1 `multiply` par text2+              (_, _ :-: _) -> text1 `multiply` par text2+              _ -> text1 `multiply` text2++      prop "subtraction with compound arguments" $ \(Compound e1) (Compound e2) ->+        let text1 = toHaskell e1+            text2 = toHaskell e2+            minus x y = x <> " - " <> y+            t = toHaskell $ e1 :-: e2+         in t `shouldBe` case (e1, e2) of+              (_, _ :+: _) -> text1 `minus` par text2+              (_, _ :-: _) -> text1 `minus` par text2+              _ -> text1 `minus` text2++      prop "operators with compound arguments" $ \(Compound e1) (Compound e2) ->+        forAll (elements [Divide, Power]) $ \op ->+          let text1 = toHaskell e1+              text2 = toHaskell e2+              optext = getBinaryFunctionText op+              t = toHaskell (BinaryApply op e1 e2)+           in t `shouldBe` par text1 <> " " <> optext <> " " <> par text2++  -- The UnaryFunction constructors have the same spelling as their corresponding function name.+  describe "correct unary function text" $ do+    mapM_+      ( \func ->+          it ("for " <> show func) $+            getUnaryFunctionText func `shouldBe` toLower (showt func)+      )+      [minBound .. maxBound]++  describe "correct binary function text" $ do+    it "for Add" $ getBinaryFunctionText Add `shouldBe` "+"+    it "for Multiply" $ getBinaryFunctionText Multiply `shouldBe` "*"+    it "for Subtract" $ getBinaryFunctionText Subtract `shouldBe` "-"+    it "for Divide" $ getBinaryFunctionText Divide `shouldBe` "/"+    it "for Power" $ getBinaryFunctionText Power `shouldBe` "**"+    it "for LogBase" $ getBinaryFunctionText LogBase `shouldBe` "logBase"++-- | Surrounds the given text with parentheses.+par :: Text -> Text+par s = "(" <> s <> ")"
+ test/Symtegration/Symbolic/LaTeXSpec.hs view
@@ -0,0 +1,117 @@+-- |+-- Description: Tests for Symtegration.Symbolic.LaTeX.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Symbolic.LaTeXSpec (spec) where++import Symtegration.Symbolic.Arbitrary ()+import Symtegration.Symbolic.LaTeX+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck++spec :: Spec+spec = parallel $ describe "toLaTeX" $ do+  -- Full-fledged property-based tests would be nice, but for now,+  -- check only for totality property and test with specific examples instead.+  prop "is total" $ \e -> total (toLaTeX e)++  describe "addition" $ do+    it "-(1 + a)" $ toLaTeX (-(1 + "a")) `shouldBe` "-\\left(1 + a\\right)"++    it "1 + 5" $ toLaTeX (1 + 5) `shouldBe` "1 + 5"++    it "1 + 2 + 3" $ toLaTeX (1 + 2 + 3) `shouldBe` "1 + 2 + 3"++    it "x * y + u * v" $ toLaTeX ("x" * "y" + "u" * "v") `shouldBe` "x y + u v"++    it "x + u * v" $ toLaTeX ("x" + "u" * "v") `shouldBe` "x + u v"++    it "(x - y) + (u - v)" $ toLaTeX (("x" - "y") + ("u" - "v")) `shouldBe` "x - y + u - v"++    it "(-1) + (-4)" $ toLaTeX ((-1) + (-4)) `shouldBe` "-1 - 4"++    it "(-x) + y" $ toLaTeX ((-"x") + "y") `shouldBe` "-x + y"++    it "(-x) + (-y)" $ toLaTeX ((-"x") + (-"y")) `shouldBe` "-x - y"++    it "sin x + cos y" $ toLaTeX (sin "x" + cos "y") `shouldBe` "\\sin x + \\cos y"++    it "10 + (-4)" $ toLaTeX (10 + (-4)) `shouldBe` "10 - 4"++    it "x + (-y)" $ toLaTeX ("x" + (-"y")) `shouldBe` "x - y"++    it "x + 3 * x ** 2 + 2" $ toLaTeX ("x" + 3 * "x" ** 2 + 2) `shouldBe` "x + 3 x^{2} + 2"++  describe "subtraction" $ do+    it "x - y" $ toLaTeX ("x" - "y") `shouldBe` "x - y"++    it "(x + y) - z" $ toLaTeX (("x" + "y") - "z") `shouldBe` "x + y - z"++    it "x - (y + z)" $ toLaTeX ("x" - ("y" + "z")) `shouldBe` "x - \\left(y + z\\right)"++    it "x - (y - z)" $ toLaTeX ("x" - ("y" - "z")) `shouldBe` "x - \\left(y - z\\right)"++    it "x - (y * z)" $ toLaTeX ("x" - ("y" * "z")) `shouldBe` "x - y z"++    it "(-x) - y" $ toLaTeX ((-"x") - "y") `shouldBe` "-x - y"++    it "x - 3 * x ** 2 + 2" $ toLaTeX ("x" - 3 * "x" ** 2 + 2) `shouldBe` "x - 3 x^{2} + 2"++  describe "multiplication" $ do+    it "2 * 5" $ toLaTeX (2 * 5) `shouldBe` "2 \\times 5"++    it "2 * 3 * 4 * 6" $ toLaTeX (2 * 3 * 4 * 6) `shouldBe` "2 \\times 3 \\times 4 \\times 6"++    it "(2 * 3) * (4 * 6)" $ toLaTeX ((2 * 3) * (4 * 6)) `shouldBe` "2 \\times 3 \\times 4 \\times 6"++    it "2 * (4 * 6)" $ toLaTeX (2 * (4 * 6)) `shouldBe` "2 \\times 4 \\times 6"++    it "sin x * 3" $ toLaTeX (sin "x" * 3) `shouldBe` "\\sin x \\times 3"++    it "abs x * y" $ toLaTeX (abs "x" * "y") `shouldBe` "\\left\\lvert x \\right\\rvert y"++    it "signum x * y" $ toLaTeX (signum "x" * "y") `shouldBe` "\\mathrm{signum}\\left(x\\right) y"++    it "exp x * y" $ toLaTeX (exp "x" * "y") `shouldBe` "e^{x} y"++    it "sin x * y" $ toLaTeX (sin "x" * "y") `shouldBe` "\\left(\\sin x\\right) y"++    it "(-2) * (-5)" $ toLaTeX ((-2) * (-5)) `shouldBe` "\\left(-2\\right) \\left(-5\\right)"++    it "sin x * cos y" $ toLaTeX (sin "x" * cos "y") `shouldBe` "\\sin x \\cos y"++    it "4 * sin x" $ toLaTeX (4 * sin "x") `shouldBe` "4 \\sin x"++    it "x * y" $ toLaTeX ("x" * "y") `shouldBe` "x y"++    it "x * y ** z" $ toLaTeX ("x" * "y" ** "z") `shouldBe` "x y^{z}"++    it "log x * y" $ toLaTeX (log "x" * "y") `shouldBe` "\\left(\\log x\\right) y"++    it "logBase x y * z" $ toLaTeX (logBase "x" "y" * "z") `shouldBe` "\\left(\\log_{x}y\\right) z"++  describe "negation" $ do+    it "-a" $ toLaTeX (-"a") `shouldBe` "-a"++    it "-19" $ toLaTeX (negate 19) `shouldBe` "-19"++    it "-x" $ toLaTeX (negate "x") `shouldBe` "-x"++    it "-(-x)" $ toLaTeX (negate (negate "x")) `shouldBe` "-\\left(-x\\right)"++    it "-(x + y)" $ toLaTeX (negate ("x" + "y")) `shouldBe` "-\\left(x + y\\right)"++    it "-(x - y)" $ toLaTeX (negate ("x" - "y")) `shouldBe` "-\\left(x - y\\right)"++    it "-(sin x)" $ toLaTeX (negate (sin "x")) `shouldBe` "-\\sin x"++    it "-(x * sin x)" $ toLaTeX (negate ("x" * sin "x")) `shouldBe` "-x \\sin x"++    it "x * negate (sin x)" $ toLaTeX ("x" * negate (sin "x")) `shouldBe` "x \\left(-\\sin x\\right)"++  describe "unary function arguments" $ do+    it "cos (log x)" $ toLaTeX (cos (log "x")) `shouldBe` "\\cos \\left(\\log x\\right)"++    it "tan (pi * x)" $ toLaTeX (tan (pi * "x")) `shouldBe` "\\tan \\left(\\pi x\\right)"
+ test/Symtegration/Symbolic/Simplify/AlgebraicRingOrderSpec.hs view
@@ -0,0 +1,28 @@+-- |+-- Description: Tests for Symtegration.Symbolic.Simplify.AlgebraicRingOrder.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Symbolic.Simplify.AlgebraicRingOrderSpec (spec) where++import Data.Map qualified as Map+import Data.Text (Text)+import Symtegration.Symbolic+import Symtegration.Symbolic.Arbitrary+import Symtegration.Symbolic.Simplify.AlgebraicRingOrder+import Symtegration.Symbolic.Simplify.Properties+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck++spec :: Spec+spec = parallel $ do+  describe "simplify" $ do+    modifyMaxSuccess (* 100) $+      prop "maintains semantics" $+        equivalentProperty' order++equivalentProperty' :: (Text -> Expression -> Expression) -> Complete -> Property+equivalentProperty' f (Complete e m) = do+  forAll (elements $ Map.keys m) $ \v ->+    equivalentProperty (f v) (Complete e m)
+ test/Symtegration/Symbolic/Simplify/FractionSpec.hs view
@@ -0,0 +1,18 @@+-- |+-- Description: Tests Symtegration.Symbolic.Simplify.Fraction+-- Copyright: Copyright 2025 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Symbolic.Simplify.FractionSpec (spec) where++import Symtegration.Symbolic.Simplify.Fraction+import Symtegration.Symbolic.Simplify.Properties+import Test.Hspec+import Test.Hspec.QuickCheck++spec :: Spec+spec = parallel $ do+  describe "simplify" $ do+    modifyMaxSuccess (* 100) $+      prop "maintains semantics" $+        equivalentProperty simplify
+ test/Symtegration/Symbolic/Simplify/NumericFoldingSpec.hs view
@@ -0,0 +1,64 @@+-- |+-- Description: Tests Symtegration.Symbolic.Simplify.NumericFolding+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Symbolic.Simplify.NumericFoldingSpec (spec) where++import Data.Text (unpack)+import Symtegration.Symbolic+import Symtegration.Symbolic.Arbitrary+import Symtegration.Symbolic.Haskell+import Symtegration.Symbolic.Simplify.NumericFolding+import Symtegration.Symbolic.Simplify.Properties+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck++spec :: Spec+spec = parallel $ do+  describe "simplify" $ do+    modifyMaxSuccess (* 100) $+      prop "maintains semantics" $+        equivalentProperty simplify++    modifyMaxSuccess (* 100) $+      prop "folds to simple numeric expressions" $+        forAll genNumeric $ \e ->+          simpleNumeric e ==>+            let e' = simplify e+             in counterexample ("e = " <> unpack (toHaskell e)) $+                  counterexample ("simplify e = " <> unpack (toHaskell e')) $+                    e' `shouldSatisfy` simpleNumeric++-- Numeric folding should be able to fold arithmetic on numbers+-- to either an integer or a fraction.  The exception is if there+-- is a divide by zero somewhere, which we intentionally leave alone.+simpleNumeric :: Expression -> Bool+simpleNumeric (Number _) = True+simpleNumeric (Number _ :/: Number _) = True+simpleNumeric x = hasDivideByZero x+  where+    hasDivideByZero (_ :/: 0) = True+    hasDivideByZero (UnaryApply _ u) = hasDivideByZero u+    hasDivideByZero (BinaryApply _ u v) = hasDivideByZero u || hasDivideByZero v+    hasDivideByZero _ = False++-- | Generate arbitrary expression involving no symbols and which are+-- guaranteed to reduce exactly to a simple numeric term.+genNumeric :: Gen Expression+genNumeric = sized $ \case+  0 -> arbitraryNumber+  n ->+    frequency+      [ (1, arbitraryNumber),+        (1, resize (max 0 (n - 1)) $ UnaryApply Negate <$> genNumeric),+        ( 1,+          resize (max 0 (n - 1)) $+            BinaryApply Power <$> genNumeric <*> (Number <$> arbitrarySizedNatural `suchThat` (/= 0))+        ),+        ( 10,+          resize (n `div` 2) $+            BinaryApply <$> elements [Add, Multiply, Subtract, Divide] <*> genNumeric <*> genNumeric+        )+      ]
+ test/Symtegration/Symbolic/Simplify/Properties.hs view
@@ -0,0 +1,23 @@+-- |+-- Description: Provides a QuickCheck property checking that simplification does not change an expression's semantics.+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Symbolic.Simplify.Properties (equivalentProperty) where++import Symtegration.FiniteDouble+import Symtegration.Symbolic+import Symtegration.Symbolic.Arbitrary+import Symtegration.Symbolic.Haskell+import Test.Hspec+import Test.QuickCheck++equivalentProperty :: (Expression -> Expression) -> Complete -> Property+equivalentProperty simplify (Complete e m) =+  let e' = simplify e+      v = evaluate e (assign m)+      v' = evaluate e' (assign m)+   in counterexample ("e = " <> show (toHaskell e)) $+        counterexample ("simplify e = " <> show (toHaskell e')) $+          maybe False isFinite v && maybe False isFinite v' ==>+            fmap Near v' `shouldBe` fmap Near v
+ test/Symtegration/Symbolic/Simplify/SymbolicFoldingSpec.hs view
@@ -0,0 +1,18 @@+-- |+-- Description: Tests Symtegration.Symbolic.Simplify.Symbolic+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Symbolic.Simplify.SymbolicFoldingSpec (spec) where++import Symtegration.Symbolic.Simplify.Properties+import Symtegration.Symbolic.Simplify.SymbolicFolding+import Test.Hspec+import Test.Hspec.QuickCheck++spec :: Spec+spec = parallel $ do+  describe "simplify" $ do+    modifyMaxSuccess (* 100) $+      prop "maintains semantics" $+        equivalentProperty simplify
+ test/Symtegration/Symbolic/Simplify/TidySpec.hs view
@@ -0,0 +1,69 @@+-- |+-- Description: Tests Symtegration.Symbolic.Simplify.Tidy.+-- Copyright: Copyright 2025 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.Symbolic.Simplify.TidySpec (spec) where++import Symtegration.Symbolic+import Symtegration.Symbolic.Arbitrary+import Symtegration.Symbolic.Simplify.Properties+import Symtegration.Symbolic.Simplify.Tidy+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck++spec :: Spec+spec = parallel $ do+  describe "tidy" $ do+    modifyMaxSuccess (* 100) $+      prop "maintains semantics" $+        equivalentProperty tidy++    prop "x + negate y" $+      forAll arbitrarySymbol $ \x ->+        forAll arbitrarySymbol $ \y ->+          tidy (x + negate y) `shouldBe` x - y++    prop "x + (-1) * y" $+      forAll arbitrarySymbol $ \x ->+        forAll arbitrarySymbol $ \y ->+          tidy (x + Number (-1) * y) `shouldBe` x - y++    prop "x + (-n) * y" $ \(Positive n) ->+      n > 1 ==>+        forAll arbitrarySymbol $ \x ->+          forAll arbitrarySymbol $ \y ->+            tidy (x + Number (-n) * y) `shouldBe` x - Number n * y++    prop "(-n) / m" $ \(Positive n) (Positive m) ->+      tidy (Number (-n) / Number m) `shouldBe` negate (Number n / Number m)++    prop "(-n) / x" $ \(Positive n) ->+      forAll arbitrarySymbol $ \x ->+        tidy (Number (-n) / x) `shouldBe` negate (Number n / x)++    prop "(-x) * y" $+      forAll arbitrarySymbol $ \x ->+        forAll arbitrarySymbol $ \y ->+          tidy ((-x) * y) `shouldBe` negate (x * y)++    prop "x * (-y)" $+      forAll arbitrarySymbol $ \x ->+        forAll arbitrarySymbol $ \y ->+          tidy (x * (-y)) `shouldBe` negate (x * y)++    prop "(-x) * (-y)" $+      forAll arbitrarySymbol $ \x ->+        forAll arbitrarySymbol $ \y ->+          tidy ((-x) * (-y)) `shouldBe` x * y++    prop "x + ((-y) + z)" $+      forAll arbitrarySymbol $ \x ->+        forAll arbitrarySymbol $ \y ->+          forAll arbitrarySymbol $ \z ->+            tidy (x + ((-y) + z)) `shouldBe` x - y + z++    prop "(x + y) ** (1/2)" $+      forAll arbitrarySymbol $ \x ->+        tidy (x ** (1 / 2)) `shouldBe` sqrt x
+ test/Symtegration/SymbolicSpec.hs view
@@ -0,0 +1,218 @@+-- |+-- Description: Tests for Symtegration.SymbolicSpec+-- Copyright: Copyright 2024 Yoo Chung+-- License: Apache-2.0+-- Maintainer: dev@chungyc.org+module Symtegration.SymbolicSpec (spec) where++import Data.Map qualified as Map+import Data.Maybe (isJust)+import Data.Ratio (denominator, numerator)+import Data.String (fromString)+import Data.Text (Text)+import Data.Text qualified as Text+import Symtegration.FiniteDouble+import Symtegration.Symbolic+import Symtegration.Symbolic.Arbitrary+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck++-- | Same as 'evaluate', except specialized to 'FiniteDouble'.+evaluate' :: Expression -> (Text -> Maybe FiniteDouble) -> Maybe FiniteDouble+evaluate' = evaluate++spec :: Spec+spec = parallel $ do+  describe "Expression from" $ modifyMaxSuccess (`div` 10) $ do+    describe "IsString" $ do+      prop "fromString" $ \(SymbolText s) ->+        fromString (Text.unpack s) `shouldBe` Symbol s++    describe "Num" $ do+      prop "+" $ \(Simple x) (Simple y) ->+        x + y `shouldBe` BinaryApply Add x y++      prop "-" $ \(Simple x) (Simple y) ->+        x - y `shouldBe` BinaryApply Subtract x y++      prop "*" $ \(Simple x) (Simple y) ->+        x * y `shouldBe` BinaryApply Multiply x y++      prop "negate" $ \(Simple x) ->+        negate x `shouldBe` UnaryApply Negate x++      prop "abs" $ \(Simple x) ->+        abs x `shouldBe` UnaryApply Abs x++      prop "signum" $ \(Simple x) ->+        signum x `shouldBe` UnaryApply Signum x++      prop "fromInteger" $ \n ->+        fromInteger n `shouldBe` Number n++    describe "Fractional" $ do+      prop "/" $ \(Simple x) (Simple y) ->+        x / y `shouldBe` BinaryApply Divide x y++      prop "recip" $ \(Simple x) ->+        recip x `shouldBe` BinaryApply Divide 1 x++      prop "fromRational" $ \x ->+        let n = fromInteger $ numerator x+            d = fromInteger $ denominator x+         in case d of+              1 -> label "integer" $ fromRational x `shouldBe` n+              _ -> label "fraction" $ fromRational x `shouldBe` BinaryApply Divide n d++    describe "Floating" $ do+      prop "pi" $ pi `shouldBe` Symbol "pi"++      prop "exp" $ \(Simple x) ->+        exp x `shouldBe` UnaryApply Exp x++      prop "log" $ \(Simple x) ->+        log x `shouldBe` UnaryApply Log x++      prop "sqrt" $ \(Simple x) ->+        sqrt x `shouldBe` UnaryApply Sqrt x++      prop "**" $ \(Simple x) (Simple y) ->+        x ** y `shouldBe` BinaryApply Power x y++      prop "logBase" $ \(Simple x) (Simple y) ->+        logBase x y `shouldBe` BinaryApply LogBase x y++      prop "sin" $ \(Simple x) ->+        sin x `shouldBe` UnaryApply Sin x++      prop "cos" $ \(Simple x) ->+        cos x `shouldBe` UnaryApply Cos x++      prop "tan" $ \(Simple x) ->+        tan x `shouldBe` UnaryApply Tan x++      prop "asin" $ \(Simple x) ->+        asin x `shouldBe` UnaryApply Asin x++      prop "acos" $ \(Simple x) ->+        acos x `shouldBe` UnaryApply Acos x++      prop "atan" $ \(Simple x) ->+        atan x `shouldBe` UnaryApply Atan x++      prop "sinh" $ \(Simple x) ->+        sinh x `shouldBe` UnaryApply Sinh x++      prop "cosh" $ \(Simple x) ->+        cosh x `shouldBe` UnaryApply Cosh x++      prop "tanh" $ \(Simple x) ->+        tanh x `shouldBe` UnaryApply Tanh x++      prop "asinh" $ \(Simple x) ->+        asinh x `shouldBe` UnaryApply Asinh x++      prop "acosh" $ \(Simple x) ->+        acosh x `shouldBe` UnaryApply Acosh x++      prop "atanh" $ \(Simple x) ->+        atanh x `shouldBe` UnaryApply Atanh x++  describe "substitute" $ do+    prop "for number" $ \n (SymbolMap m) ->+      substitute (Number n) (assign m) `shouldBe` Number n++    prop "for unmapped symbol" $ \(SymbolText s) ->+      substitute (Symbol s) (const Nothing) `shouldBe` Symbol s++    prop "for mapped symbol" $ \(SymbolText s) e (SymbolMap m) ->+      let m' = Map.insert s e m+       in substitute (Symbol s) (assign m') `shouldBe` e++    prop "for unary function" $ \func e (SymbolMap m) ->+      substitute (UnaryApply func e) (assign m)+        `shouldBe` UnaryApply func (substitute e (assign m))++    prop "for binary function" $ \func x y (SymbolMap m) ->+      substitute (BinaryApply func x y) (assign m)+        `shouldBe` BinaryApply func (substitute x $ assign m) (substitute y $ assign m)++  describe "Expression exactly evaluates as" $ do+    prop "number" $ \n (SymbolMap m) ->+      evaluate' (Number n) (assign m) `shouldBe` Just (fromInteger n)++    prop "symbol" $ \(SymbolText s) x ->+      evaluate' (Symbol s) (\s' -> if s' == s then Just x else Nothing) `shouldBe` Just x++    prop "unary function" $ \(Complete e m) func ->+      fmap Exact (evaluate' (UnaryApply func e) (assign m))+        `shouldBe` fmap (Exact . getUnaryFunction func) (evaluate' e (assign m))++    prop "binary function" $ \(Complete e1 m1) (Complete e2 m2) func ->+      let m = Map.union m1 m2+          f = getBinaryFunction func+       in fmap Exact (evaluate' (BinaryApply func e1 e2) (assign m))+            `shouldBe` fmap Exact (f <$> evaluate' e1 (assign m) <*> evaluate' e2 (assign m))++    prop "nothing" $ \(Complete e m) ->+      not (Map.null m) ==> evaluate' e (const Nothing) `shouldBe` Nothing++  describe "Expression fractionally evaluates as" $ do+    prop "number" $ \n ->+      fractionalEvaluate (Number n) (const Nothing) `shouldBe` Just (fromInteger n :: Rational)++    prop "symbol" $ \(SymbolText s) x ->+      fractionalEvaluate (Symbol s) (const $ Just x) `shouldBe` Just (x :: Rational)++    prop "similar to evaluate" $ \(Complete e m) ->+      let v = fractionalEvaluate e (fmap toRational . assign m)+          v' = evaluate e (assign m)+       in maybe False isFinite v' && isJust v ==>+            Near . FiniteDouble . fromRational <$> v `shouldBe` Near <$> v'++  describe "unary functions are correctly mapped for" $ do+    mapM_+      ( \(func, f) -> prop (show func) $ \x ->+          Exact (getUnaryFunction func x) `shouldBe` Exact (f x)+      )+      ( [ (Negate, negate),+          (Abs, abs),+          (Signum, signum),+          (Exp, exp),+          (Log, log),+          (Sqrt, sqrt),+          (Sin, sin),+          (Cos, cos),+          (Tan, tan),+          (Asin, asin),+          (Acos, acos),+          (Atan, atan),+          (Sinh, sinh),+          (Cosh, cosh),+          (Tanh, tanh),+          (Asinh, asinh),+          (Acosh, acosh),+          (Atanh, atanh)+        ] ::+          [(UnaryFunction, FiniteDouble -> FiniteDouble)]+      )++  describe "binary functions are correctly mapped for" $ do+    mapM_+      ( \(func, f) -> prop (show func) $+          \x y -> Exact (getBinaryFunction func x y) `shouldBe` Exact (f x y)+      )+      ( [ (Add, (+)),+          (Multiply, (*)),+          (Subtract, (-)),+          (Divide, (/)),+          (Power, (**)),+          (LogBase, logBase)+        ] ::+          [(BinaryFunction, FiniteDouble -> FiniteDouble -> FiniteDouble)]+      )++  describe "show" $ do+    prop "has inverse with read" $ \e ->+      read (show e) `shouldBe` (e :: Expression)