{-
Fractional exponentiation without exp and log.

Usage: $0 x a n

computes x^(a/2^n).  a and n must be integers.

Computation of y = x^a done by augmenting the Double type to allow
arbitrarily large exponents, then y^(1/2^n) is computed with repeated
square root.

Copyright 2022 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, PackageImports #-}
module Main(main) where {
import System.Environment(getArgs);
import Data.Function((&));

main :: IO();
main = getArgs >>= \case{
[x,a,n] -> binaryexp (R (read x) 0) (read a) & iterate mysqrt & (!! (read n)) & tofloat & print;
_ -> undefined;
};

type Ii = Integer;
data R = R Double Ii deriving (Show);

instance Num R where {
(*) (R x a) (R y b) = reduce $ R (x*y) (a+b);
fromInteger i = reduce $ R (fromInteger i) 0;
(+) = undefined;
abs = undefined;
signum = undefined;
negate = undefined;
};

-- exponentiation by squaring
binaryexp :: forall a . (Num a) => a -> Ii -> a;
binaryexp _ 0 = 1;
binaryexp x 1 = x;
binaryexp _ negative | negative < 0 = error "negative exponent";
binaryexp x b = let {half ::Ii = div b 2
; s ::a = binaryexp x half
} in if half*2 == b
then s*s
else s*s*x;

reduce :: R -> R;
reduce x@(R m e) = if m>2
then reduce $ R (m/2) (e+1)
else x;

mysqrt :: R -> R;
mysqrt (R m e) = case divMod e 2 of {
(q,0) -> R (sqrt m) q;
_ -> mysqrt $ R (m*2) (e-1);
};

tofloat :: R -> Double;
tofloat (R m e) | e<0 = m/fromInteger(2^negate e)
| True = m * fromInteger(2^e);
} --end