diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,4 @@ +# 0.1.0.0 + +Initial version. + diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,202 @@ + + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright [yyyy] [name of copyright owner] + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/numeric-logarithms.cabal b/numeric-logarithms.cabal new file mode 100644 --- /dev/null +++ b/numeric-logarithms.cabal @@ -0,0 +1,57 @@ +cabal-version: 1.12 + +-- This file has been generated from package.yaml by hpack version 0.33.0. +-- +-- see: https://github.com/sol/hpack +-- +-- hash: 055716fd894ce8e71c2be1db7a3f5a30b6ef588291635e8f57429d608c19eb71 + +name: numeric-logarithms +version: 0.1.0.0 +synopsis: Efficient rounded log2 algorithms on integral and rational types. +description: This provides a user-facing API for the GHC primitive 'integerLog2#', in the + form of functions computing the floor, ceiling, or other rounded form of the + base-2 logarithm of 'Integral' or 'Rational' types. +category: Numeric +homepage: https://github.com/google/hs-numeric-logarithms#readme +bug-reports: https://github.com/google/hs-numeric-logarithms/issues +author: Andrew Pritchard +maintainer: Andrew Pritchard +copyright: 2018-2021 Google LLC +license: Apache-2.0 +license-file: LICENSE +build-type: Simple +extra-source-files: + CHANGELOG.md + +source-repository head + type: git + location: https://github.com/google/hs-numeric-logarithms + +library + exposed-modules: + Numeric.Logarithms + other-modules: + Paths_numeric_logarithms + hs-source-dirs: + src + build-depends: + base >=4.9 && <4.16 + , integer-gmp >=0.5 && <1.1 + default-language: Haskell2010 + +test-suite Logarithms-test + type: exitcode-stdio-1.0 + main-is: LogarithmsTest.hs + other-modules: + Paths_numeric_logarithms + hs-source-dirs: + test + build-depends: + QuickCheck + , base >=4.9 && <4.16 + , integer-gmp >=0.5 && <1.1 + , numeric-logarithms + , test-framework + , test-framework-quickcheck2 + default-language: Haskell2010 diff --git a/src/Numeric/Logarithms.hs b/src/Numeric/Logarithms.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Logarithms.hs @@ -0,0 +1,170 @@ +-- Copyright 2021 Google LLC +-- +-- Licensed under the Apache License, Version 2.0 (the "License"); +-- you may not use this file except in compliance with the License. +-- You may obtain a copy of the License at +-- +-- https://www.apache.org/licenses/LICENSE-2.0 +-- +-- Unless required by applicable law or agreed to in writing, software +-- distributed under the License is distributed on an "AS IS" BASIS, +-- WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +-- See the License for the specific language governing permissions and +-- limitations under the License. + +-- | Rounded base-2 logarithms of 'Integral' and 'Real' types. + +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE LambdaCase #-} +{-# LANGUAGE MagicHash #-} + +module Numeric.Logarithms + ( -- * Real Logarithms + log2Floor, log2Ceiling + -- ** Advanced + , log2Approx, log2With + -- * Integer Logarithms + , ilog2Floor, ilog2Ceiling + -- ** Advanced + , ilog2Approx, ilog2With + ) where + +import Data.Bits (bit) +import Data.Ratio (denominator, numerator) +import GHC.Exts (Int(I#)) +import GHC.Integer (shiftLInteger) +import GHC.Integer.Logarithms (integerLog2#) +import GHC.Stack (HasCallStack) + +integerLog2Floor :: Integer -> Int +integerLog2Floor x = I# (integerLog2# x) + +-- Derivation of log2 algorithms for Rational: +-- +-- floor(log2(r)) +-- = floor(log2(x/y)) with x, y the numerator and denominator of r. +-- = floor(log2((2^n * (1 + k/2^n)) / (2^m * (1 + j/2^m)))) +-- with x = 2^n + k, k < 2^n, n an integer +-- y = 2^m + j, j < 2^m, m an integer +-- = floor(log2(2^n/2^m) + log2((1+k/2^n)/(1+j/2^m))) +-- = floor(n - m + log2((1+k/2^n)/(1+j/2^m))) +-- = n - m + floor(log2((1+k/2^n)/(1+j/2^m))) +-- = n - m + if (1+k/2^n) < (1+j/2^m) then -1 else 0 +-- because 1+k/2^n and 1+j/2^m are both in [1, 2), +-- so their quotient is in (1/2, 2) +-- so log2 of their quotient is in (-1, 1) +-- so the floor of log2 of their quotient is -1 or 0: +-- -1 if the quotient is less than 1; 0 if the quotient is >= 1 +-- = n - m + if k/2^n < j/2^m then -1 else 0 +-- = n - m + if k*2^m < j*2^n then -1 else 0 +-- = n - m + if (k< Integer -> (Int, Integer, Integer) +splitLog2Unchecked num den = (n - m, num', den') + where + n = integerLog2Floor num + m = integerLog2Floor den + -- Note: the Bits instance for Integer doesn't define unsafeShiftL or shiftL, + -- so it defaults to testing the sign of the argument and dispatching to left + -- or right shifts. So, we call what its definition should have been instead. + shl x (I# s) = shiftLInteger x s + (num', den') = if m > n + then (num `shl` (m-n), den) + else (num, den `shl` (n-m)) + +{-# INLINE log2Unchecked #-} +log2Unchecked :: Integer -> Integer -> (Int, Ordering) +log2Unchecked num den = + let (lg, num', den') = splitLog2Unchecked num den + in (lg, compare num' den') + +-- The underscore-suffixed variants here are inlined versions so that GHC will +-- implement the exported functions by big chunks of optimized code, but not +-- try to include that code in the module interface. + +{-# INLINE log2Approx_ #-} +log2Approx_ :: (HasCallStack, Real a) => a -> (Int, Ordering) +log2Approx_ x = + let !xr = toRational x + !num = numerator xr + !den = denominator xr + in if num <= 0 + then error "log2With_: x <= 0" + else log2Unchecked num den + +{-# INLINE log2With_ #-} +log2With_ :: (HasCallStack, Real a) => (Int -> Ordering -> Int) -> a -> Int +log2With_ adj x = + let !(lg, cmp) = log2Approx_ x + in adj lg cmp + +-- | Returns an approximate base-2 logarithm of the argument. +-- +-- The returned @Int@ is one of the two nearest integers to the exact result, +-- and the returned @Ordering@ tells whether the exact result is greater than, +-- equal to, or less than the @Int@. `LT` means the exact result is less than +-- the rounded result. +-- +-- This effectively gives results like "between 7 and 8" and leaves it up to +-- the caller to decide how (or whether) to round the result. See also +-- `log2With` for a version that will always round, but leaves the particular +-- rounding strategy up to the caller. +log2Approx :: (HasCallStack, Real a) => a -> (Int, Ordering) +log2Approx = log2Approx_ + +-- | Returns the base-2 logarithm with custom rounding. +-- +-- The first parameter is a rounding adjustment function: given the approximate +-- result and an 'Ordering' indicating its relation to the exact result, return +-- the rounded result. +-- +-- This could be useful if you want something other than floor or ceil, e.g. +-- round-towards-0, round-towards-even, etc. +log2With :: (HasCallStack, Real a) => (Int -> Ordering -> Int) -> a -> Int +log2With = log2With_ + +-- | Returns the floor of the base-2 logarithm of a 'Real' argument. +log2Floor :: (HasCallStack, Real a) => a -> Int +log2Floor = log2With_ (\x -> \case LT -> x - 1; _ -> x) + +-- | Returns the ceiling of the base-2 logarithm of a 'Real' argument. +log2Ceiling :: (HasCallStack, Real a) => a -> Int +log2Ceiling = log2With_ (\x -> \case GT -> x + 1; _ -> x) + +{-# INLINE withPositiveInteger #-} +withPositiveInteger :: HasCallStack => (Integer -> r) -> Integer -> r +withPositiveInteger f x = if x <= 0 + then error "withPositiveInteger: x <= 0" + else f x + +-- | Returns the floor of the base-2 logarithm of an 'Integral' argument. +{-# INLINABLE ilog2Floor #-} +ilog2Floor :: (HasCallStack, Integral a) => a -> Int +ilog2Floor x = withPositiveInteger integerLog2Floor (toInteger x) + +-- | Returns the ceiling of the base-2 logarithm of an 'Integral' argument. +{-# INLINABLE ilog2Ceiling #-} +ilog2Ceiling :: (HasCallStack, Integral a) => a -> Int +ilog2Ceiling = ilog2With (\x -> \case GT -> x + 1; _ -> x) + +-- | Returns the approximate base-2 logarithm of an 'Integral' argument. +{-# INLINABLE ilog2Approx #-} +ilog2Approx :: (HasCallStack, Integral a) => a -> (Int, Ordering) +ilog2Approx x = withPositiveInteger + (\xi -> + let lg = integerLog2Floor xi + in (lg, if xi > bit lg then GT else EQ)) + (toInteger x) + +-- | Returns the base-2 logarithm of an 'Integral' with custom rounding. +{-# INLINABLE ilog2With #-} +ilog2With :: (HasCallStack, Integral a) => (Int -> Ordering -> Int) -> a -> Int +ilog2With adj x = + let !(lg, cmp) = ilog2Approx x + in adj lg cmp diff --git a/test/LogarithmsTest.hs b/test/LogarithmsTest.hs new file mode 100644 --- /dev/null +++ b/test/LogarithmsTest.hs @@ -0,0 +1,108 @@ +-- Copyright 2021 Google LLC +-- +-- Licensed under the Apache License, Version 2.0 (the "License"); +-- you may not use this file except in compliance with the License. +-- You may obtain a copy of the License at +-- +-- http://www.apache.org/licenses/LICENSE-2.0 +-- +-- Unless required by applicable law or agreed to in writing, software +-- distributed under the License is distributed on an "AS IS" BASIS, +-- WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +-- See the License for the specific language governing permissions and +-- limitations under the License. + +{-# LANGUAGE AllowAmbiguousTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeApplications #-} + +module Main(main) where + +import Test.Framework.Providers.QuickCheck2 (testProperty) +import Test.QuickCheck + ( (===), (.&&.), Arbitrary(..), Property + , forAll, suchThat, withMaxSuccess + ) + +import Test.Framework(defaultMain) + +import Numeric.Logarithms + ( log2Floor, log2Ceiling, log2Approx + , ilog2Floor, ilog2Ceiling + ) + +exp2 :: Int -> Rational +exp2 x | x < 0 = recip (exp2 (-x)) +exp2 x = 2 ^ x + +prop_log2FloorCorrect :: forall a. (Arbitrary a, Real a, Show a) => Property +prop_log2FloorCorrect = + withMaxSuccess 10000 $ + forAll (arbitrary @a `suchThat` (>0)) $ \x -> + let xr = toRational x + lg = log2Floor x + lo = exp2 lg + hi = exp2 (lg+1) + in (lo <= xr) === True .&&. (xr < hi) === True + +prop_log2CeilingCorrect :: forall a. (Arbitrary a, Real a, Show a) => Property +prop_log2CeilingCorrect = + withMaxSuccess 10000 $ + forAll (arbitrary @a `suchThat` (>0)) $ \x -> + let xr = toRational x + lg = log2Ceiling x + lo = exp2 (lg-1) + hi = exp2 lg + in (lo < xr) === True .&&. (xr <= hi) === True + +prop_log2ApproxCorrect :: forall a. (Arbitrary a, Real a, Show a) => Property +prop_log2ApproxCorrect = + withMaxSuccess 10000 $ + forAll (arbitrary @a `suchThat` (>0)) $ \x -> + let xr = toRational x + (lg, o) = log2Approx x + lt = exp2 (lg-1) + eq = exp2 lg + gt = exp2 (lg+1) + in case o of + LT -> (lt < xr) === True .&&. (xr < eq) === True + EQ -> xr === eq + GT -> (eq < xr) === True .&&. (xr < gt) === True + +prop_ilog2FloorAgrees :: forall a. (Arbitrary a, Integral a, Show a) => Property +prop_ilog2FloorAgrees = + withMaxSuccess 10000 $ + forAll (arbitrary @a `suchThat` (>0)) $ \x -> + ilog2Floor x === log2Floor x + +prop_ilog2CeilingAgrees :: forall a. (Arbitrary a, Integral a, Show a) => Property +prop_ilog2CeilingAgrees = + withMaxSuccess 10000 $ + forAll (arbitrary @a `suchThat` (>0)) $ \x -> + ilog2Ceiling x === log2Ceiling x + +main :: IO () +main = defaultMain + [ testProperty "log2Floor @Int" (prop_log2FloorCorrect @Int) + , testProperty "log2Floor @Integer" (prop_log2FloorCorrect @Integer) + , testProperty "log2Floor @Double" (prop_log2FloorCorrect @Double) + , testProperty "log2Floor @Rational" (prop_log2FloorCorrect @Rational) + + , testProperty "log2Ceiling @Int" (prop_log2CeilingCorrect @Int) + , testProperty "log2Ceiling @Integer" (prop_log2CeilingCorrect @Integer) + , testProperty "log2Ceiling @Double" (prop_log2CeilingCorrect @Double) + , testProperty "log2Ceiling @Rational" (prop_log2CeilingCorrect @Rational) + + , testProperty "log2Approx @Int" (prop_log2ApproxCorrect @Int) + , testProperty "log2Approx @Integer" (prop_log2ApproxCorrect @Integer) + , testProperty "log2Approx @Double" (prop_log2ApproxCorrect @Double) + , testProperty "log2Approx @Rational" (prop_log2ApproxCorrect @Rational) + + , testProperty "ilog2Floor @Int" (prop_ilog2FloorAgrees @Int) + , testProperty "ilog2Floor @Integer" (prop_ilog2FloorAgrees @Integer) + , testProperty "ilog2Floor @Word" (prop_ilog2FloorAgrees @Word) + + , testProperty "ilog2Ceiling @Int" (prop_ilog2CeilingAgrees @Int) + , testProperty "ilog2Ceiling @Integer" (prop_ilog2CeilingAgrees @Integer) + , testProperty "ilog2Ceiling @Word" (prop_ilog2CeilingAgrees @Word) + ]