let mq= 1000;;
let print_neg= false;;
let print_inv= false;;

open Array;;
open Printf;;

type vl= B|G of (int*((int*(int array))array));;
type tc= C|E of int;;

let ca=init mq(fun q->make(2*(q+1))C);;

let coef zr i p=match p with (m,a)->if((i<m)||(i>=(length a)+m))then zr else a.(i-m);;

let sum sm zr p1 p2=let ((m1,a1),(m2,a2))=(p1,p2) in let m=min m1 m2 in
(m,init(max(length a1+m1)(length a2+m2)-m)(fun i->sm(coef zr(i+m)p1)(coef zr(i+m)p2)));;

let prod pr sm p1 p2=let ((m1,a1),(m2,a2))=(p1,p2) in fold_left(sm)
(0,[||])(mapi(fun i c->(i+m1+m2,map(fun b->pr b c)a2))a1);;

let sign p=if p>0 then 1 else -1;;
let m2 p=2*(abs(p) mod 2)-1;;

let rv p=(1+p/2)*(m2 p);;
let iv p=abs(p)*2+(sign(p)-3)/2;;

let gd p q=let d=q/abs(p) in 
sign(p)*(if d mod 2==0 then d else d+1);;

let np p q=sign(p)*(q-(gd p q)*p);;

let sp k=match k with 0->(0,[||])|_->if k>0 then(1,init (k-1)(fun i->if((i mod 2)==0)then -1 else 0))
else(k+1,init(-k-1)(fun i->if((i mod 2)==0)then 1 else 0));;

let rec cpr p q=
if p=0 then false else
if abs(p)=1 then true else
cpr (q mod p)p;;

let isgood p q=(p*q mod 2==0)&&cpr p q;;

let gi p0 q0=
let q=q0*(2-(q0 mod 2)) in
let p=if p0>=0 then p0 else q+p0 in
let rec f np=
let nnp=(p*np) mod q in 
if nnp==1 then np else f nnp in
let r=f 1 in
if q0 mod 2==0 then
if r<q0 then r else r-q
else
if r mod 2==0 then r else r-q0
;;


let gmod p q=(if p>0 then p else q+p)mod q;;

for i=0 to mq-1 do(
let q=i+1 in
for j=0 to 2*q-1 do(
let p=rv j in
if isgood p (i+1) then(
let np=gi p q in
let inp=iv np in
if inp<j then ca.(i).(j)<-E inp;
))done;)done;;

let prpoly ch p=let (k,h)=p in
for i=0 to length h-1 do(
if h.(i)<>(0,[||]) then(
fprintf ch "%d"(i+k);fprintf ch "  ";
let (nk,nh)=h.(i) in fprintf ch "%d"nk;fprintf ch "  ";
for j=0 to length nh-1 do(
if(j mod 2==0)then(
fprintf ch "%d"(nh.(j));fprintf ch " ";);
)done;
fprintf ch "\n";)
)done;;

let tb=Array.init mq(fun q->(init(2*(q+1))(fun _->ref B)));;

let (lb,rb)=(ref 0,ref 0);;

let ut1 p=
let (k,h)=p in
lb:=0;
rb:=length h-1;
while(!lb<length h)&&(h.(!lb)==0) do(lb:=!lb+1)done;
while(!rb>=0)&&(h.(!rb)==0) do(rb:=!rb-1)done;
if lb<=rb then
(k+(!lb),init(!rb-(!lb)+1)(fun i->h.(i+(!lb))))
else (0,[||])
;;

let ut2 p=
let (k,h)=p in (k,map ut1 h);;


for i=0 to mq-1 do(
let q=i+1 in 
for j=0 to 2*q-1 do(
let p=rv j in
if isgood p q then(
let b=(gd p q)*(m2 p) in
let spb=sp b in
if abs(p)==1 then
tb.(i).(j):=(G(-1,[|(b-1,[|-1;0;1|]);(0,[||]);spb|]))
else((
match ca.(i).(j) with
E k->tb.(i).(j):=!(tb.(i).(k))
|C->
let(p1,q1)=(np p q,abs(p)) in
match !(tb.(q1-1).(iv p1)) with
G h1->
let h2=(if abs(p1)==1 then
(0,[|(0,[|1|])|])
else(
let(p2,q2)=(np p1 q1,abs(p1)) in
match !(tb.(q2-1).(iv p2)) with
G h->h))in
let npoly=((sum(sum(+)0)(0,[||]))(prod(prod( * )(sum(+)0))(sum(sum(+)0)(0,[||]))(0,[|(b,[|1|])|])h2)(prod(prod( * )(sum(+)0))(sum(sum(+)0)(0,[||]))(1,[|spb|])h1))in
tb.(i).(j):=G(ut2 npoly);
))))done;
)done;;


print_string (String.concat""["\n";"homfly";(string_of_int mq);".txt";"\n"]);;

let ch=open_out (String.concat""["homfly";(string_of_int mq);".txt"]) in

for i=0 to mq-1 do(
let q=i+1 in 
for j=0 to 2*q-1 do(
let p=rv j in
if isgood p q then(
if ((ca.(i).(j)==C)||print_inv) then(
if ((p<0)||print_neg)then(
fprintf ch "%d"((-1)*p);fprintf ch "/";fprintf ch"%d" q;fprintf ch "\n";
match !(tb.(i).(iv(p*(if q mod 2==0 then -1 else 1)))) with 
B->print_string"!";
|G poly->prpoly ch(poly);
))))done;)done;

close_out ch;;