import sys import re # # Input Format: # ============= # # Specify the exponents: # s = (, , ..., ) # Specify any number of canonical tuples as follows. # A nonzero Jordan tuple with eigenvalue x of size m: # J(m,x) # A zero Jordan tuple of size m at position k (1,...,p): # N(m,k) # A right singular tuple of size m at positions k_1 and k_2: # R(m,k_1,k_2) # A left singular tuple of size m at positions k_1 and k_2: # L(m,k_1,k_2) # # # Internal Datastructures: # ======================== # # Internally, p is an int, s is a list of ints. # The canonical structure is represented as a list of tuples of the following forms: # ('J', m, x) # ('N', m, k) # ('R', m, k_1, k_2) # ('L', m, k_1, k_2) # A placement type is represented by numerating the blocks row by row. # Example: # (I, B_1, B_2, I) # (C_1, I, C_2, I) # is represented as 'IBBICICI' # This representation is unique and unambiguous. # A sign pattern is represented as a string of - and +. # Example: # (+1, +1, -1, +1) # is represented as '++-+' # A tuple pair is represented as a two letter string. # Examples: # LN, LR, RJ, ... # The placement_signs_dict dictionary maps from (placement,signs) to (equation,variant) # The equation_tuples_dict dictionary maps from (equation,tuples) to (interaction expression) # # # placement_signs_dict # placement_signs_dict = { # Jordan tuples ('BIIC','++') : ('I', 1), ('BIIC','+-') : ('III', 1), ('BIIC','-+') : ('III', 1), ('BIIC','--') : ('I', 1) , ('IBCI','++') : ('I', 1), ('IBCI','+-') : ('III', 1), ('IBCI','-+') : ('III', 1), ('IBCI','--') : ('I', 1) , # Jordan and singular tuples ('BIIICC','+++') : ('I', 3), ('BIIICC','++-') : ('II', 1), ('BIIICC','+-+') : ('II', 1), ('BIIICC','+--') : ('III', 1), ('BIIICC','-++') : ('III', 1), ('BIIICC','-+-') : ('II', 1), ('BIIICC','--+') : ('II', 1), ('BIIICC','---') : ('I', 1) , ('IBBCII','+++') : ('I', 2), ('IBBCII','++-') : ('IV', 1), ('IBBCII','+-+') : ('IV', 1), ('IBBCII','+--') : ('III', 2), ('IBBCII','-++') : ('III', 2), ('IBBCII','-+-') : ('IV', 1), ('IBBCII','--+') : ('IV', 1), ('IBBCII','---') : ('I', 1) , ('IBICIC','+++') : ('I', 1), ('IBICIC','++-') : ('II', 1), ('IBICIC','+-+') : ('III', 1), ('IBICIC','+--') : ('II', 1), ('IBICIC','-++') : ('II', 1), ('IBICIC','-+-') : ('III', 1), ('IBICIC','--+') : ('II', 1), ('IBICIC','---') : ('I', 3) , ('BIBICI','+++') : ('I', 1), ('BIBICI','++-') : ('IV', 1), ('BIBICI','+-+') : ('III', 2), ('BIBICI','+--') : ('IV', 1), ('BIBICI','-++') : ('IV', 1), ('BIBICI','-+-') : ('III', 2), ('BIBICI','--+') : ('IV', 1), ('BIBICI','---') : ('I', 2) , ('IIBCCI','+++') : ('I', 3), ('IIBCCI','++-') : ('III', 3), ('IIBCCI','+-+') : ('II', 1), ('IIBCCI','+--') : ('II', 1), ('IIBCCI','-++') : ('II', 1), ('IIBCCI','-+-') : ('II', 1), ('IIBCCI','--+') : ('III', 1), ('IIBCCI','---') : ('I', 1) , ('BBIIIC','+++') : ('I', 2), ('BBIIIC','++-') : ('III', 2), ('BBIIIC','+-+') : ('IV', 1), ('BBIIIC','+--') : ('IV', 1), ('BBIIIC','-++') : ('IV', 1), ('BBIIIC','-+-') : ('IV', 1), ('BBIIIC','--+') : ('III', 1), ('BBIIIC','---') : ('I', 1) , # singular tuples ('BIBIICIC','++++') : ('VIII', 1), ('BIBIICIC','+++-') : ('II', 1), ('BIBIICIC','++-+') : ('IV', 3), ('BIBIICIC','++--') : ('V', 1), ('BIBIICIC','+-++') : ('II', 2), ('BIBIICIC','+-+-') : ('III', 2), ('BIBIICIC','+--+') : ('V', 1), ('BIBIICIC','+---') : ('IV', 3), ('BIBIICIC','-+++') : ('IV', 1), ('BIBIICIC','-++-') : ('V', 1), ('BIBIICIC','-+-+') : ('III', 1), ('BIBIICIC','-+--') : ('II', 1), ('BIBIICIC','--++') : ('V', 1), ('BIBIICIC','--+-') : ('IV', 1), ('BIBIICIC','---+') : ('II', 2), ('BIBIICIC','----') : ('VIII', 3) , ('IBIBCICI','++++') : ('VIII', 3), ('IBIBCICI','+++-') : ('IV', 1), ('IBIBCICI','++-+') : ('II', 2), ('IBIBCICI','++--') : ('V', 1), ('IBIBCICI','+-++') : ('IV', 3), ('IBIBCICI','+-+-') : ('III', 4), ('IBIBCICI','+--+') : ('V', 1), ('IBIBCICI','+---') : ('II', 2), ('IBIBCICI','-+++') : ('II', 1), ('IBIBCICI','-++-') : ('V', 1), ('IBIBCICI','-+-+') : ('III', 2), ('IBIBCICI','-+--') : ('IV', 1), ('IBIBCICI','--++') : ('V', 1), ('IBIBCICI','--+-') : ('II', 1), ('IBIBCICI','---+') : ('IV', 3), ('IBIBCICI','----') : ('VIII', 1) , ('BBIIIICC','++++') : ('I', 4), ('BBIIIICC','+++-') : ('II', 2), ('BBIIIICC','++-+') : ('II', 2), ('BBIIIICC','++--') : ('III', 2), ('BBIIIICC','+-++') : ('IV', 3), ('BBIIIICC','+-+-') : ('V', 1), ('BBIIIICC','+--+') : ('VI', 1), ('BBIIIICC','+---') : ('IV', 1), ('BBIIIICC','-+++') : ('IV', 3), ('BBIIIICC','-++-') : ('VII', 1), ('BBIIIICC','-+-+') : ('V', 1), ('BBIIIICC','-+--') : ('IV', 1), ('BBIIIICC','--++') : ('III', 1), ('BBIIIICC','--+-') : ('II', 1), ('BBIIIICC','---+') : ('II', 1), ('BBIIIICC','----') : ('I', 1) , ('IIBBCCII','++++') : ('I', 4), ('IIBBCCII','+++-') : ('IV', 3), ('IIBBCCII','++-+') : ('IV', 3), ('IIBBCCII','++--') : ('III', 4), ('IIBBCCII','+-++') : ('II', 2), ('IIBBCCII','+-+-') : ('V', 1), ('IIBBCCII','+--+') : ('VII', 1), ('IIBBCCII','+---') : ('II', 1), ('IIBBCCII','-+++') : ('II', 2), ('IIBBCCII','-++-') : ('VI', 1), ('IIBBCCII','-+-+') : ('V', 1), ('IIBBCCII','-+--') : ('II', 1), ('IIBBCCII','--++') : ('III', 2), ('IIBBCCII','--+-') : ('IV', 1), ('IIBBCCII','---+') : ('IV', 1), ('IIBBCCII','----') : ('I', 1) , ('BIIBICCI','++++') : ('I', 3), ('BIIBICCI','+++-') : ('IV', 3), ('BIIBICCI','++-+') : ('II', 1), ('BIIBICCI','++--') : ('VII', 1), ('BIIBICCI','+-++') : ('II', 1), ('BIIBICCI','+-+-') : ('V', 1), ('BIIBICCI','+--+') : ('III', 2), ('BIIBICCI','+---') : ('IV', 1), ('BIIBICCI','-+++') : ('IV', 3), ('BIIBICCI','-++-') : ('III', 4), ('BIIBICCI','-+-+') : ('V', 1), ('BIIBICCI','-+--') : ('II', 2), ('BIIBICCI','--++') : ('VI', 1), ('BIIBICCI','--+-') : ('II', 2), ('BIIBICCI','---+') : ('IV', 1), ('BIIBICCI','----') : ('I', 2) , ('IBBICIIC','++++') : ('I', 2), ('IBBICIIC','+++-') : ('II', 2), ('IBBICIIC','++-+') : ('IV', 1), ('IBBICIIC','++--') : ('VI', 1), ('IBBICIIC','+-++') : ('IV', 1), ('IBBICIIC','+-+-') : ('V', 1), ('IBBICIIC','+--+') : ('III', 1), ('IBBICIIC','+---') : ('II', 1), ('IBBICIIC','-+++') : ('II', 2), ('IBBICIIC','-++-') : ('III', 2), ('IBBICIIC','-+-+') : ('V', 1), ('IBBICIIC','-+--') : ('IV', 3), ('IBBICIIC','--++') : ('VII', 1), ('IBBICIIC','--+-') : ('IV', 3), ('IBBICIIC','---+') : ('II', 1), ('IBBICIIC','----') : ('I', 3) } # # equation_tuples_dict # equation_tuples_dict = { # Jordan tuples ('I','JJ') : 'min(l,m)', ('I','JN') : '0', ('I','NJ') : '0', ('I','NN') : 'min(l,m)' , #('III','JJ') : '0', (impossible) ('III','JN') : '0', ('III','NJ') : '0', ('III','NN') : '0' , # Jordan and singular tuples ('I','JR') : '0', ('I','NR') : 'min(l,m)', ('I','JL') : '0', ('I','NL') : 'min(l,m+1)', ('I','RJ') : '0', ('I','RN') : 'min(l,m)', ('I','LJ') : '0', ('I','LN') : 'min(l+1,m)' , ('II','JR') : '0', ('II','NR') : '0', ('II','JL') : 'l', ('II','NL') : 'l' , ('III','JR') : '0', ('III','NR') : '0', ('III','JL') : 'l', ('III','NL') : 'l', ('III','RJ') : '0', ('III','RN') : '0', ('III','LJ') : '0', ('III','LN') : '0' , ('IV','RJ') : 'm', ('IV','RN') : 'm', ('IV','LJ') : '0', ('IV','LN') : '0' , # singular tuples ('I','RR') : 'min(l,m)', ('I','RL') : 'min(l,m+1)', ('I','LR') : 'min(l+1,m)', ('I','LL') : '1+min(l,m)' , ('II','RR') : '0', ('II','RL') : 'l', ('II','LR') : '0', ('II','LL') : 'l+1' , ('III','RR') : '0', ('III','RL') : '0', ('III','LR') : '0', ('III','LL') : '0' , ('IV','RR') : 'm', ('IV','RL') : 'm+1', ('IV','LR') : '0', ('IV','LL') : '0' , ('V','RR') : 'max(0,m-l)', ('V','RL') : 'l+m+1', ('V','LR') : '0', ('V','LL') : 'max(0,l-m)' , ('VI','RR') : 'max(0,m-l-1)', ('VI','RL') : 'l+m+2', ('VI','LR') : '0', ('VI','LL') : 'max(0,l-m-1)' , ('VII','RR') : 'max(0,m-l+1)', ('VII','RL') : 'l+m', ('VII','LR') : '0', ('VII','LL') : 'max(0,l-m+1)' , ('VIII','RR') : '1+min(l,m)', ('VIII','RL') : 'min(l,m)', ('VIII','LR') : 'min(l,m)', ('VIII','LL') : '1+min(l,m)' } # # Return the dimensions of the given block at a certain position # def get_one_block_size(s, blkdesc, pos): blktype = blkdesc[0] right = (blktype == 'R') left = (blktype == 'L') if blktype == 'J' or blktype == 'N': m,n = blkdesc[1], blkdesc[1] elif left or right: k_1,k_2 = blkdesc[2], blkdesc[3] s1 = s[k_1-1] s2 = s[k_2-1] samesign = (s1 == s2) m0 = blkdesc[1] m1 = blkdesc[1]+1 # determine the possible dimensions if right: F_dim = (m0,m1) if samesign: G_dim = (m1,m0) else: G_dim = (m0,m1) if s1 == 1: outside_dim = (m1,m1) between_dim = (m0,m0) else: outside_dim = (m0,m0) between_dim = (m1,m1) elif left: F_dim = (m1,m0) if samesign: G_dim = (m0,m1) else: G_dim = (m1,m0) if s1 == 1: outside_dim = (m0,m0) between_dim = (m1,m1) else: outside_dim = (m1,m1) between_dim = (m0,m0) # determine where pos is in relation to F and G outside = (pos < k_1 or pos > k_2) F = (pos == k_1) between = (pos > k_1 and pos < k_2) G = (pos == k_2) # pick the appropriate dimension if outside: m,n = outside_dim elif F: m,n = F_dim elif between: m,n = between_dim elif G: m,n = G_dim return (m,n) # # Return all the dimensions of the given tuple # def get_all_block_sizes(s, blkdesc): p = len(s) n = [0]*p for i in range(p): mi,ni = get_one_block_size(s, blkdesc, i+1) if s[i] == 1: n[i] += ni else: n[i] += mi return n # # Get nontrivials positions # def get_nontrivial_positions(B): if B[0] == 'J': return [1] elif B[0] == 'N': return [B[2]] else: return [B[2],B[3]] # # Return the interaction between B and C # def get_interaction(s, B, C): p = len(s) regular = ('J','N') singular = ('R','L') Bmod = list(B) Cmod = list(C) smod = list(s) pmod = p # extract nontrivial positions Bpos = get_nontrivial_positions(Bmod) Cpos = get_nontrivial_positions(Cmod) # remove trivial positions i = 1 while i <= pmod: if i not in Bpos and i not in Cpos: # remove position i from the list of signs del smod[i-1] # update positions to the right of i for li in (Bpos, Cpos): for j in range(len(li)): if li[j] > i: li[j] -= 1 i -= 1 # dont advance pmod -= 1 i += 1 # remove shared positions i = 1 while i <= pmod: if i in Bpos and i in Cpos: # expand position i if smod[i-1] == 1: # move C one step to the right smod.insert(i,+1) for j in range(len(Cpos)): if Cpos[j] >= i: Cpos[j] += 1 for j in range(len(Bpos)): if Bpos[j] > i: Bpos[j] += 1 else: # move B one step to the right smod.insert(i,-1) for j in range(len(Bpos)): if Bpos[j] >= i: Bpos[j] += 1 for j in range(len(Cpos)): if Cpos[j] > i: Cpos[j] += 1 pmod += 1 i += 1 # create modified tuples if Bmod[0] != 'J': Bmod[2] = Bpos[0] if Bmod[0] in singular: Bmod[3] = Bpos[1] if Cmod[0] != 'J': Cmod[2] = Cpos[0] if Cmod[0] in singular: Cmod[3] = Cpos[1] # form sign pattern signs = ''.join([{1:'+',-1:'-'}[x] for x in smod]) # form placement pattern v1 = [{True:'B',False:'I'}[x in Bpos] for x in range(1,pmod+1)] v2 = [{True:'C',False:'I'}[x in Cpos] for x in range(1,pmod+1)] placement = ''.join(v1) + ''.join(v2) # (i) loopup equation type and variant equation, variant = placement_signs_dict[(placement,signs)] # (ii) modify tuple description tuple_desc = [Bmod[0],Cmod[0]] tuple_org = ''.join(tuple_desc) if variant in (2,4): tuple_desc[0] = {'R':'L','L':'R'}[tuple_desc[0]] if variant in (3,4): tuple_desc[1] = {'R':'L','L':'R'}[tuple_desc[1]] tuple_desc = ''.join(tuple_desc) # (iii) lookup interaction expression interaction_expr = equation_tuples_dict[(equation,tuple_desc)] # (iv) evaluate interaction expression l = Bmod[1] m = Cmod[1] if tuple_desc == 'JJ' and Bmod[2] != Cmod[2]: # Jordan tuples belonging to mutually different eigenvalues have no interaction interaction_expr = '0' d = eval(interaction_expr, {'min':min, 'max':max, 'l':l, 'm':m}) print ' [%s] Equation and variant: %s:%d' % (signs, equation, variant) print ' [%s] Tuple & modified tuple: %s --> %s' % (''.join(v1), tuple_org, tuple_desc) print ' [%s] Interaction: %s = %d' % (''.join(v2), interaction_expr, d) return d # # Return the lack of dimension # def get_lack_of_dimension(n): x = 0 for i0 in range(len(n)): i1 = (i0+1) % len(n) x += (n[i1] - n[i0])**2 return - x // 2 # # Return definitions of MATLAB functions that construct canonical blocks # def matlab_auxilliary_functions(): li = """ function A = J(m,x) A = gallery('jordbloc', m, x); function A = F(m) A = eye(m+1); A = A(1:m,:); function A = G(m) A = J(m+1,0); A = A(1:m,:); """.splitlines() return li def describe_block(B): if B[0] in ('J','N'): return '%s(%s,%s)' % (B[0], B[1], B[2]) else: return '%s(%s,%s,%s)' % (B[0], B[1], B[2], B[3]) def header1(s): print print s print '='*len(s) if __name__ == '__main__': # regexp to find text within parentheses re_paren = re.compile('\((.*)\)') def get_paren(s): return re_paren.search(s).groups()[0] def get_args(s): ps = get_paren(s) return [x.strip() for x in ps.split(',')] # exponents s = None # blocks blocks = [] # read input for line in sys.stdin.readlines(): line = line.strip() # parse line if line.find('=') != -1: # contains '=' lhs, rhs = line.split('=') lhs, rhs = lhs.strip(), rhs.strip() if lhs == 's': # set exponents s = [int(x) for x in get_args(rhs)] else: # does not contain '=' blktype = line[0] args = get_args(line) if blktype == 'J': m = int(args[0]) x = args[1] args = [m, x] elif blktype == 'N': m = int(args[0]) k = int(args[1]) args = [m, k] elif blktype == 'R': m = int(args[0]) k_1 = int(args[1]) k_2 = int(args[2]) args = [m, k_1, k_2] elif blktype == 'L': m = int(args[0]) k_1 = int(args[1]) k_2 = int(args[2]) args = [m, k_1, k_2] blkdesc = tuple([blktype] + args) blocks.append(blkdesc) # detect invalid input err = False p = len(s) if s is None: print 'Error: you must specify s' err = True else: if len(s) < 1: print 'Error: must specify at least one exponent' err = True if len([x for x in s if x not in (1,-1)]) > 0: print 'Error: specify only +1 or -1 as exponents' err = True for blkdesc in blocks: if blkdesc[0] == 'J': if blkdesc[1] < 1: print 'Error: Size of Jordan block (J) must be at least one' err = True elif blkdesc[0] == 'N': if blkdesc[1] < 1: print 'Error: Size of nilpotent Jordan block (N) must be at least one' err = True if blkdesc[2] < 1 or blkdesc[2] > p: print 'Error: Nilpotent Jordan block (N) has an invalid position (must be 1,...,%d)' % p err = True elif blkdesc[0] == 'R': if blkdesc[1] < 0: print 'Error: Size of right singular block (R) must be at least zero' err = True if blkdesc[2] < 1 or blkdesc[2] > p or blkdesc[3] < 1 or blkdesc[3] > p: print 'Error: Right singular block (R) has an invalid position (must be 1,...,%d)' % p err = True if blkdesc[3] <= blkdesc[2]: print 'Error: Right singular block must end after it starts' err = True elif blkdesc[0] == 'L': if blkdesc[1] < 0: print 'Error: Size of left singular block (L) must be at least zero' err = True if blkdesc[2] < 1 or blkdesc[2] > p or blkdesc[3] < 1 or blkdesc[3] > p: print 'Error: Left singular block (L) has an invalid position (must be 1,...,%d)' % p err = True if blkdesc[3] <= blkdesc[2]: print 'Error: Left singular block must end after it starts' err = True if err: sys.exit(1) # determine dimensions n n = [0]*p for blkdesc in blocks: nn = get_all_block_sizes(s, blkdesc) n = [x+y for x,y in zip(n,nn)] header1('Length of tuples:') print 'p = %d' % p header1('Dimension tuple:') print 'n = (' + ', '.join([str(x) for x in n]) + ')' header1('Sign tuple:') print 's = (' + ', '.join(['%+d' % x for x in s]) + ')' header1('Canonical tuples:') print print ' ' + ' '.join(['%+2d' % x for x in s]) for i,blkdesc in enumerate(blocks): print '#%-3d %10s: ' % (i+1, describe_block(blkdesc)), if blkdesc[0] == 'J': lo = 1 else: lo = blkdesc[2] if blkdesc[0] in ('J','N'): hi = lo a = blkdesc[0] b = a else: hi = blkdesc[3] s1 = s[lo-1] s2 = s[hi-1] if s1 == s2: if blkdesc[0] == 'R': a = 'F' b = 'GT' else: a = 'FT' b = 'G' else: if blkdesc[0] == 'R': a = 'F' b = 'G' else: a = 'FT' b = 'GT' sys.stdout.write(' I ' * (lo-1)) sys.stdout.write(' %-2s' % a) sys.stdout.write(' I ' * (hi-lo-1)) if hi > lo: sys.stdout.write(' %-2s' % b) print ' I ' * (p-hi) # compute interaction header1('Pairwise interactions') d = 0 for i,B in enumerate(blocks): for j,C in enumerate(blocks): # compute interaction between block i and block j (NOT the combined interaction) print print 'Tuple #%d, %s, interacting with tuple #%d, %s' % (i+1, describe_block(B), j+1, describe_block(C)) d += get_interaction(s, B, C) # compute cointeraction c = d + get_lack_of_dimension(n) header1('Summary') print 'Kernel dimension: %3d' % (d) print 'Codimension: %3d' % (c) # output matlab code that computes the codimension directly from the Kroneckerized Sylvester operator header1('MATLAB code') print 'function [c,d,Z] = compute_codimension()' # compute size of Kroneckerized system matrix Z mz = 0 nz = 0 for i0 in range(p): i1 = (i0+1)%p nk = n[i0] nk1 = n[i1] mz += nk*nk1 nz += nk*nk print '%% System matrix Z is %d x %d' % (mz, nz) print 'Z = zeros(%d, %d);' % (mz, nz) # build Ak for k in range(p): k1 = (k+1)%p nk = n[k] nk1 = n[k1] sk = s[k] li = [] for B in blocks: Bpos = get_nontrivial_positions(B) mblk, nblk = get_one_block_size(s, B, k+1) if (k+1) not in Bpos: # identity blk = 'eye(%d)' % (mblk) elif B[0] == 'J': # Jordan block blk = 'J(%d,%f)' % (B[1], float(B[2])) elif B[0] == 'N': # Nilpotent Jordan block blk = 'J(%d,0)' % (B[1]) elif B[0] in ('R','L'): # L/R s1 = s[Bpos[0]-1] s2 = s[Bpos[1]-1] if (k+1) == Bpos[0]: # F transpose = '' if B[0] == 'L': transpose = "'" blk = 'F(%d)%s' % (B[1],transpose) else: # G transpose = '' if (B[0] == 'R' and s1 == s2) or (B[0] == 'L' and s1 != s2): transpose = "'" blk = 'G(%d)%s' % (B[1],transpose) li.append(blk) print 'A%d = blkdiag(%s);' % (k+1,','.join(li)) # assemble Z iz1 = 1 jz1 = 1 for k in range(p): # fill in the first block row of Z k1 = (k+1)%p nk = n[k] nk1 = n[k1] sk = s[k] mblk1, nblk1 = nk*nk1, nk*nk mblk2, nblk2 = nk*nk1, nk1*nk1 iz2 = iz1 if nz == 0: jz2 = jz1+nblk1-1 + 1 else: jz2 = ((jz1+nblk1-1) % nz) + 1 # second block (first superdiagonal) if sk == 1: print "Z(%d:%d,%d:%d) = - kron(A%d',eye(%d));" % (iz2, iz2+mblk2-1, jz2, jz2+nblk2-1, k+1, nk1) else: print "Z(%d:%d,%d:%d) = kron(eye(%d),A%d);" % (iz2, iz2+mblk2-1, jz2, jz2+nblk2-1, nk1, k+1) # first block (main diagonal) if sk == 1: print "Z(%d:%d,%d:%d) = kron(eye(%d),A%d);" % (iz1, iz1+mblk1-1, jz1, jz1+nblk1-1, nk, k+1) else: print "Z(%d:%d,%d:%d) = - kron(A%d',eye(%d));" % (iz1, iz1+mblk1-1, jz1, jz1+nblk1-1, k+1, nk) # advance if mz == 0: iz1 = (iz1 + mblk1 - 1) + 1 else: iz1 = ((iz1 + mblk1 - 1) % mz) + 1 if nz == 0: jz1 = (jz1 + nblk1 - 1) + 1 else: jz1 = ((jz1 + nblk1 - 1) % nz) + 1 # compute kernel dimension print 'd = size(Z,2) - rank(Z);' # compute codimension lack = get_lack_of_dimension(n) print 'c = d + (%d);' % lack # auxilliary functions print '\n'.join(matlab_auxilliary_functions())