program polys; (* Program written by H.R.Morton and H.B.Short, University of Liverpool, May 1985. Amended most recently, March 1997, minor modification Nov 1997. Further modification May 2014, to alter input for braid *) const maxl = 160; globalnfact = 362880; (* largest n_factorial used *) globaln = 9; (* largest number of strings used *) globaldn = 8; (* = globaln - 1 *) jmax = 130; type row1 = array[1..globaldn] of 0..globaln; row = array[-1..maxl] of integer; ptr1 = ^term; term = record coefft : integer; expon : 0..maxl; nexterm : ptr1 end; label start; var s, g, k, i, j, n, dn, nfact, len, r, sig, tr : integer; algcrt, count : integer; mistake : boolean; zerovect : row; z : row; conw : array[1..9] of row; braidcomment : array[1..9] of packed array [1..160] of char; braid : array[0..maxl] of -9..9; braidstr, anstr,filename : packed array [1..160] of char; m : array[1..globalnfact] of row1; grow : row1; factor : array[0..globaln] of 1..globalnfact; ch,ans : char; c : array[1..globalnfact] of ptr1; outfile : text; w, t : array[0..maxl, 0..globaln] of integer; results : array[1..9, 0..maxl,0..globaln] of integer; alex : array[1..9, -jmax..jmax] of integer; jopoly : array[1..9, -maxl..jmax] of integer; algcr : array[1..9] of integer; example : array[1..9] of packed array [1..160] of char; function max(x, y : integer):integer; var m1 : integer; begin m1 := y; if x > y then m1 := x; max := m1; end; function bin(v, v1 : integer) : integer; var j,denom,num: integer; begin denom := 1; num := 1; for j := 1 to v do num := num*j; for j := 1 to v1 do denom := -denom*j; for j := 1 to v-v1 do denom := denom*j; bin:= num div denom ; end; function binreal(v, v1 : integer) : double; var j, vmax,vmin : integer; denom,num : double ; begin denom:=1; num:=1; vmax:=max(v1,v-v1);if vmax=v1 then vmin:=v-v1 else vmin:=v1; for j := vmax+1 to v do num := num*j; for j := 1 to vmin do denom := denom*j; if odd(v1) then num:=-num; binreal := num / denom; end; procedure addmat(i : integer); var j: integer; begin j := grow[i]+1; if (j>i) then begin grow[i] := 0; if i+1 < n then addmat(i+1); end else grow[i] := j; end; procedure extract(g : integer); var d : integer; t1 : ptr1; begin t1 := c[g]; z := zerovect; d := -1; while (t1 <> nil) do begin d := t1^.expon; z[d] := t1^.coefft; t1 := t1^.nexterm; end; z[-1] := d; end; procedure alexander(x : integer); var i, j, k, t : integer; alexx : array[-jmax..jmax] of double; begin for i := -jmax to jmax do alexx[i] := 0; for i := 0 to maxl do begin t := conw[x,i]; k := i-n+1; if t<>0 then for j := 0 to k do alexx[2*j-k] := alexx[2*j-k] + binreal(k,k-j)*t; end; for i := -jmax to jmax do if abs(alexx[i]) < maxint then alex[x,i] := round(alexx[i]) else writeln(' Alexander coeffts greater than maxint'); end; procedure stack(t : ptr1); var t1 : ptr1; begin while t <> nil do begin t1 := t; t := t^.nexterm; dispose(t1); end; end; procedure insert(q : integer); var i, d : integer; t1,t2 : ptr1; begin d := -1; new(t1); c[q] := t1; t2 := t1; for i := 0 to z[-1] do if z[i] <> 0 then begin t1^.coefft := z[i]; d := i; t1^.expon := i; t2^.nexterm := t1; t2 := t1; new(t1); end; if d = -1 then c[q] := nil; dispose(t1); t2^.nexterm := nil; end; procedure multiply(a,b : integer); var j, k, r1, g, h, deg : integer; x, y : row; t1, t2 : ptr1; begin r1 := abs(r); for g := a to b do begin if r1 > 1 then k :=m [g,r1-1] else k := 0; j := k - m[g,r1]; if j >= 0 then begin h := (g+(j+1)*factor[r1]) - j*factor[r1-1]; t1 := c[g]; t2 := c[h]; if (t1 <> nil) or (t2 <> nil) then begin extract(g); x := z; extract(h); y := z; deg := max(x[-1],y[-1]); if r>0 then begin c[g] := t2; stack(t1); z[0] := x[0]; for j := 0 to deg do z[j+1] := x[j+1] + y[j]; z[-1] := max(deg,y[-1]+1); insert(h); end else begin c[h] := t1; stack(t2); z[0] := y[0]; for j := 0 to deg do z[j+1] := y[j+1] - x[j]; z[-1] := max(deg,x[-1]+1); insert(g); end; end; end; end; end; procedure wrapup; var i, j, k, s, a, b, g, rr, f, ff, gg : integer; x : row; begin for s := n-1 downto 1 do begin b := 0; f := factor[s]; for i := 0 to n-s-1 do begin b := b + factor[s+i]; a := b+1-f; for rr := 1 to s-1 do begin k := (s-rr+1) * f; r := rr; multiply(a+k,b+k); for g := a+k to b+k do if c[g] <> nil then begin extract(g-f); x := z; extract(g); z[-1] := max(x[-1],z[-1]); for j := 0 to z[-1] do z[j] := z[j]+x[j]; insert(g-f); end; end; ff := f-factor[s+i]; for g := a to b do if (i > 0) and (c[g] <> nil) then begin extract(g); x := z; gg := g+ff; extract(gg); z[-1] := max(x[-1],z[-1]); for j := 0 to z[-1] do z[j] := z[j] + x[j]; insert(gg); end; end; writeln('wrapup done for s = ',s); end; f := 0; for j := 0 to n-1 do begin f := f + factor[j]; extract(f); for i := 0 to maxl do t[i,j] := z[i]; end; end; procedure init; var i, j : integer; begin dn:=n-1; nfact:=factor[n]; for i := 0 to maxl do for j := 0 to n do begin t[i,j] := 0; w[i,j] := 0; end; for i:=-1 to maxl do zerovect[i]:=0; (* needed on SGI as basic initialisation doesn't work *) new(c[1]); c[1]^.coefft := 1; c[1]^.expon := 0; c[1]^.nexterm := nil; for i := 2 to nfact do c[i] := nil; end; procedure jones; var i,j,s,k,u,tr:integer; jones:array[-maxl..jmax] of double; begin for i:=-maxl to jmax do jones[i]:=0; for j:=0 to len do for s:=0 to n-1 do begin tr:=t[j,s]; k:=n-s-1; sig:=-1; if tr<>0 then begin for i:=0 to k do sig:=-sig; for i:=0 to k do for u:=0 to j do begin r:=2*i-2*u+k+j; if (r>-maxl) and (r0 then begin g:=(k+algcr[x]*2-2*(n-1)); if fn=1 then writeln(outfile,t,' (',g:4,')') else writeln(output,t,' (',g:4,')'); end; end; end; procedure clearscreen; begin writeln(output,' press return'); readln(input,anstr); end; procedure screenout(x:integer); var j,k,u:integer; line:packed array [1..80] of char; begin clearscreen; line:=' ___________________________________________________________'; writeln(output,'braid : ',example[x]);writeln(output); writeln(output,'algebraic cross no:',algcr[x]:3);writeln(output); for u:=0 to n-1 do write(output,algcr[x]-n+1+2*u:10); writeln(output); writeln(output,line); writeln(output); for k:=0 to maxl do begin j:=0; s:=1; for u:=0 to n-1 do begin tr:=results[x,k,u]; if not(tr=0) then begin j:=1; write(output,tr:10*s); s:=1; end else s:=s+1; end; s:=10*(s-1)+6; if j=1 then writeln(output,'|':s,k-n+1:3); end;writeln(output); write(output,' for Conway poly '); clearscreen; writeln(output,'Conway polynomial'); for k:=0 to maxl do if conw[x,k]<>0 then writeln(output,conw[x,k],' (',k-n+1:3,')'); writeln(output,'Alexander Polynomial'); for k:=-jmax to jmax do if alex[x,k]<>0 then writeln(output,alex[x,k],' (',k:3,')'); writeln(output); writeln(output,'Jones Polynomial'); jonesprint(x,0); writeln(output); end; procedure fileout(x:integer); var j,k,u:integer; line:packed array [1..80] of char; begin writeln(outfile);writeln(outfile,'*** ',example[x],' ***'); writeln('outfile'); line:=' _______________________________________________________'; writeln(outfile,anstr); for u:=0 to n-1 do write(outfile,algcr[x]-n+1+2*u:10); writeln(outfile);writeln(outfile,line); for k:=0 to maxl do begin j:=0; s:=1; for u:=0 to n-1 do begin tr:=results[x,k,u]; if not(tr=0) then begin j:=1; write(outfile,tr:10*s); s:=1; end else s:=s+1; end; if j=1 then writeln(outfile,'|':8*s,k-n+1:3); end;writeln(outfile); writeln(outfile,'Conway Polynomial'); for k:=0 to maxl do if conw[x,k]<>0 then writeln(outfile,conw[x,k],' (',k-n+1:3,')'); writeln(outfile,' Alexander Polynomial'); for k:=-jmax to jmax do if alex[x,k]<>0 then writeln(outfile,alex[x,k],' (',k:3,')'); writeln(outfile);writeln(outfile,'Jones Polynomial'); jonesprint(x,1); end; procedure compare(x,y:integer); var i,j,k,ax,ay,tr1,tr2,diff,difa,difb:integer; trj1,trj2,ja,jb:integer; ch:char; flag:boolean; begin clearscreen; writeln(output,x:2,' : ',example[x]); writeln(output,y:2,' : ',example[y]); ax:=algcr[x]-n+1; ay:=algcr[y]-n+1; flag:=false; difa:=0; difb:=0; diff:=ax-ay; writeln(output,'----------',x:2,' ------------',y:2); if diff>0 then difb:=diff else difa:=-diff;diff:=abs(diff); if trunc(diff/2)=diff/2 then begin for j:=-diff to n-1+diff do for i:=0 to maxl do begin ja:=j+difa; jb:=j+difb; tr1:=0; tr2:=0; if (ja>=0) and (ja<= n-1) then tr1:=results[x,i,ja]; if (jb>=0) and (jb<= n-1) then tr2:=results[y,i,jb]; k:= ax+2*(j+difa); if tr1=tr2 then ch:=' ' else begin flag:=true;ch:='*'; end; if not((tr1=0) and (tr2=0)) then writeln(output,tr1,' ',tr2,' ',ch,' (',i-n+1:3,',',k:3,')'); end; if ((flag=false) and (diff=0)) then writeln(output,' !! polynomials are identical !'); writeln(output,'Conway polynomials are '); flag:=false; for i:=0 to maxl do begin tr1:=conw[x,i];tr2:=conw[y,i]; if (tr1<>tr2) then begin ch:='*'; flag:=true end else ch:=' '; if not((tr1=0) and (tr2=0)) then writeln(output,tr1,' ',tr2,' ',ch); end; if (flag=false) then writeln(output,' !! Conway polynomials identical !!'); writeln(output,'Jones polynomials are '); flag:=false; j:=2*algcr[x] -2*(n-1); for i:=-maxl to jmax do begin tr1:=i+2*difa; tr2:=i+2*difb; trj1:=0; trj2:=0; if (tr1>-maxl) and (tr1-maxl) and (tr2trj2) then begin ch:='*'; flag:=true end else ch:=' '; if not((trj1=0) and (trj2=0)) then writeln(output,trj1,' ',trj2,' ',ch,' ',i+j+diff); end; if (flag=false) then writeln(output,' !! Jones polynomials identical !!'); end else writeln(output,'braids not comparable as different number of components'); end; procedure menu; var i, j, k : integer; menu : array[1..4] of packed array [1..80] of char; begin menu[1] := ' review information on one example'; menu[2] := ' send information of one example to file'; menu[3] := ' compare information for two braids'; menu[4] := ' exit from menu'; repeat writeln(' ref number braid'); for i := 1 to 9 do if example[i] <> '*' then begin writeln(i:8, ' ', example[i],' : ',braidcomment[i]); end; writeln; writeln(' Choose one of the following options :'); for i := 1 to 4 do writeln(i:2, ' : ' ,menu[i]); readln(i); k := 0; if (i <> 1) and (i <> 2) and (i <> 3) then i := 4; case i of 1: begin writeln(' give ref number'); writeln(' ref number braid'); for i := 1 to 9 do if example[i] <> '*' then begin writeln(i:8, ' ', example[i],' : ',braidcomment[i]); end; readln(j); screenout(j); end; 3: begin writeln(' give ref numbers (separated by space)'); writeln(' ref number braid'); for i := 1 to 9 do if example[i] <> '*' then begin writeln(i:8, ' ', example[i],' : ',braidcomment[i]); end; readln(j,k); compare(j,k); end; 2: begin writeln('The output will be sent to the file called: ', filename); writeln('This will overwrite any information in this file'); writeln('If you want to write to a different file type y'); writeln('Otherwise press RETURN twice to use the file ', filename); readln(ch); if ch='y' then begin writeln('Give filename now'); readln(filename); writeln('The output will be sent to the file called ', filename); end; //rewrite(outfile,filename); writeln(' ref number braid'); for i := 1 to 9 do if example[i] <> '*' then begin writeln(i:8, ' ', example[i],' : ',braidcomment[i]); end; writeln(' give number of example (as in displayed list) '); readln(j); writeln(' give comment if required'); readln(anstr); fileout(j); end; 4: k := -100; end; until k = -100; end; (* ******************* MAIN PROGRAM STARTS HERE ***************** *) begin; count := 0; filename:='braidout'; //rewrite(outfile,filename); n:=globaln;nfact:=globalnfact; (* At this point in the original IBM program the start of the stack was recorded using mark(top_of_stackmarker); *) for i := 1 to n-1 do grow[i] := 0; for g := 1 to nfact do begin for i := 1 to n-1 do m[g][i]:= grow[i] ; addmat(1); end; for i := 1 to 9 do example[i] := '*'; factor[0] := 1; for i := 1 to n do factor[i] := i * factor[i-1]; repeat writeln; start: writeln('Give braid to be calculated, on at most', globaln:3,' strings'); writeln('Use 1,2,.. for the standard braid generators and -1,-2,... for their inverses.'); writeln; writeln('Write the braid word without any separators, and finish with 0'); writeln('for example 1-232-3-10'); writeln; writeln('Type * to exit the program'); readln(braidstr); (* writeln('and second part'); readln(anstr); braidstr := braidstr + anstr; *) (* Removed for SGI implementation as I cannot find how to adjoin two strings *) if braidstr[1] <> '*' then begin writeln('The number of strings to be used in this braid is currently',n:3); writeln('If you want to change this number type y'); writeln('else press any key, and return, to continue with' ,n:3, ' strings'); ans:='N'; readln(ans); while ans='y' do begin writeln('Give the number of strings in the braid'); readln(n); if (n>globaln) or (n<2) then writeln('The number of strings should be >1 and at most',globaln:3) else ans:='N'; end; init; (* This next is an alternative means to find the length of the input word. *) i:=1; while braidstr[i]<>'0' do begin i:=i+1;end; len:=i-1; (* len := length(braidstr); *) count := count + 1; if count = 10 then count := 1; example[count] := braidstr; sig := 1; algcrt := 0; i := 0; mistake := false; for j := 1 to len do begin ch := braidstr[j]; if ch = '-' then sig := -1 else begin i := i+1; braid[i] := sig*(ord(ch) - ord('0')); algcrt := algcrt + sig; sig := 1; if abs(braid[i]) > n-1 then mistake := true; end; end; len := i; algcr[count] := algcrt; if (len > maxl-15) then writeln('DANGER..braid too long!!'); if mistake = true then begin writeln('wrong elements in braid'); goto start; end; writeln('The polynomial of the closed braid is being calculated'); for i := 1 to len do begin r := braid[i]; multiply(1,nfact); end; wrapup; for j := 0 to maxl do for s := 0 to n-1 do begin tr := t[j,s]; k := n-s-1; if tr <> 0 then for i := 0 to k do w[j+s,i] := w[j+s,i]+tr* bin(k,i); end; for i := 0 to maxl do begin r := 0; for j := 0 to n-1 do begin r := r+w[i,j]; results[count,i,j] := w[i,j]; end; conw[count,i] := r; end; alexander(count); jones; screenout(count); writeln('do you want to add a comment to this braid (y/n)?'); readln(anstr); braidcomment[count]:=''; if (anstr = 'Y') or (anstr = 'y') then begin writeln('type comment (one line) now'); readln(anstr); braidcomment[count]:=anstr; end; writeln('to consult menu answer M '); readln(anstr); if (anstr = 'm') or (anstr = 'M') then menu; end; until braidstr[1] = '*'; end.