--
-- base data types for use in the SOME project
--
-- John Ringland 2005/08/27

include FDData.e  -- includes SOME_Utilities.e

-- PNData is a sequence of 1,2,4 or 8 FDDatas
global type PNData(object d)
    integer len
    if sequence(d) then
        len = length(d)
        if len > 0 then
            for i = 1 to length(d) do
                if not FDData(d[i]) then
                    return 0
                end if
            end for
        else
            return 0
        end if
    else
        return 0
    end if
    return 1
end type

global function smSetNull(PNData d)
    for i = 1 to length(d) do
        d[i][FDIND_] = -1 -- set all FDDatas to undefined
    end for
    return d
end function

-- real => {FDData}
global type SMReal(object r)
    return PNData(r) and length(r)=1
end type

-- complex => {FDData, FDData}
global type SMComplex(object r)
    return PNData(r) and length(r)=2
end type

-- quaternion => {FDData, FDData, FDData, FDData}
global type SMQuaternion(object r)
    return PNData(r) and length(r)=4
end type

-- octonion => {FDData x 8}
global type SMOctonion(object r)
    return PNData(r) and length(r)=8
end type

-- RootVariable, q^2 = empirical value
global type SMRootVariable(object q)
    return PNData(q)
end type

-- tPower returns the transcendent power of the PNData
-- equivalent to z.z* i.e. conjugated
global function tPower(PNData d)
    FDData result result = fdZero
    for i = 1 to length(d) do
        --result += power(d[i],2)
        result = fdAdd(result,result,fdMultiply(d[i],d[i],d[i]))
    end for
    return {result} -- an SMReal
end function

-- ePower returns the empirical value of the PNData
-- equivalent to z.z i.e. NOT conjugated
global function ePower(PNData d)
    FDData result result = fdZero
    result = fdAdd(result,result,fdMultiply(d[1],d[1],d[1]))
    for i = 2 to length(d) do
        result = fdSubtract(result,result,fdMultiply(d[i],d[i],d[i]))
    end for
    return {result} -- an SMReal
end function

-- tAssign assigns a transcendent power to the PNData
-- without effecting the vectors orientation its power is defined
global function tAssign(PNData d, SMReal val)
    integer len
    FDData tp, multiplier
    FDValue fdv
    fdv = fdGet(val[1])
    len = length(d)
    tp = first(tPower(d))
    multiplier = fdZero

    if equal(tp,fdZero) then  -- all power is shared equally amongst the components
        -- not really a multiplier here, just a temp variable
        multiplier = fdSet(multiplier,{sqrt(fdv) / len})
        for i = 1 to len do
            d[i] = multiplier
        end for
        return d
    end if
    if fdv >= 0 then
        multiplier = fdSet(multiplier,sqrt(fdv))
        multiplier = fdDivide(multiplier,multiplier,tp)
        for i = 1 to len do
            d[i] = fdMultiply(d[i],d[i],multiplier)
        end for
        return d
    else -- val[1] < 0 then
        return smSetNull(d)
    end if
end function

-- eAssign assigns an empirical value to the PNData
-- which effects the vectors orientation
global function eAssign(PNData d, SMReal val)
    integer len
    FDValue fdv
    len = length(d)
    fdv = fdGet(val[1])
    if fdv >= 0 then -- all power to the real component
        d[1] = fdSet(d[1],sqrt(fdv))
        for i = 2 to len do
            d[i] = fdSet(d[i],0)
        end for
        return d
    else -- val[1] < 0 so power is distributed amongst the imaginary components
        d[1] = fdSet(d[1],0)
        for i = 2 to len do
            d[i] = fdSet(d[i],sqrt(fdv)/len)
        end for
        return d
    end if
end function

global function conjugate(PNData d)
    for i = 2 to length(d) do
        d[i] = fdMultiply(d[i],d[i],fdMinusOne)
    end for
    return d
end function

global function smAdd(PNData d, PNData d1, PNData d2)
    integer len len = length(d)
    if length(d1) = len and length(d2) = len then -- they can be added
        for i = 1 to len do
            d[i] = fdAdd(d[i],d1[i],d2[i])
        end for
        return d
    end if
    return smSetNull(d) -- incompatible data types
end function

global function smSubtract(PNData d, PNData d1, PNData d2)
    integer len len = length(d)
    if length(d1) = len and length(d2) = len then -- they can be subtracted
        for i = 1 to len do
            d[i] = fdSubtract(d[i],d1[i],d2[i])
        end for
        return d
    end if
    return smSetNull(d) -- incompatible data types
end function

global function smMultiply(PNData d, PNData d1, PNData d2)
    integer len, len1, len2
    len = length(d)
    len1 = length(d1)
    len2 = length(d2)
    if len1=1 and len2=len then
        for i = 1 to len do
            d[i] = fdMultiply(d[i],d1[1],d2[i])
        end for
        return d
    elsif len2=1 and len1=len then
        for i = 1 to len do
            d[i] = fdMultiply(d[i],d1[i],d2[1])
        end for
        return d
    elsif len1 = len2 and len2 = len then -- they can be multiplied out
        if len = 2 then -- they are complex and the elements are FDDatas
            return {fdSubtract(d1[1],fdMultiply(d1[1],d1[1],d2[1]),fdMultiply(d1[2],d1[2],d2[2])),
                    fdAdd(d1[1],fdMultiply(d1[1],d1[1],d2[2]),fdMultiply(d1[2],d1[2],d2[1]))}
        elsif len = 4 then -- they are quaternion
            -- not implemented yet
        elsif len = 8 then -- they are octonion
            -- not implemented yet
        end if
    end if
    return smSetNull(d)
end function

-- need to implement this one ...
global function smDivide(PNData d, PNData d1, PNData d2)
    d1 = d1
    d2 = d2
    return smSetNull(d) -- not implemented yet
end function


procedure testPNData()
    FDData v0, v1
    SMComplex c1, c2, c3, c4,c5,c6,tmp
    --atom b, min, max
    
    v0 = fdZero
    v1 = fdOne
    trace(1)
    ?fdGet(v0)
    ?fdGet(v1)
    c1 = {v1,v0}
    c2 = {v0,v1}
    c3 = {v1,v1}
    tmp = smAdd(c1,c1,c2)
    ?{fdGet(tmp[1]),fdGet(tmp[2])}
    tmp = smSubtract(c1,c1,c2)
    ?{fdGet(tmp[1]),fdGet(tmp[2])}
    tmp = smMultiply(c1,c1,c2)
    ?{fdGet(tmp[1]),fdGet(tmp[2])}
    tmp = smDivide(c1,c1,c2)
    ?{fdGet(tmp[1]),fdGet(tmp[2])}
    c4 = c1
    c5 = c2
    c6 = c3
    ?c1
    ?tPower(c1)
    ?fdGet(first(tPower(c1)))
    ?ePower(c1)
    ?fdGet(first(ePower(c1)))
    c1 = smMultiply(c1,c1,c1)
    ?c1
    ?{fdGet(c1[1]),fdGet(c1[2])}
    ?tPower(c2)
    ?fdGet(first(tPower(c2)))
    ?ePower(c2)
    ?fdGet(first(ePower(c2)))
    c2 = smMultiply(c2,c2,c2)
    ?c2
    ?{fdGet(c2[1]),fdGet(c2[2])}
    ?tPower(c3)
    ?fdGet(first(tPower(c3)))
    ?ePower(c3)
    ?fdGet(first(ePower(c3)))
    c3 = smMultiply(c3,c3,c3)
    ?c3
    ?{fdGet(c3[1]),fdGet(c3[2])}
    trace(1)
    ?{fdGet(c4[1]),fdGet(c4[2])}
    ?{fdGet(c5[1]),fdGet(c5[2])}
    c5 = conjugate(c5)
    ?{fdGet(c5[1]),fdGet(c5[2])}
    ?{fdGet(c6[1]),fdGet(c6[2])}
    ?fdGet(first(tPower(c6)))
    ?fdGet(first(ePower(c6)))
    c1 = smMultiply(c6,c6, conjugate(c6))
    ?{fdGet(c1[1]),fdGet(c1[2])}
    
end procedure

--testPNData()
-- I think this is working, 
-- so I can now add subtract and multiply, 
-- real or complex, finite discrete data.