diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -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 diff --git a/Flint2-Examples.cabal b/Flint2-Examples.cabal new file mode 100644 --- /dev/null +++ b/Flint2-Examples.cabal @@ -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 diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. + + + Copyright (C) + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 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. + + , 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. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,171 @@ +![examples of complex_plot](https://github.com/monien/Flint2-Examples/raw/main/docs/out.png) + +# 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. + +

+ +

+ +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 diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/app/bernoulli/Main.hs b/app/bernoulli/Main.hs new file mode 100644 --- /dev/null +++ b/app/bernoulli/Main.hs @@ -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") diff --git a/app/binet/Main.hs b/app/binet/Main.hs new file mode 100644 --- /dev/null +++ b/app/binet/Main.hs @@ -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") + diff --git a/app/class_poly/Main.hs b/app/class_poly/Main.hs new file mode 100644 --- /dev/null +++ b/app/class_poly/Main.hs @@ -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." + diff --git a/app/complex_plot/ColorFunction.hs b/app/complex_plot/ColorFunction.hs new file mode 100644 --- /dev/null +++ b/app/complex_plot/ColorFunction.hs @@ -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 diff --git a/app/complex_plot/Functions.hs b/app/complex_plot/Functions.hs new file mode 100644 --- /dev/null +++ b/app/complex_plot/Functions.hs @@ -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 + + + + diff --git a/app/complex_plot/Main.hs b/app/complex_plot/Main.hs new file mode 100644 --- /dev/null +++ b/app/complex_plot/Main.hs @@ -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") diff --git a/app/crt/Main.hs b/app/crt/Main.hs new file mode 100644 --- /dev/null +++ b/app/crt/Main.hs @@ -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 ++ " ." + _ -> putStrLn $ "usage: " ++ prog ++ " ." + +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 + + diff --git a/app/delta_qexp/Main.hs b/app/delta_qexp/Main.hs new file mode 100644 --- /dev/null +++ b/app/delta_qexp/Main.hs @@ -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" diff --git a/app/dft/Main.hs b/app/dft/Main.hs new file mode 100644 --- /dev/null +++ b/app/dft/Main.hs @@ -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 + diff --git a/app/elementary/Main.hs b/app/elementary/Main.hs new file mode 100644 --- /dev/null +++ b/app/elementary/Main.hs @@ -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 + + diff --git a/app/expression/Main.hs b/app/expression/Main.hs new file mode 100644 --- /dev/null +++ b/app/expression/Main.hs @@ -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 + diff --git a/app/factor_integer/Main.hs b/app/factor_integer/Main.hs new file mode 100644 --- /dev/null +++ b/app/factor_integer/Main.hs @@ -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") diff --git a/app/fmpq_poly/Main.hs b/app/fmpq_poly/Main.hs new file mode 100644 --- /dev/null +++ b/app/fmpq_poly/Main.hs @@ -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 + diff --git a/app/fmpz_mod_poly/Main.hs b/app/fmpz_mod_poly/Main.hs new file mode 100644 --- /dev/null +++ b/app/fmpz_mod_poly/Main.hs @@ -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" diff --git a/app/fmpz_mpoly_factor/Main.hs b/app/fmpz_mpoly_factor/Main.hs new file mode 100644 --- /dev/null +++ b/app/fmpz_mpoly_factor/Main.hs @@ -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") diff --git a/app/fmpz_poly_factor_zassenhaus/Main.hs b/app/fmpz_poly_factor_zassenhaus/Main.hs new file mode 100644 --- /dev/null +++ b/app/fmpz_poly_factor_zassenhaus/Main.hs @@ -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" diff --git a/app/fpwrap/Main.hs b/app/fpwrap/Main.hs new file mode 100644 --- /dev/null +++ b/app/fpwrap/Main.hs @@ -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) diff --git a/app/fq_poly/Main.hs b/app/fq_poly/Main.hs new file mode 100644 --- /dev/null +++ b/app/fq_poly/Main.hs @@ -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 + + + + + diff --git a/app/function_benchmark/Main.hs b/app/function_benchmark/Main.hs new file mode 100644 --- /dev/null +++ b/app/function_benchmark/Main.hs @@ -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 + diff --git a/app/hilbert_matrix/Main.hs b/app/hilbert_matrix/Main.hs new file mode 100644 --- /dev/null +++ b/app/hilbert_matrix/Main.hs @@ -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 () + diff --git a/app/hilbert_matrix_ca/Main.hs b/app/hilbert_matrix_ca/Main.hs new file mode 100644 --- /dev/null +++ b/app/hilbert_matrix_ca/Main.hs @@ -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") diff --git a/app/integrals/Integrands.hs b/app/integrals/Integrands.hs new file mode 100644 --- /dev/null +++ b/app/integrals/Integrands.hs @@ -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 () diff --git a/app/integrals/Main.hs b/app/integrals/Main.hs new file mode 100644 --- /dev/null +++ b/app/integrals/Main.hs @@ -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 + diff --git a/app/keiper_li/Main.hs b/app/keiper_li/Main.hs new file mode 100644 --- /dev/null +++ b/app/keiper_li/Main.hs @@ -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 + + + diff --git a/app/l_values/Main.hs b/app/l_values/Main.hs new file mode 100644 --- /dev/null +++ b/app/l_values/Main.hs @@ -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") + + diff --git a/app/logistic/Main.hs b/app/logistic/Main.hs new file mode 100644 --- /dev/null +++ b/app/logistic/Main.hs @@ -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" + + diff --git a/app/machin/Main.hs b/app/machin/Main.hs new file mode 100644 --- /dev/null +++ b/app/machin/Main.hs @@ -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)]) + ] diff --git a/app/multi_crt/Main.hs b/app/multi_crt/Main.hs new file mode 100644 --- /dev/null +++ b/app/multi_crt/Main.hs @@ -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") diff --git a/app/padic/Main.hs b/app/padic/Main.hs new file mode 100644 --- /dev/null +++ b/app/padic/Main.hs @@ -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" + + + + diff --git a/app/partitions/Main.hs b/app/partitions/Main.hs new file mode 100644 --- /dev/null +++ b/app/partitions/Main.hs @@ -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 ++ " " + _ -> putStrLn $ "usage: " ++ prog ++ " " + +run n = do + result <- newFmpz + withFmpz result $ \p -> arith_number_of_partitions p n + putStrLn $ "p(" ++ show n ++ ") = " ++ show result + diff --git a/app/pi_digits/Main.hs b/app/pi_digits/Main.hs new file mode 100644 --- /dev/null +++ b/app/pi_digits/Main.hs @@ -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") + diff --git a/app/poly_roots/Main.hs b/app/poly_roots/Main.hs new file mode 100644 --- /dev/null +++ b/app/poly_roots/Main.hs @@ -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 Easy polynomial 1 + 2x + ... + (n+1)x^n" + , "t Chebyshev polynomial T_n" + , "u Chebyshev polynomial U_n" + , "p Legendre polynomial P_n" + , "c Cyclotomic polynomial Phi_n" + , "s Swinnerton-Dyer polynomial S_n" + , "b Bernoulli polynomial B_n" + , "w Wilkinson polynomial W_n" + , "e Taylor series of exp(x) truncated to degree n" + , "m The Mignotte-like polynomial x^n + (100x+1)^m, n > m" + , "coeffs 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 diff --git a/app/primegen/Main.hs b/app/primegen/Main.hs new file mode 100644 --- /dev/null +++ b/app/primegen/Main.hs @@ -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") diff --git a/app/qadic/Main.hs b/app/qadic/Main.hs new file mode 100644 --- /dev/null +++ b/app/qadic/Main.hs @@ -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" + + + + diff --git a/app/radix/Main.hs b/app/radix/Main.hs new file mode 100644 --- /dev/null +++ b/app/radix/Main.hs @@ -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 () + + + + + + + diff --git a/app/real_roots/FFI.hsc b/app/real_roots/FFI.hsc new file mode 100644 --- /dev/null +++ b/app/real_roots/FFI.hsc @@ -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 + diff --git a/app/real_roots/Functions.hs b/app/real_roots/Functions.hs new file mode 100644 --- /dev/null +++ b/app/real_roots/Functions.hs @@ -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) diff --git a/app/real_roots/Main.hs b/app/real_roots/Main.hs new file mode 100644 --- /dev/null +++ b/app/real_roots/Main.hs @@ -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 () diff --git a/app/real_roots/Run.hs b/app/real_roots/Run.hs new file mode 100644 --- /dev/null +++ b/app/real_roots/Run.hs @@ -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" diff --git a/app/real_roots/z_param.h b/app/real_roots/z_param.h new file mode 100644 --- /dev/null +++ b/app/real_roots/z_param.h @@ -0,0 +1,16 @@ +#ifndef Z_PARAM_H_ +#define Z_PARAM_H_ + +#include +#include + +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_ diff --git a/app/stirling_matrix/Main.hs b/app/stirling_matrix/Main.hs new file mode 100644 --- /dev/null +++ b/app/stirling_matrix/Main.hs @@ -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 ++ " " + _ -> putStrLn $ "usage: " ++ prog ++ " " + +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 + + diff --git a/app/swinnerton_dyer_poly/Main.hs b/app/swinnerton_dyer_poly/Main.hs new file mode 100644 --- /dev/null +++ b/app/swinnerton_dyer_poly/Main.hs @@ -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" + diff --git a/app/taylor_integrals/Main.hs b/app/taylor_integrals/Main.hs new file mode 100644 --- /dev/null +++ b/app/taylor_integrals/Main.hs @@ -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 + diff --git a/app/zeta_zeros/Main.hs b/app/zeta_zeros/Main.hs new file mode 100644 --- /dev/null +++ b/app/zeta_zeros/Main.hs @@ -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 + + diff --git a/docs/mj.png b/docs/mj.png new file mode 100644 Binary files /dev/null and b/docs/mj.png differ diff --git a/docs/out.png b/docs/out.png new file mode 100644 Binary files /dev/null and b/docs/out.png differ diff --git a/src/Lib.hs b/src/Lib.hs new file mode 100644 --- /dev/null +++ b/src/Lib.hs @@ -0,0 +1,6 @@ +module Lib + ( someFunc + ) where + +someFunc :: IO () +someFunc = putStrLn "someFunc" diff --git a/test/Spec.hs b/test/Spec.hs new file mode 100644 --- /dev/null +++ b/test/Spec.hs @@ -0,0 +1,2 @@ +main :: IO () +main = putStrLn "Test suite not yet implemented"