#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# @Author: Gillett Hernandez
# @Date:   2016-07-14 17:06:40
# @Last Modified by:   Gillett Hernandez
# @Last Modified time: 2017-12-16 14:25:17

# from math import log
import time
from itertools import chain, count
import sys
import os
import fractions
import decimal
# from typing import *
if sys.hexversion > 0x03000000:
    xrange = range

from operator import mul
from functools import reduce
import math

prod = lambda L: reduce(mul, L, 1)

def timed(f):
    def inner(*args, **kwargs):
        c1 = time.clock()
        res = f(*args, **kwargs)
        c2 = time.clock()
        print(c2-c1)
        return res
    inner.__name__ = f.__name__
    inner.__doc__ = f.__doc__
    return inner


def digit_sum(n):
    return sum(int(c) for c in str(n))


def phi(n):
    if n % 2 == 0:
        if (n//2) % 2 == 0:
            return 2 * phi(n//2)
        else:
            return phi(n//2)
    from fractions import gcd
    return sum([1]+[1 for i in range(2, n+1) if gcd(i, n) == 1])


def phi_from_pf(pf):
    n = 1
    for p, k in pf.items():
        # n *= p**k
        # but phi(n) = n * prod(1-1/pi))
        n *= p**k - p**(k-1)
    return n


def prime_sieve(limit=0):
    starting = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]
    if limit < starting[-1]:
        return [p for p in starting if p < limit]
    tests = set(range(starting[-1]+2, limit, 2))
    if limit != 0:
        for p in starting:
            for i in range(p, 1+limit//p):
                tests.discard(p*i)
        slimit = 1+int(limit**0.5)
        for pc in range(starting[-1]+2, slimit, 2):
            if pc not in tests:
                continue
            for i in range(pc, 1+limit//pc):
                tests.discard(i*pc)
    R = list(tests) # type: List[int]
    R.sort()
    return starting + R

def primes(less_than):
    n=less_than
    generator = lambda: filter(lambda e: e % 5 != 0, chain([2], range(3, n, 2)))
    S = {i for i in generator()}
    for pc in generator():
        if pc in S:
            for i in range(pc, 1+n//pc):
                S.discard(i*pc)
    return list(sorted(list(S)))


def div(a, b):
    """divides a/b repeatedly, returning the number of times it could be divided and the leftover amount
>>> div(33, 3)
(1, 11)"""
    c = 0
    while a % b == 0:
        a //= b
        c += 1
    return c, a


def list_to_limit(L, x):
    for el in L:
        if el > x:
            break
        else:
            yield el


def get_primes(limit=None, number=None, filepath="../primes_to_10000000.txt"):
    """gets *limit* primes or the primes from 2 to number, whichever comes first\nhard limit of primes up to 100_000_000, about 2_700_000 primes"""
    if limit is None and number is None:
        raise RuntimeError("Invalid Call")

    if limit is not None:
        assert isinstance(limit, int)
    if number is not None:
        assert isinstance(number, int)
    if limit is not None and number is not None:
        pass

    if not os.path.exists(filepath):
        return None

    primes = []
    c = 0
    with open(filepath, 'r') as fd:
        for line in fd:
            line = line.strip()
            if line == '':
                continue
            p = int(line)
            if limit is not None and p > limit:
                break
            c += 1
            if number is not None and c > number:
                break
            primes.append(p)
        else:
            print("Warning: only primes below {} are obtainable by thing function. function called with limit={}".format("10_000_000", limit))
            return primes
    # if the loop goes to completion, then there werent enough primes
    if number is not None and number - c > 0:
        # idk not enough primes read but they were read.
        raise RuntimeError("somehow not enough primes were read")
        # another problem, this needs to be a generator or something.

    return primes


# print(len(get_primes(number=100)))


def prime_factorization(n, primes=[]):
    if primes != []:
        return pf_with_primes_given(n, primes)
    factors = {}
    LIMIT = int(n**0.5) + 1
    for p in get_primes(limit=LIMIT):
        k = 0
        if n % p == 0:
            k, n = div(n, p)
            factors[p] = k
        if n == 1:
            break
    if n != 1:
        factors[n] = 1
    return factors


def pf_with_primes_given(n, primes):
    factors = {}
    LIMIT = int(n**0.5) + 1
    assert primes[-1] >= LIMIT, (primes[-1], LIMIT, n)
    for p in primes:
        k = 0
        if n % p == 0:
            k, n = div(n, p)
            factors[p] = k
        if n == 1:
            break
    if n != 1:
        factors[n] = 1
    return factors


def n_from_factorization(pf):
    P = 1
    for k, v in pf.items():
        P *= k**v
    return P


def sum_of_divisors(n, func=lambda x: x):
    s = 1
    sqrtn = 1+int(n**(0.5))
    for i in range(2, sqrtn):
        if n % i == 0 and n // i != i:
            s += func(i)
            s += func(n//i)
        elif n % i == 0 and n // i == i:
            s += func(i)
    return s


def divisor_list(n):
    ds = []
    def mask(d):
        nonlocal ds
        ds.append(d)
        return 1
    # sum_of_divisors(n, lambda d: 1 if ds.append(d) is None else 1)
    sum_of_divisors(n, mask)
    ds.sort()
    return ds


def divisor_count(n):
    """extremely similar to divisors(n) but this one does NOT keep track of found divisors"""
    ct = sum_of_divisors(n, lambda d: 1)
    return ct


def __divisor_list_pf(pf, prev=1, limit=float("inf")):
    """works recursively to yield all divisors of n=prod(p**f for p,f in pf.items())"""
    p = next(iter(pf))
    f = pf[p]
    for i in range(0,f+1):
        # print(p,i)
        new = prev * p**i
        if i != 0:
            if new > limit:
                break
            yield new
        copy = {}
        copy.update(pf)
        del copy[p]
        yield from __divisor_list_pf(copy, prev=new, limit=limit)


def divisor_list_pf(pf, limit=None):
    yield 1
    if limit is not None:
        # limit = limited
        limit = int(n_from_factorization(pf)**0.5) if limit == True else limit
        assert isinstance(limit, int)
    yield from __divisor_list_pf(pf, limit=limit)

def divisor_pairs_pf(pf):
    n = n_from_factorization(pf)
    for d in divisor_list_pf(pf, limit=int(n**0.5)):
        yield (d, n//d)


def divisor_pairs(n, primes=None):
    if primes is not None:
        pf = prime_factorization(n, primes)
        yield from divisor_pairs_pf(pf)
    else:
        for d in divisor_list(n):
            if d > n//d:
                break
            yield (d, n//d)


def divisor_0(pf):
    """divisor function with k = 0 based off of factorization"""
    P = 1
    for p,f in pf.items():
        P *= (f+1)
    return P


def divisor_1(pf):
    """divisor function with k = 1"""
    P = 1
    for p, f in pf.items():
        P *= (pow(p, f+1) - 1) // (p-1)
    return P


def is_prime(n, primes=[]):
    if primes != []:
        return is_prime_w_primes_given(n, primes)
    """not fast because there is no caching in the get_primes call"""
    LIMIT = int(n**0.5) + 1
    if n == 1:
        return False
    try:
        primes = get_primes(limit=LIMIT)
    except FileNotFoundError:
        primes = prime_sieve(limit)
    for p in get_primes(limit=LIMIT):
        if n % p == 0:
            return False
    return True


def is_prime_w_primes_given(n, primes):
    LIMIT = int(n**0.5) + 1
    assert primes[-1] >= LIMIT, (n, primes[-1], LIMIT)
    if n == 1:
        return False
    for p in primes:
        if p > LIMIT:
            return True
        if n == p:
            return True
        try:
            if n % p == 0:
                return False
        except ZeroDivisionError:
            print(locals())
            raise
    return True


# print(is_prime(1231237))


def is_permutation(a, b):
    return sorted(str(a)) == sorted(str(b))


def is_permutation2(a, b):
    return (a % 9) == (b % 9)


def lag(L):
    assert isinstance(L, list)
    iL = iter(L)
    v1, v2 = next(iL), next(iL)
    yield (v1, v2)
    for v in iL:
        v1, v2 = v2, v
        yield (v1, v2)


def replace_repeated(S, ruleset, N=0):
    i = 0
    while i < N or N == 0:
        Sprev = S
        for match, replacement in ruleset.items():
            S = S.replace(match, replacement)
        if Sprev == S:
            return S
        i += 1


def tee(*args):
    for arg in args:
        print(arg,)
    print('')
    return tuple(args)


def modified_prime_sieve(n, lowerlimit):
    """a modified prime sieve that works off of primes already generated"""
    NLIMIT = lowerlimit  # 10 million
    assert n < NLIMIT**2

    N = set(range(NLIMIT+1, n, 2))
    # work by having the multiples be ranged and calculated to be in the N range
    primes = get_primes(limit=NLIMIT)
    for p in primes:
        lrange = -1+NLIMIT//p
        if lrange < 2:
            lrange = 2
        urange = 1+n//p
        for i in range(lrange, urange):
            N.discard(i*p)
    return list(sorted(N))

bind = lambda instance, func, asname=None: setattr(instance, asname or func.__name__, func.__get__(instance))
testable = (lambda f: (lambda g: bind(f, g)))


def large_prime_sieve(n, sentinel=False):
    """uses modified_prime_sieve"""
    return []
    limit = l = int(n**0.5+1)
    diff = n-limit
    factor = n//limit
    chunk = diff//factor
    print(chunk+l, l, chunk)
    print(factor*chunk+limit)
    actually_all_primes_below_n = None
    if not sentinel:
        actually_all_primes_below_n = []
        print("Warning: this function can output to file instead of saving to memory. Pass sentinel=True to confirm modifying 'primes_to_{limit}.txt'".format(limit=n))
        def write_to_file(chunk):
            nonlocal actually_all_primes_below_n
            actually_all_primes_below_n.extend(chunk)
    else:
        def write_to_file(chunk):
            with open("primes_to_{limit}.txt".format(limit=n), 'a') as fd:
                for p in chunk:
                    fd.write('\n'+str(p))
    for i in range(1, factor+1):
        # print(i)
        # break
        l_chunk = modified_prime_sieve(l + chunk, l)
        l = i*chunk+limit
        write_to_file(l_chunk)
        # print("successfully wrote")
    return actually_all_primes_below_n

def basic_large_prime_sieve(n, debug=False):
    """generates primes up to n using bytearrays instead of a set object or normal array\nworks up to n = <free memory size in bytes> * 8\nfor instance at 100MB memory usage, n could be ~838 million"""
    limit = l = int(n**0.5+1)
    allprimes = []
    candidates = bytearray(1+l//8)
    primes = []
    startprimes = [2,3,5,7,11,13,17,19,23]
    primes.extend(startprimes)
    for p in startprimes:
        for i in range(p, l//p):
            m = i*p
            offset, idx = divmod(m, 8)
            # bmask sets the bit cooresponding to m to 1, marking it as composite
            bmask = 1 << idx
            # print(bmask)
            candidates[offset] |= bmask
    for offset, byte in enumerate(candidates):
        for idx in range(8):
            # if the bit at 8*offset+idx is 0
            # then 8*offset + idx is prime
            if ((byte>>idx) & 1) ^ 1:
                p = (offset << 3) + idx
                if p == 0 or p == 1:
                    continue
                primes.append(p)
                for i in range(p, l//p):
                    m = i*p
                    g_offset, g_idx = divmod(m, 8)
                    # bmask sets the bit cooresponding to m to 1, marking it as composite
                    bmask = 1 << g_idx
                    candidates[g_offset] |= bmask
    #print(primes)
    # and then generate all the rest using a sliding window type iteration
    # need to use offset to coorespondingly adjust primes and bmasks and stuff
    # for offset in range(l, n+l+1, l):
        # offset is the lower bound of the next frame of computation
    # for i in range(len(candidates)):
        # candidates[i] = 0
    candidates = bytearray(1+n//8)
    # reset candidates to all 0s
    # marks all composites by fixing p then multiplying by i in [p, n//p]
    for p in primes:
        for i in range(p, 1+n//p):
            m = i*p
            # print(m)
            # generate a composite
            if m > n:
                break
            # g_offset and idx are where in the n-space the numbers are
            g_offset, idx = divmod(m, 8)
            # print( g_offset, idx)
            # subtract offset
            # g_offset -= offset//8
            # bmask sets the bit cooresponding to m to 1, marking it as composite
            bmask = 1 << idx
            # print(bmask)
            candidates[g_offset] |= bmask
    for g_offset, byte in enumerate(candidates):
        for idx in range(8):
            # if the bit at 8*offset+idx is 0
            # then 8*offset + idx is prime
            # `<<3` <=> `*8`
            if (g_offset << 3) + idx > n:
                break
            if ((byte>>idx) & 1) ^ 1:
                p = (g_offset << 3) + idx
                if p == 0 or p == 1:
                    # 0 and 1 will not be marked composite because they arent multiples of any two primes
                    # so just skip them
                    if debug: print(g_offset, idx, bin(byte)[2:], bin(candidates[g_offset+1])[2:])
                    continue
                allprimes.append(p)
    # print(candidates)
    return allprimes

def really_large_prime_sieve(n, debug=False):
    """uses much less memory when calculating primes than large_prime_sieve.\n at n=1_000_000_000, memory usage to actually generate them is sqrt(n)//8 bytes, or 3.9kb,\n while memory usage to keep them all in memory is ~n/log(n)*8, or 368MB at n=1_000_000_000"""
    # principles:
    # A solution to these problems is offered by segmented sieves, where only portions of the range are sieved at a time. These have been known since the 1970s, and work as follows:
    # Divide the range 2 through n into segments of some size Δ ≤ √n.
    # Find the primes up to Δ, using the "regular" sieve.
    # For each Δ-sized block from √n + 1 to n, set up a Boolean array of size Δ. Eliminate the multiples of each prime p ≤ √n found in step 2.
    # If Δ is chosen to be √n, the space complexity of the algorithm is O(√n), while the time complexity is the same as that of the regular sieve.
    # multiples of the form p*2, p*3, etc do not need to be eliminated as they would already be eliminated when eliminating multiples of 2, 3, 5, etc
    # only multiples of the form p*p, p*(p+1),... need to be eliminated
    limit = l = 1+int(n**0.5)
    # print(l)
    allprimes = []
    # note: candidates may be slightly oversized
    candidates = bytearray(1+l//8)
    primes = []
    startprimes = [2,3,5,7,11,13,17,19,23]
    # primes.extend(startprimes)
    for p in startprimes:
        for i in range(p, 1+l//p):
            m = i*p
            offset, idx = divmod(m, 8)
            # bmask sets the bit cooresponding to m to 1, marking it as composite
            bmask = 1 << idx
            # print(bmask)
            candidates[offset] |= bmask
    for offset, byte in enumerate(candidates):
        for idx in range(8):
            # if the bit at 8*offset+idx is 0
            # then 8*offset + idx is prime
            if ((byte>>idx) & 1) ^ 1:
                p = (offset << 3) + idx
                if p == 0 or p == 1:
                    continue
                if p > l:
                    break
                primes.append(p)
                for i in range(p, l//p):
                    m = i*p
                    g_offset, g_idx = divmod(m, 8)
                    # bmask sets the bit cooresponding to m to 1, marking it as composite
                    bmask = 1 << g_idx
                    candidates[g_offset] |= bmask
    allprimes.extend(primes)

    # and then generate all the rest using a sliding window type iteration
    # need to use offset to coorespondingly adjust primes and bmasks and stuff
    # for offset in range(l, n+l+1, l):
        # offset is the lower bound of the next frame of computation
    for i in range(1, l):
        # gehhh
        for _ in range(len(candidates)):
            candidates[_] = 0
        lower = i*l
        upper = i*l + l
        # print("lower limit is {}, upper is {}".format(lower,upper))
        # reset candidates to all 0s
        # marks all composites by fixing p then multiplying by i in [lower//p, upper//p]
        # l_offset and l_idx are where in the n space we need to move to
        l_offset, l_idx = divmod(lower, 8)
        if debug: print(lower, upper)
        for p in primes:
            if debug: print(p*(lower//p-1), p*(1+upper//p))
            for k in range(lower//p-1, 2+upper//p):
                if k <= 1:
                    continue
                m = k*p
                if m < lower:
                    continue
                if m > upper:
                    break
                # print(m)
                # generate a composite
                if m > n:
                    break
                # g_offset and idx are where in the n-space the numbers are
                g_offset, idx = divmod(m, 8)
                # subtract offset because g_offset and idx are only settable if the bytearray spans the entirety of [0..n]
                # but the bytearray only goes to [0..sqrt(n)] and offsets need to be added accordingly
                g_offset -= l_offset
                idx -= l_idx
                if idx < 0:
                    g_offset -= 1
                    idx += 8
                # bmask sets the bit cooresponding to m to 1, marking it as composite
                bmask = 1 << idx
                # print(bmask)
                candidates[g_offset] |= bmask
                # print("good", locals())
        for g_offset, byte in enumerate(candidates):
            for idx in range(8):
                # if the bit at 8*offset+idx is 0
                # then 8*offset + idx is prime
                p = (g_offset << 3) + idx + lower
                if p > upper or p > n:
                    break
                if ((byte>>idx) & 1) ^ 1:
                    # print("freaking oops", locals())
                    if p == 0 or p == 1:
                        continue
                    allprimes.append(p)
        # print(candidates)
    return allprimes

def tlps_gen(n):
    # principles:
    # A solution to these problems is offered by segmented sieves, where only portions of the range are sieved at a time. These have been known since the 1970s, and work as follows:
    # Divide the range 2 through n into segments of some size Δ ≤ √n.
    # Find the primes up to Δ, using the "regular" sieve.
    # For each Δ-sized block from √n + 1 to n, set up a Boolean array of size Δ. Eliminate the multiples of each prime p ≤ √n found in step 2.
    # If Δ is chosen to be √n, the space complexity of the algorithm is O(√n), while the time complexity is the same as that of the regular sieve.
    # multiples of the form p*2, p*3, etc do not need to be eliminated as they would already be eliminated when eliminating multiples of 2, 3, 5, etc
    limit = l = int(n**0.5+1)
    allprimes = []
    candidates = bytearray(1+l//8)
    if True:
    # if l < 10000000:
        # primes = get_primes(limit=l)
    # if primes is []:
        # if its unfeasable, memory-wise, to generate the primes using the usual method
        # generate them using the 8x less memory method
        primes = []
        startprimes = [2,3,5,7,11,13,17,19,23]
        primes.extend(startprimes)
        for p in startprimes:
            for i in range(p, l//p):
                m = i*p
                offset, idx = divmod(m, 8)
                # bmask sets the bit cooresponding to m to 1, marking it as composite
                bmask = 1 << idx
                # print(bmask)
                candidates[offset] |= bmask
        for offset, byte in enumerate(candidates):
            for idx in range(8):
                # if the bit at 8*offset+idx is 0
                # then 8*offset + idx is prime
                if ((byte>>idx) & 1) ^ 1:
                    p = (offset << 3) + idx
                    if p == 0 or p == 1:
                        continue
                    primes.append(p)
                    for i in range(p, l//p):
                        m = i*p
                        g_offset, g_idx = divmod(m, 8)
                        # bmask sets the bit cooresponding to m to 1, marking it as composite
                        bmask = 1 << g_idx
                        candidates[g_offset] |= bmask
        for p in primes:
            yield p
    #print(primes)
    # and then generate all the rest using a sliding window type iteration
    # need to use offset to coorespondingly adjust primes and bmasks and stuff
    # for offset in range(l, n+l+1, l):
        # offset is the lower bound of the next frame of computation
    # for i in range(len(candidates)):
        # candidates[i] = 0
    candidates = bytearray(1+n//8)
    # reset candidates to all 0s
    raise NotImplementedError
    for p in primes:
        for i in range(p, 1+n//p):
            m = i*p
            # print(m)
            # generate a composite
            if m > n:
                break
            # g_offset and idx are where in the n-space the numbers are
            g_offset, idx = divmod(m, 8)
            # print( g_offset, idx)
            # subtract offset
            # g_offset -= offset//8
            # bmask sets the bit cooresponding to m to 1, marking it as composite
            bmask = 1 << idx
            # print(bmask)
            candidates[g_offset] |= bmask
    for g_offset, byte in enumerate(candidates):
        for idx in range(8):
            # if the bit at 8*offset+idx is 0
            # then 8*offset + idx is prime
            if (g_offset << 3) + idx > n:
                break
            if ((byte>>idx) & 1) ^ 1:
                p = (g_offset << 3) + idx
                yield p
    # print(candidates)


def matmul(A, B):
    # A = n x m matrix
    # B = m x p matrix
    # len(A) would be its # of rows
    # len(B[0]) would be its # of columns
    if not isinstance(A[0], (list, tuple)):
        A = [[v] for v in A]
    lA = len(A[0])
    lB = len(B)
    m = len(A)
    assert m == len(B[0])
    new = [[0 for _ in range(lA)] for _ in range(lB)]
    # print(new)
    for i in range(lA):
        for j in range(lB):
            new[j][i] = sum(A[k][i]*B[j][k] for k in range(m))
    return new

def flatten(L):
    l = []
    for sub in L:
        if isinstance(sub, (list, tuple)):
            l.extend(sub)
        else:
            l.append(sub)
    return l

A = [[1, -2, 2],
     [2, -1, 2],
     [2, -2, 3]]


B = [[1, 2, 2],
     [2, 1, 2],
     [2, 2, 3]]

C = [[-1, 2, 2],
     [-2, 1, 2],
     [-2, 2, 3]]


def primitive_generator_depth_first_w_condition(root=[3,4,5], cond=lambda a,b,c: sum([a,b,c]) < 1500000):
    """starting at a root in the form [[a], [b], [c]], and a condition with form lambda g: sum(v[0] for v in g) < 150000\ntraverse depth first the tree of primitive pythagorean triples"""
    if cond(*root) is False:
        return
    # assert root[0]**2 + root[1]**2 == root[2]**2
    yield root
    # print(root)
    Am = flatten(matmul(root, A))
    # yield from primitive_generator_depth_first_w_condition(Am, cond=cond)
    for rt in primitive_generator_depth_first_w_condition(Am, cond=cond):
        yield rt
    Bm = flatten(matmul(root, B))
    # yield from primitive_generator_depth_first_w_condition(Bm, cond=cond)
    for rt in primitive_generator_depth_first_w_condition(Bm, cond=cond):
        yield rt
    Cm = flatten(matmul(root, C))
    # yield from primitive_generator_depth_first_w_condition(Cm, cond=cond)
    for rt in primitive_generator_depth_first_w_condition(Cm, cond=cond):
        yield rt

def R(I, rule="2*a[i-1]+i"):
    offset = I[0][0]
    a = [] # type: List[int]
    for i in range(offset-1):
        a.append(None)
    for k, v in I:
        a.append(v)
    i = offset + 1
    while True:
        v = eval(rule)
        a.append(v)
        i += 1
        yield v

def lag(iterable):
    iterator = iter(iterable)
    last = next(iterator)
    for e in iterator:
        yield last, e
        last = e


def egcd(a,b):
    s, o_s = 0, 1
    t, o_t = 1, 0
    r, o_r = b, a
    while r != 0:
        q = o_r // r
        o_r, r = r, o_r - q*r
        o_s, s = s, o_s - q*s
        o_t, t = t, o_t - q*t
    return (o_s, o_t), o_r, (t, s)


def to_base(base):
    def f(n):
        digits = []
        while (n >= 1):
            digits.append(n % base)
            n //= base
        return digits[::-1]
    return f


def isint(n):
    return int(n)==n


def is_square(n):
    return int(n**0.5)**2 == n


def lagged_fibonacci_gen():
    sk = []
    l1 = lambda k: (100003 - 200003*k + 300007*k**3) % 1000000 - 500000
    l2 = lambda sk, k: (sk[-24] + sk[-55] + 1000000) % 1000000 - 500000
    for k in range(1,55+1):
        n = l1(k)
        sk.append(n)
        # yield n
    for k in range(56, 4000000+1):
        n = l2(sk, k)
        # for i in range(len(sk)-1):
        #     sk[i] = sk[i+1]
        # # print(sk)
        # sk[-1] = n
        # yield n
        sk.append(n)
    return sk


class Graph:
    def __init__(self, vertices, edges):
        assert isinstance(vertices, list)
        assert isinstance(edges, dict)
        # self.__G = {u: {v: w for v, w in edges[u]} for u in vertices}
        self.__G = {u: {v: edges[(u, v)] for v in vertices if (u, v) in edges} for u in vertices}
        self.vertices = vertices

    def edge_weight(self, u, v):
        return self.__G[u][v]

    def path_weight(self, path):
        return sum(self.edge_weight(u, v) for u, v in path)

    def shortest_path_weight(self, u, v):
        if not self.connected(u, v):
            return float("inf")
        raise NotImplementedError

    def initialize_single_source(self, s):
        self.d = {}
        self.p = {}
        for v in self.__G.keys():
            self.d[v] = float("inf")
            self.p[v] = None
        self.d[s] = 0

    def relax(self, u, v):
        w = self.edge_weight(u, v)
        if self.d[v] > self.d[u] + w:
            self.d[v] = self.d[u] + w
            self.p[v] = u

    def dijkstra(self):
        S = set([])
        # Q = PriorityQueue()
        # for v in self.__G.keys():
            # Q[v] = self.d[v]
        Q = list(sorted(list(self.__G.keys()), key=lambda el: self.d[el]))
        while len(Q) != 0:
            # u = Q.smallest()
            u = min(Q, key=lambda el: self.d[el])
            Q.remove(u)
            S |= {u}
            for v in self.__G[u].keys():
                self.relax(u, v)




class ContinuedFraction:
    def __init__(self, start, repeat, debug=False):
        self.start = start
        self._l = WrapList(repeat)
        self.cng = self.convergent_gen

    @classmethod
    def from_(cls, d):
        """returns the continued fraction representation of √d"""
        start, repeat, length = analyze_triplets(d)
        return cls(start, list(repeat))

    def truncate(self,i):
        return [self[k] for k in range(i)]

    def convergent_gen(self, n=None):
        # h_n=a_n*h_(n-1)+h_(n-2)
        # k_n=a_n*k_(n-1)+k_(n-2)
        # n  :-2|-1| 0| 1| 2| 3| 4| 5| 6| 7| 8|
        # a_n:  |  | 2| 1| 2| 1| 1| 4| 1| 1| 6|
        # ------+------------------------------
        # h_n: 0| 1| 2| 3| 8|  |  |  |  |  |  |
        # k_n: 1| 0| 1| 1| 3|  |  |  |  |  |  |
        hs = [0, 1]
        ks = [1, 0]
        i = 0

        an = self[i]


        hs[-2], hs[-1] = hs[-1], an*hs[-1]+hs[-2]
        ks[-2], ks[-1] = ks[-1], an*ks[-1]+ks[-2]

        yield fractions.Fraction(hs[-1], ks[-1])

        i += 1
        while True:
            if n is not None and i >= n:
                return
            an = self[i]
            hs[-2], hs[-1] = hs[-1], an*hs[-1]+hs[-2]
            ks[-2], ks[-1] = ks[-1], an*ks[-1]+ks[-2]
            yield fractions.Fraction(hs[-1], ks[-1])

        # hs.append(an*hs[-1]+hs[-2])
        # ks.append(an*ks[-1]+ks[-2])

        # yield fractions.Fraction(hs[-1], ks[-1])

        # i += 1
        # while True:
        #     if n is not None and i >= n:
        #         raise StopIteration
        #     an = self[i]
        #     hs.append(an*hs[-1]+hs[-2])
        #     ks.append(an*ks[-1]+ks[-2])
        #     yield fractions.Fraction(hs[-1], ks[-1])

            i+=1
        if self.debug:
            print('|'.join([str(elem).rjust(5) for elem in hs]))
            print('|'.join([str(elem).rjust(5) for elem in ks]))

    def convergent(self, n, debug=False):
        # wasteful
        return list(self.cng(n+1))[-1]

    def convergent_list(self, n):
        r = fractions.Fraction(self[n], 1)
        l = [r]
        for elem in reversed(self[1:n-1]):
            r = elem+1/r
            l.append(self.start+1/r)
        return l

    def __getitem__(self, key):
        if type(key) == int:
            if key < 0:
                raise IndexError
            if key == 0:
                return self.start
            else:
                return self._l[key-1]
        elif type(key) == slice:
            if key.step is not None:
                raise NotImplementedError
            return [self[i] for i in range(key.start,key.stop)]


class WrapList:
    def __init__(self,l):
        self.l = l

    def __getitem__(self,key):
        ll = len(self.l)
        if (type(self.l[key%ll]).__name__=='function'):
            return self.l[key%ll](key)
        else:
            return self.l[key%ll]

    def __len__(self):
        return len(self.l)


def sconvergents(x):
    if int(x**0.5)**2 == x:
        return
    m = 0
    d = 1
    a = a0 = math.floor(x**0.5)
    triplets = [(m, d, a)]
    while True:
        m = d*a - m
        d = (x-m**2)/d
        a = math.floor((a0+m)/d)
        if (m, d, a) in triplets:
            return triplets + [(m, d, a)]
        else:
            triplets.append((m, d, a))


def analyze_triplets(x):
    triplets = sconvergents(x)
    if triplets is None:
        return
    last = triplets[-1] # last occurence of the repeated meme
    index = triplets.index(last) # first occurence
    an = list(zip(*triplets))[2]
    return an[:index][0], list(an[index:-1]), len(triplets)-index-1


def error(p,q,target):
    # p,q = fraction.numerator, fraction.denominator
    assert isinstance(target, decimal.Decimal)
    return (p - q*target)/(q)

def construct(L):
    f = fractions.Fraction(1,L[-1])
    for e in reversed(L[:-1]):
        f = 1/(f + e)
    return 1/f



def construct2(L):
    # h_n=a_n*h_(n-1)+h_(n-2)
    # k_n=a_n*k_(n-1)+k_(n-2)
    # n  :-2|-1| 0| 1| 2| 3| 4| 5| 6| 7| 8|
    # a_n:  |  | 2| 1| 2| 1| 1| 4| 1| 1| 6|
    # ------+------------------------------
    # h_n: 0| 1| 2| 3| 8|  |  |  |  |  |  |
    # k_n: 1| 0| 1| 1| 3|  |  |  |  |  |  |
    hs = [0, 1]
    ks = [1, 0]
    i = 0

    an = L[i]

    i = 1


    hs[-2], hs[-1] = hs[-1], an*hs[-1]+hs[-2]
    ks[-2], ks[-1] = ks[-1], an*ks[-1]+ks[-2]

    yield fractions.Fraction(hs[-1], ks[-1])

    while i < len(L):
        an = L[i]
        hs[-2], hs[-1] = hs[-1], an*hs[-1]+hs[-2]
        ks[-2], ks[-1] = ks[-1], an*ks[-1]+ks[-2]
        yield fractions.Fraction(hs[-1], ks[-1])
        i += 1


def best_rational_approximation(target,cf,d):
    # l = []
    best_error = float("inf")
    hs = [0, 1]
    ks = [1, 0]
    best_pair = (0,0)
    hs[-2], hs[-1] = hs[-1], cf[0]*hs[-1]+hs[-2]
    ks[-2], ks[-1] = ks[-1], cf[0]*ks[-1]+ks[-2]
    for i in count(1):
        limit = cf[i]
        # l2 = 1+limit//2
        newh, newk = sum(hs), sum(ks)
        for k in range(1,limit+1):
            # semi = fractions.Fraction(newh,newk)
            if newk > d:
                # print(newh, newk)
                return best_pair

            if k < limit//2:
                newh += hs[-1]
                newk += ks[-1]
                continue

            es = abs(error(newh, newk, target))
            if es < best_error:
                best_error = es
                best_pair = (newh,newk)
                # print(newh, newk)
                # l.append((newh, newk))
            newh += hs[-1]
            newk += ks[-1]
        newh -= hs[-1]
        newk -= ks[-1]
        hs[-2], hs[-1] = hs[-1], newh
        ks[-2], ks[-1] = ks[-1], newk
        # hs[-1] = newh
        # ks[-1] = newk
    return l[-1]


def solves_pell_pos(D, x, y):
    return x**2 - D*y**2 == 1


def solves_pell_neg(D, x, y):
    return x**2 - D*y**2 == -1


def solve_pell(D, filter=solves_pell_pos, terminate=lambda i, x, y, c: False):
    """solves basic pell equation of the form x^2-dy^2=(+-)1 by searching through convergents of D. filter determines whether a particular solution is yielded,\n and terminate is a lambda that takes (i) the index of the current convergent, (x,y) the actual convergent, and c, the count of successful solutions by `filter`'s criteria"""
    cf = ContinuedFraction.from_(D)
    c = 0
    for i in count(1):
        conv = cf.convergent(i)
        x,y = conv.numerator, conv.denominator
        if filter(D, x, y):
            c += 1
            yield (x, y, i)
        if terminate(i,x,y,c):
            return


def solve_pell_like(a,b,c,d,e,f):
    """solves pell-like equation of the form ax^2 + bxy + cy^2 + dx + ey + f = 0"""
    D = b**2 - 4*a*c
    E = b*d - 2*a*e
    F = d**2 - 4*a*f
    N = E**2 - D * F
    print(D,E,F,N)
    # solve X^2 - D Y^2 = N
    # X = Dy + E
    # y = (X - E)/D
    # Y = 2ax + by + d
    # x = (Y - by - d)/2a

    solution_iter = solve_pell(D, filter=lambda D,x,y: x**2 - D*y**2 == N)
    for X, Y, i in solution_iter:
        y = (X - E) / D
        if not isint(y):
            continue
        x = (Y - b*y - d) / (2*a)
        if not isint(x):
            continue
        assert a*x**2 + b*x*y + c*y**2 + d*x + e*y + f == 0
        yield (x,y,i)


def nfib(k):
    last = []
    for i in range(1,k):
        last.append(2**i)
        yield last[-1]
    last.append(2**k-1)
    yield last[-1]
    while True:
        new = sum(last)
        last.pop(0)
        last.append(new)
        yield new



if __name__ == '__main__':
    import os
    nlimit1 = 10000000   # 10 million
    nlimit2 = 100000000 # 100 million
    if not os.path.exists("../primes_to_{}.txt".format(nlimit1)):
        with open("../primes_to_{}.txt".format(nlimit1), "w") as fd:
            print("primes to 10 million")
            primes = tlps1(nlimit1)
            for p in primes:
                fd.write(str(p))
                fd.write("\n") 

    # if not os.path.exists("../primes_to_{}.txt".format(nlimit2)):
    #     with open("../primes_to_{}.txt".format(nlimit1), "w") as fd:
    #         print("primes to 100 million")
    #         primes = tlps1(nlimit2)
    # time_test()
    # for pyg_triple in primitive_generator_depth_first_w_condition(cond=lambda g: sum(v[0] for v in g) < 1000):
    #     print(pyg_triple)
