https en wikipedia orgwikiautomaticdifferentiationautomaticdifferentiationusingdualnumbers https blog jliszka org20131024exactnumericnthderivatives html note this only works for basic functions fx where the power of x is positive differentiatelambda x x2 2 2 2 differentiatelambda x x2 x4 9 2 196830 differentiatelambda y 0 5 y 3 6 3 5 4 7605 0 differentiatelambda y y 2 4 3 0 differentiate8 8 8 traceback most recent call last valueerror differentiate requires a function as input for func differentiatelambda x x 2 1 traceback most recent call last valueerror differentiate requires a float as input for position differentiatelambda x x2 3 traceback most recent call last valueerror differentiate requires an int as input for order https en wikipedia org wiki automatic_differentiation automatic_differentiation_using_dual_numbers https blog jliszka org 2013 10 24 exact numeric nth derivatives html note this only works for basic functions f x where the power of x is positive differentiate lambda x x 2 2 2 2 differentiate lambda x x 2 x 4 9 2 196830 differentiate lambda y 0 5 y 3 6 3 5 4 7605 0 differentiate lambda y y 2 4 3 0 differentiate 8 8 8 traceback most recent call last valueerror differentiate requires a function as input for func differentiate lambda x x 2 1 traceback most recent call last valueerror differentiate requires a float as input for position differentiate lambda x x 2 3 traceback most recent call last valueerror differentiate requires an int as input for order
from math import factorial




class Dual:
    def __init__(self, real, rank):
        self.real = real
        if isinstance(rank, int):
            self.duals = [1] * rank
        else:
            self.duals = rank

    def __repr__(self):
        return (
            f"{self.real}+"
            f"{'+'.join(str(dual)+'E'+str(n+1)for n,dual in enumerate(self.duals))}"
        )

    def reduce(self):
        cur = self.duals.copy()
        while cur[-1] == 0:
            cur.pop(-1)
        return Dual(self.real, cur)

    def __add__(self, other):
        if not isinstance(other, Dual):
            return Dual(self.real + other, self.duals)
        s_dual = self.duals.copy()
        o_dual = other.duals.copy()
        if len(s_dual) > len(o_dual):
            o_dual.extend([1] * (len(s_dual) - len(o_dual)))
        elif len(s_dual) < len(o_dual):
            s_dual.extend([1] * (len(o_dual) - len(s_dual)))
        new_duals = []
        for i in range(len(s_dual)):
            new_duals.append(s_dual[i] + o_dual[i])
        return Dual(self.real + other.real, new_duals)

    __radd__ = __add__

    def __sub__(self, other):
        return self + other * -1

    def __mul__(self, other):
        if not isinstance(other, Dual):
            new_duals = []
            for i in self.duals:
                new_duals.append(i * other)
            return Dual(self.real * other, new_duals)
        new_duals = [0] * (len(self.duals) + len(other.duals) + 1)
        for i, item in enumerate(self.duals):
            for j, jtem in enumerate(other.duals):
                new_duals[i + j + 1] += item * jtem
        for k in range(len(self.duals)):
            new_duals[k] += self.duals[k] * other.real
        for index in range(len(other.duals)):
            new_duals[index] += other.duals[index] * self.real
        return Dual(self.real * other.real, new_duals)

    __rmul__ = __mul__

    def __truediv__(self, other):
        if not isinstance(other, Dual):
            new_duals = []
            for i in self.duals:
                new_duals.append(i / other)
            return Dual(self.real / other, new_duals)
        raise ValueError

    def __floordiv__(self, other):
        if not isinstance(other, Dual):
            new_duals = []
            for i in self.duals:
                new_duals.append(i // other)
            return Dual(self.real // other, new_duals)
        raise ValueError

    def __pow__(self, n):
        if n < 0 or isinstance(n, float):
            raise ValueError("power must be a positive integer")
        if n == 0:
            return 1
        if n == 1:
            return self
        x = self
        for _ in range(n - 1):
            x *= self
        return x


def differentiate(func, position, order):
    
    if not callable(func):
        raise ValueError("differentiate() requires a function as input for func")
    if not isinstance(position, (float, int)):
        raise ValueError("differentiate() requires a float as input for position")
    if not isinstance(order, int):
        raise ValueError("differentiate() requires an int as input for order")
    d = Dual(position, 1)
    result = func(d)
    if order == 0:
        return result.real
    return result.duals[order - 1] * factorial(order)


if __name__ == "__main__":
    import doctest

    doctest.testmod()

    def f(y):
        return y**2 * y**4

    print(differentiate(f, 9, 2))