{- 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; <http://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, -- ** Constructors -- 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 = 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 = fst getRemainder = 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, lhs) rhs | n < d = (0, n) -- In some degenerate cases, the result of the subsequent calculation can be more simply determined. | otherwise = Control.Arrow.first (* Data.Ratio.numerator base) $ n `quotRem` d where d, n :: I d = Data.Ratio.denominator base n = lhs + getQuotient 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 preDigits l | overflowMargin > 1 = preDigits ++ nextRow [digit] -- There's neither overflow, nor risk of impact from subsequent overflow. | overflowMargin == 1 = nextRow $ preDigits ++ [digit] -- There's no overflow, but risk of impact from subsequent overflow. | otherwise = map ((`rem` decimal) . succ) preDigits ++ nextRow [0] -- Overflow => propagate the excess to previously withheld preDigits. where results :: [QuotRem] results = init $ scanr carryAndDivide (0, undefined) l digit :: I digit = getQuotient $ head results overflowMargin :: I overflowMargin = decimal - digit nextRow :: [I] -> [I] nextRow preDigits' = mkRow series preDigits' $ map getRemainder 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 preDigits = processColumns series preDigits . zip (recip (fromIntegral decimal) : Math.Implementations.Pi.Spigot.Series.bases series) . map (* 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 decimalDigits = read . map ( Data.Char.intToDigit . fromIntegral ) . take decimalDigits . mkRow series [] . take ( Math.Implementations.Pi.Spigot.Series.nTerms series decimalDigits ) $ Math.Implementations.Pi.Spigot.Series.coefficients series