diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -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. diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -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 diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,3 @@ +import Distribution.Simple + +main = defaultMain diff --git a/docs/CODE_OF_CONDUCT.md b/docs/CODE_OF_CONDUCT.md new file mode 100644 --- /dev/null +++ b/docs/CODE_OF_CONDUCT.md @@ -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. diff --git a/docs/CONTRIBUTING.md b/docs/CONTRIBUTING.md new file mode 100644 --- /dev/null +++ b/docs/CONTRIBUTING.md @@ -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/ diff --git a/docs/SECURITY.md b/docs/SECURITY.md new file mode 100644 --- /dev/null +++ b/docs/SECURITY.md @@ -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 diff --git a/src/Symtegration.hs b/src/Symtegration.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration.hs @@ -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 diff --git a/src/Symtegration/Differentiation.hs b/src/Symtegration/Differentiation.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Differentiation.hs @@ -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 diff --git a/src/Symtegration/Integration.hs b/src/Symtegration/Integration.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Integration.hs @@ -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] diff --git a/src/Symtegration/Integration/Exponential.hs b/src/Symtegration/Integration/Exponential.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Integration/Exponential.hs @@ -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 diff --git a/src/Symtegration/Integration/Factor.hs b/src/Symtegration/Integration/Factor.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Integration/Factor.hs @@ -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 diff --git a/src/Symtegration/Integration/Parts.hs b/src/Symtegration/Integration/Parts.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Integration/Parts.hs @@ -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 diff --git a/src/Symtegration/Integration/Powers.hs b/src/Symtegration/Integration/Powers.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Integration/Powers.hs @@ -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 diff --git a/src/Symtegration/Integration/Rational.hs b/src/Symtegration/Integration/Rational.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Integration/Rational.hs @@ -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 diff --git a/src/Symtegration/Integration/Substitution.hs b/src/Symtegration/Integration/Substitution.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Integration/Substitution.hs @@ -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 diff --git a/src/Symtegration/Integration/Sum.hs b/src/Symtegration/Integration/Sum.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Integration/Sum.hs @@ -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 diff --git a/src/Symtegration/Integration/Term.hs b/src/Symtegration/Integration/Term.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Integration/Term.hs @@ -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' diff --git a/src/Symtegration/Integration/Trigonometric.hs b/src/Symtegration/Integration/Trigonometric.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Integration/Trigonometric.hs @@ -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 diff --git a/src/Symtegration/Numeric.hs b/src/Symtegration/Numeric.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Numeric.hs @@ -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 diff --git a/src/Symtegration/Polynomial.hs b/src/Symtegration/Polynomial.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Polynomial.hs @@ -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 diff --git a/src/Symtegration/Polynomial/Indexed.hs b/src/Symtegration/Polynomial/Indexed.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Polynomial/Indexed.hs @@ -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 diff --git a/src/Symtegration/Polynomial/Solve.hs b/src/Symtegration/Polynomial/Solve.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Polynomial/Solve.hs @@ -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)) + ] diff --git a/src/Symtegration/Polynomial/Symbolic.hs b/src/Symtegration/Polynomial/Symbolic.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Polynomial/Symbolic.hs @@ -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 diff --git a/src/Symtegration/Symbolic.hs b/src/Symtegration/Symbolic.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Symbolic.hs @@ -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 diff --git a/src/Symtegration/Symbolic/Haskell.hs b/src/Symtegration/Symbolic/Haskell.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Symbolic/Haskell.hs @@ -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" diff --git a/src/Symtegration/Symbolic/LaTeX.hs b/src/Symtegration/Symbolic/LaTeX.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Symbolic/LaTeX.hs @@ -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 <> "}" diff --git a/src/Symtegration/Symbolic/Simplify.hs b/src/Symtegration/Symbolic/Simplify.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Symbolic/Simplify.hs @@ -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 diff --git a/src/Symtegration/Symbolic/Simplify/AlgebraicRingOrder.hs b/src/Symtegration/Symbolic/Simplify/AlgebraicRingOrder.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Symbolic/Simplify/AlgebraicRingOrder.hs @@ -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 diff --git a/src/Symtegration/Symbolic/Simplify/Fraction.hs b/src/Symtegration/Symbolic/Simplify/Fraction.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Symbolic/Simplify/Fraction.hs @@ -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 diff --git a/src/Symtegration/Symbolic/Simplify/NumericFolding.hs b/src/Symtegration/Symbolic/Simplify/NumericFolding.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Symbolic/Simplify/NumericFolding.hs @@ -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 diff --git a/src/Symtegration/Symbolic/Simplify/SymbolicFolding.hs b/src/Symtegration/Symbolic/Simplify/SymbolicFolding.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Symbolic/Simplify/SymbolicFolding.hs @@ -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 diff --git a/src/Symtegration/Symbolic/Simplify/Tidy.hs b/src/Symtegration/Symbolic/Simplify/Tidy.hs new file mode 100644 --- /dev/null +++ b/src/Symtegration/Symbolic/Simplify/Tidy.hs @@ -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 diff --git a/symtegration.cabal b/symtegration.cabal new file mode 100644 --- /dev/null +++ b/symtegration.cabal @@ -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 diff --git a/test/Examples.hs b/test/Examples.hs new file mode 100644 --- /dev/null +++ b/test/Examples.hs @@ -0,0 +1,7 @@ +module Main (main) where + +import System.Environment (getArgs) +import Test.DocTest (mainFromCabal) + +main :: IO () +main = mainFromCabal "symtegration" =<< getArgs diff --git a/test/Spec.hs b/test/Spec.hs new file mode 100644 --- /dev/null +++ b/test/Spec.hs @@ -0,0 +1,1 @@ +{-# OPTIONS_GHC -F -pgmF hspec-discover #-} diff --git a/test/Symtegration/ErrorDouble.hs b/test/Symtegration/ErrorDouble.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/ErrorDouble.hs @@ -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 diff --git a/test/Symtegration/FiniteDouble.hs b/test/Symtegration/FiniteDouble.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/FiniteDouble.hs @@ -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 diff --git a/test/Symtegration/Integration/ExponentialSpec.hs b/test/Symtegration/Integration/ExponentialSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Integration/ExponentialSpec.hs @@ -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" diff --git a/test/Symtegration/Integration/FactorSpec.hs b/test/Symtegration/Integration/FactorSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Integration/FactorSpec.hs @@ -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" diff --git a/test/Symtegration/Integration/PartsSpec.hs b/test/Symtegration/Integration/PartsSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Integration/PartsSpec.hs @@ -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" diff --git a/test/Symtegration/Integration/PowersSpec.hs b/test/Symtegration/Integration/PowersSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Integration/PowersSpec.hs @@ -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" diff --git a/test/Symtegration/Integration/Properties.hs b/test/Symtegration/Integration/Properties.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Integration/Properties.hs @@ -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 diff --git a/test/Symtegration/Integration/RationalSpec.hs b/test/Symtegration/Integration/RationalSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Integration/RationalSpec.hs @@ -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" diff --git a/test/Symtegration/Integration/SubstitutionSpec.hs b/test/Symtegration/Integration/SubstitutionSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Integration/SubstitutionSpec.hs @@ -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" diff --git a/test/Symtegration/Integration/SumSpec.hs b/test/Symtegration/Integration/SumSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Integration/SumSpec.hs @@ -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" diff --git a/test/Symtegration/Integration/TermSpec.hs b/test/Symtegration/Integration/TermSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Integration/TermSpec.hs @@ -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" diff --git a/test/Symtegration/Integration/TrigonometricSpec.hs b/test/Symtegration/Integration/TrigonometricSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Integration/TrigonometricSpec.hs @@ -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" diff --git a/test/Symtegration/IntegrationSpec.hs b/test/Symtegration/IntegrationSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/IntegrationSpec.hs @@ -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 diff --git a/test/Symtegration/NumericSpec.hs b/test/Symtegration/NumericSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/NumericSpec.hs @@ -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 diff --git a/test/Symtegration/Polynomial/Indexed/Arbitrary.hs b/test/Symtegration/Polynomial/Indexed/Arbitrary.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Polynomial/Indexed/Arbitrary.hs @@ -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] diff --git a/test/Symtegration/Polynomial/IndexedSpec.hs b/test/Symtegration/Polynomial/IndexedSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Polynomial/IndexedSpec.hs @@ -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]] diff --git a/test/Symtegration/Polynomial/SolveSpec.hs b/test/Symtegration/Polynomial/SolveSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Polynomial/SolveSpec.hs @@ -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 diff --git a/test/Symtegration/Polynomial/SymbolicSpec.hs b/test/Symtegration/Polynomial/SymbolicSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Polynomial/SymbolicSpec.hs @@ -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 diff --git a/test/Symtegration/PolynomialSpec.hs b/test/Symtegration/PolynomialSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/PolynomialSpec.hs @@ -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 diff --git a/test/Symtegration/Symbolic/Arbitrary.hs b/test/Symtegration/Symbolic/Arbitrary.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Symbolic/Arbitrary.hs @@ -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 diff --git a/test/Symtegration/Symbolic/HaskellSpec.hs b/test/Symtegration/Symbolic/HaskellSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Symbolic/HaskellSpec.hs @@ -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 <> ")" diff --git a/test/Symtegration/Symbolic/LaTeXSpec.hs b/test/Symtegration/Symbolic/LaTeXSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Symbolic/LaTeXSpec.hs @@ -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)" diff --git a/test/Symtegration/Symbolic/Simplify/AlgebraicRingOrderSpec.hs b/test/Symtegration/Symbolic/Simplify/AlgebraicRingOrderSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Symbolic/Simplify/AlgebraicRingOrderSpec.hs @@ -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) diff --git a/test/Symtegration/Symbolic/Simplify/FractionSpec.hs b/test/Symtegration/Symbolic/Simplify/FractionSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Symbolic/Simplify/FractionSpec.hs @@ -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 diff --git a/test/Symtegration/Symbolic/Simplify/NumericFoldingSpec.hs b/test/Symtegration/Symbolic/Simplify/NumericFoldingSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Symbolic/Simplify/NumericFoldingSpec.hs @@ -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 + ) + ] diff --git a/test/Symtegration/Symbolic/Simplify/Properties.hs b/test/Symtegration/Symbolic/Simplify/Properties.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Symbolic/Simplify/Properties.hs @@ -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 diff --git a/test/Symtegration/Symbolic/Simplify/SymbolicFoldingSpec.hs b/test/Symtegration/Symbolic/Simplify/SymbolicFoldingSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Symbolic/Simplify/SymbolicFoldingSpec.hs @@ -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 diff --git a/test/Symtegration/Symbolic/Simplify/TidySpec.hs b/test/Symtegration/Symbolic/Simplify/TidySpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/Symbolic/Simplify/TidySpec.hs @@ -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 diff --git a/test/Symtegration/SymbolicSpec.hs b/test/Symtegration/SymbolicSpec.hs new file mode 100644 --- /dev/null +++ b/test/Symtegration/SymbolicSpec.hs @@ -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)