{-

Calculate reciprocal square roots in arbitrary bases to arbitrary
precision.

Copyright 2016 Ken Takusagawa

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/>.
-}

{-# LANGUAGE ScopedTypeVariables, LambdaCase #-}
module Main where {
import Control.Category((>>>));
import Prelude hiding((>>));
import Data.List;
import Data.Number.Fixed(Fixed,Epsilon,dynamicEps);
import Data.Ratio((%));
import Data.Function((&));
import System.Environment(getArgs);

main :: IO();
main = getArgs >>= \case{
[base, x, numdigits] -> go (read base) (read x) (read numdigits);
_ -> error "$0 base x numdigits\ncomputes sqrt(1/x) to base numdigits"
};

safety_margin :: Integer;
safety_margin = 10;

go :: Integer -> Integer -> Integer -> IO();
go base x numdigits = do {
let { eps1 :: Rational;
 eps1 = (1%base) ^ (numdigits + safety_margin);
};
 recip_sqrt eps1 base (fromInteger x) & genericTake numdigits & map show & unwords & putStrLn;
};

-- converts a number x (0 <= x < 1) to base.
radix_convert_fractional :: forall e . Epsilon e => Integer -> Fixed e -> Maybe (Integer, Fixed e);
radix_convert_fractional base x = let
{ y :: Fixed e
; y = x * fromInteger base
; ret :: Integer
; ret = floor y} in if x == 0 then Nothing else Just (ret, y - fromInteger ret);

-- consider dropping leading zeroes, but this makes the calculation of the proper eps tricky

recip_sqrt :: Rational -> Integer -> Rational -> [Integer];
recip_sqrt eps base = dynamicEps eps (recip >>> sqrt >>> unfoldr (radix_convert_fractional base));
} --end