{-
	Copyright (C) 2011 Dr. Alistair Ward

	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 3 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, see <http://www.gnu.org/licenses/>.
-}
{- |
 [@AUTHOR@]	Dr. Alistair Ward

 [@DESCRIPTION@]

	* Implements a /spigot/-algorithm; <https://en.wikipedia.org/wiki/Spigot_algorithm>.

	* Uses the traditional algorithm, rather than the /unbounded/ algorithm described by <http://www.comlab.ox.ac.uk/jeremy.gibbons/publications/spigot.pdf>.
-}

module Factory.Math.Implementations.Pi.Spigot.Spigot(
-- * Types
-- ** Type-synonyms
--	Base,
--	Coefficients,
--	I,
--	Pi,
--	PreDigits,
--	QuotRem,
-- * Constants
	decimal,
-- * Functions
--	carryAndDivide,
--	processColumns,
	openI,
-- ** Accessors
--	getQuotient,
--	getRemainder,
-- ** Constructor
--	mkRow
) where

import qualified	Control.Arrow
import qualified	Data.Char
import qualified	Data.Ratio
import qualified	Factory.Math.Implementations.Pi.Spigot.Series	as Math.Implementations.Pi.Spigot.Series
import qualified	Factory.Math.Precision				as Math.Precision

{- |
	* The type in which all arithmetic is performed.

	* A small dynamic range, 32 bits or more, is typically adequate.
-}
type I	= Int

-- | The constant base in which we want the resulting value of /Pi/ to be expressed.
decimal :: I
decimal :: I
decimal	= I
10

-- | Coerce the polymorphic type 'Data.Ratio.Ratio' to suit the base used in our series.
type Base	= Data.Ratio.Ratio I

-- | Coerce the polymorphic type returned by 'quotRem' to our specific requirements.
type QuotRem	= (I, I)

-- Accessors.
getQuotient, getRemainder :: QuotRem -> I
getQuotient :: QuotRem -> I
getQuotient	= QuotRem -> I
forall a b. (a, b) -> a
fst
getRemainder :: QuotRem -> I
getRemainder	= QuotRem -> I
forall a b. (a, b) -> b
snd

type PreDigits		= [I]
type Pi			= [I]
type Coefficients	= [I]

{- |
	* For a digit on one row of the spigot-table, add any numerator carried from the similar calculation one column to the right.

	* Divide the result of this summation, by the denominator of the base, to get the quotient and remainder.

	* Determine the quantity to carry to the similar calculation one column to the left, by multiplying the quotient by the numerator of the base.
-}
carryAndDivide :: (Base, I) -> QuotRem -> QuotRem
carryAndDivide :: (Base, I) -> QuotRem -> QuotRem
carryAndDivide (Base
base, I
lhs) QuotRem
rhs
	| I
n I -> I -> Bool
forall a. Ord a => a -> a -> Bool
< I
d		= (I
0, I
n)	-- In some degenerate cases, the result of the subsequent calculation can be more simply determined.
	| Bool
otherwise	= (I -> I) -> QuotRem -> QuotRem
forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (b, d) (c, d)
Control.Arrow.first (I -> I -> I
forall a. Num a => a -> a -> a
* Base -> I
forall a. Ratio a -> a
Data.Ratio.numerator Base
base) (QuotRem -> QuotRem) -> QuotRem -> QuotRem
forall a b. (a -> b) -> a -> b
$ I
n I -> I -> QuotRem
forall a. Integral a => a -> a -> (a, a)
`quotRem` I
d
	where
		d, n :: I
		d :: I
d	= Base -> I
forall a. Ratio a -> a
Data.Ratio.denominator Base
base
		n :: I
n	= I
lhs I -> I -> I
forall a. Num a => a -> a -> a
+ QuotRem -> I
getQuotient QuotRem
rhs	-- Carry numerator from the column to the right and add it to the current digit.

{- |
	* Fold 'carryAndDivide', from right to left, over the columns of a row in the spigot-table, continuously checking for overflow.

	* Release any previously withheld result-digits, after any adjustment for overflow in the current result-digit.

	* Withhold the current result-digit until the risk of overflow in subsequent result-digits has been assessed.

	* Call 'mkRow'.
-}
processColumns
	:: Math.Implementations.Pi.Spigot.Series.Series I
	-> PreDigits
	-> [(Base, I)]	-- ^ Data-row.
	-> Pi
processColumns :: Series I -> PreDigits -> [(Base, I)] -> PreDigits
processColumns Series I
series PreDigits
preDigits [(Base, I)]
l
	| I
overflowMargin I -> I -> Bool
forall a. Ord a => a -> a -> Bool
> I
1	= PreDigits
preDigits PreDigits -> PreDigits -> PreDigits
forall a. [a] -> [a] -> [a]
++ PreDigits -> PreDigits
nextRow [I
digit]				-- There's neither overflow, nor risk of impact from subsequent overflow.
	| I
overflowMargin I -> I -> Bool
forall a. Eq a => a -> a -> Bool
== I
1	= PreDigits -> PreDigits
nextRow (PreDigits -> PreDigits) -> PreDigits -> PreDigits
forall a b. (a -> b) -> a -> b
$ PreDigits
preDigits PreDigits -> PreDigits -> PreDigits
forall a. [a] -> [a] -> [a]
++ [I
digit]			-- There's no overflow, but risk of impact from subsequent overflow.
	| Bool
otherwise		= (I -> I) -> PreDigits -> PreDigits
forall a b. (a -> b) -> [a] -> [b]
map ((I -> I -> I
forall a. Integral a => a -> a -> a
`rem` I
decimal) (I -> I) -> (I -> I) -> I -> I
forall b c a. (b -> c) -> (a -> b) -> a -> c
. I -> I
forall a. Enum a => a -> a
succ) PreDigits
preDigits PreDigits -> PreDigits -> PreDigits
forall a. [a] -> [a] -> [a]
++ PreDigits -> PreDigits
nextRow [I
0]	-- Overflow => propagate the excess to previously withheld preDigits.
	where
		results :: [QuotRem]
		results :: [QuotRem]
results	= [QuotRem] -> [QuotRem]
forall a. [a] -> [a]
init ([QuotRem] -> [QuotRem]) -> [QuotRem] -> [QuotRem]
forall a b. (a -> b) -> a -> b
$ ((Base, I) -> QuotRem -> QuotRem)
-> QuotRem -> [(Base, I)] -> [QuotRem]
forall a b. (a -> b -> b) -> b -> [a] -> [b]
scanr (Base, I) -> QuotRem -> QuotRem
carryAndDivide (I
0, I
forall a. HasCallStack => a
undefined) [(Base, I)]
l

		digit :: I
		digit :: I
digit	= QuotRem -> I
getQuotient (QuotRem -> I) -> QuotRem -> I
forall a b. (a -> b) -> a -> b
$ [QuotRem] -> QuotRem
forall a. [a] -> a
head [QuotRem]
results

		overflowMargin :: I
		overflowMargin :: I
overflowMargin	= I
decimal I -> I -> I
forall a. Num a => a -> a -> a
- I
digit

		nextRow :: [I] -> [I]
		nextRow :: PreDigits -> PreDigits
nextRow PreDigits
preDigits'	= Series I -> PreDigits -> PreDigits -> PreDigits
mkRow Series I
series PreDigits
preDigits' (PreDigits -> PreDigits) -> PreDigits -> PreDigits
forall a b. (a -> b) -> a -> b
$ (QuotRem -> I) -> [QuotRem] -> PreDigits
forall a b. (a -> b) -> [a] -> [b]
map QuotRem -> I
getRemainder [QuotRem]
results

{- |
	* Multiply the remainders from the previous row.

	* Zip them with the constant bases, with an addition one stuck on the front to perform the conversion to decimal, to create a new row of the spigot-table.

	* Call 'processColumns'.
-}
mkRow :: Math.Implementations.Pi.Spigot.Series.Series I -> PreDigits -> Coefficients -> Pi
mkRow :: Series I -> PreDigits -> PreDigits -> PreDigits
mkRow Series I
series PreDigits
preDigits	= Series I -> PreDigits -> [(Base, I)] -> PreDigits
processColumns Series I
series PreDigits
preDigits ([(Base, I)] -> PreDigits)
-> (PreDigits -> [(Base, I)]) -> PreDigits -> PreDigits
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Base] -> PreDigits -> [(Base, I)]
forall a b. [a] -> [b] -> [(a, b)]
zip (Base -> Base
forall a. Fractional a => a -> a
recip (I -> Base
forall a b. (Integral a, Num b) => a -> b
fromIntegral I
decimal) Base -> [Base] -> [Base]
forall a. a -> [a] -> [a]
: Series I -> [Base]
forall i. Integral i => Series i -> [Ratio i]
Math.Implementations.Pi.Spigot.Series.bases Series I
series) (PreDigits -> [(Base, I)])
-> (PreDigits -> PreDigits) -> PreDigits -> [(Base, I)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (I -> I) -> PreDigits -> PreDigits
forall a b. (a -> b) -> [a] -> [b]
map (I -> I -> I
forall a. Num a => a -> a -> a
* I
decimal)

{- |
	* Initialises a /spigot/-table with the row of 'Math.Implementations.Pi.Spigot.Series.coefficients'.

	* Ensures that the row has suffient terms to accurately generate the required number of digits.

	* Extracts only those digits which are guaranteed to be accurate.

	* CAVEAT: the result is returned as an 'Integer', i.e. without any decimal point.
-}
openI :: Math.Implementations.Pi.Spigot.Series.Series I -> Math.Precision.DecimalDigits -> Integer
openI :: Series I -> I -> Integer
openI Series I
series I
decimalDigits	= String -> Integer
forall a. Read a => String -> a
read (String -> Integer)
-> (PreDigits -> String) -> PreDigits -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (I -> Char) -> PreDigits -> String
forall a b. (a -> b) -> [a] -> [b]
map (
	I -> Char
Data.Char.intToDigit (I -> Char) -> (I -> I) -> I -> Char
forall b c a. (b -> c) -> (a -> b) -> a -> c
. I -> I
forall a b. (Integral a, Num b) => a -> b
fromIntegral
 ) (PreDigits -> String)
-> (PreDigits -> PreDigits) -> PreDigits -> String
forall b c a. (b -> c) -> (a -> b) -> a -> c
. I -> PreDigits -> PreDigits
forall a. I -> [a] -> [a]
take I
decimalDigits (PreDigits -> PreDigits)
-> (PreDigits -> PreDigits) -> PreDigits -> PreDigits
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Series I -> PreDigits -> PreDigits -> PreDigits
mkRow Series I
series [] (PreDigits -> PreDigits)
-> (PreDigits -> PreDigits) -> PreDigits -> PreDigits
forall b c a. (b -> c) -> (a -> b) -> a -> c
. I -> PreDigits -> PreDigits
forall a. I -> [a] -> [a]
take (
	Series I -> I -> I
forall i. Series i -> I -> I
Math.Implementations.Pi.Spigot.Series.nTerms Series I
series I
decimalDigits
 ) (PreDigits -> Integer) -> PreDigits -> Integer
forall a b. (a -> b) -> a -> b
$ Series I -> PreDigits
forall i. Series i -> [i]
Math.Implementations.Pi.Spigot.Series.coefficients Series I
series