Flint2-Examples (empty) → 0.1.0.0
raw patch · 53 files changed
+6438/−0 lines, 53 filesdep +Flint2dep +Flint2-Examplesdep +JuicyPixelssetup-changedbinary-added
Dependencies added: Flint2, Flint2-Examples, JuicyPixels, base, containers, filepath, gloss, gloss-juicy, mtl, optparse-applicative, time, timeit
Files
- CHANGELOG.md +11/−0
- Flint2-Examples.cabal +813/−0
- LICENSE +339/−0
- README.md +171/−0
- Setup.hs +2/−0
- app/bernoulli/Main.hs +41/−0
- app/binet/Main.hs +60/−0
- app/class_poly/Main.hs +56/−0
- app/complex_plot/ColorFunction.hs +113/−0
- app/complex_plot/Functions.hs +193/−0
- app/complex_plot/Main.hs +134/−0
- app/crt/Main.hs +53/−0
- app/delta_qexp/Main.hs +43/−0
- app/dft/Main.hs +257/−0
- app/elementary/Main.hs +244/−0
- app/expression/Main.hs +179/−0
- app/factor_integer/Main.hs +70/−0
- app/fmpq_poly/Main.hs +20/−0
- app/fmpz_mod_poly/Main.hs +14/−0
- app/fmpz_mpoly_factor/Main.hs +74/−0
- app/fmpz_poly_factor_zassenhaus/Main.hs +18/−0
- app/fpwrap/Main.hs +40/−0
- app/fq_poly/Main.hs +58/−0
- app/function_benchmark/Main.hs +83/−0
- app/hilbert_matrix/Main.hs +69/−0
- app/hilbert_matrix_ca/Main.hs +113/−0
- app/integrals/Integrands.hs +532/−0
- app/integrals/Main.hs +619/−0
- app/keiper_li/Main.hs +158/−0
- app/l_values/Main.hs +110/−0
- app/logistic/Main.hs +101/−0
- app/machin/Main.hs +88/−0
- app/multi_crt/Main.hs +57/−0
- app/padic/Main.hs +70/−0
- app/partitions/Main.hs +22/−0
- app/pi_digits/Main.hs +43/−0
- app/poly_roots/Main.hs +180/−0
- app/primegen/Main.hs +55/−0
- app/qadic/Main.hs +94/−0
- app/radix/Main.hs +130/−0
- app/real_roots/FFI.hsc +27/−0
- app/real_roots/Functions.hs +143/−0
- app/real_roots/Main.hs +168/−0
- app/real_roots/Run.hs +152/−0
- app/real_roots/z_param.h +16/−0
- app/stirling_matrix/Main.hs +31/−0
- app/swinnerton_dyer_poly/Main.hs +69/−0
- app/taylor_integrals/Main.hs +137/−0
- app/zeta_zeros/Main.hs +160/−0
- docs/mj.png binary
- docs/out.png binary
- src/Lib.hs +6/−0
- test/Spec.hs +2/−0
+ CHANGELOG.md view
@@ -0,0 +1,11 @@+# Changelog for `Flint2-Examples`++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.1.0.0 - YYYY-MM-DD
+ Flint2-Examples.cabal view
@@ -0,0 +1,813 @@+cabal-version: 1.18++-- This file has been generated from package.yaml by hpack version 0.35.2.+--+-- see: https://github.com/sol/hpack++name: Flint2-Examples+version: 0.1.0.0+synopsis: Examples for the Flint2 library+description: Examples for Flint2 Haskell bindings for the flint library for number theory.+category: Math+homepage: https://github.com/githubuser/Flint2-Examples#readme+bug-reports: https://github.com/githubuser/Flint2-Examples/issues+author: monien+maintainer: hmonien@uni-bonn.de+copyright: 2023 Hartmut Monien+license: GPL-2+license-file: LICENSE+build-type: Simple+extra-source-files:+ README.md+ CHANGELOG.md+extra-doc-files:+ docs/mj.png+ docs/out.png++source-repository head+ type: git+ location: https://github.com/githubuser/Flint2-Examples++library+ exposed-modules:+ Lib+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ src+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable bernoulli+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/bernoulli+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , time >=1.12.2 && <1.13+ , timeit ==2.0.*+ default-language: Haskell2010++executable binet+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/binet+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable class_poly+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/class_poly+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , time >=1.12.2 && <1.13+ , timeit ==2.0.*+ default-language: Haskell2010++executable complex_plot+ main-is: Main.hs+ other-modules:+ ColorFunction+ Functions+ Paths_Flint2_Examples+ hs-source-dirs:+ app/complex_plot+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , JuicyPixels >=3.3.8 && <3.4+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , filepath >=1.4.2 && <1.5+ , gloss >=1.13.2 && <1.14+ , gloss-juicy >=0.2.3 && <0.3+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable crt+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/crt+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable delta_qexp+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/delta_qexp+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable dft+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/dft+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable elementary+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/elementary+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable expression+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/expression+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable factor_integer+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/factor_integer+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable fmpq_poly+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/fmpq_poly+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable fmpz_mod_poly+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/fmpz_mod_poly+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable fmpz_mpoly_factor+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/fmpz_mpoly_factor+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , time >=1.12.2 && <1.13+ , timeit ==2.0.*+ default-language: Haskell2010++executable fmpz_poly_factor_zassenhaus+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/fmpz_poly_factor_zassenhaus+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable fpwrap+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/fpwrap+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable fq_poly+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/fq_poly+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable function_benchmark+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/function_benchmark+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable hilbert_matrix+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/hilbert_matrix+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable hilbert_matrix_ca+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/hilbert_matrix_ca+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable integrals+ main-is: Main.hs+ other-modules:+ Integrands+ Paths_Flint2_Examples+ hs-source-dirs:+ app/integrals+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable keiper_li+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/keiper_li+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable l_central+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/l_values+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable l_values+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/l_values+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable logistic+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/logistic+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable machin+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/machin+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable multi_crt+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/multi_crt+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable padic+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/padic+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable partitions+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/partitions+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable pi_digits+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/pi_digits+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable poly_roots+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/poly_roots+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable primegen+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/primegen+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable qadic+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/qadic+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable radix+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/radix+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable real_roots+ main-is: Main.hs+ other-modules:+ FFI+ Functions+ Run+ Paths_Flint2_Examples+ hs-source-dirs:+ app/real_roots+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable stirling_matrix+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/stirling_matrix+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable swinnerton_dyer_poly+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/swinnerton_dyer_poly+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable taylor_integrals+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/taylor_integrals+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++executable zeta_zeros+ main-is: Main.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ app/zeta_zeros+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010++test-suite Flint2-Examples-test+ type: exitcode-stdio-1.0+ main-is: Spec.hs+ other-modules:+ Paths_Flint2_Examples+ hs-source-dirs:+ test+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ include-dirs:+ app/real_roots+ install-includes:+ app/real_roots/z_param.h+ build-depends:+ Flint2 >=0.1.0.4 && <0.2.0.0+ , Flint2-Examples+ , base >=4.7 && <5+ , containers >=0.6.7 && <0.7+ , mtl >=2.2.2 && <2.3+ , optparse-applicative >=0.18.1 && <0.19+ , timeit ==2.0.*+ default-language: Haskell2010
+ LICENSE view
@@ -0,0 +1,339 @@+ GNU GENERAL PUBLIC LICENSE+ Version 2, June 1991++ Copyright (C) 1989, 1991 Free Software Foundation, Inc.,+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA+ Everyone is permitted to copy and distribute verbatim copies+ of this license document, but changing it is not allowed.++ Preamble++ The licenses for most software are designed to take away your+freedom to share and change it. By contrast, the GNU General Public+License is intended to guarantee your freedom to share and change free+software--to make sure the software is free for all its users. This+General Public License applies to most of the Free Software+Foundation's software and to any other program whose authors commit to+using it. (Some other Free Software Foundation software is covered by+the GNU Lesser General Public License instead.) You can apply it to+your programs, too.++ When we speak of free software, we are referring to freedom, not+price. Our General Public Licenses are designed to make sure that you+have the freedom to distribute copies of free software (and charge for+this service if you wish), that you receive source code or can get it+if you want it, that you can change the software or use pieces of it+in new free programs; and that you know you can do these things.++ To protect your rights, we need to make restrictions that forbid+anyone to deny you these rights or to ask you to surrender the rights.+These restrictions translate to certain responsibilities for you if you+distribute copies of the software, or if you modify it.++ For example, if you distribute copies of such a program, whether+gratis or for a fee, you must give the recipients all the rights that+you have. You must make sure that they, too, receive or can get the+source code. And you must show them these terms so they know their+rights.++ We protect your rights with two steps: (1) copyright the software, and+(2) offer you this license which gives you legal permission to copy,+distribute and/or modify the software.++ Also, for each author's protection and ours, we want to make certain+that everyone understands that there is no warranty for this free+software. If the software is modified by someone else and passed on, we+want its recipients to know that what they have is not the original, so+that any problems introduced by others will not reflect on the original+authors' reputations.++ Finally, any free program is threatened constantly by software+patents. We wish to avoid the danger that redistributors of a free+program will individually obtain patent licenses, in effect making the+program proprietary. To prevent this, we have made it clear that any+patent must be licensed for everyone's free use or not licensed at all.++ The precise terms and conditions for copying, distribution and+modification follow.++ GNU GENERAL PUBLIC LICENSE+ TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION++ 0. This License applies to any program or other work which contains+a notice placed by the copyright holder saying it may be distributed+under the terms of this General Public License. The "Program", below,+refers to any such program or work, and a "work based on the Program"+means either the Program or any derivative work under copyright law:+that is to say, a work containing the Program or a portion of it,+either verbatim or with modifications and/or translated into another+language. (Hereinafter, translation is included without limitation in+the term "modification".) Each licensee is addressed as "you".++Activities other than copying, distribution and modification are not+covered by this License; they are outside its scope. The act of+running the Program is not restricted, and the output from the Program+is covered only if its contents constitute a work based on the+Program (independent of having been made by running the Program).+Whether that is true depends on what the Program does.++ 1. You may copy and distribute verbatim copies of the Program's+source code as you receive it, in any medium, provided that you+conspicuously and appropriately publish on each copy an appropriate+copyright notice and disclaimer of warranty; keep intact all the+notices that refer to this License and to the absence of any warranty;+and give any other recipients of the Program a copy of this License+along with the Program.++You may charge a fee for the physical act of transferring a copy, and+you may at your option offer warranty protection in exchange for a fee.++ 2. You may modify your copy or copies of the Program or any portion+of it, thus forming a work based on the Program, and copy and+distribute such modifications or work under the terms of Section 1+above, provided that you also meet all of these conditions:++ a) You must cause the modified files to carry prominent notices+ stating that you changed the files and the date of any change.++ b) You must cause any work that you distribute or publish, that in+ whole or in part contains or is derived from the Program or any+ part thereof, to be licensed as a whole at no charge to all third+ parties under the terms of this License.++ c) If the modified program normally reads commands interactively+ when run, you must cause it, when started running for such+ interactive use in the most ordinary way, to print or display an+ announcement including an appropriate copyright notice and a+ notice that there is no warranty (or else, saying that you provide+ a warranty) and that users may redistribute the program under+ these conditions, and telling the user how to view a copy of this+ License. (Exception: if the Program itself is interactive but+ does not normally print such an announcement, your work based on+ the Program is not required to print an announcement.)++These requirements apply to the modified work as a whole. If+identifiable sections of that work are not derived from the Program,+and can be reasonably considered independent and separate works in+themselves, then this License, and its terms, do not apply to those+sections when you distribute them as separate works. But when you+distribute the same sections as part of a whole which is a work based+on the Program, the distribution of the whole must be on the terms of+this License, whose permissions for other licensees extend to the+entire whole, and thus to each and every part regardless of who wrote it.++Thus, it is not the intent of this section to claim rights or contest+your rights to work written entirely by you; rather, the intent is to+exercise the right to control the distribution of derivative or+collective works based on the Program.++In addition, mere aggregation of another work not based on the Program+with the Program (or with a work based on the Program) on a volume of+a storage or distribution medium does not bring the other work under+the scope of this License.++ 3. You may copy and distribute the Program (or a work based on it,+under Section 2) in object code or executable form under the terms of+Sections 1 and 2 above provided that you also do one of the following:++ a) Accompany it with the complete corresponding machine-readable+ source code, which must be distributed under the terms of Sections+ 1 and 2 above on a medium customarily used for software interchange; or,++ b) Accompany it with a written offer, valid for at least three+ years, to give any third party, for a charge no more than your+ cost of physically performing source distribution, a complete+ machine-readable copy of the corresponding source code, to be+ distributed under the terms of Sections 1 and 2 above on a medium+ customarily used for software interchange; or,++ c) Accompany it with the information you received as to the offer+ to distribute corresponding source code. (This alternative is+ allowed only for noncommercial distribution and only if you+ received the program in object code or executable form with such+ an offer, in accord with Subsection b above.)++The source code for a work means the preferred form of the work for+making modifications to it. For an executable work, complete source+code means all the source code for all modules it contains, plus any+associated interface definition files, plus the scripts used to+control compilation and installation of the executable. However, as a+special exception, the source code distributed need not include+anything that is normally distributed (in either source or binary+form) with the major components (compiler, kernel, and so on) of the+operating system on which the executable runs, unless that component+itself accompanies the executable.++If distribution of executable or object code is made by offering+access to copy from a designated place, then offering equivalent+access to copy the source code from the same place counts as+distribution of the source code, even though third parties are not+compelled to copy the source along with the object code.++ 4. You may not copy, modify, sublicense, or distribute the Program+except as expressly provided under this License. Any attempt+otherwise to copy, modify, sublicense or distribute the Program is+void, and will automatically terminate your rights under this License.+However, parties who have received copies, or rights, from you under+this License will not have their licenses terminated so long as such+parties remain in full compliance.++ 5. You are not required to accept this License, since you have not+signed it. However, nothing else grants you permission to modify or+distribute the Program or its derivative works. These actions are+prohibited by law if you do not accept this License. Therefore, by+modifying or distributing the Program (or any work based on the+Program), you indicate your acceptance of this License to do so, and+all its terms and conditions for copying, distributing or modifying+the Program or works based on it.++ 6. Each time you redistribute the Program (or any work based on the+Program), the recipient automatically receives a license from the+original licensor to copy, distribute or modify the Program subject to+these terms and conditions. You may not impose any further+restrictions on the recipients' exercise of the rights granted herein.+You are not responsible for enforcing compliance by third parties to+this License.++ 7. If, as a consequence of a court judgment or allegation of patent+infringement or for any other reason (not limited to patent issues),+conditions are imposed on you (whether by court order, agreement or+otherwise) that contradict the conditions of this License, they do not+excuse you from the conditions of this License. If you cannot+distribute so as to satisfy simultaneously your obligations under this+License and any other pertinent obligations, then as a consequence you+may not distribute the Program at all. For example, if a patent+license would not permit royalty-free redistribution of the Program by+all those who receive copies directly or indirectly through you, then+the only way you could satisfy both it and this License would be to+refrain entirely from distribution of the Program.++If any portion of this section is held invalid or unenforceable under+any particular circumstance, the balance of the section is intended to+apply and the section as a whole is intended to apply in other+circumstances.++It is not the purpose of this section to induce you to infringe any+patents or other property right claims or to contest validity of any+such claims; this section has the sole purpose of protecting the+integrity of the free software distribution system, which is+implemented by public license practices. Many people have made+generous contributions to the wide range of software distributed+through that system in reliance on consistent application of that+system; it is up to the author/donor to decide if he or she is willing+to distribute software through any other system and a licensee cannot+impose that choice.++This section is intended to make thoroughly clear what is believed to+be a consequence of the rest of this License.++ 8. If the distribution and/or use of the Program is restricted in+certain countries either by patents or by copyrighted interfaces, the+original copyright holder who places the Program under this License+may add an explicit geographical distribution limitation excluding+those countries, so that distribution is permitted only in or among+countries not thus excluded. In such case, this License incorporates+the limitation as if written in the body of this License.++ 9. The Free Software Foundation may publish revised and/or new versions+of the General Public License from time to time. Such new versions will+be similar in spirit to the present version, but may differ in detail to+address new problems or concerns.++Each version is given a distinguishing version number. If the Program+specifies a version number of this License which applies to it and "any+later version", you have the option of following the terms and conditions+either of that version or of any later version published by the Free+Software Foundation. If the Program does not specify a version number of+this License, you may choose any version ever published by the Free Software+Foundation.++ 10. If you wish to incorporate parts of the Program into other free+programs whose distribution conditions are different, write to the author+to ask for permission. For software which is copyrighted by the Free+Software Foundation, write to the Free Software Foundation; we sometimes+make exceptions for this. Our decision will be guided by the two goals+of preserving the free status of all derivatives of our free software and+of promoting the sharing and reuse of software generally.++ NO WARRANTY++ 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES+PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,+REPAIR OR CORRECTION.++ 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE+POSSIBILITY OF SUCH DAMAGES.++ END OF TERMS AND CONDITIONS++ How to Apply These Terms to Your New Programs++ If you develop a new program, and you want it to be of the greatest+possible use to the public, the best way to achieve this is to make it+free software which everyone can redistribute and change under these terms.++ To do so, attach the following notices to the program. It is safest+to attach them to the start of each source file to most effectively+convey the exclusion of warranty; and each file should have at least+the "copyright" line and a pointer to where the full notice is found.++ <one line to give the program's name and a brief idea of what it does.>+ Copyright (C) <year> <name of author>++ This program is free software; you can redistribute it and/or modify+ it under the terms of the GNU General Public License as published by+ the Free Software Foundation; either version 2 of the License, or+ (at your option) any later version.++ This program is distributed in the hope that it will be useful,+ but WITHOUT ANY WARRANTY; without even the implied warranty of+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the+ GNU General Public License for more details.++ You should have received a copy of the GNU General Public License along+ with this program; if not, write to the Free Software Foundation, Inc.,+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.++Also add information on how to contact you by electronic and paper mail.++If the program is interactive, make it output a short notice like this+when it starts in an interactive mode:++ Gnomovision version 69, Copyright (C) year name of author+ Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.+ This is free software, and you are welcome to redistribute it+ under certain conditions; type `show c' for details.++The hypothetical commands `show w' and `show c' should show the appropriate+parts of the General Public License. Of course, the commands you use may+be called something other than `show w' and `show c'; they could even be+mouse-clicks or menu items--whatever suits your program.++You should also get your employer (if you work as a programmer) or your+school, if any, to sign a "copyright disclaimer" for the program, if+necessary. Here is a sample; alter the names:++ Yoyodyne, Inc., hereby disclaims all copyright interest in the program+ `Gnomovision' (which makes passes at compilers) written by James Hacker.++ <signature of Ty Coon>, 1 April 1989+ Ty Coon, President of Vice++This General Public License does not permit incorporating your program into+proprietary programs. If your program is a subroutine library, you may+consider it more useful to permit linking proprietary applications with the+library. If this is what you want to do, use the GNU Lesser General+Public License instead of this License.
+ README.md view
@@ -0,0 +1,171 @@+++# Flint2-Examples++## Introduction++Examples for the [Flint2](https://hackage.haskell.org/package/Flint2)+Haskell wrapper of the [Flint](https://flintlib.org) library.++Most of these are directly "translated" from C to Haskell so they are+written using "imperative" constructs in Haskell to demonstrate the+use of the Haskell wrapper (not intended as examples of +elegant Haskell code). The structure of programs follows closely the+structure of the C code (with some exceptions:+e.g. swinnerton_dyer_poly uses a different algorithm, logistic uses a+monad transformer instead of a for loop ...)++There are basic benchmarks in some of the programs. In some cases+experimental support of multithreading is available. Checkout the+options by typing prog -h.++## Installation++- Install the C-library available from [Flint](https://flintlib.org). + There are packages available for various operating systems.++- Install the Haskell interface with++```bash+cabal install Flint2 --lib+```++- Check that your Cabal default local bin directory does not contain any+ conflicting names (see below for a list of binaries)!+- Install the examples with ++```bash+cabal install Flint2-Examples+```+'+## Source code++To study the source code download the code from Github with. More+documentation will hopefully available soon. Also check the flintlib site.++```bash+git clone https://github.com/monien/Flint2-Examples.git+```++## Quick start++A simple example would be the factorization of integers+(factor_integers). To find out what options are available use -h:++```bash+factor_integers -h+```++which prints++```bash+Factor integers.++Usage: factor_integer INTEGER [-t|--threads THREADS] [--timing]++ Factor integers.++Available options:+ INTEGER Integer given as expression (e.g. 2^64+1)+ -t,--threads THREADS number of threads+ --timing timing+ -h,--help Show this help text+```++Now typing ++```bash+factor_integer 2^256-1+```+returns+```bash+[(3,1),(5,1),(17,1),(257,1),(641,1),(65537,1),(274177,1),(6700417,1),(67280421310721,1),(59649589127497217,1),(5704689200685129054721,1)]+```++## A more advanced example: complex_plot+Try complex_plot by typing:++```bash+complex_plot+```+This will plot the phase of the *Klein invariant* in the upper half+plane. ++<p align="center">+ <img src="https://github.com/monien/Flint2-Examples/raw/main/docs/mj.png" />+</p>++The app complex_plot has many more options. Using the help option one obtains++```bash+Plotting special functions in the complex plane.++Usage: complex_plot [--xa XA] [--xb XB] [--ya YA] [--yb YB] [--width WIDTH] + [--height HEIGHT] [-c|--color-mode COLOR-MODE] + [-f|--function FUNCTION] [-o|--output IMAGE-FILE]++ plotting special functions.++Available options:+ -c,--color-mode COLOR-MODE+ possible values: 0 .. 6+ -f,--function FUNCTION possible values: agm, ai, barnesg, besseli, besselj,+ besselk, bessely, bi, digamma, ellipp, ellipsigma,+ ellipzeta, erf, fresnelc, fresnels, gamma, lgamma,+ modeta, modetaq, modj, modjq, modlambda, modlambdaq,+ zeta+ -o,--output IMAGE-FILE write output to IMAGE-FILE+ -h,--help Show this help tex+```++## Multithreading ++To use the experimental support of multithreading it is most+convenient to set the number of cores available using the environment+variable GHCRTS (e.g. for 10 cores)++```bash+GHCRTS='-N10'+export GHCRTS+```++## List of available examples++- bernoulli+- binet+- class_poly+- complex_plot+- crt+- delta_qexp+- dft+- elementary+- expression+- factor_integer+- fmpq_poly+- fmpz_mod_poly+- fmpz_mpoly_factor+- fmpz_poly_factor_zassenhaus+- fpwrap+- fq_poly+- function_benchmark+- hilbert_matrix+- hilbert_matrix_ca+- integrals+- keiper_li+- l_central+- l_values+- logistic+- machin+- multi_crt+- padic+- partitions+- pi_digits+- poly_roots+- primegen+- qadic+- radix+- real_roots+- stirling_matrix+- swinnerton_dyer_poly+- taylor_integrals+- zeta_zeros
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ app/bernoulli/Main.hs view
@@ -0,0 +1,41 @@+import Options.Applicative+import Control.Monad+import Foreign.C.Types+import Data.Time+import System.TimeIt++import Data.Number.Flint++main = timeItNamed "time"+ $ run =<< customExecParser (prefs showHelpOnEmpty) opts where+ desc = "Calculates bernoulli numbers."+ opts = info (parameters <**> helper) (+ fullDesc+ <> progDesc desc+ <> header desc)++run params@(Parameters n display num_threads) = do+ flint_set_num_threads num_threads+ x <- newFmpq+ withFmpq x $ \x -> bernoulli_fmpq_ui x n+ when (n <= 100 || display) $ print x+ +data Parameters = Parameters {+ n :: CULong+ , display :: Bool+ , num_threads :: CInt+} deriving Show++parameters :: Parser Parameters+parameters = Parameters+ <$> argument auto (+ help "n"+ <> metavar "n")+ <*> switch (+ help "display value"+ <> short 'd')+ <*> option auto (+ help "number of threads"+ <> long "threads"+ <> value 1+ <> metavar "threads")
+ app/binet/Main.hs view
@@ -0,0 +1,60 @@+import Options.Applicative+import Control.Monad+import System.TimeIt+import Foreign.C.Types++import Data.Number.Flint++main = timeItNamed "binet"+ $ run =<< customExecParser (prefs showHelpOnEmpty) opts where+ desc = "This program computes the n-th Fibonacci number using Binet’s formula"+ opts = info (parameters <**> helper) (+ fullDesc+ <> progDesc desc+ <> header desc)++run params@(Parameters n limit) = do++ ctx <- newCaCtx++ [sqrt5, phi, psi, t, u] <- replicateM 5 (newCa ctx)++ withCaCtx ctx $ \ctx -> do+ case limit of+ Just prec -> ca_ctx_set_option ctx ca_opt_prec_limit prec+ _ -> return ()+ withCa sqrt5 $ \sqrt5 -> do+ withCa phi $ \phi -> do+ withCa psi $ \psi -> do+ withCa t $ \t -> do+ withCa u $ \u -> do+ ca_sqrt_ui sqrt5 5 ctx+ ca_add_ui phi sqrt5 1 ctx+ ca_div_ui phi phi 2 ctx+ ca_ui_sub psi 1 phi ctx+ withFmpz n $ \n -> do+ ca_pow_fmpz t phi n ctx+ ca_pow_fmpz u psi n ctx+ ca_sub t t u ctx+ ca_div t t sqrt5 ctx+ ca_print t ctx; putStr "\n"++-- Parser ----------------------------------------------------------------------++data Parameters = Parameters {+ n :: Fmpz+ , limit :: Maybe CLong+} deriving Show++parameters :: Parser Parameters+parameters = Parameters+ <$> argument auto (+ help "nth power"+ <> metavar "n")+ <*> optional optionLimit++optionLimit = option auto (+ help "limiting precision"+ <> long "limit"+ <> metavar "limit")+
+ app/class_poly/Main.hs view
@@ -0,0 +1,56 @@+import System.Environment+import System.TimeIt++import Options.Applicative+import Control.Monad+import Control.Monad.State+import Text.Read (readMaybe)++import Options.Applicative++import Foreign.C.Types+import Foreign.Ptr++import Data.Number.Flint++main = timeItNamed "time"+ $ run =<< customExecParser (prefs showHelpOnEmpty) opts where+ desc = "Calculates class polynomial for negative discriminant."+ opts = info (parameters <**> helper) (+ fullDesc+ <> progDesc desc+ <> header desc)++run (Parameters d num_threads) = do+ flint_set_num_threads num_threads+ res <- newFmpzPoly+ withFmpzPoly res $ \res -> acb_modular_hilbert_class_poly res d+ when (abs d <= 100) $ print res++data Parameters = Parameters {+ d :: CLong+ , num_threads :: CInt+} deriving Show++parameters :: Parser Parameters+parameters = Parameters+ <$> argument discriminant (+ help "absolute value of (-D)"+ <> metavar "D")+ <*> option auto (+ help "number of threads"+ <> long "threads"+ <> value 1+ <> metavar "threads")++discriminant :: (Read a, Integral a, Show a) => ReadM a+discriminant = eitherReader $ \s -> do+ let d = negate $ read s+ if d < 0 then+ if d `mod` 4 == 1 || d `mod` 4 == 0 then+ Right d+ else+ Left $ "discriminant D (=" ++ show d ++ ") `mod` 4 /= 0, 1."+ else+ Left $ "discriminant D (=" ++ show d ++ ") > 0."+
+ app/complex_plot/ColorFunction.hs view
@@ -0,0 +1,113 @@+module ColorFunction where++import Data.List (findIndex)+import Data.Complex++colorFunction mode z+ | not (isInfinite r || isNaN r) = colorFunction' mode z+ | otherwise = (0.5, 0.5, 0.5)+ where r = magnitude z+ +colorFunction' 0 z = toRGB (h, l, s) where+ phi = phase z+ tmp = (phi + pi)/(2*pi) + 1/2+ h = tmp - fromIntegral (floor tmp)+ l | log (magnitude z) > 200 = 1+ | log (magnitude z) < -200 = 0+ | otherwise = 1 - 1 / (1 + (magnitude z) ** 0.2)+ s = 0.8++colorFunction' 1 z = (r, g, b) where+ phi = phase z+ h = max (min (phi/pi) 1) (-1)+ Just j = findIndex ((>h) . head) blueOrangeColors+ [ha, ra, ga, ba] = blueOrangeColors !! (j-1)+ [hb, rb, gb, bb] = blueOrangeColors !! j+ s = (h - ha) / (hb - ha)+ r = ra + (rb - ra) * s+ g = ga + (gb - ga) * s+ b = ba + (bb - ba) * s+ blueOrangeColors = [+ [-1.0, 0.0, 0.0, 0.0],+ [-0.95, 0.1, 0.2, 0.5],+ [-0.5, 0.0, 0.5, 1.0],+ [-0.05, 0.4, 0.8, 0.8],+ [ 0.0, 1.0, 1.0, 1.0],+ [ 0.05, 1.0, 0.9, 0.3],+ [ 0.5, 0.9, 0.5, 0.0],+ [ 0.95, 0.7, 0.1, 0.0],+ [ 1.0, 0.0, 0.0, 0.0],+ [ 2.0, 0.0, 0.0, 0.0]+ ]++colorFunction' mode z+ | mode == 2 = (r, g, b)+ | mode == 3 = mix 0.0 (-0.5) 0.2 0.0 0.0 (-0.1) 0.0 (-1.0) (-0.2)+ | mode == 4 = mix 0.0 (-0.5) 0.2 0.0 0.5 (-0.1) 0.0 (-0.3) (-1.0)+ | mode == 5 = mix 0.0 (-0.5) (-1.0) 0.0 (-0.1) (-0.67) 0.0 (-0.55) (-0.12)+ | mode == 6 = mix 0.86 0.0 0.13 0.57 0.19 (-0.52) 0.31 (-0.30) (-0.94)+ where+ (r1, g1, b1) = colorFunction 0 z+ (r2, g2, b2) = colorFunction 1 z+ f x y = blend x (clamp (dodge x y))+ [r, g, b] = zipWith f [r1, g1, b1] [r2, g2, b2] + mix = balance (r, g, b)+ +toRGB (h, l, s)+ | s == 0 = (l, l, l)+ | otherwise = (vv m1 m2 (h+1/3), vv m1 m2 h, vv m1 m2 (h-1/3))+ where+ m2 = if l <= 0.5 then l * (1 + s) else l + s - l*s+ m1 = 2*l - m2 + vv m1 m2 hue+ | 6*h < 1 = m1 + (m2 - m1)*h*6+ | 2*h < 1 = m2+ | 3*h < 2 = m1 + (m2 - m1)*(2/3 - h)*6+ | otherwise = m1+ where h = hue - fromIntegral (floor hue)++toHLS (r, g, b)+ | hi == lo = (0, l, 0)+ | l <= 0.5 = (h', l, d / (hi + lo))+ | otherwise = (h', l, d / (2 - hi - lo))+ where+ hi = max (max r g) b+ lo = min (min r g) b+ l = 0.5 * (lo + hi)+ d = hi -lo+ d' = if d /= 0 then d else 1+ h | r == hi = (g - b) / d'+ | g == hi = (b - r) / d' + 2+ | otherwise = (r - g) / d' + 4+ h' = h / 6+ h'' = if h' < 0 then h' + 1 else h'++clamp x = max 0 (min x 1)+blend x y = (x + y)/2+dodge a b = a / ((1 - b) + 1/256)++balanceChannel value x (shadows, midtones, highlights) = value' where+ [a, b, scale] = [1/4, 1/3, 7/10]+ shadows' = clamp ((x - b) / (-a) + 0.5) * scale * shadows+ midtones' = clamp ((x - b) / a + 0.5) *+ clamp ((x + b - 1.0) / (-a) + 0.5) * scale * midtones+ highlights' = clamp ((x + b - 1.0) / a + 0.5) * scale * highlights+ value' = clamp (value + shadows' + midtones' + highlights')++balance (r, g, b) u' u'' u''' v' v'' v''' w' w'' w''' = (r'', g'', b'') where+ [u, v, w] = [(u', u'', u'''), (v', v'', v'''), (w', w'', w''')]+ [r', g', b'] = zipWith (\c x -> balanceChannel c c x) [r, g, b] [u, v, w]+ (h, l, s ) = toHLS (r, g, b )+ (h', l', s') = toHLS (r', g', b')+ (r'', g'', b'') = toRGB (h', l, s')++hlsFunction 0 z = (h, l, s) where+ x = realPart z+ y = imagPart z+ phi = atan2 y x+ tmp = (phi + pi)/(2*pi) + 1/2+ h = tmp - fromIntegral (floor tmp)+ l | log(magnitude z) > 138 = 1+ | log(magnitude z) < -138 = 0+ | otherwise = 1 - 1 / (1 + (magnitude z) ** 0.2)+ s = 0.8
+ app/complex_plot/Functions.hs view
@@ -0,0 +1,193 @@+module Functions where++import Control.Monad+import qualified Data.Map as Map+import Foreign.Ptr (Ptr, nullPtr)+import Foreign.C.Types (CLong)+import Data.Complex+import System.IO.Unsafe++import Data.Number.Flint.Arb+import Data.Number.Flint.Arb.Arf+import Data.Number.Flint.Arb.Hypgeom+import Data.Number.Flint.Acb+import Data.Number.Flint.Acb.Hypgeom+import Data.Number.Flint.Acb.Modular+import Data.Number.Flint.Acb.Elliptic++functions = Map.fromList+ [ ("gamma" , gamma)+ , ("digamma" , digamma)+ , ("lgamma" , lgamma)+ , ("zeta" , zeta)+ , ("erf" , erf)+ , ("ai" , ai)+ , ("bi" , bi)+ , ("besselj" , besselj)+ , ("bessely" , bessely)+ , ("besseli" , besseli)+ , ("besselk" , besselk)+ , ("modj" , modj)+ , ("modjq" , modjq)+ , ("modeta" , modeta)+ , ("modetaq" , modetaq)+ , ("modlambda" , modlambda)+ , ("modlambdaq", modlambdaq)+ , ("ellipp" , ellipp)+ , ("ellipzeta" , ellipzeta)+ , ("ellipsigma", ellipsigma)+ , ("barnesg" , barnesg)+ , ("agm" , agm)+ , ("fresnels" , fresnels)+ , ("fresnelc" , fresnelc)]++gamma = acb_gamma+digamma = acb_digamma+lgamma = acb_lgamma+zeta = acb_zeta++agm = acb_agm1++ai res =+ acb_hypgeom_airy res+ (nullPtr :: Ptr CAcb)+ (nullPtr :: Ptr CAcb)+ (nullPtr :: Ptr CAcb)++bi res =+ acb_hypgeom_airy+ (nullPtr :: Ptr CAcb)+ (nullPtr :: Ptr CAcb)+ res+ (nullPtr :: Ptr CAcb)++barnesg = acb_barnes_g++besselj res z prec = do + withNewAcb $ \nu -> acb_hypgeom_bessel_j res nu z prec+ return ()++bessely res z prec = do + withNewAcb $ \nu -> acb_hypgeom_bessel_y res nu z prec+ return ()++besseli res z prec = do + withNewAcb $ \nu -> acb_hypgeom_bessel_i res nu z prec+ return ()++besselk res z prec = do + withNewAcb $ \nu -> acb_hypgeom_bessel_k res nu z prec+ return ()++erf = acb_hypgeom_erf++modj res z prec = do+ acb_modular_j res z prec+ acb_div_ui res res 1728 prec++modjq res z prec = do+ withNewAcb $ \t -> do + acb_log res z prec+ acb_const_pi t prec+ acb_div res res t prec+ acb_mul_2exp_si res res (-1)+ acb_div_onei res res+ acb_modular_j res res prec+ acb_div_ui res res 1728 prec+ return ()+ +modeta = acb_modular_eta+modetaq res z prec = do+ withNewAcb $ \t -> do + acb_log res z prec+ acb_log res z prec+ acb_const_pi t prec+ acb_div res res t prec+ acb_mul_2exp_si res res (-1)+ acb_div_onei res res+ acb_modular_eta res res prec+ return ()+ +modlambda = acb_modular_lambda+modlambdaq res z prec = do+ withNewAcb $ \t -> do + acb_log res z prec+ acb_log res z prec+ acb_const_pi t prec+ acb_div res res t prec+ acb_mul_2exp_si res res (-1)+ acb_div_onei res res+ acb_modular_lambda res res prec+ return ()+ +ellipp res z prec = do+ acb_onei res+ acb_elliptic_p res z res prec++ellipzeta res z prec = do+ acb_onei res+ acb_elliptic_zeta res z res prec++ellipsigma res z prec = do+ acb_onei res+ acb_elliptic_sigma res z res prec++fresnels res z prec = acb_hypgeom_fresnel res (nullPtr :: Ptr CAcb) z 0 prec+fresnelc res z prec = acb_hypgeom_fresnel (nullPtr :: Ptr CAcb) res z 0 prec++eval zp prec (xa, xb, xnum) (ya, yb, ynum) f i j = do+ withNewArb $ \re -> do+ withNewArb $ \im -> do+ withNewArf $ \x -> do+ withNewArf $ \xap -> do+ withNewArf $ \xbp -> do+ withNewArf $ \y -> do+ withNewArf $ \yap -> do+ withNewArf $ \ybp -> do+ arf_set_d xap (realToFrac xa)+ arf_set_d xbp (realToFrac xb)+ arf_sub x xbp xap prec arf_rnd_down+ arf_mul_ui x x (fromIntegral i) prec arf_rnd_down+ arf_div_ui x x (fromIntegral xnum-1) prec arf_rnd_down+ arf_add x x xap prec arf_rnd_down+ arb_set_arf re x+ arf_set_d yap (realToFrac ya)+ arf_set_d ybp (realToFrac yb)+ arf_sub y ybp yap prec arf_rnd_down+ arf_mul_ui y y (fromIntegral j) prec arf_rnd_down+ arf_div_ui y y (fromIntegral ynum-1) prec arf_rnd_down+ arf_add y y yap prec arf_rnd_down+ arb_set_arf im y+ acb_set_arb_arb zp re im+ f zp zp prec+ return ()++evalSafe (xa, xb, xnum) (ya, yb, ynum) f i j = unsafePerformIO $ do+ let iter prec = do+ (result, _) <- withNewAcb $ \value -> do+ eval value prec (xa, xb, xnum) (ya, yb, ynum) f i j+ finite <- acb_is_finite value+ bits <- acb_rel_accuracy_bits value+ when ((finite /= 1 || bits < 4) && prec < 500) $ do+ _ <- iter (prec + 30)+ return ()+ return value+ return result+ arb_get_d x = do+ d <- arf_get_d (arb_midref x) arf_rnd_near+ return d+ result <- iter 30+ (_, (_, (_, functionValue)))+ <- withAcb result $ \result -> do+ withNewArb $ \x -> do+ withNewArb $ \y-> do+ acb_get_real x result+ acb_get_imag y result+ dx <- arb_get_d x+ dy <- arb_get_d y+ return $ dx :+ dy+ return functionValue++ + +
+ app/complex_plot/Main.hs view
@@ -0,0 +1,134 @@+import System.FilePath++import Data.List (intercalate)+import qualified Data.Map as Map++import Control.Monad++import Options.Applicative++import Codec.Picture+import Graphics.Gloss +import Graphics.Gloss.Juicy++import ColorFunction+import Functions++main :: IO ()+main = run =<< execParser opts where+ opts = info (parameters <**> helper) + ( fullDesc+ <> progDesc "\nplotting special functions.\n"+ <> header "Plotting special functions in the complex plane.")++run :: Parameters -> IO ()+run (Parameters xa xb ya yb w h colorMode f imgFile num_threads) = do+ when (colorMode < 0 || colorMode > 6) $ do error "colorMode not available."+ case Map.lookup f functions of + Just g -> do let u i j = evalSafe (xa, xb, w) (ya, yb, h) g i (h-j)+ v i j = rgba colorMode (u i j)+ img = ImageRGBA8 (generateImage v w h)+ case imgFile of+ Just out -> saveImage out img+ _ -> return ()+ case fromDynamicImage img of+ Just picture -> display (InWindow "arb plot" (w, h) (0, 0))+ white picture+ _ -> putStrLn "could not display picture."+ _ -> putStrLn $ "function '" ++ f ++ "' not available."++rgba colorMode z = PixelRGBA8 r' g' b' alpha where+ (r, g, b) = colorFunction colorMode z+ alpha = if (r, g, b) /= (0.5, 0.5, 0.5) then 255 else 0+ [r', g', b'] = map (min 255 . floor . (*255)) [r, g, b]++saveImage out img = do+ if hasExtension out then do+ let f = case takeExtension out of+ ".bmp" -> Just saveBmpImage+ ".png" -> Just savePngImage+ ".jpg" -> Just (saveJpgImage 100)+ ".tiff" -> Just saveTiffImage+ ".hdr" -> Just saveRadianceImage+ _ -> Nothing+ case f of+ Just f -> do f out img+ return ()+ _ -> putStrLn "Could not save image: unknown file extension.\n\+ \Use bmp, jpg, png, tiff or hdr as file extension." + else+ putStrLn "Could not save image: no file extension.\n\+ \Use bmp, jpg, png, tiff or hdr as file extension." + +data Range = Range Double Double Double Double deriving Show+data Size = Size Int Int deriving Show++data Parameters = Parameters {+ xa :: Double+ , xb :: Double+ , ya :: Double+ , yb :: Double+ , with :: Int+ , height :: Int+ , colorMode :: Int+ , fun :: String+ , imgFile :: Maybe String+ , num_threads :: Int+ } deriving Show++parameters :: Parser Parameters+parameters = Parameters+ <$> option auto (+ long "xa" <>+ value (-1.0) <>+ showDefault <>+ metavar "XA")+ <*> option auto (+ long "xb" <>+ value 1.0 <>+ showDefault <>+ metavar "XB")+ <*> option auto (+ long "ya" <>+ value 0.0 <>+ showDefault <>+ metavar "YA") + <*> option auto (+ long "yb" <>+ value 2 <>+ showDefault <>+ metavar "YB")+ <*> option auto (+ long "width" <>+ -- short 'w' <>+ value 512 <>+ showDefault <>+ metavar "WIDTH")+ <*> option auto (+ long "height" <>+ -- short 'h' <>+ value 512 <>+ showDefault <>+ metavar "HEIGHT")+ <*> option auto (+ long "color-mode" <>+ short 'c' <>+ value 0 <>+ help ("possible values: 0 .. 6") <>+ metavar "COLOR-MODE")+ <*> strOption (+ short 'f' <>+ long "function" <>+ value "modj" <> + help ("possible values: " ++ intercalate ", " (Map.keys functions)) <>+ metavar "FUNCTION")+ <*> optional (strOption (+ long "output" <>+ short 'o' <>+ metavar "IMAGE-FILE" <>+ help "write output to IMAGE-FILE"))+ <*> option auto (+ help "number of threads"+ <> long "threads"+ <> value 1+ <> metavar "THREADS")
+ app/crt/Main.hs view
@@ -0,0 +1,53 @@+import System.Environment+import System.TimeIt++import Options.Applicative+import Control.Monad+import Control.Monad.State+import Text.Read (readMaybe)++import Foreign.C.Types+import Foreign.Ptr++import Data.Number.Flint++main = do+ prog <- getProgName+ args <- getArgs+ case args of+ [arg0] -> case readMaybe arg0 :: Maybe Fmpz of+ Just x-> run x+ _ -> putStrLn $ "usage: " ++ prog ++ " <integer>."+ _ -> putStrLn $ "usage: " ++ prog ++ " <integer>."++run x = do+ print x+ withFmpz x $ \x -> do+ bits <- fmpz_bits x+ withNewFmpz $ \y -> do+ fmpz_zero y+ withNewFmpz $ \prod -> do+ fmpz_one prod+ _ <- runStateT (next x (bits + 2)) (0, y, prod, 0)+ return ()+ return ()++next :: Ptr CFmpz -> CFBitCnt+ -> StateT (Int, Ptr CFmpz, Ptr CFmpz, CULong) IO ()+next x bit_bound = do+ (i, y, prod, prime) <- get+ (prod', y', prime', bits) <- liftIO $ do+ prime' <- n_nextprime prime 0+ res <- fmpz_fdiv_ui x prime'+ fmpz_CRT_ui y y prod res prime' 1+ putStr $ "residue mod " ++ show prime' ++ " = " ++ show res+ putStr $ "; reconstruction = "+ fmpz_print y+ putStr "\n"+ fmpz_mul_ui prod prod prime'+ bits <- fmpz_bits prod+ return (prod, y, prime', bits)+ put (succ i, y', prod', prime')+ when (bits < bit_bound) $ next x bit_bound+ +
+ app/delta_qexp/Main.hs view
@@ -0,0 +1,43 @@+import Options.Applicative+import Control.Monad+import Control.Monad.State++import Foreign.C.Types+import Foreign.Ptr++import Data.Number.Flint++main = run =<< customExecParser (prefs showHelpOnEmpty) opts where+ desc = "Compute nth coefficient of the q-expansion of the weight \+ \12 cusp form."+ opts = info (parameters <**> helper) (+ fullDesc+ <> progDesc desc+ <> header desc)++run params@(Parameters n) = do+ withNewFmpz $ \k -> do+ fmpz_set_si k n+ withNewFmpz $ \c -> do+ arith_ramanujan_tau c k+ putStr $ "Coefficient of q^" ++ show n ++ " is "+ fmpz_print c+ putStr "\n"+ +data Parameters = Parameters {+ n :: CLong+} deriving Show++parameters :: Parser Parameters+parameters = Parameters+ <$> argument pos (+ help "positive integer n"+ <> metavar "n")++pos :: (Read a, Integral a) => ReadM a+pos = eitherReader $ \s -> do+ let result = read s+ if result >= 0 then + Right result+ else+ Left "expected positive number"
+ app/dft/Main.hs view
@@ -0,0 +1,257 @@+import System.TimeIt++import Foreign.Ptr+import Foreign.ForeignPtr+import Foreign.C.Types+import Foreign.Marshal.Array++import Control.Monad+import Options.Applicative+import Text.ParserCombinators.ReadP hiding (option)++import Data.Char+import Data.Number.Flint++main = run =<< customExecParser (prefs showHelpOnEmpty) opts where+ opts = info (options <**> helper) (+ fullDesc+ <> progDesc "Discrete Fourier Transform for functions [0..5]."+ <> header "Calculate Discrete Fourier Transform.")++run opts@(Options n verbosity range qqbar_limit nogb timing num_threads) = do+ let desc = [ "x_k = k + 2"+ , "x_k = sqrt(k + 2)"+ , "x_k = log(k + 2)"+ , "x_k = exp(2 pi i / (k + 2))"+ , "x_k = 1 / (1 + (k + 2) pi)"+ , "x_k = 1 / (1 + sqrt(k + 2) pi)"+ ]+ Range (start, end) = range+ gb = if nogb then 0 else 1+ ctx <- newCaCtx+ withNewCa ctx $ \t -> do+ withCaCtx ctx $ \ctx -> do+ forM_ [start .. end] $ \j -> do+ putStrLn $ desc !! j+ putStr "\n"+ if timing then do+ timeIt $ benchmark_DFT n j verbosity qqbar_limit gb t ctx+ else do+ benchmark_DFT n j verbosity qqbar_limit gb t ctx++data Options = Options {+ n :: CLong+ , verbosity :: Int+ , range :: Range+ , limit :: CLong+ , nogb :: Bool+ , doTiming :: Bool+ , num_threads :: CInt+} deriving Show++options :: Parser Options+options = Options+ <$> argument pos (+ help "n"+ <> value 2+ <> metavar "n")+ <*> option pos (+ help "verbosity 0 .. "+ <> long "verbosity"+ <> short 'v'+ <> value 0+ <> metavar "verbosity")+ <*> option (rng 0 5) (+ help "select function to be transformed [0..5]"+ <> long "range"+ <> short 'r'+ <> metavar "range")+ <*> option pos (+ help "limit"+ <> long "limit"+ <> value 0 + <> metavar "limit")+ <*> switch (+ help "do not use Groebner bases"+ <> long "nogb")+ <*> switch (+ help "timing"+ <> long "timing"+ <> short 't')+ <*> option pos (+ help "number of threads"+ <> long "threads"+ <> value 1+ <> metavar "threads")+ +-- ReadM parsers ---------------------------------------------------------------++rng :: Int -> Int -> ReadM Range+rng start end = eitherReader $ \s -> do+ let result@(Range (a, b)) = read s :: Range+ if start <= a && a <= b && b <= end then+ Right result+ else+ Left $ "Could not parse range " ++ s+ +pos :: (Read a, Integral a) => ReadM a+pos = eitherReader $ \s -> do+ let result = read s+ if result >= 0 then + Right result+ else+ Left "expected positive number"++between a b = eitherReader $ \s -> do+ let result = read s+ if a <= result && result <= b then + Right result+ else+ Left $ "expected number in range [" ++ show a ++ " .. " ++ show b ++ "]."++-- instances Range -------------------------------------------------------------++newtype Range = Range (Int, Int)++instance Read Range where+ readsPrec _ = readP_to_S (parseRange <|> parseIndex)++instance Show Range where+ show (Range (a, b)) = "[" ++ show a ++ ":" ++ show b ++ "]"++parseRange :: ReadP Range +parseRange = do+ a <- read <$> munch1 isNumber+ char ':'+ b <- read <$> munch1 isNumber+ return $ Range (a, b)++parseIndex :: ReadP Range+parseIndex = do+ a <- read <$> munch1 isNumber+ return $ Range (a, a)+ +--------------------------------------------------------------------------------++benchmark_DFT n input verbose qqbar_limit gb t ctx = do++ let (.+.) x y = x `advancePtr` (fromIntegral y)+ + x <- _ca_vec_init n ctx+ x' <- _ca_vec_init n ctx+ y <- _ca_vec_init n ctx+ w <- _ca_vec_init (2*n) ctx++ -- ca_ctx_set_option ctx ca_opt_print_flags ca_print_debug+ -- ca_ctx_set_option ctx ca_opt_verbose 1+ + ca_ctx_set_option ctx ca_opt_use_groebner gb+ + when (qqbar_limit /= 0) $ do+ ca_ctx_set_option ctx ca_opt_qqbar_deg_limit qqbar_limit++ -- Construct input vector++ when (verbose > 0) $ do putStr "[x] =\n"+ + forM_ [0 .. n - 1] $ \i -> do+ let xi = x .+. i+ case input of+ 0 -> do+ ca_set_si xi (i + 2) ctx+ 1 -> do+ ca_set_si xi (i + 2) ctx+ ca_sqrt xi xi ctx+ 2 -> do+ ca_set_si xi (i + 2) ctx+ ca_log xi xi ctx+ 3 -> do+ ca_pi_i xi ctx+ ca_mul_ui xi xi 2 ctx+ ca_div_si xi xi (i + 2) ctx+ ca_exp xi xi ctx+ 4 -> do+ ca_pi xi ctx+ ca_mul_si xi xi (i + 2) ctx+ ca_add_ui xi xi 1 ctx+ ca_inv xi xi ctx+ 5 -> do+ ca_pi xi ctx+ ca_sqrt_ui w (fromIntegral i + 2) ctx+ ca_mul xi xi w ctx+ ca_add_ui xi xi 1 ctx+ ca_inv xi xi ctx + when (verbose > 0) $ do+ ca_print xi ctx+ putStr "\n"++ -- construct roots of unity++ when (verbose > 1) $ do putStrLn "\n[w] = "+ + forM_ [0 .. 2 * n - 1] $ \i -> do+ let wi = w .+. i+ if i == 0 then do+ ca_one wi ctx+ else do+ if i == 1 then do+ ca_pi_i wi ctx+ ca_mul_ui wi wi 2 ctx+ ca_div_si wi wi n ctx+ ca_exp wi wi ctx+ else do+ ca_mul wi (w .+. (i - 1)) (w .+. 1) ctx+ when (verbose > 1) $ do+ ca_print wi ctx; putStr "\n"+ + -- forward dft: x -> x'++ when (verbose > 2) $ do putStrLn "\nDFT([x]) = "++ forM_ [0 .. n - 1] $ \k -> do+ let x'k = x' .+. k+ ca_zero x'k ctx+ forM_ [0 .. n - 1] $ \j -> do+ let xj = x .+. j+ wj = w .+. (((2 * n - k) * j) `mod` (2 * n))+ ca_mul t xj wj ctx+ ca_add x'k x'k t ctx+ when (verbose > 2) $ do+ ca_print x'k ctx; putStr "\n"++ -- inverse dft++ when (verbose > 2) $ do putStrLn "\nIDFT(DFT([x])) ="++ forM_ [0 .. n - 1] $ \k -> do+ let yk = y `advancePtr` (fromIntegral k)+ ca_zero yk ctx+ forM_ [0 .. n - 1] $ \j -> do+ let x'j = x' .+. j+ wj = w .+. ((k * j) `mod` (2 * n))+ ca_mul t x'j wj ctx+ ca_add yk yk t ctx+ ca_div_si yk yk n ctx+ when (verbose > 2) $ do+ ca_print yk ctx; putStr "\n"++ when (verbose > 0) $ do putStrLn "\n[x] - IDFT(DFT([x])) ="++ forM_ [0 .. n - 1] $ \k -> do+ let xk = x .+. k+ yk = y .+. k+ ca_sub t xk yk ctx+ is_zero <- ca_check_is_zero t ctx+ + when (verbose > 0) $ do+ ca_print t ctx+ putStrLn $ " (= 0 " ++ show is_zero ++ ")"+ when (is_zero /= t_true) $ do error "Failed to prove equality!"++ when (verbose > 0) $ do putStr "\n"++ _ca_vec_clear x n ctx+ _ca_vec_clear x' n ctx+ _ca_vec_clear y n ctx+ _ca_vec_clear w (2 * n) ctx+
+ app/elementary/Main.hs view
@@ -0,0 +1,244 @@+import System.TimeIt+import Foreign.ForeignPtr++import Data.Number.Flint++main = do++ ctx <- newCaCtx++ withNewCa ctx $ \x -> do+ withNewCa ctx $ \y -> do+ withNewCa ctx $ \z -> do+ withCaCtx ctx $ \ctx -> do+ timeIt $ run x y z ctx+ return ()+ +run x y z ctx = do++ putStrLn $ ">>> Exp(Pi*I) + 1"+ ca_pi_i x ctx+ ca_exp x x ctx+ ca_add_ui x x 1 ctx+ putStr "\n"; ca_print x ctx; putStr "\n\n"++ putStrLn $ ">>> Log(-1) / (Pi*I)"+ ca_set_si x (-1) ctx+ ca_log x x ctx+ ca_pi_i y ctx+ ca_div x x y ctx+ putStr "\n"; ca_print x ctx; putStr "\n\n"++ putStrLn $ ">>> Log(-I) / (Pi*I)"+ ca_neg_i x ctx+ ca_log x x ctx+ ca_pi_i y ctx+ ca_div x x y ctx+ putStr "\n"; ca_print x ctx; putStr "\n\n"++ putStrLn $ ">>> Log(1 / 10^123) / Log(100)"+ ca_set_ui x 10 ctx+ ca_pow_ui x x 123 ctx+ ca_inv x x ctx+ ca_log x x ctx+ ca_set_ui y 100 ctx+ ca_log y y ctx+ ca_div x x y ctx+ putStr "\n"; ca_print x ctx; putStr "\n\n"++ putStrLn $ ">>> Log(1 + Sqrt(2)) / Log(3 + 2*Sqrt(2))"+ ca_sqrt_ui x 2 ctx+ ca_add_ui x x 1 ctx+ ca_log x x ctx+ ca_sqrt_ui y 2 ctx+ ca_mul_ui y y 2 ctx+ ca_add_ui y y 3 ctx+ ca_log y y ctx+ ca_div x x y ctx+ putStr "\n"; ca_print x ctx; putStr "\n\n"++ putStrLn $ ">>> Sqrt(2)*Sqrt(3) - Sqrt(6)"+ ca_sqrt_ui x 2 ctx+ ca_sqrt_ui y 3 ctx+ ca_sqrt_ui z 6 ctx+ ca_mul x x y ctx+ ca_sub x x z ctx+ putStr "\n"; ca_print x ctx; putStr "\n\n"++ putStrLn $ ">>> Exp(1+Sqrt(2)) * Exp(1-Sqrt(2)) / (Exp(1)^2)"+ ca_sqrt_ui x 2 ctx+ ca_add_ui x x 1 ctx+ ca_exp x x ctx+ ca_sqrt_ui y 2 ctx+ ca_ui_sub y 1 y ctx+ ca_exp y y ctx+ ca_one z ctx+ ca_exp z z ctx+ ca_pow_ui z z 2 ctx+ ca_mul x x y ctx+ ca_div x x z ctx+ putStr "\n"; ca_print x ctx; putStr "\n\n"++ putStrLn $ ">>> I^I - Exp(-Pi/2)"+ ca_i x ctx+ ca_pow x x x ctx+ ca_pi y ctx+ ca_div_ui y y 2 ctx+ ca_neg y y ctx+ ca_exp y y ctx+ ca_sub x x y ctx+ putStr "\n"; ca_print x ctx; putStr "\n\n"++ putStrLn $ ">>> Exp(Sqrt(3))^2 - Exp(Sqrt(12))"+ ca_sqrt_ui x 3 ctx+ ca_exp x x ctx+ ca_pow_ui x x 2 ctx+ ca_sqrt_ui y 12 ctx+ ca_exp y y ctx+ ca_sub x x y ctx+ putStr "\n"; ca_print x ctx; putStr "\n\n"++ putStrLn $ ">>> 2*Log(Pi*I) - 4*Log(Sqrt(Pi)) - Pi*I"+ ca_pi_i x ctx+ ca_log x x ctx+ ca_mul_ui x x 2 ctx+ ca_pi y ctx+ ca_sqrt y y ctx+ ca_log y y ctx+ ca_mul_ui y y 4 ctx+ ca_sub x x y ctx+ ca_pi_i y ctx+ ca_sub x x y ctx+ putStr "\n"; ca_print x ctx; putStr "\n\n"++ -- Example 1 in [BBK2014]+ putStrLn $ ">>> -I*Pi/8*Log(2/3-2*I/3)^2 + I*Pi/8*Log(2/3+2*I/3)^2 + \+ \Pi^2/12*Log(-1-I) + Pi^2/12*Log(-1+I) + Pi^2/12*Log(1/3-I/3) + \+ \Pi^2/12*Log(1/3+I/3) - Pi^2/48*Log(18)"+ ca_zero x ctx+ ca_set_qqi_si y 2 3 (-2) 3 ctx+ ca_log y y ctx+ ca_sqr y y ctx+ ca_pi_i z ctx+ ca_mul y y z ctx+ ca_div_si y y (-8) ctx+ ca_add x x y ctx++ ca_set_qqi_si y 2 3 2 3 ctx+ ca_log y y ctx+ ca_sqr y y ctx+ ca_pi_i z ctx+ ca_mul y y z ctx+ ca_div_si y y 8 ctx+ ca_add x x y ctx++ ca_set_qqi_si y (-1) 1 (-1) 1 ctx+ ca_log y y ctx+ ca_pi z ctx+ ca_sqr z z ctx+ ca_mul y y z ctx+ ca_div_si y y 12 ctx+ ca_add x x y ctx++ ca_set_qqi_si y (-1) 1 1 1 ctx+ ca_log y y ctx+ ca_pi z ctx+ ca_sqr z z ctx+ ca_mul y y z ctx+ ca_div_si y y 12 ctx+ ca_add x x y ctx++ ca_set_qqi_si y 1 3 (-1) 3 ctx+ ca_log y y ctx+ ca_pi z ctx+ ca_sqr z z ctx+ ca_mul y y z ctx+ ca_div_si y y 12 ctx+ ca_add x x y ctx++ ca_set_qqi_si y 1 3 1 3 ctx+ ca_log y y ctx+ ca_pi z ctx+ ca_sqr z z ctx+ ca_mul y y z ctx+ ca_div_si y y 12 ctx+ ca_add x x y ctx++ ca_set_ui y 18 ctx+ ca_log y y ctx+ ca_pi z ctx+ ca_sqr z z ctx+ ca_mul y y z ctx+ ca_div_si y y (-48) ctx+ ca_sub x x y ctx+ putStr "\n"; ca_print x ctx; putStr "\n\n"++ putStrLn $ ">>> Sqrt(5 + 2*Sqrt(6)) - Sqrt(2) - Sqrt(3)"+ ca_sqrt_ui x 6 ctx+ ca_mul_ui x x 2 ctx+ ca_add_ui x x 5 ctx+ ca_sqrt x x ctx+ ca_sqrt_ui y 2 ctx+ ca_sub x x y ctx+ ca_sqrt_ui y 3 ctx+ ca_sub x x y ctx+ putStr "\n"; ca_print x ctx; putStr "\n\n"++ putStrLn $ ">>> Sqrt(I) - (1+I)/Sqrt(2)"+ ca_i x ctx+ ca_sqrt x x ctx+ ca_i y ctx+ ca_add_ui y y 1 ctx+ ca_sqrt_ui z 2 ctx+ ca_div y y z ctx+ ca_sub x x y ctx+ putStr "\n"; ca_print x ctx; putStr "\n\n"++ putStrLn $ ">>> Exp(Pi*Sqrt(163)) - (640320^3 + 744)"+ ca_pi x ctx+ ca_sqrt_ui y 163 ctx+ ca_mul x x y ctx+ ca_exp x x ctx+ ca_set_ui y 640320 ctx+ ca_pow_ui y y 3 ctx+ ca_add_ui y y 744 ctx+ ca_sub x x y ctx+ putStr "\n"; ca_print x ctx; putStr "\n\n"++ -- Taken slightly tweaked from:+ -- https://reference.wolfram.com/language/ref/PossibleZeroQ.html+ putStrLn $ ">>> Erf(2*Log(Sqrt(1/2-Sqrt(2)/4))+Log(4)) - Erf(Log(2-Sqrt(2)))"+ ca_sqrt_ui x 2 ctx+ ca_div_ui x x 4 ctx+ ca_one y ctx+ ca_div_ui y y 2 ctx+ ca_sub x y x ctx+ ca_sqrt x x ctx+ ca_log x x ctx+ ca_mul_ui x x 2 ctx+ ca_set_ui y 4 ctx+ ca_log y y ctx+ ca_add x y x ctx+ ca_erf x x ctx+ ca_sqrt_ui y 2 ctx+ ca_ui_sub y 2 y ctx+ ca_log y y ctx+ ca_erf y y ctx+ ca_sub x x y ctx+ putStr "\n"; ca_print x ctx; putStr "\n\n"++ putStr "\n"++ca_set_qqi_si res a b c d ctx = do+ t <- mallocForeignPtr+ withForeignPtr t $ \t -> do+ ca_init t ctx+ ca_set_si t a ctx+ ca_div_si t t b ctx+ ca_i res ctx+ ca_mul_si res res c ctx+ ca_div_si res res d ctx+ ca_add res res t ctx+ ca_clear t ctx++
+ app/expression/Main.hs view
@@ -0,0 +1,179 @@+{-# language FlexibleContexts #-}++import System.IO.Unsafe++import Control.Monad+import Foreign.C.Types+import Foreign.C.String+import Foreign.Marshal.Alloc+import Foreign.Storable++import Data.Map (Map, (!), (!?))++import Data.Number.Flint++main = do+ let n = 12+ x <- toFexpr "x"+ expr <- polyToFexpr (laguerreL 12) x+ heading $ "Hermite polynomial of degree " ++ show n+ print expr+ heading "expanded to normal form"+ print =<< expand expr+ heading "LaTeX representation"+ print =<< latex expr+ heading "another expression"+ print =<< g+ heading "LaTeX representation"+ print =<< latex =<< g++heading x = putStrLn $ "\n\ESC[1;34m" ++ x + ++ "\n" ++ replicate (length x) '‾' ++ "\ESC[0m"+ +--------------------------------------------------------------------------------++g = do+ f <- toFexpr FEXPR_Hypergeometric2F1+ a <- toFexpr (1//2 :: Fmpq)+ b <- toFexpr (1//3 :: Fmpq)+ c <- toFexpr (1//7 :: Fmpq)+ x <- toFexpr "x"+ res <- newFexpr+ withFexpr f $ \f -> do+ withFexpr a $ \a -> do+ withFexpr b $ \b -> do+ withFexpr c $ \c -> do+ withFexpr x $ \x -> do + withFexpr res $ \res -> do+ fexpr_call4 res f a b c x+ return res++horner :: (Foldable t, Num a) => a -> t a -> a+horner x = foldr (\c sum -> c + x * sum) 0++-- polyToFexpr :: (FlintExpression (Item a)) => a -> Fexpr -> IO Fexpr+polyToFexpr poly var = do+ coeffs <- mapM toFexpr $ toList poly+ return $ horner var coeffs++-- polynomials -----------------------------------------------------------++cyclotomicC n = unsafePerformIO $ do+ poly <- newFmpzPoly+ withFmpzPoly poly $ \poly -> fmpz_poly_cyclotomic poly n+ return poly++cosMinPoly n = unsafePerformIO $ do+ poly <- newFmpzPoly+ withFmpzPoly poly $ \poly -> fmpz_poly_cos_minpoly poly n+ return poly++hermiteH n = unsafePerformIO $ do+ poly <- newFmpzPoly+ withFmpzPoly poly $ \poly -> fmpz_poly_hermite_h poly n+ return poly++chebyshevT n = unsafePerformIO $ do+ poly <- newFmpzPoly+ withFmpzPoly poly $ \poly -> fmpz_poly_chebyshev_t poly n+ return poly++chebyshevU n = unsafePerformIO $ do+ poly <- newFmpzPoly+ withFmpzPoly poly $ \poly -> fmpz_poly_chebyshev_u poly n+ return poly++legendreP n = unsafePerformIO $ do+ poly <- newFmpqPoly+ withFmpqPoly poly $ \poly -> fmpq_poly_legendre_p poly n+ return poly++laguerreL n = unsafePerformIO $ do+ poly <- newFmpqPoly+ withFmpqPoly poly $ \poly -> fmpq_poly_laguerre_l poly n+ return poly++monienH n = unsafePerformIO $ do+ poly <- newFmpqPoly+ withFmpqPoly poly $ \poly -> fmpq_poly_monien_h poly n+ return poly++-- manipulating expressions ----------------------------------------------------++expand x = do+ withFexpr x $ \ x -> do+ fexpr_expanded_normal_form x x 0+ return x++latex x = do+ ptr <- malloc :: IO (Ptr CCalciumStream)+ calcium_stream_init_str ptr+ (_, result) <- withFexpr x $ \x -> do+ fexpr_write_latex ptr x 0+ CCalciumStream _ cs len _ <- peek ptr+ s <- peekCString cs+ free cs+ return s+ free ptr+ return result++--------------------------------------------------------------------------------++-- testCalcium = do+-- ctx <- newCaCtx+-- x <- newCa ctx+-- withCa x $ \x -> do+-- withCaCtx ctx $ \ctx -> do+-- ca_one x ctx+-- ca_div_ui x x 2 ctx+-- v <- newCaVec 6 ctx+-- poly <- newCaPoly ctx+-- prod <- newCaPoly ctx+-- ext <- newCaExtFx ca_Cos x ctx+-- withCaCtx ctx $ \ctx -> do+-- withCaExt ext $ \ext -> do+-- ca_ext_print ext ctx; putStr "\n"+-- withCaPoly poly $ \poly -> do+-- withCaPoly prod $ \prod -> do+-- withCa x $ \x -> do+-- ca_poly_one prod ctx+-- ca_poly_x poly ctx+-- forM_ [0..10] $ \j -> do+-- ca_set_si x (-j) ctx+-- ca_poly_x poly ctx+-- ca_poly_set_coeff_ca poly 0 x ctx+-- ca_poly_mul prod prod poly ctx+-- ca_poly_print prod ctx; putStr "\n"+-- ca_pi_i x ctx+-- ca_print x ctx; putStr "\n"+-- ca_euler x ctx+-- ca_print x ctx; putStr "\n"+-- withCaVec v $ \v -> ca_vec_print v ctx+ +-- testCalciumStream = do+-- h <- hermiteH 11+-- w <- g+-- fileName <- newCString "calcium.out"+-- mode <- newCString "w"+-- fp <- fopen fileName mode+-- cs <- newCalciumStreamFile fp+-- withCalciumStream cs $ \cs -> do+-- withFexpr h $ \h -> do+-- fexpr_write cs h+-- fexpr_write_latex cs h 0+-- fexpr_expanded_normal_form h h 0+-- withFexpr w $ \w -> do+-- fexpr_write cs w+-- fexpr_write_latex cs w 0+-- flag <- fclose fp+-- return ()++-- foreign import ccall "stdio.h fopen"+-- fopen :: CString -> CString -> IO (Ptr CFile)++-- foreign import ccall "stdio.h fclose"+-- fclose :: Ptr CFile -> IO CInt++-- foreign import ccall "stdio.h fputs"+-- fputs :: CString -> Ptr CFile -> IO CInt+
+ app/factor_integer/Main.hs view
@@ -0,0 +1,70 @@+import System.TimeIt+import System.IO.Unsafe++import Options.Applicative+import Control.Monad+import Control.Monad.State++import Foreign.C.Types+import Foreign.C.String+import Foreign.Ptr+import Foreign.Marshal.Array+import Foreign.Storable ++import Data.Number.Flint++main = run =<< customExecParser (prefs showHelpOnEmpty) opts where+ desc = "Factor integers."+ opts = info (parameters <**> helper) (+ fullDesc+ <> progDesc desc+ <> header desc)++run params@(Parameters expression num_threads timing) = do+ flint_set_num_threads num_threads+ case parseExpression expression of+ Just n ->+ if timing then do+ timeItNamed "time for factorization" $ print $ factor n+ else do+ print $ factor n+ _ -> putStrLn "Could not parse expression."++parseExpression expression = unsafePerformIO $ do+ mctx <- newFmpzMPolyCtx 0 ord_lex+ f <- newFmpzMPoly mctx+ fac <- newFmpzMPolyFactor mctx+ (_, (_, flag)) <- withFmpzMPolyCtx mctx $ \mctx -> do+ withFmpzMPoly f $ \f -> do+ withCString expression $ \poly -> do + fmpz_mpoly_set_str_pretty f poly nullPtr mctx+ if flag == 0 then do+ n <- newFmpz+ withFmpzMPolyCtx mctx $ \mctx -> do+ withFmpzMPoly f $ \f -> do+ withFmpz n $ \n -> do+ fmpz_mpoly_get_fmpz n f mctx+ return $ Just n+ else do+ return $ Nothing+ +data Parameters = Parameters {+ expression :: String+ , num_threads :: CInt+ , timing :: Bool+} deriving Show++parameters :: Parser Parameters+parameters = Parameters+ <$> argument str (+ help "Integer given as expression (e.g. 2^64+1)"+ <> metavar "INTEGER")+ <*> option auto (+ help "number of threads"+ <> short 't'+ <> long "threads"+ <> value 1+ <> metavar "THREADS")+ <*> switch (+ help "timing"+ <> long "timing")
+ app/fmpq_poly/Main.hs view
@@ -0,0 +1,20 @@+import System.IO.Unsafe+import Data.Char+import Foreign.C.String+import Foreign.C.Types+import Foreign.Marshal.Alloc +import GHC.Read+-- import Text.ParserCombinators.ReadPrec+import Text.ParserCombinators.ReadP+import Control.Monad+import Data.Number.Flint+import Debug.Trace++main = do+ let f, g :: FmpqPoly+ f = read "2 1/2 3/5"+ g = read "4 1/3 2 3/2 -1/2"+ print f+ print g+ print $ f*g+
+ app/fmpz_mod_poly/Main.hs view
@@ -0,0 +1,14 @@+import Data.Number.Flint++main = do++ ctx <- newFmpzModCtx 7++ withNewFmpzModPoly ctx $ \x -> do+ withNewFmpzModPoly ctx $ \y -> do+ withFmpzModCtx ctx $ \ctx -> do+ fmpz_mod_poly_set_coeff_ui x 3 5 ctx+ fmpz_mod_poly_set_coeff_ui x 0 6 ctx+ fmpz_mod_poly_sqr y x ctx+ fmpz_mod_poly_print x ctx; putStr "\n"+ fmpz_mod_poly_print y ctx; putStr "\n"
+ app/fmpz_mpoly_factor/Main.hs view
@@ -0,0 +1,74 @@+import System.TimeIt++import Options.Applicative+import Control.Monad+import Control.Monad.State++import Foreign.C.Types+import Foreign.C.String+import Foreign.Ptr+import Foreign.Marshal.Array+import Foreign.Storable ++import Data.Number.Flint++main = run =<< customExecParser (prefs showHelpOnEmpty) opts where+ desc = "Factor multivariate polynomial in three variables {x, y, z}."+ opts = info (parameters <**> helper) (+ fullDesc+ <> progDesc desc+ <> header desc)++run params@(Parameters poly num_threads timing suppress) = do+ mctx <- newFmpzMPolyCtx 3 ord_lex+ f <- newFmpzMPoly mctx+ fac <- newFmpzMPolyFactor mctx+ vars <- newArray =<< traverse newCString ["x", "y", "z"]+ (_, (_, flag)) <- withFmpzMPolyCtx mctx $ \mctx -> do+ withFmpzMPoly f $ \f -> do+ withCString poly $ \poly -> do + fmpz_mpoly_set_str_pretty f poly vars mctx+ if flag /= 0 then do+ putStrLn "unable to parse polynomial."+ else do+ withFmpzMPolyCtx mctx $ \mctx -> do+ withFmpzMPoly f $ \f -> do+ withFmpzMPolyFactor fac $ \fac -> do+ if timing then do+ timeItNamed "time for factorization" $ fmpz_mpoly_factor fac f mctx+ else do+ fmpz_mpoly_factor fac f mctx+ if not suppress then do+ fmpz_mpoly_print_pretty f vars mctx+ putStr " ="+ fmpz_mpoly_factor_print_pretty fac vars mctx+ putStr "\n"+ else do+ CFmpzMPolyFactor _ _ _ _ num _ <- peek fac+ putStrLn $ "polynomial has " ++ show num ++ " irreducible factors."+ return ()++data Parameters = Parameters {+ poly :: String+ , num_threads :: Int+ , timing :: Bool+ , suppress :: Bool+} deriving Show++parameters :: Parser Parameters+parameters = Parameters+ <$> argument str (+ help "multivariate polynomial in x, y and z."+ <> metavar "POLYNOMIAL")+ <*> option auto (+ help "number of threads"+ <> short 't'+ <> long "threads"+ <> value 1+ <> metavar "THREADS")+ <*> switch (+ help "timing"+ <> long "timing")+ <*> switch (+ help "suppress output"+ <> long "suppress")
+ app/fmpz_poly_factor_zassenhaus/Main.hs view
@@ -0,0 +1,18 @@+import Foreign.C.String++import Data.Number.Flint++main = do+ f <- newFmpzPoly+ facs <- newFmpzPolyFactor+ withFmpzPoly f $ \f -> do+ withCString s $ \s -> do+ fmpz_poly_set_str f s+ withFmpzPolyFactor facs $ \facs -> do+ fmpz_poly_factor_zassenhaus facs f+ fmpz_poly_print f+ putStr "\n"+ fmpz_poly_factor_print facs+ putStr "\n"++s = "63 1 1 1 -4 -7 -2 -6 -3 -7 18 7 25 -11 95 36 21 16 69 56 35 36 32 33 26 -26 -15 -14 -53 -96 67 72 -67 40 -79 -116 -452 -312 -260 -29 -1393 327 69 -28 -241 230 -54 -309 -125 -74 -450 -69 -3 66 -27 73 68 50 -63 -1290 372 31 -16 2"
+ app/fpwrap/Main.hs view
@@ -0,0 +1,40 @@+{-# language FlexibleInstances #-}++import System.IO.Unsafe++import Foreign.C.Types+import Foreign.Marshal.Alloc+import Foreign.Storable++import Text.Printf++import Data.Complex+import Data.Number.Flint hiding (magnitude, phase, realPart, imagPart)++main = do+ let x = 2.0 :: Double+ fx = zeta x+ z = 0.5 :+ 123.0 :: Complex Double+ fz = zeta z+ printf "zeta(%g) = %.16g\n" x fx+ printf "zeta(%g + %gi) = %.16g + %.16gi\n"+ (realPart z ) (imagPart z )+ (realPart fz) (imagPart fz)++-- implement instances for zeta for Double and Complex Double ------------------++instance Special Double where+ zeta x = unsafePerformIO $ do+ alloca $ \f -> do+ flag <- arb_fpwrap_double_zeta f (realToFrac x) 0+ result <- peek f+ return $ realToFrac $ result++instance Special (Complex Double) where+ zeta z = unsafePerformIO $ do+ alloca $ \f -> do+ alloca $ \zp -> do+ poke zp ((realToFrac.realPart) z :+ (realToFrac.imagPart) z)+ flag <- arb_fpwrap_cdouble_zeta f zp 0+ result <- peek f+ return $ ((realToFrac.realPart) result :+ (realToFrac.imagPart) result)
+ app/fq_poly/Main.hs view
@@ -0,0 +1,58 @@+import System.CPUTime+import Text.Printf+import Control.Monad+import Data.Number.Flint++main = do++ putStrLn "Polynomial multiplication over GF(q)\n\+ \------------------------------------\n\+ \\ntimes measured per function application\n"++ state <- newFRandState++ withFRandState state $ \state -> do+ + putStrLn "1) Two length-10,000 polynomials over GF(3^2)"+ + ctx <- newFqCtxConway 3 2 "X"+ fqPolyTest state ctx 10000 1 100 100 ++ putStrLn "2) Two length-500 polynomials over GF(3^263)"+ + ctx <- newFqCtxConway 3 263 "X"+ fqPolyTest state ctx 500 1 1 100++ putStrLn "3) Two length-5 polynomials over GF(109987^4)"+ + ctx <- newFqCtxConway 109987 4 "X"+ fqPolyTest state ctx 5 10000 10000 10000++-- test and timing -------------------------------------------------------------++fqPolyTest state ctx len c r k = do+ putStr "\n"+ withNewFqPoly ctx $ \f -> do+ withNewFqPoly ctx $ \g -> do+ withNewFqPoly ctx $ \h -> do+ withFqCtx ctx $ \ctx -> do++ fq_poly_randtest g state len ctx+ fq_poly_randtest h state len ctx++ timeItNamed "Classical" (fq_poly_mul_classical f g h ctx) c+ timeItNamed "Reorder" (fq_poly_mul_reorder f g h ctx) r+ timeItNamed "KS" (fq_poly_mul_KS f g h ctx) k+ putStr "\n"+ +timeItNamed s f n = do+ t0 <- getCPUTime+ replicateM n f+ t1 <- getCPUTime+ let dt = fromIntegral (t1 - t0) / 10^9 / (fromIntegral n) :: Double+ printf "%s: %.6f ms\n" s dt+++++
+ app/function_benchmark/Main.hs view
@@ -0,0 +1,83 @@+import Text.Printf+import System.CPUTime+import Foreign.C.Types+import Control.Monad++import Data.Number.Flint++main = do+ withNewArb $ \x -> do+ withNewArb $ \y -> do+ withNewArb $ \res -> do+ forM_ funs $ \(fun, desc, limit) -> do+ printf "\n%s\n" desc+ forM_ (takeWhile (<=limit) $ map (10^) [1..]) $ \n -> do+ let prec = round (fromIntegral n / logBase 10 2 + 1)+ printf "%12d" (toInteger n)+ arb_sqrt_ui x 2 prec+ arb_sub_ui x x 1 prec+ + flint_set_num_threads 1+ timeItNamed " 1 thread " $ fun res x n prec++ putStr " "++ flint_set_num_threads 8+ timeItNamed " 8 threads " $ fun res x n prec++ putStr "\n"+ ++funs = + [ ( f0, "const_pi, n digits", 10^8)+ , ( f1, "const_euler, n digits", 10^7)+ , ( f2, "exp(x), n digits", 10^7)+ , ( f3, "log(x), n digits", 10^7)+ , ( f4, "sin(x), n digits", 10^7)+ , ( f5, "atan(x), n digits", 10^7)+ , ( f6, "erf(x), n digits", 10^6)+ , ( f7, "gamma(x), n digits", 10^5)+ , ( f8, "zeta(x), n digits", 10^4)+ , ( f9, "eta(ix), n digits", 10^6)+ , ( f10, "bernoulli(n)", 10^6)+ , ( f11, "partitions(n^2)", 10^7)+ ]++f0 res x n prec = arb_const_pi res prec+f1 res x n prec = arb_const_euler res prec+f2 res x n prec = arb_exp res x prec+f3 res x n prec = arb_log res x prec+f4 res x n prec = arb_sin res x prec+f5 res x n prec = arb_atan res x prec+f6 res x n prec = arb_hypgeom_erf res x prec+f7 res x n prec = arb_gamma res x prec+f8 res x n prec = arb_zeta res x prec++f9 res x n prec = do+ withNewAcb $ \z -> do+ arb_set (acb_imagref z) x+ acb_modular_eta z z prec+ arb_set res (acb_realref z)+ return ()++-- ?? +f10 res x n prec = do+ withNewFmpq $ \t -> do+ bernoulli_fmpq_ui t (fromIntegral n)+ return ()+ +-- ??+f11 res x n prec = do+ withNewFmpz $ \t -> do+ fmpz_set_si t n+ fmpz_mul t t t+ partitions_fmpz_fmpz t t 0+ return ()+ +timeItNamed s f = do+ t0 <- getCPUTime+ result <- f + t1 <- getCPUTime+ let dt = fromIntegral (t1 - t0) / 10^9 :: Double+ printf "%s: %12.4f ms" s dt+
+ app/hilbert_matrix/Main.hs view
@@ -0,0 +1,69 @@+import Options.Applicative+import Control.Monad+import Foreign.Ptr (nullPtr, castPtr)+import Foreign.C.Types+import Foreign.Marshal.Array++import Data.Number.Flint++main = run =<< customExecParser (prefs showHelpOnEmpty) opts where+ opts = info (parameters <**> helper) (+ fullDesc+ <> progDesc "calculates the determinant of the Hilbert matrix"+ <> header "Hilbert matrix determinant")++run p@(Parameters eig n) = do+ print p+ hilbertMatrix eig (fromIntegral n) 20+ +data Parameters = Parameters {+ eig :: Bool+ , n :: Int+} deriving Show++parameters :: Parser Parameters+parameters = Parameters+ <$> switch (+ long "eig"+ <> help "calculating det as a product of eigenvalues.")+ <*> argument auto (+ help "dimension of Hilbert matrix."+ <> metavar "n")+++hilbertMatrix :: Bool -> CLong -> CLong -> IO ()+hilbertMatrix eig n prec = do+ _ <- withNewArb $ \det -> do+ _ <- withNewArbMat n n $ \a -> do+ arb_mat_hilbert a prec+ if not eig then+ arb_mat_det det a prec+ else do+ _ <- withNewAcbMat n n $ \r -> do+ _ <- withNewAcbMat n n $ \c -> do+ acb_mat_set_arb_mat c a+ e <- _acb_vec_init n+ _ <- acb_mat_approx_eig_qr e nullPtr r c nullPtr 0 prec+ simple <- acb_mat_eig_simple e nullPtr nullPtr c e r prec+ if simple == 1 then do+ arb_one det+ forM_ [0 .. fromIntegral n-1] $ \j -> do+ arb_zero (castPtr (e `advancePtr` j) `advancePtr` 1)+ _acb_vec_sort_pretty e n+ acb_get_real det e+ else+ arb_indeterminate det+ _acb_vec_clear e n+ return ()+ return ()+ return ()+ zero <- arb_contains_zero det+ if zero == 1 then do+ hilbertMatrix eig n (prec + 20)+ else do+ putStrLn $ "success with prec " ++ show prec ++ " bits."+ putStr "value: "+ arb_printd det 16+ putStr "\n"+ return ()+
+ app/hilbert_matrix_ca/Main.hs view
@@ -0,0 +1,113 @@+import Options.Applicative+import Control.Monad++import System.TimeIt++import Foreign.C.Types+import Foreign.Marshal.Array++import Text.Printf++import Data.Number.Flint++main = timeIt $ run =<< customExecParser (prefs showHelpOnEmpty) opts where+ desc = "This program constructs the Hilbert matrixq as exact \+ \algebraic numbers, and verifies the exact trace and \+ \determinant formulas."++ opts = info (options <**> helper) (+ fullDesc+ <> progDesc desc+ <> header desc)++run params@(Options n useQQbar useVieta) = do+ print params+ if useQQbar then do+ withNewFmpqMat n n $ \mat -> do+ withNewQQbar $ \trace -> do+ withNewQQbar $ \det -> do+ + eig <- _qqbar_vec_init n+ + fmpq_mat_hilbert_matrix mat+ qqbar_eigenvalues_fmpq_mat eig mat 0+ + putStrLn "Trace:"+ forM_ [0 .. fromIntegral n - 1] $ \i -> do+ qqbar_add trace trace (eig `advancePtr` i)+ deg <- qqbar_degree trace+ putStr $ show i ++ " " ++ show n ++ ": "; print deg+ qqbar_print trace; putStr "\n"+ + putStrLn "Determinant:"+ qqbar_one det+ forM_ [0 .. fromIntegral n - 1] $ \i -> do+ qqbar_mul det det (eig `advancePtr` i)+ deg <- fromIntegral <$> qqbar_degree det + putStr $ show i ++ " " ++ show n ++ ": "; print deg+ qqbar_print det; putStr "\n"+ + _qqbar_vec_clear eig n+ return ()+ else do+ ctx <- newCaCtx+ withCaCtx ctx $ \ctx -> do+ -- Verification requires high-degree algebraics.+ ca_ctx_set_option ctx ca_opt_qqbar_deg_limit 10000+ if useVieta then do + ca_ctx_set_option ctx ca_opt_vieta_limit n+ else do+ ca_ctx_set_option ctx ca_opt_vieta_limit 0+ withNewCaMat n n ctx $ \mat -> do+ withNewCa ctx $ \trace -> do+ withNewCa ctx $ \det -> do+ withNewCa ctx $ \t -> do+ withNewCaVec 0 ctx $ \eig -> do + withCaCtx ctx $ \ctx -> do + + ca_mat_hilbert mat ctx++ allocaArray (fromIntegral n) $ \mul -> do+ ca_mat_eigenvalues eig mul mat ctx++ -- note: in general, we should use the multiplicities, but+ -- we happen to know that the eigenvalues are simple here++ ca_mat_trace trace mat ctx+ ca_mat_det det mat ctx++ putStrLn "Trace:"+ forM_ [0 .. fromIntegral n - 1] $ \i -> do+ ca_add t t (ca_vec_entry_ptr eig i) ctx+ ca_print trace ctx; putStr "\n"+ ca_print t ctx; putStr "\n"+ res <- ca_check_equal trace t ctx+ putStr $ "Equal: " ++ show res ++ "\n\n"++ putStrLn "Det:"+ ca_one t ctx+ forM_ [0 .. fromIntegral n - 1] $ \i -> do+ ca_mul t t (ca_vec_entry_ptr eig i) ctx+ ca_print det ctx; putStr "\n"+ ca_print t ctx; putStr "\n"+ res <- ca_check_equal det t ctx+ putStr $ "Equal: " ++ show res ++ "\n\n"+ putStr "\n."+ +data Options = Options {+ n :: CLong+ , useQQbar :: Bool+ , useVieta :: Bool+} deriving Show++options :: Parser Options+options = Options+ <$> argument auto (+ help "n, the dimension of the Hilbert matrix"+ <> metavar "n")+ <*> switch (+ help "use QQbar arithmetic."+ <> long "qqbar")+ <*> switch (+ help "use Vieta formula."+ <> long "vieta")
+ app/integrals/Integrands.hs view
@@ -0,0 +1,532 @@+module Integrands (+ description+, makeFunPtr+-- * Functions+, f_airy_ai+, f_atanderiv+, f_circle+, f_elliptic_p_laurent_n+, f_erf_bent+, f_essing+, f_essing2+, f_exp+, f_exp_airy+, f_factorial1000+, f_floor+, f_gamma+, f_gaussian+, f_gaussian_twist+, f_helfgott+, f_horror+, f_lambertw+, f_log_div1p+, f_log_div1p_transformed+, f_max_sin_cos+, f_monster+, f_rgamma+, f_rsqrt+, f_rump+, f_scaled_bessel+, f_sech+, f_sech3+, f_sin+, f_sin_cos_frac+, f_sin_near_essing+, f_sin_plus_small+, f_spike+, f_sqrt+, f_zeta+, f_zeta_frac+) where++import Foreign.C.Types+import Foreign.Ptr+import Foreign.Storable+import Foreign.Marshal.Array++import Control.Monad++import Data.Number.Flint++foreign import ccall safe "wrapper"+ makeFunPtr :: CAcbCalcFunc -> IO (FunPtr CAcbCalcFunc)++--------------------------------------------------------------------------------++-- f(z) = Ai(z)+f_airy_ai res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ acb_hypgeom_airy res nullPtr nullPtr nullPtr z prec+ return 0++-- f(z) = 1/(1+z^2)+f_atanderiv res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ acb_mul res z z prec+ acb_add_ui res res 1 prec+ acb_inv res res prec+ return 0++-- f(z) = sqrt(1-z^2)+f_circle res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ acb_one res + acb_submul res z z prec + acb_real_sqrtpos res res (if order /= 0 then 1 else 0) prec+ return 0++f_elliptic_p_laurent_n res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ n <- peek (castPtr param) :: IO CLong+ withNewAcb $ \tau -> do+ acb_onei tau+ acb_modular_elliptic_p res z tau prec+ acb_pow_si tau z (-n-1) prec+ acb_mul res res tau prec+ return 0++-- f(z) = erf(z/sqrt(0.0002)*0.5 +1.5)*exp(-z)+-- example provided by Silviu-Ioan Filip+f_erf_bent res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ withNewAcb $ \t -> do+ acb_set_ui t 1250+ acb_sqrt t t prec+ acb_mul t t z prec+ acb_set_d res 1.5+ acb_add res res t prec+ acb_hypgeom_erf res res prec++ acb_neg t z+ acb_exp t t prec+ acb_mul res res t prec+ return 0++-- f(z) = sin(1/z)+-- Assume z on real interval+f_essing res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ isReal <- (==1) <$> acb_is_real z+ containsZero <- (==1) <$> arb_contains_zero (acb_realref z)+ if order == 0 && isReal && containsZero then do+ acb_zero res+ mag_one (arb_radref (acb_realref res))+ else do+ acb_inv res z prec+ acb_sin res res prec+ return 0++-- f(z) = z*sin(1/z)+-- Assume z on real interval+f_essing2 res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ isReal <- (==1) <$> acb_is_real z+ containsZero <- (==1) <$> arb_contains_zero (acb_realref z)+ if order == 0 && isReal && containsZero then do+ acb_zero res+ mag_one (arb_radref (acb_realref res))+ else do+ acb_inv res z prec+ acb_sin res res prec+ acb_mul res res z prec+ return 0++-- f(z) = exp(z)+f_exp res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ acb_exp res z prec+ return 0++-- f(z) = exp(-z) Ai(-z)+f_exp_airy res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ withNewAcb $ \t -> do+ acb_neg t z+ acb_hypgeom_airy res nullPtr nullPtr nullPtr t prec+ acb_exp t t prec+ acb_mul res res t prec+ return 0++-- f(z) = exp(-z)*z^1000+f_factorial1000 res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ withNewAcb $ \t -> do+ acb_pow_ui t z 1000 prec+ acb_neg res z+ acb_exp res res prec+ acb_mul res res t prec+ return 0++-- f(z) = floor(z)+f_floor res z param order prec = do+ when (order > 1) $ error "f_floor: Would be needed for Taylor method."+ acb_real_floor res z (if order /= 0 then 1 else 0) prec+ return 0++-- f(z) = gamma(z)+f_gamma res z param order prec = do+ when (order > 1) $ error "f_floor: Would be needed for Taylor method."+ acb_gamma res z prec+ return 0++-- f(z) = exp(-z^2)+f_gaussian res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ acb_mul z z z prec+ acb_neg z z+ acb_exp res z prec+ return 0++-- f(z) = exp(-z^2+iz)+f_gaussian_twist res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ acb_mul_onei res z+ acb_submul res z z prec+ acb_exp res res prec+ return 0++-- | /f_helfgott/ /res/ /z/ /param/ /order/ /prec/+--+-- f(z) = |z^4 + 10z^3 + 19z^2 - 6z - 6| exp(z)+-- (for real z)+-- Helfgott's integral on MathOverflow+f_helfgott res z param order prec = do+ when (order > 1) $ error "f_helfgott: Would be needed for Taylor method."+ acb_add_si res z 10 prec+ acb_mul res res z prec+ acb_add_si res res 19 prec+ acb_mul res res z prec+ acb_add_si res res (-6) prec+ acb_mul res res z prec+ acb_add_si res res (-6) prec+ acb_real_abs res res (if order /= 0 then 1 else 0) prec+ isFinite <- (==1) <$> acb_is_finite res+ when isFinite $ do+ withNewAcb $ \t -> do+ acb_exp t z prec+ acb_mul res res t prec+ return ()+ return 0++f_horror res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ withNewAcb $ \s -> do+ withNewAcb $ \t -> do+ acb_real_floor res z (if order /= 0 then 1 else 0) prec+ isFinite <- (==1) <$> acb_is_finite res+ when isFinite $ do+ acb_sub res z res prec+ acb_set_d t 0.5+ acb_sub res res t prec+ acb_sin_cos s t z prec+ acb_real_max s s t (if order /= 0 then 1 else 0) prec+ acb_mul res res s prec+ return 0++f_lambertw res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ bits <- acb_rel_accuracy_bits z+ let prec' = min prec (bits + 10)+ withNewAcb $ \t -> do+ when (order /= 0 ) $ do+ arb_const_e (acb_realref t) prec'+ acb_inv t t prec'+ acb_add t t z prec'+ containsZero <- (==1) <$> arb_contains_zero (acb_imagref t)+ nonPositive <- (==1) <$> arb_contains_nonpositive (acb_realref t)+ when (containsZero && nonPositive) $ acb_indeterminate t+ return ()+ isFinite <- (==1) <$> acb_is_finite t+ if isFinite then do+ withNewFmpz $ \k -> do+ acb_lambertw res z k 0 prec'+ return ()+ else do+ acb_indeterminate res+ return ()+ return 0++-- f(z) = -log(z) / (1 + z)+f_log_div1p res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ withNewAcb $ \t -> do+ acb_add_ui t z 1 prec+ acb_log res z prec+ acb_div res res t prec+ acb_neg res res+ return 0++-- f(z) = z exp(-z) / (1 + exp(-z))+f_log_div1p_transformed res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ withNewAcb $ \t -> do+ acb_neg t z+ acb_exp t t prec+ acb_add_ui res t 1 prec+ acb_div res t res prec+ acb_mul res res z prec+ return 0++-- f(z) = max(sin(z), cos(z))+f_max_sin_cos res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ withNewAcb $ \s -> do+ withNewAcb $ \c -> do+ acb_sin_cos s c z prec+ acb_real_max res s c (if order /= 0 then 1 else 0) prec+ return 0++f_monster res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ withNewAcb $ \t -> do+ acb_exp t z prec+ acb_real_floor res t (if order /= 0 then 1 else 0) prec+ isFinite <- (==1) <$> acb_is_finite res+ when isFinite $ do+ acb_sub res t res prec+ acb_add t t z prec+ acb_sin t t prec+ acb_mul res res t prec+ return 0++-- f(z) = rgamma(z)+f_rgamma res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ acb_rgamma res z prec+ return 0++-- f(z) = rsqrt(z)+f_rsqrt res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ acb_rsqrt_analytic res z (if order /= 0 then 1 else 0) prec+ return 0++-- f(z) = sin(z + exp(z)) -- Rump's oscillatory example+f_rump res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ acb_exp res z prec+ acb_add res res z prec+ acb_sin res res prec+ return 0++-- f(z) = exp(-z) (I_0(z/k))^k, from Bruno Salvy+f_scaled_bessel res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ k <- peek (castPtr param)+ withNewAcb $ \nu -> do+ acb_init nu+ acb_div_ui res z k prec+ acb_hypgeom_bessel_i_scaled res nu res prec+ acb_pow_ui res res k prec+ acb_clear nu+ return 0++-- f(z) = sech(z)+f_sech res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ acb_sech res z prec+ return 0++-- f(z) = sech^3(z)+f_sech3 res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ acb_sech res z prec+ acb_cube res res prec+ return 0++-- f(z) = sin(z)+f_sin res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ acb_sin res z prec+ return 0++-- f(z) = z sin(z) / (1 + cos(z)^2)+f_sin_cos_frac res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ withNewAcb $ \s -> do+ withNewAcb $ \c -> do+ acb_sin_cos s c z prec+ acb_mul c c c prec+ acb_add_ui c c 1 prec+ acb_mul s s z prec+ acb_div res s c prec+ return 0++-- f(z) = sin((1/1000 + (1-z)^2)^(-3/2)), example from+-- Mioara Jolde's thesis (suggested by Nicolas Brisebarre)+f_sin_near_essing res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ withNewAcb $ \t -> do+ withNewAcb $ \u -> do+ acb_sub_ui t z 1 prec+ acb_neg t t+ acb_mul t t t prec+ acb_one u+ acb_div_ui u u 1000 prec+ acb_add t t u prec+ acb_set_d u (-1.5)+ acb_pow_analytic t t u (if order /= 0 then 1 else 0) prec+ acb_sin res t prec+ return 0++-- f(z) = sin(z) + exp(-200-z^2)+f_sin_plus_small res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ withNewAcb $ \t -> do+ acb_mul t z z prec+ acb_add_ui t t 200 prec+ acb_neg t t+ acb_exp t t prec+ acb_sin res z prec+ acb_add res res t prec+ return 0++-- f(z) = sech(10(x-0.2))^2 + sech(100(x-0.4))^4 + sech(1000(x-0.6))^6+f_spike res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ withNewAcb $ \a -> do+ withNewAcb $ \b -> do+ withNewAcb $ \c -> do+ acb_mul_ui a z 10 prec+ acb_sub_ui a a 2 prec+ acb_sech a a prec+ acb_pow_ui a a 2 prec++ acb_mul_ui b z 100 prec+ acb_sub_ui b b 40 prec+ acb_sech b b prec+ acb_pow_ui b b 4 prec++ acb_mul_ui c z 1000 prec+ acb_sub_ui c c 600 prec+ acb_sech c c prec+ acb_pow_ui c c 6 prec++ acb_add res a b prec+ acb_add res res c prec+ return 0++-- f(z) = sqrt(z)+f_sqrt res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ acb_sqrt_analytic res z (if order /= 0 then 1 else 0) prec+ return 0++-- f(z) = zeta(z)+f_zeta res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ acb_zeta res z prec+ return 0++-- f(z) = zeta'(z) / zeta(z)+f_zeta_frac res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ t <- _acb_vec_init 2+ acb_dirichlet_zeta_jet t z 0 2 prec+ acb_div res (t `advancePtr` 1) t prec+ _acb_vec_clear t 2+ return 0++-- examples --------------------------------------------------------------------++description = + [ "int_0^100 sin(x) dx"+ , "4 int_0^1 1/(1+x^2) dx"+ , "2 int_0^{inf} 1/(1+x^2) dx (using domain truncation)"+ , "4 int_0^1 sqrt(1-x^2) dx"+ , "int_0^8 sin(x+exp(x)) dx"+ , "int_1^101 floor(x) dx"+ , "int_0^1 |x^4+10x^3+19x^2-6x-6| exp(x) dx"+ , "1/(2 pi i) int zeta(s) ds (closed path around s = 1)"+ , "int_0^1 sin(1/x) dx (slow convergence, use -heap and/or -tol)"+ , "int_0^1 x sin(1/x) dx (slow convergence, use -heap and/or -tol)"+ , "int_0^10000 x^1000 exp(-x) dx"+ , "int_1^{1+1000i} gamma(x) dx"+ , "int_{-10}^{10} sin(x) + exp(-200-x^2) dx"+ , "int_{-1020}^{-1010} exp(x) dx (use -tol 0 for relative error)"+ , "int_0^{inf} exp(-x^2) dx (using domain truncation)"+ , "int_0^1 sech(10(x-0.2))^2 + sech(100(x-0.4))^4 + sech(1000(x-0.6))^6 dx"+ , "int_0^8 (exp(x)-floor(exp(x))) sin(x+exp(x)) dx (use higher -eval)"+ , "int_0^{inf} sech(x) dx (using domain truncation)"+ , "int_0^{inf} sech^3(x) dx (using domain truncation)"+ , "int_0^1 -log(x)/(1+x) dx (using domain truncation)"+ , "int_0^{inf} x exp(-x)/(1+exp(-x)) dx (using domain truncation)"+ , "int_C wp(x)/x^(11) dx (contour for 10th Laurent coefficient of Weierstrass p-function)"+ , "N(1000) = count zeros with 0 < t <= 1000 of zeta(s) using argument principle"+ , "int_0^{1000} W_0(x) dx"+ , "int_0^pi max(sin(x), cos(x)) dx"+ , "int_{-1}^1 erf(x/sqrt(0.0002)*0.5+1.5)*exp(-x) dx"+ , "int_{-10}^10 Ai(x) dx"+ , "int_0^10 (x-floor(x)-1/2) max(sin(x),cos(x)) dx"+ , "int_{-1-i}^{-1+i} sqrt(x) dx"+ , "int_0^{inf} exp(-x^2+ix) dx (using domain truncation)"+ , "int_0^{inf} exp(-x) Ai(-x) dx (using domain truncation)"+ , "int_0^pi x sin(x) / (1 + cos(x)^2) dx"+ , "int_0^3 sin(0.001 + (1-x)^2)^(-3/2)) dx (slow convergence, use higher -eval)"+ , "int_0^{inf} exp(-x) I_0(x/3)^3 dx (using domain truncation)"+ , "int_0^{inf} exp(-x) I_0(x/15)^{15} dx (using domain truncation)"+ , "int_{-1-i}^{-1+i} 1/sqrt(x) dx"+ , "int_0^{inf} 1/gamma(x) dx (using domain truncation)"+ ]++--------------------------------------------------------------------------------++functions :: [CAcbCalcFunc]+functions =+ [+ f_airy_ai+ , f_atanderiv+ , f_circle+ , f_elliptic_p_laurent_n+ , f_erf_bent+ , f_essing+ , f_essing2+ , f_exp+ , f_exp_airy+ , f_factorial1000+ , f_floor+ , f_gamma+ , f_gaussian+ , f_gaussian_twist+ , f_helfgott+ , f_horror+ , f_lambertw+ , f_log_div1p+ , f_log_div1p_transformed+ , f_max_sin_cos+ , f_monster+ , f_rgamma+ , f_rsqrt+ , f_rump+ , f_scaled_bessel+ , f_sech+ , f_sech3+ , f_sin+ , f_sin_cos_frac+ , f_sin_near_essing+ , f_sin_plus_small+ , f_spike+ , f_sqrt+ , f_zeta+ , f_zeta_frac+ ]++lift :: (Ptr CAcb -> Ptr CAcb -> Ptr () -> CLong -> CLong -> IO ())+ -> Ptr CAcb -> Ptr CAcb -> Ptr () -> CLong -> CLong -> IO CInt+lift f res z param order prec = do+ when (order > 1) $ error "order > 1 would be needed for Taylor method."+ f res z param order prec+ return 0++testFunction f x = do+ withNewAcb $ \res -> do+ withNewAcb $ \t -> do+ acb_set_d t x+ putStr"testFunction: arg = "+ acb_printn t 16 arb_str_no_radius+ putStr "\n"+ flag <- f res t nullPtr 0 1024+ putStr "testFunction: res = "+ acb_printn res 16 arb_str_no_radius+ putStr "\n"+ return ()
+ app/integrals/Main.hs view
@@ -0,0 +1,619 @@+import GHC.Conc++import System.IO.Unsafe+import System.CPUTime++import Options.Applicative++import Control.Monad+import Control.Applicative ((<|>))++import Foreign.C.Types+import Foreign.C.String+import Foreign.Ptr+import Foreign.Marshal.Alloc (free)+import Foreign.Marshal.Utils (with)+import Foreign.Storable++import Text.ParserCombinators.ReadP hiding (option)+import Text.Read (readMaybe)+import Text.Printf++import Data.Bits+import Data.Char+import Data.List (intercalate)+import Data.Number.Flint++import Integrands++main = run =<< execParser opts where+ hDesc = "Calculate integrals using acb_calculate."+ desc = "Calculate integrals from list in range 0:"+ ++ show (length description) ++ ". "+ ++ "For a list of implemented integrals use --list."+ opts = info (parameters <**> helper) (+ fullDesc+ <> progDesc desc+ <> header hDesc)++run params = do+ if list params then do+ putStrLn "List of implemented integrals:\n"+ mapM_ (uncurry (printf " %2d %s\n")) $ zip [0 :: Int ..] description+ else do+ calc params+ +calc params@(Parameters list (Range (start, end)) prec opt_goal tol twice+ heap verbose deg eval depth num_threads) = do+ numProcessors <- getNumProcessors+ if num_threads <= fromIntegral numProcessors then do + flint_set_num_threads num_threads+ else do+ let msg = "number of threads (=" ++ show num_threads ++ ") greater than "+ ++ "number of processors (=" ++ show numProcessors ++ ").\n\n"+ ++ "choose [1.." ++ show numProcessors ++ "] threads."+ error msg+ let use_heap = if heap then 1 else 0+ goal = if opt_goal == 0 then prec else opt_goal+ opts <- newAcbCalcIntegrateOpt_ deg eval depth use_heap verbose+ print opts+ withAcbCalcIntegrateOpt opts $ \opts -> do+ withMag tol $ \tol -> do+ flag <- mag_is_zero tol+ when (flag == 1) $ do mag_set_ui_2exp_si tol 1 (-prec)+ withNewAcb $ \a -> do+ withNewAcb $ \b -> do+ withNewAcb $ \s -> do+ withNewAcb $ \t -> do+ forM_ [start .. end] $ \j -> do+ let desc = description !! j+ startTime <- getCPUTime+ replicateM_ (if twice then 2 else 1) $ do+ case j of+ 0 -> do+ f <- makeFunPtr f_sin+ acb_set_si a 0+ acb_set_si b 100+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ 1 -> do+ f <- makeFunPtr f_atanderiv+ acb_set_si a 0+ acb_set_si b 1+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ acb_mul_2exp_si s s 2+ return ()+ 2 -> do+ f <- makeFunPtr f_atanderiv+ acb_set_si a 0+ acb_set_si b 1+ acb_mul_2exp_si b b goal+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ arb_add_error_2exp_si (acb_realref s) (-goal)+ acb_mul_2exp_si s s 1+ return ()+ 3 -> do+ f <- makeFunPtr f_circle+ acb_set_si a 0+ acb_set_si b 1+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ acb_mul_2exp_si s s 2+ return ()+ 4 -> do+ f <- makeFunPtr f_rump+ acb_set_si a 0+ acb_set_si b 8+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ 5 -> do+ f <- makeFunPtr f_floor+ acb_set_si a 1+ acb_set_si b 101+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ 6 -> do+ f <- makeFunPtr f_helfgott+ acb_set_si a 0+ acb_set_si b 1+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ 7 -> do+ f <- makeFunPtr f_zeta+ let points = [(-1, -1), (2, -1), (2, 1), (-1, 1)]+ path = zip points $ tail $ cycle points+ acb_zero s+ forM_ path $ \((x1, y1), (x2, y2)) -> do + acb_set_si_si a x1 y1+ acb_set_si_si b x2 y2+ acb_calc_integrate t f nullPtr a b goal tol opts prec+ acb_add s s t prec+ acb_const_pi t prec+ acb_div s s t prec+ acb_mul_2exp_si s s (-1)+ acb_div_onei s s+ return()+ 8 -> do+ f <- makeFunPtr f_essing+ acb_set_si a 0+ acb_set_si b 1+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ 9 -> do+ f <- makeFunPtr f_essing2+ acb_set_si a 0+ acb_set_si b 1+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ 10 -> do+ f <- makeFunPtr f_factorial1000+ acb_set_si a 0+ acb_set_si b 10000+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ 11 -> do+ f <- makeFunPtr f_gamma+ acb_set_si_si a 1 0+ acb_set_si_si b 1 1000+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ 12 -> do+ f <- makeFunPtr f_sin_plus_small+ acb_set_si a (-10)+ acb_set_si b 10+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ 13 -> do+ f <- makeFunPtr f_exp+ acb_set_si a (-1020)+ acb_set_si b (-1010)+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ 14 -> do+ f <- makeFunPtr f_gaussian+ acb_zero a+ acb_set_si b $ fromIntegral+ $ ceiling + $ sqrt (fromIntegral goal * log 2) + 1+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ acb_mul b b b prec+ acb_neg b b+ acb_exp b b prec+ arb_add_error (acb_realref s) (acb_realref b)+ return ()+ 15 -> do+ f <- makeFunPtr f_spike+ acb_zero a+ acb_one b + acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ 16 -> do+ f <- makeFunPtr f_monster+ acb_zero a+ acb_set_si b 8+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ 17 -> do+ f <- makeFunPtr f_sech+ acb_set_d a 0+ acb_set_d b $ fromIntegral+ $ ceiling + $ fromIntegral goal * log 2 + 1+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ acb_neg b b+ acb_exp b b prec+ acb_mul_2exp_si b b 1+ arb_add_error (acb_realref s) (acb_realref b)+ return ()+ 18 -> do+ f <- makeFunPtr f_sech3+ acb_zero a+ acb_set_d b $ fromIntegral+ $ ceiling + $ fromIntegral goal * log 2 / 3 + 2+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ acb_neg b b + acb_mul_ui b b 3 prec+ acb_exp b b prec+ acb_mul_2exp_si b b 3+ acb_div_ui b b 3 prec+ arb_add_error (acb_realref s) (acb_realref b)+ return ()+ 19 -> do+ when (goal < 0) $ do error "goal < 0"+ -- error bound 2^-N (1+N) when truncated at 2^-N+ let bitCount x = finiteBitSize x - countLeadingZeros x+ n = goal + fromIntegral (bitCount goal)+ acb_one a+ acb_mul_2exp_si a a (-n)+ acb_one b+ f <- makeFunPtr f_log_div1p+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ acb_set_si b (n+1)+ acb_mul_2exp_si b b (-n)+ arb_add_error (acb_realref s) (acb_realref b)+ return ()+ 20 -> do + f <- makeFunPtr f_log_div1p_transformed+ when (goal < 0) $ do error "goal < 0"+ -- error bound (N+1) exp(-N) when truncated at N+ let bitCount x = finiteBitSize x - countLeadingZeros x+ n = goal + fromIntegral (bitCount goal)+ acb_zero a+ acb_set_si b n+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ acb_neg b b+ acb_exp b b prec+ acb_mul_si b b (n+1) prec+ arb_add_error (acb_realref s) (acb_realref b)+ return ()+ 21 -> do+ f <- makeFunPtr f_elliptic_p_laurent_n+ let n = 10 ::CLong+ points = [(1, -1), (1, 1), (-1, 1), (-1, -1)]+ path = zip points $ tail $ cycle points+ acb_zero s+ with n $ \np -> do+ let p = castPtr np :: Ptr ()+ forM_ path $ \((x1, y1), (x2, y2)) -> do+ acb_set_si_si a x1 y1+ acb_set_si_si b x2 y2+ acb_div_ui a a 2 prec+ acb_div_ui b b 2 prec+ acb_calc_integrate t f p a b goal tol opts prec+ acb_add s s t prec+ acb_const_pi t prec+ acb_div s s t prec+ acb_mul_2exp_si s s (-1)+ acb_div_onei s s+ return ()+ 22 -> do+ f <- makeFunPtr f_zeta_frac+ let n = 1000 :: CLong+ acb_zero s++ acb_set_si_si a 100 0+ acb_set_si_si b 100 n+ acb_calc_integrate t f nullPtr a b goal tol opts prec+ acb_add s s t prec++ acb_set_si_si a 100 n+ acb_set_si_si b 1 (2*n)+ acb_div_ui b b 2 prec+ acb_calc_integrate t f nullPtr a b goal tol opts prec+ acb_add s s t prec++ acb_div_onei s s+ arb_zero (acb_imagref s)++ acb_set_si t n+ acb_dirichlet_hardy_theta t t nullPtr nullPtr 1 prec+ acb_add s s t prec++ acb_const_pi t prec+ acb_div s s t prec+ acb_add_ui s s 1 prec+ return ()+ 23 -> do+ acb_set_si a 0+ acb_set_si b 1000+ f <- makeFunPtr f_lambertw + acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ 24 -> do+ acb_zero a+ acb_const_pi b prec+ f <- makeFunPtr f_max_sin_cos+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ 25 -> do+ acb_set_si a (-1)+ acb_set_si b 1+ f <- makeFunPtr f_erf_bent+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ 26 -> do+ acb_set_si a (-10)+ acb_set_si b 10+ f <- makeFunPtr f_airy_ai+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ 27 -> do+ acb_set_si a 0+ acb_set_si b 10+ f <- makeFunPtr f_horror+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ 28 -> do+ acb_set_si_si a (-1) (-1)+ acb_set_si_si b (-1) 1+ f <- makeFunPtr f_sqrt+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ 29 -> do+ f <- makeFunPtr f_gaussian_twist+ acb_zero a+ acb_set_si b $ ceiling+ $ sqrt (fromIntegral goal * log 2) + 1+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ acb_mul b b b prec+ acb_neg b b+ acb_exp b b prec+ arb_add_error (acb_realref s) (acb_realref b)+ arb_add_error (acb_imagref s) (acb_realref b)+ return ()+ 30 -> do+ f <- makeFunPtr f_exp_airy+ acb_zero a+ acb_set_si b $ ceiling+ $ fromIntegral goal * log 2 + 1+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ acb_mul b b b prec+ acb_neg b b+ acb_exp b b prec+ acb_mul_2exp_si b b 1+ arb_add_error (acb_realref s) (acb_realref b)+ return ()+ 31 -> do+ f <- makeFunPtr f_sin_cos_frac+ acb_zero a+ acb_const_pi b prec+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ 32 -> do+ f <- makeFunPtr f_sin_near_essing+ acb_zero a+ acb_set_ui b 3+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ 33 -> do+ f <- makeFunPtr f_scaled_bessel+ acb_zero a+ acb_zero b+ let k = 3+ scaled_bessel_select_N (acb_realref b) k prec+ with k $ \kp -> do+ let p = castPtr kp :: Ptr ()+ acb_calc_integrate s f p a b goal tol opts prec+ scaled_bessel_tail_bound (acb_realref a) k+ (acb_realref b) prec+ arb_add_error (acb_realref s) (acb_realref a)+ return ()+ 34 -> do+ f <- makeFunPtr f_scaled_bessel+ acb_zero a+ acb_zero b+ let k = 15+ scaled_bessel_select_N (acb_realref b) k prec+ with k $ \kp -> do+ let p = castPtr kp :: Ptr ()+ acb_calc_integrate s f p a b goal tol opts prec+ scaled_bessel_tail_bound (acb_realref a) k+ (acb_realref b) prec+ arb_add_error (acb_realref s) (acb_realref a)+ return ()+ 35 -> do+ f <- makeFunPtr f_rsqrt+ acb_set_si_si a (-1) (-1)+ acb_set_si_si b (-1) 1+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ 36 -> do+ f <- makeFunPtr f_rgamma+ when (goal < 0) $ do error "goal < 0."+ acb_zero a+ acb_set_si b (4 + (goal + 1) `div` 2)+ acb_calc_integrate s f nullPtr a b goal tol opts prec+ return ()+ _ -> do+ error "case does not exists."+ return ()+ endTime <- getCPUTime+ putStrLn $ "I" ++ show j ++ " = " ++ desc+ let diff = fromIntegral (endTime - startTime) / 10^12 :: Double+ printf "computation time: %0.3f sec\n" diff+ let digits = round (0.333 * fromIntegral prec) :: CLong + acb_printn s digits arb_str_none+ putStr "\n\n"+ return ()++-- auxiliary routines ----------------------------------------------------------++scaled_bessel_tail_bound b k n prec = do+ arb_const_pi b prec+ arb_mul b b n prec+ arb_ui_div b k b prec+ arb_sqrt b b prec+ arb_pow_ui b b k prec+ arb_mul b b n prec+ arb_mul_ui b b 2 prec+ arb_div_ui b b (k-2) prec++scaled_bessel_select_N n k prec = do+ let f = log (fromIntegral k / pi) / log 2 :: CDouble+ c e = fromIntegral e * fromIntegral k / 2 - f - fromIntegral e+ < fromIntegral prec + 5+ e = last $ takeWhile c [1, 2 ..]+ arb_one n+ arb_mul_2exp_si n n e++-- Parser ----------------------------------------------------------------------++data Parameters = Parameters {+ list :: Bool+ , range :: Range+ , prec :: CLong+ , goal :: CLong+ , tol :: Mag+ , twice :: Bool+ , heap :: Bool+ , verbosity :: CInt+ , deg :: CLong+ , eval :: CLong+ , depth :: CLong+ , num_threads :: CInt+} deriving Show+ +parameters :: Parser Parameters+parameters = Parameters+ <$> switch (+ help "show list of implemented integrals."+ <> long "list"+ <> short 'l')+ <*> option rng (+ help "range (for example --range 2:6 or --range 7)" + <> long "range"+ <> short 'r'+ <> value (Range (0, length description - 1))+ <> metavar "range")+ <*> option pos (+ help "precision in bits (default p = 64)"+ <> long "prec"+ <> short 'p'+ <> value 64+ <> metavar "p")+ <*> option pos (+ help "approximate relative accuracy goal (p)"+ <> long "goal"+ <> value 0+ <> metavar "goal")+ <*> option mag (+ help "approximate absolute accuracy goal (default 2^-p)"+ <> long "tol"+ <> value (read "0")+ <> metavar "abstol")+ <*> switch (+ help "run twice (to see overhead of computing nodes)"+ <> long "twice")+ <*> switch (+ help "use heap for subinterval queue"+ <> long "heap")+ <*> option pos (+ help "verbosity level"+ <> short 'v'+ <> long "verbosity"+ <> value 0+ <> metavar "verbosity")+ <*> option pos (+ help "use quadrature degree up to n"+ <> long "deg"+ <> value 0+ <> metavar "degree")+ <*> option pos (+ help "limit number of function evaluations to n"+ <> long "eval"+ <> value 0+ <> metavar "eval")+ <*> option pos (+ help "limit subinterval queue size to n"+ <> long "depth"+ <> value 0+ <> metavar "depth")+ <*> option auto (+ help "number of threads"+ <> long "threads"+ <> value 1 + <> metavar "threads")++-- ReadM parsers ---------------------------------------------------------------++rng :: ReadM Range+rng = eitherReader $ \s -> do+ let result@(Range (a, b)) = read s :: Range+ if 0 <= a && a <= b && b < length description then+ Right result+ else+ Left "could not parse range"+ +pos :: (Read a, Integral a) => ReadM a+pos = eitherReader $ \s -> do+ let result = read s+ if result >= 0 then + Right result+ else+ Left "expected positive number"++between a b = eitherReader $ \s -> do+ let result = read s+ if a <= result && result <= b then + Right result+ else+ Left $ "expected number in range [" ++ show a ++ " .. " ++ show b ++ "]."++mag = eitherReader $ \s -> do+ case readMaybe s of+ Just result -> Right result+ _ -> Left $ "parsing " ++ show s ++ " failed."++-- instances Mag ---------------------------------------------------------------++instance Read Mag where+ readsPrec _ = readP_to_S (mkMag 10 <$> parseArb)++-- instances Range -------------------------------------------------------------++newtype Range = Range (Int, Int)++instance Read Range where+ readsPrec _ = readP_to_S (parseRange <|> parseIndex)++instance Show Range where+ show (Range (a, b)) = "[" ++ show a ++ ":" ++ show b ++ "]"+ +-- parsers --------------------------------------------------------------------++parseArb = do+ (res, _) <- gather $ choice+ [ char '[' *> parseFloat *> pm *> parseFloat <* char ']'+ , char '[' *> pm *> parseFloat <* char ']'+ , char '[' *> parseFloat <* char ']'+ , parseFloat+ ]+ return res+ where pm = skipSpaces *> string "+/-" <* many1 (char ' ')++mkMag prec s = unsafePerformIO $ do+ (result, flag) <- withNewMag $ \m -> do+ withNewArb $ \x -> do+ withCString s $ \s -> do+ arb_set_str x s prec+ arb_get_mag m x+ return result+ +parseFloat = do+ munch (\x -> x == '+' || x == '-')+ choice [nan, inf, num *> e, num]+ where+ nan = string "nan"+ inf = string "inf"+ num = munch1 isNumber *> munch (== '.') *> munch isNumber+ e = do+ char 'e'+ munch (\x -> x == '+' || x == '-')+ munch1 isNumber++parseRange :: ReadP Range +parseRange = do+ a <- read <$> munch1 isNumber+ char ':'+ b <- read <$> munch1 isNumber+ return $ Range (a, b)++parseIndex :: ReadP Range+parseIndex = do+ a <- read <$> munch1 isNumber+ return $ Range (a, a)+ +--------------------------------------------------------------------------------++instance Show AcbCalcIntegrateOpt where+ show x = unsafePerformIO $ do+ (_, result) <- withAcbCalcIntegrateOpt x $ \x -> do+ CAcbCalcIntegrateOpt deg eval depth use_heap verbosity <- peek x+ return $ "options:"+ ++ " deg=" ++ show deg+ ++ " eval=" ++ show eval+ ++ " depth=" ++ show depth+ ++ " heap=" ++ show use_heap+ ++ " verbosity=" ++ show verbosity+ return result+
+ app/keiper_li/Main.hs view
@@ -0,0 +1,158 @@+import System.IO+import System.TimeIt++import Control.Monad +import Options.Applicative++import Foreign.C.Types+import Foreign.C.String+import Foreign.Marshal.Alloc (free) +import Foreign.Marshal.Array (advancePtr)++import Data.Number.Flint++main = timeItNamed "time"+ $ run =<< customExecParser (prefs showHelpOnEmpty) opts where+ desc = "Calculates keipers series."+ opts = info (parameters <**> helper) (+ fullDesc+ <> progDesc desc+ <> header desc)++run params@(Parameters n prec num_threads outFile) = do+ print params+ + flint_set_num_threads num_threads+ + let len = n + 1+ workingPrecision = case prec of+ Just p -> p+ Nothing -> round (1.1 * fromIntegral len + 50)+ print workingPrecision+ z <- _arb_vec_init len++ keiper_li_series z len workingPrecision++ case outFile of+ Just fileName -> do+ mode <- newCString "w"+ out_file_name <- newCString fileName+ fp <- fopen out_file_name mode+ withNewFmpz $ \man -> do+ withNewFmpz $ \exp -> do+ withNewArf $ \t -> do+ forM_ [0 .. fromIntegral len - 1] $ \j -> do+ + arf_get_fmpz_2exp man exp (arb_midref (z `advancePtr` j))+ + withCString (show j ++ " ") $ \s -> do fputs s fp+ fmpz_fprint fp man+ withCString " " $ \s -> fputs s fp+ fmpz_fprint fp exp+ withCString " +/- " $ \s -> fputs s fp+ + arf_set_mag t (arb_radref (z `advancePtr` j))+ arf_get_fmpz_2exp man exp t++ fmpz_fprint fp man+ withCString " " $ \s -> fputs s fp+ fmpz_fprint fp exp+ withCString "\n" $ \s -> fputs s fp+ fclose fp+ mapM_ free [mode, out_file_name]+ return ()+ Nothing -> do + forM_ [0 .. len - 1] $ \j -> do+ putStr $ show j ++ ": "+ arb_printd (z `advancePtr` (fromIntegral j)) 50+ putStr "\n"++-- c file ----------------------------------------------------------------------++foreign import ccall "stdio.h fopen"+ fopen :: CString -> CString -> IO (Ptr CFile)++foreign import ccall "stdio.h fclose"+ fclose :: Ptr CFile -> IO CInt++foreign import ccall "stdio.h fputs"+ fputs :: CString -> Ptr CFile -> IO CInt++data Parameters = Parameters {+ n :: CLong+ , prec :: Maybe CLong+ , num_threads :: CInt+ , outFile :: Maybe String+} deriving Show++parameters :: Parser Parameters+parameters = Parameters+ <$> argument auto (+ help "n"+ <> metavar "n")+ <*> optional optionPrecision+ <*> option auto (+ help "number of threads"+ <> long "threads"+ <> value 1+ <> metavar "threads")+ <*> optional optionOutput++optionPrecision = option auto (+ help "precision"+ <> long "prec"+ <> short 'p'+ <> metavar "precision")++optionOutput = strOption (+ help "output file"+ <> short 'o'+ <> metavar "write output to file")++--------------------------------------------------------------------------------++keiper_li_series z len prec = do++ t <- _arb_vec_init len+ u <- _arb_vec_init len+ v <- _arb_vec_init len++ timeItNamed "zeta: " $ do + arb_zero t+ arb_one (next t)+ arb_one u+ _arb_poly_zeta_series v t 2 u 0 len prec+ _arb_vec_neg v v len++ timeItNamed "log: " $ do+ _arb_poly_log_series t v len len prec++ -- add log(gamma(1+s/2))+ timeItNamed "gamma: " $ do+ arb_one u+ arb_one (next u)+ arb_mul_2exp_si (next u) (next u) (-1)+ _arb_poly_lgamma_series v u 2 len prec+ _arb_vec_add t t v len prec++ -- subtract 0.5 s log(pi)+ arb_const_pi u prec+ arb_log u u prec+ arb_mul_2exp_si u u (-1)+ arb_sub (next t) (next t) u prec++ -- add log(1-s)+ arb_one u+ arb_set_si (next u) (-1)+ _arb_poly_log_series v u 2 len prec+ _arb_vec_add t t v len prec++ timeItNamed "binomial transform: " $ do+ arb_set z t+ _arb_vec_neg (next t) (next t) (len - 1)+ _arb_poly_binomial_transform (next z) (next t) (len - 1) (len - 1) prec+ +next x = x `advancePtr` 1++ +
+ app/l_values/Main.hs view
@@ -0,0 +1,110 @@+import System.TimeIt++import Foreign.Ptr (Ptr, nullPtr, plusPtr)+import Foreign.C.Types+import Foreign.C.String+import Foreign.Storable+import Foreign.Marshal.Array (advancePtr)++import Options.Applicative+import Control.Monad++import Data.Number.Flint++main = timeItNamed "time" $ run =<< execParser opts where+ opts = info (parameters <**> helper) (+ fullDesc+ <> progDesc "Print value of Dirichlet L-function at s = x + i y."+ <> header "Print value of Dirichlet L-function at s = x + i y.")++run p@(Parameters character n xs ys prec zfun deflate len) = do+ case gcd character n of+ 1 -> do let deflate' = if deflate then 1 else 0+ lvalue character n xs ys prec zfun deflate' len+ _ -> putStrLn "need gcd(q, n) == 1 to define a character."+ +lvalue character n xs ys prec zfun deflate len = do+ z <- _acb_vec_init len+ p <- forM [0..len-1] $ \j -> return $ z `advancePtr` (fromIntegral j)+ g <- newDirichletGroup character+ withNewDirichletChar g $ \chi -> do+ withDirichletGroup g $ \g -> do+ putStr $ "character chi(" ++ show character ++ ", " ++ show n ++ ") is "+ dirichlet_char_log chi g n+ primitive <- dirichlet_char_is_primitive g chi+ if primitive == 1 then do+ putStrLn "primitive."+ withNewAcb $ \s -> do+ withNewArb $ \x -> do+ withNewArb $ \y -> do+ withCString xs $ \xs -> arb_set_str x xs prec+ withCString ys $ \ys -> arb_set_str y ys prec+ acb_set_arb_arb s x y+ if zfun then do+ acb_dirichlet_hardy_z z s g chi len prec+ else do+ acb_dirichlet_l_jet z s g chi deflate len prec+ forM_ (zip [0..] p) $ \(j, p) -> do+ let f = if zfun then "Z" else "L"+ prec' = fromIntegral $ floor $ fromIntegral prec / 3 + 1+ case j of+ 0 -> putStr $ f ++ "(s) "+ 1 -> putStr $ f ++ "'(s) "+ _ -> putStrLn $ "[x^" ++ show j ++ "] " ++ f+ acb_printn p prec' arb_str_none+ putStr "\n"+ return ()+ else+ putStrLn "not primitive."+ _acb_vec_clear z (fromIntegral len)++data Parameters = Parameters {+ character :: CULong,+ n :: CULong,+ xs :: String,+ ys :: String,+ prec :: CLong,+ zfun :: Bool,+ deflate :: Bool,+ len :: CLong+} deriving Show++-- default values give L(chi, 1/2)+-- +-- see http://www.lmfdb.org/L/Character/Dirichlet/12/11/+parameters :: Parser Parameters+parameters = Parameters+ <$> option auto (+ long "character" <>+ value 12 <>+ help "Dirichlet character" <>+ metavar "character")+ <*> option auto (+ short 'n' <>+ value 11 <>+ metavar "n")+ <*> strOption (+ short 'x' <>+ value "0.5" <>+ metavar "x")+ <*> strOption (+ short 'y' <>+ value "0.0" <>+ metavar "y")+ <*> option auto (+ short 'p' <>+ long "precision" <>+ value 1024 <>+ metavar "prec")+ <*> switch (+ short 'z')+ <*> switch (+ short 'd' <>+ long "deflate")+ <*> option auto (+ long "len" <>+ short 'l' <>+ value 1 <>+ metavar "len")++
+ app/logistic/Main.hs view
@@ -0,0 +1,101 @@+import Options.Applicative++import Control.Monad+import Control.Monad.State++import Foreign.Ptr+import Foreign.ForeignPtr+import Foreign.C.Types+import Foreign.C.String+import Foreign.Marshal.Array++import Data.Number.Flint++main = run =<< customExecParser (prefs showHelpOnEmpty) opts where+ desc = "Compute nth iterate of the logistic map x_{n+1} = r x_n (1-x_n)."+ opts = info (parameters <**> helper) (+ fullDesc+ <> progDesc desc+ <> header desc)++run params@(Parameters n xs rs digits) = do+ print params+ let goal = round (fromIntegral digits / logBase 10 2 + 3) :: CLong+ withNewArb $ \x -> do+ withNewArb $ \r -> do+ withNewArb $ \s -> do+ withNewArb $ \t -> do+ _ <- runStateT (next n goal) (0, xs, rs, x, r, s, t, 64)+ putStr $ "x_" ++ show n ++ " = "+ arb_printn x digits arb_str_none+ putStr "\n"++next :: CLong -> CLong+ -> StateT (CLong, String, String, Ptr CArb, Ptr CArb, Ptr CArb, Ptr CArb, CLong) IO ()+next n goal = do+ (i, xs, rs, x, r, s, t, prec) <- get+ when (i == 0) $ liftIO $ do+ putStr $ "Trying precision " ++ show prec ++ " bits ... "+ getArb x xs prec+ getArb r rs prec+ if i < n then do + success <- liftIO $ do+ arb_sub_ui t x 1 prec+ arb_neg t t+ arb_mul x x t prec+ arb_mul x x r prec+ p <- arb_rel_accuracy_bits x+ return $ p >= goal+ if success then do+ put (i + 1, xs, rs, x, r, s, t, prec)+ else do+ liftIO $ putStrLn $ "ran out of precision at step " ++ show i+ put (0, xs, rs, x, r, s, t, 2 * prec)+ next n goal+ else do+ liftIO $ putStrLn "success!"+ +getArb x s prec = do+ withCString s $ \cs -> do+ flag <- arb_set_str x cs prec+ is_finite <- arb_is_finite x+ when (flag /= 0 || is_finite /= 1) $ do error $ "could no parse " ++ s+ +data Parameters = Parameters {+ n :: CLong+ , x0 :: String+ , r :: String+ , digits :: CLong+} deriving Show++parameters :: Parser Parameters+parameters = Parameters+ <$> argument pos (+ help "nth iterate of the logistic map."+ <> metavar "n")+ <*> strOption (+ help "starting point."+ <> long "x0"+ <> value "0.5"+ <> metavar "x0")+ <*> strOption (+ help "r parameter of logistic map."+ <> long "r"+ <> value "3.75"+ <> metavar "r")+ <*> option auto (+ help "number of digits."+ <> long "digits"+ <> short 'd'+ <> value 10+ <> metavar "digits")++pos :: (Read a, Integral a) => ReadM a+pos = eitherReader $ \s -> do+ let result = read s+ if result >= 0 then + Right result+ else+ Left "expected positive number"++
+ app/machin/Main.hs view
@@ -0,0 +1,88 @@+import System.TimeIt++import Control.Monad+import Data.List (intercalate)+import Foreign.ForeignPtr+import Foreign.C.Types++import Data.Number.Flint++-- This program checks several variations of Machin’s formula ------------------++main = timeItNamed "machin" $ do+ checkFormula m_formulas simple_ca_atan_p_q "atan"+ checkFormula h_formulas simple_ca_atanh_p_q "atanh"+ +checkFormula formulas f fName = do+ ctx <- newCaCtx+ withNewCa ctx $ \x -> do+ withNewCa ctx $ \y -> do+ withCaCtx ctx $ \ctx -> do+ withNewFmpq $ \arg -> do+ forM_ formulas $ \(rhs, coeffs) -> do+ -- evaluate formula -> x+ ca_zero x ctx+ s <- forM coeffs $ \(c, q) -> do+ f y 1 q ctx+ ca_mul_si y y c ctx+ ca_add x x y ctx+ return $ "(" ++ show c ++ ")*" ++ fName ++ "(1/" ++ show q ++ ")"+ putStr $ "[" ++ intercalate "+" s ++ "] - "+ -- evaluate result -> y+ let (p, q) = rhs+ fmpq_set_si arg p (fromIntegral q)+ putStr $ fName ++ "("; fmpq_print arg; putStr ") = "+ f y p q ctx+ -- formula - result+ ca_sub x x y ctx+ ca_print x ctx; putStr "\n"++simple_ca_atan_p_q res p q ctx = do+ ca_set_si res p ctx+ ca_div_si res res q ctx+ ca_atan res res ctx++simple_ca_atanh_p_q res p q ctx = do+ ca_set_si res p ctx+ ca_div_si res res q ctx+ simple_ca_atanh res res ctx++-- valid for -1 < x < 1+simple_ca_atanh res x ctx = do+ t <- mallocForeignPtr+ u <- mallocForeignPtr+ withForeignPtr t $ \t -> do+ withForeignPtr u $ \u -> do+ ca_init t ctx+ ca_init u ctx+ ca_add_si t x 1 ctx+ ca_sub_si u x 1 ctx+ ca_neg u u ctx+ ca_div res t u ctx+ ca_log res res ctx+ ca_div_si res res 2 ctx+ ca_clear t ctx+ ca_clear u ctx++m_formulas :: [((CLong, CLong), [(CLong, CLong)])]+m_formulas = + [ ((1, 1), [(1, 1)])+ , ((1, 1), [(1, 2), (1, 3)])+ , ((1, 1), [(2, 2), (-1, 7)])+ , ((1, 1), [(2, 3), (1, 7)])+ , ((1, 1), [(4, 5), (-1, 239)])+ , ((1, 1), [(1, 2), (1, 5), (1, 8)])+ , ((1, 1), [(1, 3), (1, 4), (1, 7), (1, 13)])+ , ((1, 1), [(12, 49), (32, 57), (-5, 239), (12, 110443)])+ ]++h_formulas :: [((CLong, CLong), [(CLong, CLong)])]+h_formulas =+ [ ((3, 5), [(14, 31), (10, 49), (6, 161)])+ , ((4, 5), [(22, 31), (16, 49), (10, 161)])+ , ((12, 13), [(32, 31), (24, 49), (14, 161)])+ , ((3, 5), [(144, 251), (54, 449), (-38, 4801), (62, 8749)])+ , ((4, 5), [(228, 251), (86, 449), (-60, 4801), (98, 8749)])+ , ((12, 13), [(334, 251), (126, 449), (-88, 4801), (144, 8749)])+ , ((24, 25), [(404, 251), (152, 449), (-106, 4801), (174, 8749)])+ ]
+ app/multi_crt/Main.hs view
@@ -0,0 +1,57 @@+{-# language ScopedTypeVariables #-}+import Options.Applicative++import Control.Monad+import Foreign.C.Types+import Foreign.Ptr +import Foreign.Marshal.Array+import Foreign.Storable++import Data.Number.Flint++main = run =<< customExecParser (prefs showHelpOnEmpty) opts where+ desc = "Reconstruct integer using the chinese remainder theorem."+ hdesc = "Fast tree version of the integer Chinese Remainder code."+ opts = info (parameters <**> helper) (+ fullDesc+ <> progDesc desc+ <> header desc)++run params@(Parameters n num_primes) = do+ print params+ primes <- mapM n_nth_prime [1..fromIntegral num_primes]+ withArray primes $ \primes -> do+ comb <- newFmpzComb (castPtr primes) (fromIntegral num_primes)+ withFmpzComb comb $ \comb -> do+ comb_temp <- newFmpzCombTemp comb+ withFmpzCombTemp comb_temp $ \comb_temp -> do+ withFmpz n $ \x -> do+ withNewFmpz $ \y -> do+ allocaArray num_primes $ \(residues :: Ptr CLong) -> do+ -- Reduce modulo all primes + fmpz_multi_mod_ui (castPtr residues) x comb comb_temp+ -- Reconstruct+ fmpz_multi_CRT_ui y (castPtr residues) comb comb_temp 1+ forM_ [0 .. fromIntegral num_primes - 1] $ \i -> do+ p <- peek (primes `advancePtr` i)+ r <- peek (residues `advancePtr` i)+ putStrLn $ "residue mod " ++ show p ++ " = " ++ show r+ putStr "reconstruction = "+ fmpz_print y+ putStr "\n"+ +data Parameters = Parameters {+ n :: Fmpz+ , num_primes :: Int+} deriving Show++parameters :: Parser Parameters+parameters = Parameters+ <$> argument auto (+ help "n to be reconstructed"+ <> metavar "n")+ <*> option auto (+ help "number of primes [2, 3, ...] to use"+ <> long "np"+ <> value 1+ <> metavar "num_primes")
+ app/padic/Main.hs view
@@ -0,0 +1,70 @@+import Control.Monad+import Foreign.Ptr +import Foreign.C.String+import Foreign.Marshal.Alloc++import Data.Number.Flint++main = do+ putStrLn "Output:\n"++ -- Case 1++ putStrLn "Positive integer: x = 127 mod 7^10"++ forM_ [padic_terse, padic_series, padic_val_unit] $ \printMode -> do+ withNewPadicCtx 7 8 12 printMode $ \ctx -> do+ withNewPadic $ \x -> do+ padic_init2 x 10+ padic_set_ui x 127 ctx+ putStr "print: "; padic_print x ctx; putStr "\n"+ cstr <- padic_get_str nullPtr x ctx+ str <- peekCString cstr+ free cstr+ putStrLn $ "getStr: " ++ str++ -- Case 2++ putStrLn "Positive integer larger than p^N: x = 1057 mod 2^10"++ forM_ [padic_terse, padic_series, padic_val_unit] $ \printMode -> do+ withNewPadicCtx 2 8 12 printMode $ \ctx -> do+ withNewPadic $ \x -> do+ padic_init2 x 10+ padic_set_ui x 1057 ctx+ putStr "print: "; padic_print x ctx; putStr "\n"+ cstr <- padic_get_str nullPtr x ctx+ str <- peekCString cstr+ free cstr+ putStrLn $ "getStr: " ++ str++ -- Case 3++ putStrLn "Negative integer: x = -127 mod 3^10"++ forM_ [padic_terse, padic_val_unit] $ \printMode -> do+ withNewPadicCtx 2 8 12 printMode $ \ctx -> do+ withNewPadic $ \x -> do+ padic_init2 x 10+ padic_set_ui x (-127) ctx+ putStr "print: "; padic_print x ctx; putStr "\n"+ cstr <- padic_get_str nullPtr x ctx+ str <- peekCString cstr+ free cstr+ putStrLn $ "getStr: " ++ str++ -- Log++ putStrLn "Log of 7380996 mod 5^20"+ + withNewPadicCtx 5 10 25 padic_series $ \ctx -> do+ withNewPadic $ \x -> do+ withNewPadic $ \y -> do+ padic_set_ui x 7380996 ctx+ padic_log y x ctx+ putStr "x = "; padic_print x ctx; putStr "\n"+ putStr "y = "; padic_print y ctx; putStr "\n"++ + +
+ app/partitions/Main.hs view
@@ -0,0 +1,22 @@+import System.Environment+import Text.Read++import Foreign.C.Types++import Data.Number.Flint++main = do+ prog <- getProgName+ args <- getArgs+ case args of+ [arg0] -> do+ case (readMaybe arg0 :: Maybe CULong) of+ Just n -> run n+ _ -> putStrLn $ "usage: " ++ prog ++ " <integer>"+ _ -> putStrLn $ "usage: " ++ prog ++ " <integer>"++run n = do+ result <- newFmpz+ withFmpz result $ \p -> arith_number_of_partitions p n+ putStrLn $ "p(" ++ show n ++ ") = " ++ show result+
+ app/pi_digits/Main.hs view
@@ -0,0 +1,43 @@+import System.TimeIt++import Foreign.Ptr (nullPtr)+import Foreign.C.Types++import Options.Applicative+import Control.Monad++import Data.Number.Flint.Arb++main = timeItNamed "time"+ $ run =<< customExecParser (prefs showHelpOnEmpty) opts where+ opts = info (parameters <**> helper) (+ fullDesc+ <> progDesc "Calculates digits of pi."+ <> header "Calculates digits of pi.")++run (Parameters digits condense) = do+ let prec = floor $ fromIntegral digits / logBase 10 2 + 5+ flag = foldr (+) arb_str_none $+ replicate condense arb_str_condense+ putStrLn $ "computing pi with a precision of " ++ show prec ++ " bits ..."+ withNewArb $ \x -> do+ arb_const_pi x prec+ arb_printn x (fromIntegral digits) flag+ putStr "\n"+ +data Parameters = Parameters {+ digits :: Int,+ condense :: Int+} deriving Show++parameters :: Parser Parameters+parameters = Parameters+ <$> argument auto (+ help "number of digits to calculate." <>+ metavar "digits")+ <*> argument auto (+ value 20 <>+ help "condense" <>+ showDefault <>+ metavar "condense")+
+ app/poly_roots/Main.hs view
@@ -0,0 +1,180 @@+import System.TimeIt++import Foreign.Ptr (Ptr, nullPtr, plusPtr)+import Foreign.C.Types+import Foreign.C.String+import Foreign.Storable+import Foreign.Marshal.Array++import Options.Applicative+import Options.Applicative.Help.Pretty++import Control.Monad+import Data.IORef++import Data.Number.Flint++-- main = putStrLn "done"++main = timeItNamed "time"+ $ run =<< customExecParser (prefs showHelpOnEmpty) opts where+ opts = info (parameters <**> helper) (+ fullDesc+ <> progDesc "Isolates all the complex roots of a polynomial with \+ \integer coefficients.."+ <> header "Roots of polynomials")++run params@(Parameters digitsRefine digitsPrint polys) = do+ let prec = case digitsRefine of+ Just digits -> round $ fromIntegral digits / logBase 10 2 + 2+ Nothing -> 16+ print params+ if not (null polys) then do+ pols <- product <$> mapM getPolynomial polys+ forM_ (factor pols) $ \(poly, e) -> do+ putStrLn $ "polynomial expression " ++ show poly + ++ " with multiplicity " ++ show e + withFmpzPoly poly $ \p -> do+ n <- fmpz_poly_degree p+ putStrLn $ "polynomial degree: " ++ show n+ z <- _acb_vec_init n+ arb_fmpz_poly_complex_roots z p arb_fmpz_poly_roots_verbose prec+ case digitsPrint of+ Just digits -> do+ forM_ [0..fromIntegral n-1] $ \j -> do+ acb_printn (z.+.j) digits arb_str_no_radius+ putStr "\n"+ Nothing -> return ()+ _acb_vec_clear z n+ return ()+ else do+ putStrLn $ "no polynomial expression given."++(.+.) x y = x `advancePtr` (fromIntegral y)++getPolynomial :: Polynomial -> IO FmpzPoly+getPolynomial x =+ case x of+ A n -> easyA n+ T n -> chebyshevT n+ U n -> chebyshevT n+ P n -> legendreP' n+ C n -> cyclotomicC n+ S n -> swinnertonDyerS n+ B n -> bernoulliB' n++data Polynomial+ = A CULong+ | T CULong+ | U CULong+ | P CULong+ | C CULong+ | S CULong+ | B CULong+ | W CULong+ | E CULong+ | M CULong+ | Coeffs [Integer]+ deriving Show++-- Parser ----------------------------------------------------------------------++data Parameters = Parameters {+ digitsRefine :: Maybe CLong+, digitsPrint :: Maybe CLong+, poly :: [Polynomial]+} deriving Show++parameters :: Parser Parameters+parameters = Parameters+ <$> optional ( option auto (+ help "the roots are refined to a relative tolerance\+ \better than 10^(-d). By default, the roots are \+ \only computed to sufficient\+ \accuracy to isolate them. The refinement is not\+ \currently done efficiently"+ <> long "refine"+ <> short 'r'+ <> metavar "RELATIVE-TOLERANCE"))+ <*> optional ( option auto (+ help "computed roots are printed to d decimals. \+ \By default, the roots are not printed."+ <> short 'd'+ <> long "print"+ <> metavar "DIGITS") )+ <*> many polynomial++polynomial :: Parser Polynomial+polynomial = subparser cmds where+ cmds = mconcat $ zipWith mkCmd desc [a, t, u, p, phi, s, b, w, e, m, coeffs]+ mkCmd x y = command (takeWhile (/=' ') x) (info y (progDesc x))++desc =+ [ "a <n> Easy polynomial 1 + 2x + ... + (n+1)x^n"+ , "t <n> Chebyshev polynomial T_n"+ , "u <n> Chebyshev polynomial U_n"+ , "p <n> Legendre polynomial P_n"+ , "c <n> Cyclotomic polynomial Phi_n"+ , "s <n> Swinnerton-Dyer polynomial S_n"+ , "b <n> Bernoulli polynomial B_n"+ , "w <n> Wilkinson polynomial W_n"+ , "e <n> Taylor series of exp(x) truncated to degree n"+ , "m <n> <m> The Mignotte-like polynomial x^n + (100x+1)^m, n > m"+ , "coeffs <c0 c1 ... cn> c0 + c1 x + ... + cn x^n"+ ]++a, t, u, p, phi, s, b, w, e, m, coeffs :: Parser Polynomial++a = A <$> argument auto (metavar "n")+t = T <$> argument auto (metavar "n")+u = U <$> argument auto (metavar "n")+p = P <$> argument auto (metavar "n")+phi = C <$> argument auto (metavar "n")+s = S <$> argument auto (metavar "n")+b = B <$> argument auto (metavar "n")+w = W <$> argument auto (metavar "n")+e = E <$> argument auto (metavar "n")+m = M <$> argument auto (metavar "n")+coeffs = Coeffs <$> many (argument auto (metavar "coeffs"))++-- polynomials -----------------------------------------------------------------++-- Legendre polynomial (denominator removed)+legendreP' n = do+ (poly, _) <- withNewFmpzPoly $ \p -> do+ withNewFmpqPoly $ \h -> do+ fmpq_poly_legendre_p h n+ fmpq_poly_get_numerator p h+ return poly++chebyshevT n = do+ (poly, _) <- withNewFmpzPoly $ \p -> do+ fmpz_poly_chebyshev_t p n+ return poly++chebyshevU n = do+ (poly, _) <- withNewFmpzPoly $ \p -> do+ fmpz_poly_chebyshev_u p n+ return poly++cyclotomicC n = do+ (poly, _) <- withNewFmpzPoly $ \p -> do+ fmpz_poly_cyclotomic p n+ return poly++swinnertonDyerS n = do+ (poly, _) <- withNewFmpzPoly $ \p -> do+ fmpz_poly_swinnerton_dyer p n+ return poly++-- Bernoulli polynomial (denominator removed)+bernoulliB' n = do+ (poly, _) <- withNewFmpzPoly $ \p -> do+ withNewFmpqPoly $ \h -> do+ arith_bernoulli_polynomial h n+ fmpq_poly_get_numerator p h+ return poly++easyA n = do+ let poly = fromList (map fromIntegral [1..n]) :: FmpzPoly+ return poly
+ app/primegen/Main.hs view
@@ -0,0 +1,55 @@+import Options.Applicative+import Control.Monad+import Control.Monad.State++import Foreign.C.Types+import Foreign.Ptr++import Data.Number.Flint++main = run =<< customExecParser (prefs showHelpOnEmpty) opts where+ opts = info (parameters <**> helper) (+ fullDesc+ <> progDesc "Generate primes."+ <> header "Generate primes.")++run params@(Parameters n count) = do+ prime_iter <- newNPrimes+ if count then do+ (_, (prime_pi, _)) <- runStateT (countPrimes params) (0, prime_iter)+ print prime_pi+ else do+ _ <- runStateT (printPrimes params) prime_iter+ return ()+ +printPrimes :: Parameters -> StateT NPrimes IO ()+printPrimes params@(Parameters n count) = do+ prime_iter <- get+ (_, prime) <- liftIO $ withNPrimes prime_iter n_primes_next+ when (prime < n) $ do+ liftIO $ print prime+ put prime_iter+ printPrimes params++countPrimes :: Parameters -> StateT (Int, NPrimes) IO ()+countPrimes params@(Parameters n count) = do+ (nprimes, prime_iter) <- get+ (_, prime) <- liftIO $ withNPrimes prime_iter n_primes_next+ when (prime < n) $ do+ put (succ nprimes, prime_iter)+ countPrimes params+ +data Parameters = Parameters {+ n :: CULong+ , count :: Bool+} deriving Show++parameters :: Parser Parameters+parameters = Parameters+ <$> argument auto (+ help "n"+ <> metavar "n")+ <*> switch (+ help "count primes"+ <> short 'c'+ <> long "count")
+ app/qadic/Main.hs view
@@ -0,0 +1,94 @@+import Foreign.Storable+import Foreign.Marshal.Array (advancePtr)++import Data.Number.Flint++main = do++ let n = 5++ putStrLn "Compute a power and a sum"+ + withNewQadicCtxConway 3 2 0 n "a" padic_series $ \ctx -> do+ CQadicCtx pctx _ _ _ _ <- peek ctx+ withNewQadicWithPrec n $ \a -> do+ withNewQadicWithPrec n $ \b -> do+ withNewQadicWithPrec n $ \c -> do++ withFmpzPoly (fromList [1, 2]) $ \poly -> do+ padic_poly_set_fmpz_poly a poly pctx+ qadic_print_pretty a ctx; putStr "\n"+ withFmpz 4 $ \four -> qadic_pow a a four ctx+ padic_poly_set_ui b 3249 pctx+ qadic_add c a b ctx+ + qadic_print_pretty a ctx; putStr "\n"+ qadic_print_pretty b ctx; putStr "\n"+ qadic_print_pretty c ctx; putStr "\n"+ putStr "\n"+ + putStrLn "Compute a Teichmuller lift"+ + withFmpzPoly (fromList [1, 2]) $ \poly -> do+ padic_poly_set_fmpz_poly a poly pctx+ qadic_teichmuller b a ctx+ withFmpz 9 $ \nine -> qadic_pow c b nine ctx+ + qadic_print_pretty a ctx; putStr "\n"+ qadic_print_pretty b ctx; putStr "\n"+ qadic_print_pretty c ctx; putStr "\n"+ putStr "\n"+ + putStrLn "Compute an inverse"++ qadic_set a b+ qadic_inv b a ctx+ qadic_mul c a b ctx++ qadic_print_pretty a ctx; putStr "\n"+ qadic_print_pretty b ctx; putStr "\n"+ qadic_print_pretty c ctx; putStr "\n"+ putStr "\n"++ ------------------------------------------------------------------------------++ putStrLn "Compute a Frobenius image"++ withNewQadicCtxConway 3 2 0 n "X" padic_terse $ \ctx -> do+ CQadicCtx pctx _ _ _ _ <- peek ctx+ withNewQadicWithPrec n $ \a -> do+ withNewQadicWithPrec n $ \b -> do+ + withFmpzPoly (fromList [78, 45]) $ \poly -> do+ padic_poly_set_fmpz_poly a poly pctx++ qadic_frobenius b a 1 ctx+ + putStr "a = "; qadic_print_pretty a ctx; putStr "\n"+ putStr "b = "; qadic_print_pretty b ctx; putStr "\n"+ putStrLn "Context:"; qadic_ctx_print ctx; putStr "\n"+ putStr "\n"+ + ------------------------------------------------------------------------------++ putStrLn "Compute a square root"++ withNewQadicCtxConway 2 3 0 n "X" padic_series $ \ctx -> do+ CQadicCtx pctx _ _ _ _ <- peek ctx+ withNewQadicWithPrec n $ \a -> do+ withNewQadicWithPrec n $ \b -> do++ withFmpzPoly (fromList [1, 3, 1]) $ \poly -> do+ padic_poly_set_fmpz_poly a poly pctx++ ans <- qadic_sqrt b a ctx++ putStr "a = "; qadic_print_pretty a ctx; putStr "\n"+ putStr "b = "; qadic_print_pretty b ctx; putStr "\n"+ putStrLn $ "ans = " ++ show ans+ putStrLn "Context:"; qadic_ctx_print ctx; putStr "\n"+ putStr "\n"+ ++ +
+ app/radix/Main.hs view
@@ -0,0 +1,130 @@+import System.CPUTime+import Foreign.Marshal.Alloc+import Foreign.Marshal.Array+import Foreign.C.Types+import Foreign.Storable+import Control.Monad+import Text.Printf++import Data.Number.Flint++main = do++ let n = 12376+ n_poly = (fromIntegral n) `div` 26+ m = 17^26 :: Fmpz++ state <- newFRandState+ + ctx <- newFmpzModCtx m+ + r <- newFmpzModPoly ctx+ t <- newFmpzModPoly ctx++ makePolynomials n m state r t ctx+ radixInfo n m r t ctx++ radix <- timeItNamed "precomputation" $ do+ newFmpzModPolyRadix r (n + 1) ctx++ b <- _fmpz_mod_poly_vec_init (fromIntegral n + 1) ctx++ timeItNamed "conversion" $ do+ withFmpzModPoly t $ \t -> do+ withFmpzModCtx ctx $ \ctx -> do+ withFmpzModPolyRadix radix $ \radix -> do+ fmpz_mod_poly_radix b t radix ctx+ -- forM_ [0..10] $ \j -> do+ -- p <- peek (b `advancePtr` j)+ -- fmpz_mod_poly_print p ctx+ -- putStr "\n"++ _fmpz_mod_poly_vec_clear b (n_poly + 1) ctx+ +radixInfo n m r t ctx = do+ withFmpzModPoly r $ \r -> do+ withFmpzModPoly t $ \t -> do+ withFmpzModCtx ctx $ \ctx -> do++ deg <- fmpz_mod_poly_degree r ctx + polDeg <- fmpz_mod_poly_degree t ctx+ (_, bits) <- withFmpz m fmpz_bits++ putStrLn "Radix conversion\n\+ \-----------------"+ putStrLn $ " Degree of the radix: " ++ show deg+ putStrLn $ " Bit size of the modulus: " ++ show bits+ putStrLn $ " Degree of the input: " ++ show polDeg+ putStr "\n"+ return ()++makePolynomials n m state r t ctx = do+ withFmpzModPoly r $ \r -> do+ withFmpzModPoly t $ \t -> do+ withNewFmpzModPoly ctx $ \u -> do+ withNewFmpzModPoly ctx $ \v -> do+ withFmpzModCtx ctx $ \ctx -> do+ + fmpz_mod_poly_set_coeff_ui u 3 5 ctx+ fmpz_mod_poly_set_coeff_ui u 4 4 ctx++ fmpz_mod_poly_set_coeff_ui v 0 1 ctx+ fmpz_mod_poly_set_coeff_ui v 2 1 ctx+ fmpz_mod_poly_set_coeff_ui v 3 5 ctx+ fmpz_mod_poly_set_coeff_ui v 4 1 ctx+ fmpz_mod_poly_set_coeff_ui v 5 5 ctx+ fmpz_mod_poly_set_coeff_ui v 8 8 ctx+ fmpz_mod_poly_set_coeff_ui v 9 8 ctx+ fmpz_mod_poly_set_coeff_ui v 10 5 ctx+ fmpz_mod_poly_set_coeff_ui v 12 6 ctx+ fmpz_mod_poly_set_coeff_ui v 13 1 ctx++ fmpz_mod_poly_pow r u 3 ctx++ withNewFmpz $ \a -> do+ fmpz_set_ui a 4+ fmpz_mod_poly_scalar_mul_fmpz r r a ctx+ fmpz_set_ui a 27+ fmpz_mod_poly_scalar_mul_fmpz t t a ctx++ fmpz_mod_poly_add r r t ctx++ withFRandState state $ \state -> do+ fmpz_mod_poly_randtest t state (fromIntegral n + 1) ctx+ return ()++timeItNamed :: String -> IO a -> IO a+timeItNamed s f = do+ t0 <- getCPUTime+ result <- f + t1 <- getCPUTime+ let dt = fromIntegral (t1 - t0) / 10^9 :: Double+ printf "%s: %.6f ms\n" s dt+ return result++_fmpz_mod_poly_vec_init :: Int -> FmpzModCtx -> IO (Ptr (Ptr CFmpzModPoly))+_fmpz_mod_poly_vec_init n ctx = do+ vec <- mallocArray n :: IO (Ptr (Ptr CFmpzModPoly))+ result <- forM [0 .. fromIntegral n - 1] $ \j -> do+ poly <- malloc :: IO (Ptr CFmpzModPoly)+ withFmpzModCtx ctx $ \ctx -> do+ fmpz_mod_poly_init poly ctx+ return poly+ pokeArray vec result+ return vec++_fmpz_mod_poly_vec_clear :: (Ptr (Ptr CFmpzModPoly))+ -> Int -> FmpzModCtx -> IO ()+_fmpz_mod_poly_vec_clear vec n ctx = do+ withFmpzModCtx ctx $ \ctx -> do+ forM [0 .. fromIntegral n - 1] $ \j -> do+ p <- peek (vec `advancePtr` j)+ fmpz_mod_poly_clear p ctx+ return ()+ ++ ++ ++
+ app/real_roots/FFI.hsc view
@@ -0,0 +1,27 @@+module FFI where++import Data.Number.Flint++import Foreign.Ptr+import Foreign.ForeignPtr+import Foreign.Marshal.Array+import Foreign.C.Types+import Foreign.Storable++#include "z_param.h"++data ZParam = ZParam {-# UNPACK #-} !(ForeignPtr CZParam)+data CZParam = CZParam (Ptr CDirichletGroup) (Ptr CDirichletChar)++instance Storable CZParam where+ {-# INLINE sizeOf #-}+ sizeOf _ = #{size z_param_t}+ {-# INLINE alignment #-}+ alignment _ = #{alignment z_param_t}+ peek ptr = CZParam+ <$> #{peek z_param_struct, G } ptr + <*> #{peek z_param_struct, chi} ptr + poke ptr (CZParam g c) = do+ #{poke z_param_struct, G } ptr g+ #{poke z_param_struct, chi} ptr c+
+ app/real_roots/Functions.hs view
@@ -0,0 +1,143 @@+module Functions (+ z_function+ , sin_x+ , sin_x2+ , sin_1x+ , airy+ , makeFunPtr+ , F+) where++import Foreign.Ptr+import Foreign.ForeignPtr+import Foreign.Marshal.Array+import Foreign.C.Types+import Foreign.Storable++import Control.Monad++import Data.Number.Flint hiding (airy)++import FFI++foreign import ccall safe "wrapper"+ makeFunPtr :: F -> IO (FunPtr F)++type F = CArbCalcFunc+ +z_function, sin_x, sin_x2, sin_1x, airy :: F++z_function out inp params order prec = do++ CZParam g chi <- peek (castPtr params :: Ptr CZParam)++ when (g /= nullPtr) $ do+ putStr "char: "+ dirichlet_char_print g chi+ putStr "\n"+ -- putStrLn $ "g = " ++ show g+ -- putStrLn $ "chi = " ++ show chi+ + if g == nullPtr then do+ x <- _arb_vec_init 2+ arb_set x inp+ arb_one (x .+. 1)+ _arb_poly_riemann_siegel_z_series out x (min 2 order) order prec+ _arb_vec_clear x 2+ else do+ tmp <- _acb_vec_init order+ acb_set_arb tmp inp+ acb_dirichlet_hardy_z tmp tmp g chi order prec+ forM_ [0 .. order - 1] $ \k -> do+ arb_set (out .+. k) (acb_realref (tmp .+. k))+ _acb_vec_clear tmp order++ return 0+++sin_x out inp params order prec = do++ let xlen = min 2 order+ + arb_set out inp+ when (xlen > 1) $ do arb_one (out .+. 1)+ _arb_poly_sin_series out out xlen order prec++ return 0+ +sin_x2 out inp params order prec = do++ let xlen = min 2 order+ ylen = min 3 order+ + x <- _arb_vec_init xlen+ _arb_poly_mullow out x xlen x xlen ylen prec+ _arb_poly_sin_series out out ylen order prec+ _arb_vec_clear x xlen++ return 0+++sin_1x out inp params order prec = do++ let xlen = min 2 order++ x <- _arb_vec_init xlen++ arb_set x inp+ when (xlen > 1) $ do arb_one (x .+. 1)++ _arb_poly_inv_series out x xlen order prec+ _arb_poly_sin_series out out order order prec++ _arb_vec_clear x xlen++ return 0++ +airy out inp params order prec = do++ let xlen = min 2 order++ which <- peek (castPtr params) :: IO CInt++ -- putStrLn $ "airy: which = " ++ show which++ withNewAcb $ \t -> do+ withNewAcb $ \u -> do+ acb_set_arb t inp++ if xlen == 1 then do+ case which of + 0 -> acb_hypgeom_airy t nullPtr nullPtr nullPtr t prec+ 1 -> acb_hypgeom_airy nullPtr t nullPtr nullPtr t prec+ 2 -> acb_hypgeom_airy nullPtr nullPtr t nullPtr t prec+ _ -> acb_hypgeom_airy nullPtr nullPtr nullPtr t t prec+ arb_set out (acb_realref t)+ else do + if which == 0 || which == 1 then do+ acb_hypgeom_airy t u nullPtr nullPtr t prec+ else do + acb_hypgeom_airy nullPtr nullPtr t u t prec+ if which == 0 || which == 2 then do+ arb_set (out .+. 0) (acb_realref t)+ arb_set (out .+. 1) (acb_realref u)+ -- f'' z = z f z + when (xlen == 3) $ do arb_mul (out .+. 2) out inp prec+ else do+ arb_set (out .+. 0) (acb_realref u)+ arb_mul (out .+. 1) (acb_realref t) inp prec+ -- f''' z = f z + z f' z + when (xlen == 3) $ do+ arb_mul (out .+. 2) out inp prec+ arb_add (out .+. 2) (out .+. 2) (acb_realref t) prec++ return 0++-- utility functions -----------------------------------------------------------++(.+.) x y = x `advancePtr` (fromIntegral y)++next x = do+ value <- peek x+ poke x (value + 1)
+ app/real_roots/Main.hs view
@@ -0,0 +1,168 @@+{-# language OverloadedStrings #-}++import Data.Typeable++import Options.Applicative+import Options.Applicative.Help.Pretty hiding (char)+import Text.ParserCombinators.ReadP hiding (option, between, optional)++import Control.Monad+import Control.Monad.State++import Foreign.Ptr+import Foreign.ForeignPtr+import Foreign.C.Types+import Foreign.Marshal.Array+import Foreign.Marshal.Alloc+import Foreign.Storable++import Data.Char+import Data.Number.Flint hiding (airy)++import Functions+import FFI+import Run++main :: IO ()+main = run =<< customExecParser (prefs showHelpOnEmpty) opts where+ opts = info (parameters <**> helper) + ( fullDesc+ <> header "\nRoot finding (real roots).\n"+ <> progDescDoc (Just desc))+ +desc = vcat funs+funs = + [ "Examples of root finding for some examples (0..7)."+ , " 0 Z(x), Z-function (Riemann zeta or Dirichlet L-function)"+ , " 1 sin(x)"+ , " 2 sin(x^2)"+ , " 3 sin(1/x)"+ , " 4 Ai(x), Airy function"+ , " 5 Ai'(x), Airy function"+ , " 6 Bi(x), Airy function"+ , " 7 Bi'(x), Airy function"+ , "With 0, specify optional Dirichlet character with \+ \[-character q n]"+ ]++-- Parser ----------------------------------------------------------------------++parameters :: Parser Parameters+parameters = Parameters+ <$> argument (between 0 7) (+ help "function"+ <> metavar "FUNCTION")+ <*> option auto (+ help "a"+ <> long "xa"+ <> metavar "A")+ <*> option auto (+ help "b"+ <> long "xb"+ <> metavar "B")+ <*> optional optionDC + <*> option auto (+ help "refine "+ <> long "refine"+ <> value 0+ <> metavar "DIGITS")+ <*> option pos (+ help "verbose"+ <> long "verbose"+ <> short 'v'+ <> value 0+ <> metavar "VERBOSE")+ <*> option auto (+ help "maximal depth"+ <> long "maxdepth"+ <> value 30+ <> metavar "MAX-DEPTH")+ <*> option auto (+ help "maximal number of evaluations"+ <> long "maxeval"+ <> value 100000+ <> metavar "MAX-EVAL")+ <*> option auto (+ help "maximal number of zeros to be found"+ <> long "maxfound"+ <> value 100000+ <> metavar "MAX-FOUND")+ <*> option auto (+ help "precision"+ <> long "precision"+ <> short 'p'+ <> value 64+ <> metavar "PRECISION")++optionDC = option dc (+ help "Dirichlet character for L-function. Specify in the form q.r, \+ \e.g. 12.7"+ <> long "character")+ +pos :: (Read a, Integral a) => ReadM a+pos = eitherReader $ \s -> do+ let result = read s+ if result >= 0 then + Right result+ else+ Left "expected positive number"++between a b = eitherReader $ \s -> do+ let result = read s+ if a <= result && result <= b then + Right result+ else+ Left $ "expected number in range [" ++ show a ++ " .. " ++ show b ++ "]."++-- read Dirichlet Character ----------------------------------------------------++dc :: ReadM DC+dc = eitherReader $ \s -> do+ let chi@(DC p r) = read s+ if p > 0 && r >= 0 then+ Right chi+ else+ Left "expected Dirichlet character, specified as p.r"++instance Read DC where+ readsPrec _ = readP_to_S parseChar++parseChar = do+ p <- read <$> munch1 isNumber+ char '.'+ q <- read <$> munch1 isNumber+ return $ DC p q+++-- testing ---------------------------------------------------------------------++testZParam = do++ g <- newDirichletGroup 12+ c <- newDirichletChar g++ withDirichletGroup g $ \g -> do+ withDirichletChar c $ \c -> do+ print g+ print c+ dirichlet_char_print g c+ putStr "\n"+ params <- malloc :: IO (Ptr CZParam)+ poke params (CZParam g c)+ putStrLn "alloc works."+ let ptr = castPtr params :: Ptr ()+ print ptr+ CZParam u v <- peek (castPtr ptr :: Ptr CZParam)+ putStrLn "new pointer."+ print u+ print v+ dirichlet_char_print u v+ putStr "\n"+ poke params (CZParam nullPtr nullPtr)+ CZParam u v <- peek (castPtr ptr :: Ptr CZParam)+ putStrLn "new pointer."+ print u+ print v+ putStr "\n"+ free params+ return ()
+ app/real_roots/Run.hs view
@@ -0,0 +1,152 @@+module Run where++import Control.Monad+import Control.Monad.State++import Foreign.Ptr+import Foreign.ForeignPtr+import Foreign.C.Types+import Foreign.Marshal.Array+import Foreign.Marshal.Alloc+import Foreign.Storable++import Data.Char+import Data.Number.Flint hiding (airy)++import Functions+import FFI++data DC = DC CULong CULong deriving Show++data Parameters = Parameters {+ function :: Int+ , a :: CDouble+ , b :: CDouble+ , character :: Maybe DC+ , refine :: CInt+ , verbose :: Int+ , maxDepth :: CLong+ , maxEval :: CLong+ , maxFound :: CLong+ , prec :: CLong+ } deriving Show++run params@(Parameters function a b character refine verbose+ maxDepth maxEval maxFound prec) = do+ + let digits = fromIntegral (max 2 refine)+ low_prec = prec+ high_prec = round $ fromIntegral digits / logBase 10 2 + 10++ putStrLn $ "maxDepth = " ++ show maxDepth+ putStrLn $ "maxEval = " ++ show maxEval+ putStrLn $ "low_prec = " ++ show low_prec++ -- here we don't use the wrapper to keep track of allocation+ + [t, interval] <- replicateM 2 malloc; mapM_ arf_interval_init [t, interval]+ [v, w, z] <- replicateM 3 malloc; mapM_ arb_init [v, w, z]+ c <- malloc; arf_init c+ + info <- malloc+ blocks <- malloc+ params <- malloc :: IO (Ptr ())++ g <- newDirichletGroup (case character of+ Just (DC p r) -> p+ Nothing -> 1)+ gx <- newDirichletChar g+ + (_,(_, z_params)) <- withDirichletGroup g $ \g -> do+ withDirichletChar gx $ \gx -> do+ dirichlet_char_log gx g (case character of+ Just (DC p r) -> r+ Nothing -> 1)+ return (CZParam g gx)+ + fp <- makeFunPtr ([z_function, sin_x, sin_x2, sin_1x, airy] !! function)++ arf_interval_set_d_d interval a b++ -- parmeters+ + case function of+ 0 -> case character of+ Just _ -> poke (castPtr params) z_params+ _ -> poke (castPtr params) (CZParam nullPtr nullPtr)+ 1 -> return ()+ 2 -> return ()+ 3 -> return ()+ 4 -> poke (castPtr params) (0 :: CInt)+ 5 -> poke (castPtr params) (1 :: CInt)+ 6 -> poke (castPtr params) (2 :: CInt)+ 7 -> poke (castPtr params) (3 :: CInt)++ num <- arb_calc_isolate_roots blocks info fp+ (castPtr params) interval maxDepth+ maxEval maxFound low_prec+ bp <- peek blocks+ ip <- peek info++ putStrLn $ "Isolated " ++ show num ++ " roots in intervals"++ res <- forM [0.. fromIntegral num-1] $ \i -> do+ info_i <- peek (ip .+. i)+ when (info_i /= 1 && verbose > 0 ) $ do + arf_interval_printd (bp.+.i) 15; endl+ when (refine > 0) $ do+ -- bisection -------------------------------------------------------------++ res <- arb_calc_refine_root_bisect t fp params (bp.+.i) 5 low_prec+ + when (res /= arb_calc_success) $ do+ putStrLn "warning: some bisection steps failed!"+ when (verbose > 0) $ do+ putStr "after bisection 1: "; arf_interval_printd t 15; endl+ + res <- arb_calc_refine_root_bisect (bp.+.i) fp params t 5 low_prec+ when (res /= arb_calc_success) $ do+ putStrLn "warning: some bisection steps failed!"+ when (verbose > 0) $ do+ putStr "after bisection 2: "; arf_interval_printd (bp.+.i) 15; endl+ + --- Newton iteration -----------------------------------------------------+ arf_interval_get_arb v t high_prec+ arb_calc_newton_conv_factor c fp params v low_prec+ + arf_interval_get_arb w (bp.+.i) high_prec+ + res <- arb_calc_refine_root_newton z fp params w v c 10 high_prec+ when (res /= arb_calc_success) $ do+ putStrLn "warning: some newton steps failed!"+ + putStrLn $ "refined root (" ++ show i ++ "/" ++ show num ++ ")"+ arb_printn z (fromIntegral digits + 2) 0+ putStr "\n\n"+ -- count zeros+ return $ if info_i == 1 then (1, 0) else (0, 1)++ let (found, unknown) = foldr (\(x, y) (ax, ay) -> (x+ax, y+ay)) (0, 0) res+ putStr $ "Found " ++ show found ++ " roots"+ if unknown > 0 then do+ putStrLn $ ", where " ++ show unknown+ ++ "intervals contained undetected roots."+ else do+ putStr ".\n"++ mapM_ (\x -> do free x; arf_interval_clear x) [t, interval]+ mapM_ (\x -> do free x; arb_clear x) [v, w, z]+ arf_clear c; free c++ free params+ free info+ free blocks++-- utility functions -----------------------------------------------------------++arf_interval_set_d_d x a b = do+ arf_set_d (castPtr x) a+ arf_set_d (castPtr x .+. 1) b++(.+.) x y = x `advancePtr` fromIntegral y+endl = putStr "\n"
+ app/real_roots/z_param.h view
@@ -0,0 +1,16 @@+#ifndef Z_PARAM_H_+#define Z_PARAM_H_++#include <flint/flint.h>+#include <flint/dirichlet.h>++typedef struct+{+ dirichlet_group_t *G;+ dirichlet_char_t *chi;+}+z_param_struct;++typedef z_param_struct z_param_t[1];++#endif // Z_PARAM_H_
+ app/stirling_matrix/Main.hs view
@@ -0,0 +1,31 @@+import System.Environment+import Text.Read++import Foreign.C.Types++import Data.Number.Flint++main = do+ prog <- getProgName+ args <- getArgs+ case args of+ [arg0] -> do+ case (readMaybe arg0 :: Maybe CLong) of+ Just n -> run n+ _ -> putStrLn $ "usage: " ++ prog ++ " <integer>"+ _ -> putStrLn $ "usage: " ++ prog ++ " <integer>"++run n = do+ s1 <- newFmpzMat n n+ s2 <- newFmpzMat n n+ withFmpzMat s1 $ \s1 -> arith_stirling_matrix_1 s1+ withFmpzMat s2 $ \s2 -> arith_stirling_matrix_2 s2+ let p = s1*s2+ putStrLn "S1 [Stirling numbers of 1st kind]:"+ print s1+ putStrLn "S2 [Stirling numbers of 2nd kind]:"+ print s2+ putStrLn "S1 * S2:"+ print p+ +
+ app/swinnerton_dyer_poly/Main.hs view
@@ -0,0 +1,69 @@+import Control.Monad+import Data.Word+import Data.Bits++import System.TimeIt+import Options.Applicative++import Foreign.ForeignPtr+import Foreign.C.Types+import Foreign.Marshal.Array++import Data.Number.Flint++main = timeItNamed "time"+ $ run =<< customExecParser (prefs showHelpOnEmpty) opts where+ desc = "Computes the coefficients of the Swinnerton-Dyer polynomial."+ opts = info (parameters <**> helper) (+ fullDesc+ <> progDesc desc+ <> header desc)++run params@(Parameters n) = swinnerton_dyer_poly n+ +-- Parser Parameters -----------------------------------------------------------++newtype Parameters = Parameters {+ n :: CULong+ } deriving (Show, Eq)++parameters :: Parser Parameters+parameters = Parameters+ <$> argument (0 `between` 20) (+ help "n"+ <> metavar "n")++between a b = eitherReader $ \s -> do+ let result = read s+ if a <= result && result <= b then + Right result+ else+ Left $ "expected number in range [" ++ show a ++ " .. " ++ show b ++ "]."++--------------------------------------------------------------------------------++swinnerton_dyer_poly n = do+ ctx <- newCaCtx+ poly <- newCaPoly ctx+ withCaPoly poly $ \poly -> do+ withNewCaPoly ctx $ \u -> do+ withNewCaPoly ctx $ \v -> do + withNewCaPoly ctx $ \tmp -> do+ withNewCa ctx $ \w -> do+ withCaCtx ctx $ \ctx -> do+ ca_poly_x poly ctx+ forM_ [1 .. n] $ \j -> do+ p <- n_nth_prime j + ca_set_si w (fromIntegral p) ctx+ ca_sqrt w w ctx+ ca_poly_x tmp ctx+ ca_poly_set_coeff_ca tmp 0 w ctx+ ca_poly_compose u poly tmp ctx+ ca_neg w w ctx+ ca_poly_set_coeff_ca tmp 0 w ctx+ ca_poly_compose v poly tmp ctx+ ca_poly_mul poly u v ctx+ withCaCtx ctx $ \ctx -> do+ withCaPoly poly $ \poly -> do+ ca_poly_print poly ctx; putStr "\n"+
+ app/taylor_integrals/Main.hs view
@@ -0,0 +1,137 @@+import Options.Applicative+import Control.Monad++import Foreign.Ptr+import Foreign.ForeignPtr+import Foreign.C.Types+import Foreign.Marshal.Array++import Data.Number.Flint++main = run =<< customExecParser (prefs showHelpOnEmpty) opts where+ desc = "Compute integrals using d decimal digits of precision."+ opts = info (parameters <**> helper) (+ fullDesc+ <> progDesc desc+ <> header desc)++run params@(Parameters digits) = do+ print params+ let goal = round (fromIntegral digits / logBase 10 2)+ prec = round (1.1 * fromIntegral goal)+ withNewAcb $ \r -> do+ withNewAcb $ \s -> do+ withNewAcb $ \a -> do+ withNewAcb $ \b -> do+ withNewArf $ \inr -> do+ withNewArf $ \outr -> do+ -- Sin integrals+ putStrLn $ replicate 64 '-'+ putStrLn "Integral of sin(t) from 0 to 100."+ arf_set_d inr 0.125+ arf_set_d outr 1.0+ acb_set_si a 0+ acb_set_si b 100+ f <- makeFunPtr sinx+ acb_calc_integrate_taylor r f nullPtr a b inr outr goal prec+ putStrLn "RESULT:"+ acb_printn r digits 0; putStr "\n"+ -- Elliptic integral+ putStrLn $ replicate 64 '-'+ putStrLn "Elliptic integral F(phi, m) = integral of \+ \1/sqrt(1 - m*sin(t)^2)"+ arf_set_d inr 0.125+ arf_set_d outr 1.0+ acb_set_si a 0+ acb_set_si b 6+ f <- makeFunPtr elliptic+ acb_calc_integrate_taylor r f nullPtr a b inr outr goal prec+ acb_set_si a 6+ arb_set_si (acb_realref b) 6+ arb_set_si (acb_imagref b) 6+ acb_calc_integrate_taylor s f nullPtr a b inr outr goal prec+ acb_add r r s prec+ putStrLn "RESULT:"+ acb_printn r digits 0; putStr "\n"+ -- Bessel integral+ putStrLn $ replicate 64 '-'+ putStrLn "Bessel function J_n(z) = (1/pi) * integral of \+ \cos(t*n - z*sin(t))"+ arf_set_d inr 0.1+ arf_set_d outr 0.5+ let prec' = 3*prec+ acb_set_si a 0+ acb_const_pi b prec'+ f <- makeFunPtr bessel+ acb_calc_integrate_taylor r f nullPtr a b inr outr prec prec'+ acb_div r r b prec+ putStrLn "RESULT:"+ acb_printn r digits 0; putStr "\n"+ + + +data Parameters = Parameters {+ digits :: CLong +} deriving Show++parameters :: Parser Parameters+parameters = Parameters+ <$> argument auto (+ help "compute integrals using d decimal digits of precision."+ <> metavar "d")++--------------------------------------------------------------------------------++foreign import ccall safe "wrapper"+ makeFunPtr :: CAcbCalcFunc -> IO (FunPtr CAcbCalcFunc)++sinx :: Ptr CAcb -> Ptr CAcb -> Ptr () -> CLong -> CLong -> IO CInt+sinx out inp params order prec = do+ let xlen = min 2 order+ acb_set out inp+ when (xlen > 1) $ do acb_one (out `advancePtr` 1)+ _acb_poly_sin_series out out xlen order prec+ return 0++elliptic :: Ptr CAcb -> Ptr CAcb -> Ptr () -> CLong -> CLong -> IO CInt+elliptic out inp params order prec = do+ t <- _acb_vec_init order+ acb_set t inp+ when (order > 1) $ do acb_one (t `advancePtr` 1)+ _acb_poly_sin_series t t (min 2 order) order prec+ _acb_poly_mullow out t order t order order prec+ _acb_vec_scalar_mul_2exp_si t out order (-1)+ acb_sub_ui t t 1 prec+ _acb_vec_neg t t order+ _acb_poly_rsqrt_series out t order order prec+ _acb_vec_clear t order+ return 0++bessel :: Ptr CAcb -> Ptr CAcb -> Ptr () -> CLong -> CLong -> IO CInt+bessel out inp params order prec = do++ t <- _acb_vec_init order++ withNewAcb $ \z -> do+ acb_set t inp+ when (order > 1) $ do acb_one (t `advancePtr` 1)++ let n = 10+ arb_set_si (acb_realref z) 20+ arb_set_si (acb_imagref z) 10++ -- z sin t+ _acb_poly_sin_series out t (min 2 order) order prec+ _acb_vec_scalar_mul out out order z prec++ -- t n+ _acb_vec_scalar_mul_ui t t (min 2 order) n prec++ _acb_poly_sub out t (min 2 order) out order prec++ _acb_poly_cos_series out out order order prec++ _acb_vec_clear t order++ return 0+
+ app/zeta_zeros/Main.hs view
@@ -0,0 +1,160 @@+import System.IO.Unsafe+import Options.Applicative+import Control.Monad+import Control.Monad.State+import Foreign.C.Types+import Foreign.Marshal.Array ++import Data.Number.Flint++main = run =<< execParser opts where+ desc = "Reports the imaginary parts of consecutive nontrivial zeros \+ \of the Riemann zeta function starting with the nth zero."+ opts = info (options <**> helper) (+ fullDesc+ <> progDesc desc+ <> header desc)++run opts@(Options n_start count accuracy platt verbosity num_threads) = do+ when (verbosity > 0) $ do print opts+ if platt && n_start < 10000 then do+ putStrLn "This implementation of the platt algorithm \+ \is not valid\n below the 10000th zero.\n"+ else do+ flint_set_num_threads num_threads+ let (prec, digs) = case accuracy of+ Just (Precision p) -> (p, d2p p)+ Just (Digits d) -> (p2d d, d)+ Nothing -> (pDefault, p2d pDefault)+ where pDefault = 64 + ceiling (logBase 2 (fromIntegral n_start))+ d2p p = round (fromIntegral p * log 2 + 1)+ p2d d = round (fromIntegral d / log 2 + 3)+ requested = fromIntegral $ min count 30000 + usePlatt = platt || (requested > 100 && n_start > 10^11)+ workingPrecision = if platt then 2*prec else prec+ digits = p2d workingPrecision+ p <- _arb_vec_init requested+ let n = fromInteger n_start :: Fmpz+ -- _ <- runStateT (calculate p platt requested prec digits) (n, 0, 0, 0)+ -- return ()+ withFmpz n $ \n -> do+ if not platt then do+ acb_dirichlet_hardy_z_zeros p n requested prec+ print_zeros p n_start requested digits+ else do + found <- acb_dirichlet_platt_local_hardy_z_zeros p n requested prec+ if ( found > 0 ) then do+ print_zeros p n_start found digits+ else do+ putStrLn "Failed to find some zero.\nIncrease precision.\n"+ _arb_vec_clear p $ fromIntegral requested++-- calculate :: Ptr CArb -> Bool -> Integer -> CLong -> CLong+-- -> StateT (Fmpz, Integer, Integer, Integer) IO ()+-- calculate p platt requested prec digits = do+-- (n, iter, count, num_old) <- get+-- let num = if count + num_old > requested then requested-count-1 else 2*num_old+-- withFmpz n $ \n -> do +-- if not platt then do+-- acb_dirichlet_hardy_z_zeros p n num prec+-- print_zeros p n_start num digits+-- fmpz_add_ui n n (fromIntegral num)+-- put (n, succ iter, count + num, num)+-- else do+-- found <- acb_dirichlet_platt_local_hardy_z_zeros p n num prec+-- print_zeros p n_start found digits+-- if ( found > 0 ) then do+-- print_zeros p n_start found digits+-- else do+-- error "Failed to find some zero.\nIncrease precision.\n"+-- fmpz_add_ui n n (fromIntegral found)+-- put (n, succ iter, count + found, num)+-- when (count < requested) $ do calculate p platt requested prec digits++ +data Options = Options {+ n_start :: Integer+, count :: Integer+, accuracy :: Maybe Accuracy+, platt :: Bool+, verbosity :: Int+, num_threads :: CInt+} deriving Show++data Accuracy = Precision CLong | Digits CLong deriving Show++-- option parser ---------------------------------------------------------------++options :: Parser Options+options = Options+ <$> option pos (+ help "integer n > 0. start from nth zero."+ <> short 'n'+ <> value 1+ <> metavar "n")+ <*> option pos (+ help "number of zeros to calculate (<30000)."+ <> long "count"+ <> short 'c'+ <> value 30000+ <> metavar "count")+ <*> optional optionAccuracy+ <*> switch (+ help "use platt algorithm."+ <> showDefault+ <> long "platt")+ <*> option pos (+ help "verbosity."+ <> long "verbosity"+ <> short 'v'+ <> value 0+ <> metavar "verbosity")+ <*> option pos (+ help "number of threads."+ <> long "threads"+ <> short 't'+ <> value 1+ <> metavar "threads")++optionAccuracy = optionPrecision <|> optionDigits++optionPrecision :: Parser Accuracy+optionPrecision = Precision <$> option pos (+ help "precision."+ <> long "prec"+ <> short 'p'+ <> metavar "precision")++optionDigits :: Parser Accuracy+optionDigits = Digits <$> option pos (+ help "number of digits."+ <> long "digits"+ <> short 'd'+ <> metavar "digits")+ +pos :: (Read a, Integral a) => ReadM a+pos = eitherReader $ \s -> do+ let result = read s+ if result >= 0 then + Right result+ else+ Left "expected positive number"++--------------------------------------------------------------------------------++print_zeros p n_start len digits = do+ forM_ [0 .. fromIntegral len - 1] $ \j -> do+ putStr $ show (fromIntegral n_start + j) ++ "\t"+ arb_printn (p `advancePtr` j) digits arb_str_no_radius+ putStr "\n"+-- print_zeros p n len digits = do+-- withNewFmpz $ \k -> do+-- fmpz_set k n+-- forM_ [0 .. len-1] $ \i -> do+-- fmpz_print k+-- putStr "\t"+-- arb_printn (p `advancePtr` i) digits arb_str_no_radius+-- putStr "\n"+-- fmpz_add_ui k k 1++
+ docs/mj.png view
binary file changed (absent → 515566 bytes)
+ docs/out.png view
binary file changed (absent → 248693 bytes)
+ src/Lib.hs view
@@ -0,0 +1,6 @@+module Lib+ ( someFunc+ ) where++someFunc :: IO ()+someFunc = putStrLn "someFunc"
+ test/Spec.hs view
@@ -0,0 +1,2 @@+main :: IO ()+main = putStrLn "Test suite not yet implemented"