"""Sinusoids of multiples of angles, via polynomials.

See S and C for details, or
http://www.chaos.org.uk/~eddy/math/multiangle.html
for theory.

$Id: multiangle.py,v 1.11 2007/03/24 22:42:21 eddy Exp $
"""

from polynomial import Polynomial

class LazySeq:
    def __init__(self, first=Polynomial(1)): self.__seq = [ first ]

    def __growto(self, key):
        val = self.growto(key)

        if __debug__:
            i = val.rank
            while i >= 0:
                assert self.check(val.coefficient(i))
                i = i-1

        if key >= len(self.__seq):
            self.__seq = self.__seq + [ None ] * (1 + key - len(self.__seq))

        self.__seq[key] = val

        return val

    def check(self, val): return val == long(val)

    def __getitem__(self, key):
        try: ans = self.__seq[key]
        except IndexError: return self.__growto(key)
        else:
            if ans is None:
                return self.__growto(key)

        return ans

    def __getslice__(self, lo, hi, step=1):
        ans = []
        if step > 0:
            while lo < hi:
                ans.append(self[lo])
                lo = lo + step
        elif step < 0:
            while lo > hi:
                ans.append(self[lo])
                lo = lo + step
        elif lo != hi:
            raise ValueError("Sequence with zero step is a very bad idea !")

        return ans

z = Polynomial(0, 1)

class Middle (LazySeq):
    def growto(self, key,
               down={0: z-1, 1: 1-z}, up=z+1, two=z+2, slab=z*(z+2),
               checks=(lambda p, x=z: (p(-x)*p(x-2)).unafter(x*(2-x)),
                       lambda p, x=z: ((p(x*(x+2)-2)*(x+1)+2)/(x+3))**.5,
                       lambda p, x=z: ((4+(x-1)*p(-x-2)**2)/(x+3))**.5)):
        assert K is self
        n, r = divmod(key, 3)
        m = n + r - 1
        assert key == 2*n + m + 1
        this, that = self[n], self[m]
        ans = .5*(down[(n+m)%2]*this(-slab)*that(-two) + up*this(slab-2)*that)
        assert not filter(lambda p, a=ans: p(a) != a, checks)
        return ans

K = Middle()
del Middle

term = 1-z*z

class seqCos (LazySeq):
    """Sequence of polynomials describing cos(n.t) in terms of cos(t)

    C[n](cos(t)) = cos(n.t)\n"""
    def growto(self, key, kate=1-4*term, term=-term, x=z):
        n, r = divmod(key, 2)
        if r: ans = x * K[n](kate)
        else: ans = self[n](term)

        assert ans.rank == key
        return ans

class seqSin (LazySeq):
    """Sequence of polynomials describing sin(n.t)/sin(t) in terms of cos(t)

    sin(t)*S[n](cos(t)) == sin(n.t)\n"""

    def growto(self, key, stem=4*term-3, term=term, x=z):
        n, r = divmod(key, 2)
        if r: ans = 2 * x * S[n](-term)
        else: ans = K[n](stem) * {0: 1, 1: -1}[n%2]

        assert ans.rank == key
        return ans

# C[n] and S[n+1] each have n roots between -1 and +1, symmetrically placed about 0
C, S = seqCos(), seqSin()
del seqCos, seqSin, LazySeq
del term, z, Polynomial