'''This code block sets the parameters for evaluating'''
P=32 #precision of L-function estimate
crit=9/4 #the real value of the critical line
delta=2 #delta, as given by Yoshida, far enough that crit+delta is in the region of absolute convergence of the series
Lerror=100

'''This code block defines the L function being considered and the associated gamma factors used to complete it'''
L=zeta
Gam(s)=pi^(-s/2)*gamma(s/2)

'''This code just defines some functions that are needed to compute variation of argument'''

def argbranch(z,t):
    '''this gives the branch of arg of z in the interval (t-pi,t+pi), might screw up on an edge case'''
    x=float(arg(z))
    if abs(x-t)<pi:
        return x
    else:
        n=floor(abs(x-t)/(2*pi))
        if x<t:
            return float(x+2*pi*(n+1))
        else:
            return float(x-2*pi*(n+1))

def argvarl(lcalcfile, path, n):
    '''lcalcfile is the name of an lcalc file and should be input as a string in quotation marks.  path is a list of complex points.  This function uses lcalc to compute the variation of the argument of f (divided by pi) on the straight line segments in CC connecting the points in path, estimated using n subdivisions each.  Should be used for the L function part of the variation of argument calculation'''
    
    # Make list of values of L using lcalc.  Format will be list of pairs (s, L(s))
    values=[]
    for i in range(0, len(path)-1):
        rawvaluesstr=lcalc('--value-line-segment -x'+str(path[i].real())+' -y'+str(path[i].imag())+' -X'+str(path[i+1].real())+' -Y'+str(path[i+1].imag())+' --number-samples='+str(n)+' -F'+str(lcalcfile)+'.txt')
        rawvalues0=rawvaluesstr.split()
        for j in range(len(rawvalues0)/4):
            values.append((float(rawvalues0[4*j])+I*float(rawvalues0[4*j+1]), float(rawvalues0[4*j+2])+I*float(rawvalues0[4*j+3])))
    lastvals=lcalc('-v -x'+str(path[len(path)-1].real())+' -y'+str(path[len(path)-1].imag())+' -F'+str(lcalcfile)+'.txt').split()
    values.append((path[len(path)-1],float(lastvals[0])+I*float(lastvals[1])))

    # Compute variation of the argument
    start=arg(values[0][1])
    t=start
    for k in range(len(values)):
        m=argbranch(values[k][1], t)
        t=m
    return t-start

def count_zeroes(lcalcfile, bleft, tright, n):
    '''Returns number of zeroes of function defined by lcalcfile.txt in the box with bottom left corner bleft and top right corner tright.  Precision is given by n, the number of steps taken along the path'''
    
    
    def delta_L(n):
        '''Estimates variation of argument for the L function'''
        return argvarl(lcalcfile, [bleft, bleft+(tright-bleft).real(), tright, tright-(tright-bleft).real(), bleft], n)
        
    return ((delta_L(n))/(2*pi)).n(P)

%%time
count_zeroes('g1_lcalcfile', .25+13.5*I, 1+15*I, 10)

def split_box(bleft, tright, hv):
    '''Cuts the box with bottom left corner bleft and top right corner tright into two boxes horizontally or vertically depending on hv.  hv will switch after cutting so the next split will go the opposite direction.  Always returns the bottom or left box.'''
    box=[]
    
    if hv==0: # 0 corresponds to cutting horizontally
        new_bleft=bleft
        new_tright=tright.real()+((bleft+tright)*1.0/2).imag()*I
    else: # cut vertically
        new_bleft=bleft
        new_tright=((bleft+tright)/2.0).real()+tright.imag()*I
        
    box=[new_bleft, new_tright, hv]    
    return box

def find_heights(bleft, tright, lcalcfile, prec):
    '''Finds heights of zeroes of function defined in lcalcfile.txt with mesh equal to prec in the box with bottom left corner bleft and top right corner tright.  Prints heights and returns heights in a list.  Estimates heights by middle (vertically) of horizontal strips.'''
    
    heights=[]
    actual_find(bleft, tright, lcalcfile, prec, heights)
    return heights

def actual_find(bleft, tright, lcalcfile, prec, heights):
    n=round(max(tright.imag()-bleft.imag(), tright.real()-bleft.real())*100)
    
    if count_zeroes(lcalcfile, bleft, tright, n)>0:
        if tright.imag()-bleft.imag()<prec:
        #    print (bleft.imag()+tright.imag())/2
            heights.append((bleft.imag()+tright.imag())/2)
    
        else:
            bottombox=split_box(bleft, tright, 0)
            topbox=[bottombox[1].imag()*I+bottombox[0].real(), bottombox[1].real()+(2*bottombox[1].imag()-bottombox[0].imag())*I, 0]
            actual_find(bottombox[0], bottombox[1], lcalcfile, prec, heights)
            actual_find(topbox[0], topbox[1], lcalcfile, prec, heights)

find_heights(.25+10*I, .75+30*I, 'g1_lcalcfile', .01)

def find_real(bleft, tright, lcalcfile, prec):
    '''Finds real part of zeroes of function defined by lcalcfile.txt in a box with bottom left corner bleft and top right corner tright.  Mesh is given by prec.  Estimate of real part is given in the middle of the vertical strip.'''
    
    reals=[]
    actual_reals(bleft, tright, lcalcfile, prec, reals)
    return reals
    
def actual_reals(bleft, tright, lcalcfile, prec, reals):
    n=round(max(tright.real()-bleft.real(), tright.imag()-bleft.imag())*100)
    
    if count_zeroes(lcalcfile, bleft, tright, n)>0:
        if tright.real()-bleft.real() < prec:
         #   print (bleft.real()+tright.real())/2
            reals.append((tright.real()+bleft.real())/2)
            
        else:
            leftbox=split_box(bleft, tright, 1)
            rightbox=[leftbox[1].real()+leftbox[0].imag()*I, leftbox[1].real()*2-leftbox[0].real()+leftbox[1].imag()*I, 1]
            actual_reals(leftbox[0], leftbox[1], lcalcfile, prec, reals)
            actual_reals(rightbox[0], rightbox[1], lcalcfile, prec, reals)

find_real(.1+14*I, .9+30*I, 'g1_lcalcfile', .1)

def find_zeroes(bleft, tright, lcalcfile, prec):
    '''Finds zeroes of function defined by lcalcfile.txt in the box with bottom left corner bleft and top right corner tright. prec sets the size of the mesh.  Estimate gives the middle of the box containing a zero.'''
    
    zeroes=[]
    heights=find_heights(bleft, tright, lcalcfile, prec)
    for i in range(len(heights)):
        temp_bleft=bleft.real()+(heights[i]-prec/2)*I
        temp_tright=tright.real()+(heights[i]+prec/2)*I
        realparts=find_real(temp_bleft, temp_tright, lcalcfile, prec)
        
        for j in range(len(realparts)):
            print heights[i]*I+realparts[j]
            zeroes.append(heights[i]*I+realparts[j])
            
    return zeroes

find_zeroes(10*I, 5+30*I, 'g1_lcalcfile', .1)

find_heights(10*I, 1.5 + 20*I, 'g1_lcalcfile', 0.3)

find_heights(.49+8*I,2+150*I, 'g1_lcalcfile', 0.1)

find_heights(.49+8*I,2+300*I, 'g1_lcalcfile', 0.1)

find_heights(.49+150*I,2+300*I, 'g1_lcalcfile', 0.1)