module Data.Realnumber where
import Data.Bits
import GHC.Real
import Debug.Trace
--debug = flip trace


-- Ein einfaches Tern\"ares Stellenwertsystem zur Basis 2
-- known as "signed-digit number", see wikipedia or "redundant binary representation"

-- Ziffer -1 (dargestellt mit "_")
-- Ziffer 0  (dargestellt mit "0")
-- Ziffer 1  (dargestellt mit "1")

-- endless sequence of figure values -1,0,1
type FigureSequence = [Int]

-- Value of the most significant figure
figureValue :: FigureSequence -> Int
figureValue (-1:_) = -1
figureValue ( 0:_) =  0
figureValue ( 1:_) =  1

-- return values: -3, -2, -1, 0, 1, 2, 3
twoFigureValue :: FigureSequence -> Int
twoFigureValue fs = 2 * (figureValue fs) + (figureValue . nextFig) fs

nFigureValue :: Integer -> FigureSequence -> Integer
nFigureValue n fs = acc (n, fs, 0)
	where
	acc (0,   fs, v) = v
	acc (i,-1:fs, v) = acc (i-1, fs, 2*v -1)
	acc (i, 0:fs, v) = acc (i-1, fs, 2*v)
	acc (i, 1:fs, v) = acc (i-1, fs, 2*v +1)

showFigureSequence :: FigureSequence -> [Char]
-- show never terminates!
showFigureSequence (-1:fs) = "_" ++ (showFigureSequence fs)
showFigureSequence ( 0:fs) = "0" ++ (showFigureSequence fs)
showFigureSequence ( 1:fs) = "1" ++ (showFigureSequence fs)

-- Sequence with zeros: Shall be unique
zeroS :: FigureSequence
zeroS = 0:zeroS

-- forbidden:
-- seq1 =  1:seq1 (equals [ 1, 1, 1,...])
-- seq2 = -1:seq2 (equals [-1,-1,-1,...])
-- otherwise zero == R 0 (-1:seq1) == R 0 (1:seq2)

nextFig :: FigureSequence -> FigureSequence
nextFig = tail

uberNextFig :: FigureSequence -> FigureSequence
uberNextFig = nextFig . nextFig

negateSequence :: FigureSequence -> FigureSequence
negateSequence (-1:fs) =  1:(negateSequence fs)
negateSequence ( 0:fs) =  0:(negateSequence fs)
negateSequence ( 1:fs) = -1:(negateSequence fs)

data RealNumber = R Integer FigureSequence

instance Show RealNumber where
	show (R i fs) = "2^" ++ (show i) ++ " * 0." ++ (showFigureSequence fs)

-----------------------------
-- sequencePlus: A simple adding algorithm.

-- Add the most significant figures of the two sequences
-- return values: 0, 1, -1, 2, -2
sumWithOneFigure :: (FigureSequence, FigureSequence) -> Int
sumWithOneFigure (xs, ys) = figureValue xs + figureValue ys

-- the return sequence has one additional carry figure first!
-- use like:  2^n * fs1 + 2^n * fs2 = 2^(n+1) * sequencePlus (fs1, fs2)
sequencePlus :: (FigureSequence, FigureSequence) -> FigureSequence
sequencePlus (fs1, fs2) = plusAcc (nextFig fs1, nextFig fs2, sumWithOneFigure (fs1, fs2))
	where
	plusAcc (fs1, fs2, c)
		-- case 1: (s-2, s+2) inside of (-8,  0)
		| -6 <= s && s <= -2 = -1:(plusAcc (nfs1, nfs2, s + 4))
		-- case 2: (s-2, s+2) inside of (-4, -4)
		| -2 <= s && s <=  2 =  0:(plusAcc (nfs1, nfs2, s + 0))
		-- case 3: (s-2, s+2) inside of ( 0,  8)
		|  2 <= s && s <=  6 =  1:(plusAcc (nfs1, nfs2, s - 4))
		where
		s = sumWithOneFigure (fs1, fs2) + 2 * c
		nfs1 = nextFig fs1
		nfs2 = nextFig fs2

-----------------------------
-- sequenceProduct: A simple product algorithm

-- use like: 2^n * ,fs1 * 2^m * ,fs2 = 2^(n+m) * ,sequenceProduct (fs1, fs2)
sequenceProduct :: (FigureSequence, FigureSequence) -> FigureSequence

sequenceProduct (fs1, fs2)
	| -1 <= v1 && v1 <= 1
		= 0: sequenceProduct (v1 : r1, fs2)
	| -1 <= v2 && v2 <= 1 -- a * b = b * a
		= sequenceProduct (fs2, fs1)

	| v1 < 0 && v2 < 0 -- (-a) * (-b) =   a * b
		= sequenceProduct (negateSequence fs1, negateSequence fs2)
	| v1 < 0 || v2 < 0 -- (-a) *   b  = -(a * b)
		= negateSequence (sequenceProduct (negateSequence fs1, fs2))

	| v1 == 3 && v2 == 2 -- a * b = b * a
		= sequenceProduct (fs2, fs1)

	| v1 == 3 && v2 == 3
		= 1: sequencePlus (sequenceProduct (1:r1, 1:r2),
		                   sequencePlus (r1, r2))
	--  If the following code fragment is called recursive,
	--  then sequencePlus would never give any figure,
	--  because sequencePlus requires two figures
	--  of the recursive call "sequenceProduct (r1, fs2)"
	--  to compute at least one figure.
	--example:
	-- > let fs1 = (1: -1: fs1 :: FigureSequence)
	-- > let fs2 = (1:  0: fs2 :: FigureSequence)
	-- > sequenceProduct (fs1, fs2)
	-- [0,1
	--  after 1 is printed, no more figure will appear
	--buggy code: 
--	| v1 == 2 && v2 == 3
--		= 1: sequencePlus (sequenceProduct (r1, fs2),
--		                   -1: r2)
	--better code:
	| v1 == 2 && v2 == 3
		= 1: sequencePlus (sequencePlus(r1, r2),
		                   -1: sequenceProduct(r1, 1:r2))
	| v1 == 2 && v2 == 2
		= 1: (forceSimplify . sequencePlus)
                             (-1: sequencePlus (r1, r2),
		               0:0: sequenceProduct (r1, r2))
	where
	r1 = uberNextFig fs1
	r2 = uberNextFig fs2
	v1 = twoFigureValue fs1
	v2 = twoFigureValue fs2

-----------------------------
-- absolute value
absSequence :: FigureSequence -> FigureSequence
absSequence ( 0:fs) = 0: (absSequence fs)
absSequence ( 1:fs) = 1: fs
absSequence (-1:fs) = 1: (negateSequence fs)

-----------------------------
-- signum
-- caution:
-- (signum zeroS) never returns!
sequenceSignum :: FigureSequence -> FigureSequence
sequenceSignum ( 0: fs) = sequenceSignum fs
sequenceSignum (-1: _) = -1: zeroS
sequenceSignum ( 1: _) =  1: zeroS

-----------------------------
-- compareZero
compareZero :: RealNumber -> Ordering
compareZero (R e (-1:fs)) = LT
compareZero (R e ( 1:fs)) = GT
compareZero (R e ( 0:fs)) = compareZero (R (e-1) fs)

-----------------------------
-- maximum
maxReal :: RealNumber -> RealNumber -> RealNumber
maxReal (R e1 fs1) (R e2 fs2)
	| e1 > e2 = maxReal (R e1 fs1) (R e1 (headingZeros (e1-e2) fs2))
	| e2 > e1 = maxReal (R e2 (headingZeros (e2-e1) fs1)) (R e2 fs2)
	| otherwise = R e1 (maxSequence 0 fs1 fs2)

maxSequence :: Int -> FigureSequence -> FigureSequence -> FigureSequence
maxSequence u (f1:fs1) (f2:fs2)
	| u >  1  = f1:fs1
	| u >  0  = f1:maxSequence (2*u +f1-f2) fs1 fs2
	| u < -1  = f2:fs2
	| u <  0  = f2:maxSequence (2*u +f1-f2) fs1 fs2
	-- u == 0
	| f1 >=f2 = f1:maxSequence (     f1-f2) fs1 fs2
	| f2 >=f1 = f2:maxSequence (     f1-f2) fs1 fs2

minReal :: RealNumber -> RealNumber -> RealNumber
minReal x y = -(maxReal (-x) (-y))

-----------------------------
-- integer reciprocal
einsdurch :: Integer -> FigureSequence
einsdurch 0 = error "division by zero"
einsdurch m = helper (1, m)
	where
	helper (n, m)
		| 2*m*n < -m*m = -1: (helper (2*(n+m), m))
		| 2*m*n >  m*m =  1: (helper (2*(n-m), m))
		| otherwise    =  0: (helper (2*n,     m))

-----------------------------
-- helper functions
msb :: Integer -> Int
msb i = msbH (abs i, -1)
	where
	msbH :: (Integer, Int) -> Int
	msbH (i, erg)
		| i == 0  = erg
		|otherwise= msbH (shiftR i 1, erg+1) 

headingZeros :: Integer -> FigureSequence -> FigureSequence
headingZeros 0 fs = fs
headingZeros n fs = headingZeros (n-1) ( 0: fs)

-- try to cut up the first figure
-- like 2^1 * 0.01 == 2^0 * 0.1
--      2^-1 * 0.1_ == 2^-2 * 0.1
--      2^5 * 0._1 == 2^4 * 0._
-- useful, for example, if we don't need the carry bit of the sum of two numbers
simplify :: RealNumber -> RealNumber
simplify (R e ( 0: fs))     = R (e-1) fs
simplify (R e ( 1: (-1: fs))) = R (e-1) ( 1: fs)
simplify (R e (-1: ( 1: fs))) = R (e-1) (-1: fs)

-- We cannot cut the first figure applying obove rules, fall back to original x
simplify x = x

-- Warning: This function does not return on input zero
simplifyAll :: RealNumber -> RealNumber
simplifyAll (R e ( 0: fs))     = simplifyAll (R (e-1) fs)
simplifyAll (R e ( 1: (-1: fs))) = simplifyAll (R (e-1) ( 1: fs))
simplifyAll (R e (-1: ( 1: fs))) = simplifyAll (R (e-1) (-1: fs))
simplifyAll x = x

-- Z^(-1)
forceSimplify :: FigureSequence -> FigureSequence
forceSimplify ( 0: fs)     = fs
forceSimplify ( 1: (-1: fs)) =  1: fs
forceSimplify (-1: ( 1: fs)) = -1: fs
forceSimplify fs = error (show (take 5 fs))

forceNSimplify :: Integer -> FigureSequence -> FigureSequence
forceNSimplify 0 fs = fs
forceNSimplify n fs = forceNSimplify (n-1) (forceSimplify fs)

-- shift
shiftReal :: RealNumber -> Integer -> RealNumber
shiftReal (R e fs) s = R (e+s) fs

-----------------------------
-- instances for RealNumber
instance Num RealNumber where
	negate (R e fs) = R e (negateSequence fs)

	-- FIXME: support integers with more than 2^32 (2^64) bits (limited by testBit :: a -> "Int" -> Bool)
	fromInteger i
		| i < 0   = negate (fromInteger (-i))
		|otherwise= R (toInteger (msb i +1)) (helper (msb i))
		where
		helper v
			| v < 0     = zeroS
			|testBit i v=  1: (helper (v-1))
			| otherwise =  0: (helper (v-1))
	R e1 fs1 + R e2 fs2
		| e1  > e2 = (R e2 fs2) + (R e1 fs1)
		| e1 <= e2 = simplify (R (e2+1) fsNew)
		where
		fsNew = sequencePlus (headingZeros (e2-e1) fs1, fs2)
	R e1 fs1 * R e2 fs2
		= R (e2+e1) (sequenceProduct (fs1, fs2))
	abs (R e fs) = R e (absSequence fs)

	-- caution:
	-- if x == zero, signum never returns.
	signum (R e fs) = R 1 (sequenceSignum fs)

instance Eq RealNumber where
	-- co-semi-entscheidbar
	x == y = compareZero (x-y) == EQ

instance Ord RealNumber where
	compare x y = compareZero (x-y)

-----------------------------
-- approximations

-- [Int] is an ending(!) sequence of figure values -1,0,1
data ApproxReal = A Integer [Int]

takeI :: Integer -> [a] -> [a]
takeI _ [] = []
takeI i (x:xs)
	| i <= 0    = []
	| otherwise = x: takeI (i-1) xs

approx :: RealNumber -> Integer -> ApproxReal
approx (R e fs) i
	| e >  i = A e (takeI (e-i) fs)
	| e <= i = A i []

type Intervall = (RealNumber, RealNumber)

intervallToApprox :: Intervall -> ApproxReal
intervallToApprox (R e1 (f1':fs1'), R e2 (f2':fs2'))
	| e1 <  e2 = intervallToApprox (R e2 (headingZeros (e2-e1) (f1':fs1')), R e2 (f2':fs2'))
	| e1 >  e2 = intervallToApprox (R e1 (f1':fs1'), R e1 (headingZeros (e1-e2) (f2':fs2')))
	| e1 == e2 = A e1 (helper f1' fs1' f2' fs2')
	where
	helper u1 (f1:fs1) u2 (f2:fs2)
		| 2*u2 +f2 +1 < 2*u1 +f1 || u2 < -1 || u1 > 1 = error "in intervallToApprox: (a, b) with b > a is forbidden!"

		| abs (2*u1 +f1) < 2 && abs (2*u2 +f2) < 2
		               =  0 : helper (2*u1 +f1   ) fs1 (2*u2 +f2   ) fs2

		| 2*u2 +f2 < 0 && 2*u1 +f1 < 0 = -1 : helper (2*u1 +f1 +2) fs1 (2*u2 +f2 +2) fs2
		| 2*u1 +f1 > 0 && 2*u1 +f1 > 0 =  1 : helper (2*u1 +f1 -2) fs1 (2*u2 +f2 -2) fs2


		| otherwise    = []

instance Show ApproxReal where
	show (A e fs) = "2^" ++ (show e) ++ " * 0." ++ (showEndingSequence fs)
		where
		showEndingSequence (-1:fs) = "_" ++ (showEndingSequence fs)
		showEndingSequence ( 0:fs) = "0" ++ (showEndingSequence fs)
		showEndingSequence ( 1:fs) = "1" ++ (showEndingSequence fs)
		showEndingSequence [] = "&c"

type ApproxRealSchachtelung = [ApproxReal]

realFromAppRealSchach :: ApproxRealSchachtelung -> RealNumber
realFromAppRealSchach (A e fs : ars)
	= R e (fs ++ moreFigures fs ars)
	where
	moreFigures fs_fix (A e_next fs_next : ars)
		-- next better approximation has (e_next - e) more most significant figures of value 0: cut them
		| e_next > e = moreFigures fs_fix (A e (forceNSimplify (e_next -e) fs_next) : ars)
		-- next better approximation has less figures
		| e_next < e = moreFigures fs_fix (A e (headingZeros (e -e_next) fs_next) :ars)
		| otherwise  = fs_next_rest ++ (moreFigures (fs_fix ++ fs_next_rest) ars)
		where
		fs_next_rest = getBest 0 fs_fix fs_next

		getBest u [] (f2:fs2)
			| abs (2*u +f2) > 1 = error "approximations given to realFromAppRealSchach are not nested!"
			| otherwise      = (2*u +f2):fs2
		getBest _ []  [] = []
		getBest u fs1 [] = [] -- `debug` ("unusual: getBest " ++ show u ++ " " ++ show fs1 ++ " []")
		getBest u (f1:fs1) (f2:fs2) = getBest (2*u +f2-f1) fs1 fs2

type IntervallSchachtelung = [Intervall]
approxRealSchachtelung :: IntervallSchachtelung -> ApproxRealSchachtelung
approxRealSchachtelung = map intervallToApprox

realFromSchach :: IntervallSchachtelung -> RealNumber
realFromSchach = realFromAppRealSchach . approxRealSchachtelung

-----------------------------
-- reciprocal
intervallEinsdurch :: RealNumber -> IntervallSchachtelung
intervallEinsdurch xs = map einsdurchSequence [0..]
	where
	e :: Integer
	fs :: FigureSequence
	R e fs = simplifyAll xs
	einsdurchSequence :: Integer -> (RealNumber, RealNumber)
	einsdurchSequence m = (R (-e+m+3) (einsdurch v1), R (-e+m+3) (einsdurch v_))
		where
		v_ = nFigureValue (m+2) fs - 1
		v1 = nFigureValue (m+2) fs + 1

-----------------------------
-- fractional RealNumber
instance Fractional RealNumber where
	fromRational (x :% y) = fromInteger x * R 1 (einsdurch y)

	recip = realFromSchach . intervallEinsdurch

-----------------------------
-- often used constants for testing
zero = R 0 zeroS
one = R 1 ( 1: zeroS)

-- tests
{-
let show50 x = approx x (-50)

show50 $ one
show50 $ zero
show50 $ one + zero
show50 $ one + one
show50 $ - one
show50 $ one - one
show50 $ one * one
show50 $ one * (-one)
show50 $ (one+one+one) * (-one-one-one)
show50 $ (one+one+one) * (one+one+one)
show50 (fromRational (2 :% 7) :: RealNumber)
show50 $ one * one * one *one * one
show50 $ simplifyAll (one * one * one *one * one)

-- not halting:
show50 $ simplifyAll zero
show50 $ one / zero
-}