#!/usr/bin/python3
# -*- coding: utf-8 -*-
"""Another quick “reactive” calculator.

This time I’m using RPN to minimize parsing effort.  The input looks
like line noise, so the output includes infix formulas to check it.

Here’s an example session, showing how it builds up a stack of
expressions and re-evaluates them after every input.  First, we set up
the quadratic formula:

    Calculator.  Type ? for help.
    b 2^4a×c×-.5^d!b _d-2a*/b _d+2a*/
    (-b - d) / (2 * a)=(lacking b) (-b + d) / (2 * a)=(lacking b)
    .l
    d ((b ^ 2 - 4 × a × c) ^ 0.5): (lacking b)

Now we have defined d as the quadratic formula’s discriminant, and we
have a stack with the two roots from the quadratic formula on it, but
this is all expressed in terms of the coefficients a, b, and c, which
we haven’t defined yet.  Let’s solve x² - 4 = 0:

    1a!0b!4_c!
    (-b - d) / (2 * a)=-2.0 (-b + d) / (2 * a)=2.0

So, that's correct; x² - 4 has zeros at ±2.  Let’s compute √2:

    2_c!
    (-b - d) / (2 * a)=-1.4142135623730951 (-b + d) / (2 * a)=1.4142135623730951

Cool.  Also we could compute √2 more directly:

    2 .5**
    (-b - d) / (2 * a)=-1.4142135623730951 (-b + d) / (2 * a)=1.4142135623730951
    2 ** 0.5=1.4142135623730951

We can discard this result with ;:

    ;
    (-b - d) / (2 * a)=-1.4142135623730951 (-b + d) / (2 * a)=1.4142135623730951

How about the solutions of x² + 3x = 0, which should be 0 and
-3?

    3b!0c!
    (-b - d) / (2 * a)=-3.0 (-b + d) / (2 * a)=0.0

If we list the variables now we can see they’re all defined:

    .l
    a (1): 1
    b (3): 3
    c (0): 0
    d ((b ^ 2 - 4 × a × c) ^ 0.5): 3.0

After removing `b` the discriminant is no longer defined; if we press
Enter again we can see that this undefined the solutions as well:

    .r b
    .l
    a (1): 1
    c (0): 0
    d ((b ^ 2 - 4 × a × c) ^ 0.5): (lacking b)

    (-b - d) / (2 * a)=(lacking b) (-b + d) / (2 * a)=(lacking b)

But then we can set it again, and change c, and see the solutions:

    5b!
    (-b - d) / (2 * a)=-5.0 (-b + d) / (2 * a)=0.0
    1c!
    (-b - d) / (2 * a)=-4.7912878474779195 (-b + d) / (2 * a)=-0.20871215252208009

Complex numbers are supported for both output and input:

    0b!
    (-b - d) / (2 * a)=(-6.123031769111886e-17-1j)
    (-b + d) / (2 * a)=(6.123031769111886e-17+1j)
    5jb!
    (-b - d) / (2 * a)=(-1.6486767597993924e-16-5.192582403567252j)
    (-b + d) / (2 * a)=(1.6486767597993924e-16+0.19258240356725187j)
    ?b
    b ((5 × j)): 5j

There’s a “@” operation which can inline the current formula of a
variable:

    b _d@-2a*/ b _d@+2a*/
    (-b - d) / (2 * a)=(-1.6486767597993924e-16-5.192582403567252j)
    (-b + d) / (2 * a)=(1.6486767597993924e-16+0.19258240356725187j)
    (-b - (b ^ 2 - 4 × a × c) ^ 0.5) / (2 * a)=(-1.6486767597993924e-16-5.192582403567252j)
    (-b + (b ^ 2 - 4 × a × c) ^ 0.5) / (2 * a)=(1.6486767597993924e-16+0.19258240356725187j)

Now let’s check that these two solutions are actually correct, and use
;; to discard the other two redundant items from the stack:

    s1!s2!;;
    .l
    a (1): 1
    b ((5 × j)): 5j
    c (1): 1
    d ((b ^ 2 - 4 × a × c) ^ 0.5): (3.297353519598785e-16+5.385164807134504j)
    s1 ((-b + (b ^ 2 - 4 × a × c) ^ 0.5) / (2 * a)): (1.6486767597993924e-16+0.19258240356725187j)
    s2 ((-b - (b ^ 2 - 4 × a × c) ^ 0.5) / (2 * a)): (-1.6486767597993924e-16-5.192582403567252j)

Okay, so let’s check what happens when x=s₁:

    a s1 2** * b s1 * + c +
    (7.771561172376096e-16+8.878396065212234e-16j) (a * s1 ** 2 + b * s1 + c)

That seems like a reasonable rounding error to 0; now let’s store that
in z₁, and do the same thing for the other way of calculating what we
hope is zero:

    z1!
    a s2 2** * b s2 * + c +
    a * s2 ** 2 + b * s2 + c=(-3.552713678800501e-15+8.878396065212234e-16j)
    z2!

    .l
    a (1): 1
    b ((5 × j)): 5j
    c (1): 1
    d ((b ^ 2 - 4 × a × c) ^ 0.5): (3.297353519598785e-16+5.385164807134504j)
    s1 ((-b + (b ^ 2 - 4 × a × c) ^ 0.5) / (2 * a)): (1.6486767597993924e-16+0.19258240356725187j)
    s2 ((-b - (b ^ 2 - 4 × a × c) ^ 0.5) / (2 * a)): (-1.6486767597993924e-16-5.192582403567252j)
    z1 (a * s1 ** 2 + b * s1 + c): (7.771561172376096e-16+8.878396065212234e-16j)
    z2 (a * s2 ** 2 + b * s2 + c): (-3.552713678800501e-15+8.878396065212234e-16j)

If we put z1 and z2 on the stack they’ll get redisplayed whenever we
change one of the coefficients:

    z1 z2
    z1=(7.771561172376096e-16+8.878396065212234e-16j)
    z2=(-3.552713678800501e-15+8.878396065212234e-16j)

As before, we can replace z2 with its definition with @:

    @
    z1=(7.771561172376096e-16+8.878396065212234e-16j)
    a * s2 ** 2 + b * s2 + c=(-3.552713678800501e-15+8.878396065212234e-16j)

Changing a, b, and c:

    .51a!
    z1=(-8.881784197001252e-16+1.623203716046916e-15j)
    a * s2 ** 2 + b * s2 + c=1.623203716046916e-15j
    0b!
    z1=(2.220446049250313e-16+1.2246063538223773e-16j)
    a * s2 ** 2 + b * s2 + c=(2.220446049250313e-16+1.2246063538223773e-16j)
    0a!
    error evaling Var(name='z1')
    Traceback (most recent call last):
    ...
    ZeroDivisionError: float division by zero
    1a!
    z1=1.2246063538223773e-16j a * s2 ** 2 + b * s2 + c=1.2246063538223773e-16j
    387c!
    z1=(5.684341886080802e-14+4.7392265892925994e-14j)
    a * s2 ** 2 + b * s2 + c=(5.684341886080802e-14+4.7392265892925994e-14j)

Seems to work!

The @ operation applied to a formula breaks it apart, and applied to a
constant removes it from the stack, like ;, so by repeated @ you can
find out what a formula is doing:

    ;@@
    a * s1 ** 2 + b * s1=(-386.99999999999994+4.7392265892925994e-14j) c=387
    @
    a * s1 ** 2 + b * s1=(-386.99999999999994+4.7392265892925994e-14j) 387=387
    @
    a * s1 ** 2 + b * s1=(-386.99999999999994+4.7392265892925994e-14j)
    @
    a * s1 ** 2=(-386.99999999999994+4.7392265892925994e-14j) b * s1=0j
    @
    a * s1 ** 2=(-386.99999999999994+4.7392265892925994e-14j) b=0
    s1=(1.2045421322489794e-15+19.672315572906j)
    @
    a * s1 ** 2=(-386.99999999999994+4.7392265892925994e-14j) b=0
    (-b + (b ^ 2 - 4 × a × c) ^ 0.5) / (2 * a)=(1.2045421322489794e-15+19.672315572906j)
    ;
    a * s1 ** 2=(-386.99999999999994+4.7392265892925994e-14j) b=0
    @
    a * s1 ** 2=(-386.99999999999994+4.7392265892925994e-14j) 0=0
    @
    a * s1 ** 2=(-386.99999999999994+4.7392265892925994e-14j)
    @
    a=1 s1 ** 2=(-386.99999999999994+4.7392265892925994e-14j)
    @
    a=1 s1=(1.2045421322489794e-15+19.672315572906j) 2=2

You can thus also use @ to edit the definitions of variables:

    ;;;3y+
    3 + y=(lacking y)
    4y!
    3 + y=7
    x!
    ?x
    x (3 + y): 7
    x@@-x!
    ?x
    x (3 - y): -1
    x@5+x!
    x
    x=4
    @
    3 - y + 5=4
    ;x@_x!
    ?x
    x (-(3 - y + 5)): -4

Here’s another session where we calculate the resistance of a 100-watt
lightbulb at two different line voltages, also verifying numerically
that we didn’t screw up the algebra:

    Calculator.  Type ? for help.
    # how many ohms is a 100-watt bulb? P = V²/R, P/V² = 1/R, R = V²/P?
    V 2** P /
    V ** 2 / P=(lacking V)
    R!
    V 2** R/
    V ** 2 / R=(lacking V)
    240V!
    V ** 2 / R=(lacking P)
    100P!
    V ** 2 / R=100.0
    ?R
    R (V ** 2 / P): 576.0
    120V!
    V ** 2 / R=100.0
    ?R
    R (V ** 2 / P): 144.0

Trying to solve it in a dumber way:

    Calculator.  Type ? for help.
    V R/I!  V I*P!
    .l
    I (V / R): (lacking V)
    P (V * I): (lacking V)
    240V!
    .l
    I (V / R): (lacking R)
    P (V * I): (lacking R)
    V (240): 240
    100R!  # try a value
    .l
    I (V / R): 2.4
    P (V * I): 576.0
    R (100): 100
    V (240): 240
    P 200R!  # let’s add P to the stack so we don't have to .l
    P=288.0
    400R!
    P=144.0
    600R!
    P=96.0
    550R!
    P=104.72727272727272
    575R!
    P=100.17391304347825

We can even do a little bit of manual algebraic manipulation with @:

    @
    V * I=100.17391304347825
    @
    V=240 I=0.41739130434782606
    @
    V=240 V / R=0.41739130434782606
    *
    V * (V / R)=100.17391304347825

At this point we could probably go for the direct solution:

    V V * 100 / R!
    V * (V / R)=100.0
    ?R
    R (V * V / 100): 576.0

There’s also a “goal seek” operator “<>” which can be used to solve the
same problem as follows:

    Calculator.  Type ? for help.
    V R/I!  V I*P!  240V!  100R!  P 100 - R <>
    576.0=576.0
    R!P
    P=100.0
    I
    P=100.0 I=0.4166666666666667

This "<>" operator, applied to P - 100 and R, manipulates R to get P -
100 to be zero.  It works stochastically using the method of secants,
so it may not always succeed, but it succeeds surprisingly often.
Consider the cubic (x + 1)(x - 2)(x + 3), which works out to x³ + x² -
2x² + 3x² - 6x + 3x - 2x - 6 = x³ + 2x² - 5x - 6.  It’s capable of
finding a solution to this cubic; it happens to find -3 if started
from 0:

    0x! x 3**2 x 2** * + 5 x * - 6-
    x ** 3 + 2 * x ** 2 - 5 * x - 6=-6
    y!
    y x <>
    -3.0=-3.0

But given the starting point of 1, it randomly finds all three
solutions:

    1x! y x <> y x <> y x <> y x <> y x <> y x <> y x <> y x <> 
    -3.0=-3.0 -1.0=-1.0 2.0=2.0 2.0=2.0 2.0=2.0 -1.0=-1.0 2.0=2.0 2.0=2.0

If we modify the cubic to use + 6 instead of - 6, it can still solve
it on the real line:

    y@
    -3.0=-3.0 -1.0=-1.0 2.0=2.0 2.0=2.0 2.0=2.0 -1.0=-1.0 2.0=2.0 2.0=2.0
    x ** 3 + 2 * x ** 2 - 5 * x - 6=-8
    @+
    -3.0=-3.0 -1.0=-1.0 2.0=2.0 2.0=2.0 2.0=2.0 -1.0=-1.0 2.0=2.0 2.0=2.0
    x ** 3 + 2 * x ** 2 - 5 * x + 6=4
    y!
    -3.0=-3.0 -1.0=-1.0 2.0=2.0 2.0=2.0 2.0=2.0 -1.0=-1.0 2.0=2.0 2.0=2.0
    ;;;;;;;;;
            ↑ lacking an operand
    y x<>
    -3.7563213575867147=-3.7563213575867147
    x!y
    y=3.552713678800501e-15

It can even find a cube root of -3 if started with an imaginary
initial guess:

    x 3** 3+
    y=3.552713678800501e-15 x ** 3 + 3=-50.00150707085777
    y!
    y=-50.00150707085777
    1jx!
    y=(3-1j)
    y x<>
    y=(3-1j)
    (0.7211247851537043-1.2490247664834064j)=(0.7211247851537043-1.2490247664834064j)
    x!
    y=0j

"""
from __future__ import division, print_function
import operator, re, collections, readline, traceback, pprint, random
namedtuple = collections.namedtuple

try:
    raw_input
except NameError:
    raw_input = input

if __name__ == '__main__':
    # enable cgitb for embedded unit tests
    import cgitb
    cgitb.enable(format='text')

tokenizer = re.compile(r'''
      (?P<num> \d+(?:\.\d*)? | \.\d+)    \s*
    | (?P<unop> [¬~_j;@] ) \s*
    | (?P<op> \*\* | <> | [-:^&|+/*×÷↑%!] ) \s*
    | (?P<id> \w+) \s*
    | (?P<comment> \#.*)
    | (?P<err> .)
''', re.VERBOSE)

def tokenize(line):
    # Following the example in https://docs.python.org/3/library/re.html
    stripped = line.lstrip()
    eaten = len(line) - len(stripped)
    for mo in tokenizer.finditer(stripped):
        yield (mo.lastgroup,
               mo.group(mo.lastgroup),
               mo.start(mo.lastgroup) + eaten)

def ok(a, b):
    assert a == b, (a, b)

ok(list(tokenize(' 3 .4 + 5.1*2.÷x10/ \n')),
   [('num', '3', 1),
    ('num', '.4', 3),
    ('op', '+', 6),
    ('num', '5.1', 8),
    ('op', '*', 11),
    ('num', '2.', 12),
    ('op', '÷', 14),
    ('id', 'x10', 15),
    ('op', '/', 18)])

class Expr:
    pass

class Const(namedtuple('Const', ('val,')), Expr):
    def eval(self, vars):
        return self.val

    def show(self, prec_l, prec_r):
        return str(self.val)

    def decompose(self, vars):
        return []

ops = {
    '**': operator.pow,
    '^': operator.pow,
    '↑': operator.pow,
    '+': operator.add,
    '-': operator.sub,
    '/': operator.truediv,
    '÷': operator.truediv,
    '*': operator.mul,
    '×': operator.mul,
    '&': operator.and_,
    '|': operator.or_,
    '%': operator.mod,
    # XXX min, max, relu, gcd? conditional? log, exp?
}

precedence = {'|': 1, '&': 2, '+': 3, '-': 3,
              '*': 4, '×': 4, '/': 4, '÷': 4, '%': 4,
              '**': 5, '^': 5, '↑': 5,
}

class Binop(namedtuple('Binop', ('op', 'a', 'b')), Expr):
    def eval(self, vars):
        return ops[self.op](self.a.eval(vars), self.b.eval(vars))

    def show(self, prec_l, prec_r):
        p = precedence[self.op]
        parenthesize = not (prec_l < p >= prec_r)
        a = self.a.show(0 if parenthesize else prec_l, p)
        b = self.b.show(p, 0 if parenthesize else prec_r)
        s = '%s %s %s' % (a, self.op, b)
        return '(%s)' % s if parenthesize else s

    def decompose(self, vars):
        return [self.a, self.b]

unops = {
    '¬': operator.invert,
    '~': operator.invert,
    '_': operator.neg,
    'j': lambda num: num * 1j,
}

class Unop(namedtuple('Unop', ('rator', 'rand')), Expr):
    def eval(self, vars):
        return unops[self.rator](self.rand.eval(vars))

    def show(self, prec_l, prec_r):
        rand = self.rand.show(100, 100)
        if self.rator in '¬~':
            return self.rator + rand
        if self.rator == '_':
            return '-' + rand
        if self.rator == 'j':
            return '(%s × j)' % rand

    def decompose(self, vars):
        return [self.rand]

class Var(namedtuple('Var', ('name',)), Expr):
    def eval(self, vars):
        return vars[self.name].eval(vars)

    def show(self, prec_l, prec_r):
        return self.name

    def decompose(self, vars):
        return [vars[self.name]]

class CantRun(namedtuple('CantRun', ('pos', 'why'))):
    pass

def run(tokens, stack, vars):
    for ttype, val, col in tokens:
        if ttype == 'num':
            stack.append(Const(float(val) if '.' in val else int(val)))

        elif ttype == 'unop':
            try:
                rand = stack.pop()
            except IndexError:
                return CantRun(col, 'lacking an operand')

            if val == '@':  # decompose variables into definition, etc.
                stack.extend(rand.decompose(vars))
            elif val != ';':  # drop
                stack.append(Unop(val, rand))

        elif ttype == 'op':
            operands = []
            try:
                b = stack.pop()
                operands.append(b)
                a = stack.pop()
            except IndexError:
                stack.extend(operands)
                return CantRun(col, 'lacking operands')

            if val == '!':
                if not hasattr(b, 'name'):
                    return CantRun(col, '%s isn’t a variable ' % b.show(0, 0))
                vars[b.name] = a
            elif val == '<>':  # goal seek
                stack.append(Const(find_zero(b, a, vars)))
            else:
                stack.append(Binop(val, a, b))

        elif ttype == 'id':
            stack.append(Var(val))

        elif ttype == 'err':
            return CantRun(col, 'unknown token')

        elif ttype == 'comment':
            pass

        else:
            return CantRun(col, 'unknown token type')

def test_vars():
    vars = {}
    stack = []
    run(tokenize('3a!4b!5a*6b*+c!'), stack, vars)
    for tree in stack:
        tree.eval(vars)
    ok(vars['c'].eval(vars), 15+24)

test_vars()

def find_zero(design_var, objective, vars):
    "Try to find a zero using the method of secants."
    var = design_var.name
    original = vars[var]
    original_val = original.eval(vars)
    original_obj = objective.eval(vars)

    try:
        best, best_obj = None, None
        for i in range(100):
            a = original_val
            a_obj = original_obj
            b = (random.random() * 4 - 2) * a + 1

            for j in range(100):
                vars[var] = Const(b)
                b_obj = objective.eval(vars)
                if b_obj == 0:
                    return b

                # find slope of objective function between two points
                if a_obj == b_obj or a == b:
                    break
                slope = (b_obj - a_obj) / (b - a)
                new_value = b - b_obj / slope
                #print(new_value)

                a, b, a_obj = b, new_value, b_obj

            if best is None or abs(a_obj) < abs(best_obj):
                best, best_obj = a, a_obj

        return best

    finally:
        vars[var] = original

def repl():
    print("Calculator.  Type ? for help.")
    vars = {}
    stack = []

    while True:
        prompt = ''
        try:
            line = raw_input(prompt).strip()
        except EOFError:
            return
        except KeyboardInterrupt:
            print("^C")
            continue
        
        try:
            if line == '.l':
                for var in sorted(vars):
                    ast = vars[var]
                    print('%s (%s):' % (var, ast.show(0, 0)), try_eval(ast, vars))
                continue

            elif line.startswith('?'):
                var = line[1:].strip()
                if not var:
                    help()
                    continue

                ast = vars.get(var)
                if not ast:
                    print('no var %r' % var)
                    continue

                print('%s (%s):' % (var, ast.show(0, 0)), try_eval(ast, vars))
                continue

            elif line.startswith('.r'):
                var = line[2:].strip()
                try:
                    del vars[var]
                except KeyError:
                    print('no var %r' % var)
                continue

            err = run(tokenize(line), stack, vars)
            if err:
                print(' ' * (len(prompt) + err.pos) + '↑', err.why)

            col = 0
            for tree in stack:
                # too bad pprint doesn’t work usefully on NamedTuples
                #pprint.pprint(tree)
                s = '%s=%s' % (tree.show(0, 0), try_eval(tree, vars))
                spc = ' ' if col > 0 else ''
                if col + len(spc) + len(s) < 80:
                    print(spc + s, end='')
                    col += len(spc + s)
                else:
                    print()
                    print(s, end='')
                    col = len(s)

            if col > 0:
                print()

        except Exception:
            traceback.print_exc()

def help():
    print("""Commands:
.l — list definitions of existing variables
?foo — list definition of variable 'foo'
.r foo — remove variable 'foo'
^D — quit

Anything else is interpreted as an RPN numerical expression.

Operators supported:
    **, ^, or ↑: exponentiation  %: modulo               _: negation (3_ is -3)
    +: addition                  * or ×: multiplication  ¬ or ~: bitwise NOT
    -: subtraction               &: bitwise AND          ;: discard a value
    / or ÷: division             |: bitwise OR
    j: imaginary (3j is 3 i)     @: decompose top of stack
    <>: find a zero; x 100 - y <> tries to find the y where x - 100 = 0

To store a formula into a variable “foo”, use “foo!”.  This will
change the value of any variable or calculation that depends
on foo.  “foo” can’t begin with “j” or “_”.

For example, to set d to √(b² - 4ac), you could say b 2^4a×
c×-.5^d!, or to set a=1, b=0, c=-4, 1a!0b!4_c!.
""")

def try_eval(ast, vars):
    try:
        return str(ast.eval(vars))
    except KeyError as e:
        return '(lacking %s)' % e.args
    except:
        print("error evaling", repr(ast))
        raise

if __name__ == '__main__':
    repl()