cexp

#!/usr/bin/env python3


import math
import cmath
import numpy


DEFAULT_P0 = -0.07969529411773515
DEFAULT_PSTEP = -73.63139197777613
DEFAULT_TOTAL = 16384

def double_exp(p0=DEFAULT_P0,
               pstep=DEFAULT_PSTEP,
               total=DEFAULT_TOTAL):
    exp = cmath.exp
    for count in range(total):
        yield exp(1j * (p0 + count * pstep))


def double_mult(p0=DEFAULT_P0,
               pstep=DEFAULT_PSTEP,
               total=DEFAULT_TOTAL):
    exp = cmath.exp
    ep0, estep = exp(1j * p0), exp(1j * pstep)
    for _ in range(total):
        yield ep0
        ep0 *= estep


def float_exp(p0=DEFAULT_P0,
              pstep=DEFAULT_PSTEP,
              total=DEFAULT_TOTAL):
    exp = numpy.exp
    c64 = numpy.complex64
    f32 = numpy.float32
    p0, pstep = c64(1j * p0), c64(1j * pstep)
    for count in range(total):
        count = f32(count)
        yield exp(p0 + count * pstep)


def float_mult(p0=DEFAULT_P0,
               pstep=DEFAULT_PSTEP,
               total=DEFAULT_TOTAL):
    exp = cmath.exp
    c64 = numpy.complex64
    ep0, epstep = c64(exp(1j * p0)), c64(exp(1j * pstep))
    for _ in range(total):
        yield ep0
        ep0 *= estep


def float_mod(p0=DEFAULT_P0,
              pstep=DEFAULT_PSTEP,
              total=DEFAULT_TOTAL):
    exp = numpy.exp
    c64 = numpy.complex64
    f32 = numpy.float32
    fmod = math.fmod
    twopi = 2. * math.pi
    p0, pstep = fmod(p0, twopi), fmod(pstep, twopi)
    p0, pstep = c64(1j * p0), c64(1j * pstep)
    for count in range(total):
        count = f32(count)
        yield exp(p0 + count * pstep)


def main():
    generators = (double_exp(), double_mult(), float_exp(), float_mult(),
                  float_mod())
    arrays = [numpy.fromiter(generator, dtype=numpy.complex64)
              for generator in generators]
    combined = numpy.vstack(arrays)
    rel_errors = abs(combined - combined[0]) / abs(combined[0])
    print(numpy.amax(rel_errors, 1))


if __name__ == '__main__':
    main()