\\
\\ rm2(n): divide n by 2 if n is even
{rm2(n,t,u,v)=
if(mod(n,2)==mod(0,2),n/2,n)}
\\
\\ modd(n): the odd part of n; divide out any power of 2
{modd(n,t,u,v)=
while(mod(n,2)==mod(0,2),n=rm2(n));
n}
\\
\\ n32(n,l): compute iterations of odd(3*n+1) with output formatted in
\\    vectors of length l (default 10) when 0 < l <100, otherwise no
\\    output, returning the number of iterations to 1, *** assuming
\\    this number is finite, as conjectured.
{n32(n,l,t,u,v,w,x)=
if(l<0,l=100,);
if(l==0,l=10,);
v=vector(l,z,0);
u=1;
v[u]=n;
x=0;
t=n;
while(u >= l,\
  if(l < 100,print(v),);\
  v=vector(l,z,0);\
  u=0);
while(t > 1,\
  u=u+1;\
  x=x+1;\
  t=modd(3*t+1);\
  v[u]=t;\
  while(u >= l,\
    if(l < 100,print(v),);\
    v=vector(l,z,0);\
    u=0);\
  );
if((l < 100) && (u > 0), w=vector(u,j,v[j]);print(w),);
x}
\\
\\ n32ta(m) = greatest number of iterations of modd(3*x+1) to 1
\\    for all j from 1 to m
\\    print the location and the max no. and return the max no.
{n32ta(m,u,j,k,x)=
while(m>0,\
u=0;j=1;k=0;\
while(j<=m,x=n32(j,-1);if(x>u,u=x;k=j;,);j=j+1);\
print(k," ",u);\
m=0;\
);
u}
\\
\\ n32t(m) = greatest number of iterations of modd(3*x+1) to 1
\\    for all ODD j from 1 to m
\\    print the location and the max no. and return the max no.
{n32t(m,u,j,k,x)=
while(m>0,\
u=0;j=1;k=0;\
while(j<=m,x=n32(j,-1);if(x>u,u=x;k=j;,);j=j+2);\
print(k," ",u);\
m=0;\
);
u}
\\
\\ n32tai(m,n) = greatest number of iterations of modd(3*x+1) to 1
\\    for all j from m to n
\\    print the location and the max no. and return the max no.
{n32tai(m,n,u,j,k,x)=
while(m>0,\
u=0;j=m;k=0;\
while(j<=n,x=n32(j,-1);if(x>u,u=x;k=j;,);j=j+1);\
print(k," ",u);\
m=0;\
);
u}
\\
\\ n32ti(m,n) = greatest number of iterations of modd(3*x+1) to 1
\\    for all ODD j from m to n
\\    print the location and the max no. and return the max no.
{n32ti(m,n,u,j,k,x)=
while(m>0,\
u=0;if(mod(m,2)==mod(0,2),j=m+1,j=m);k=0;\
while(j<=n,x=n32(j,-1);if(x>u,u=x;k=j;,);j=j+2);\
print(k," ",u);\
m=0;\
);
u}
\\
\\ n32m(m,n,l): compute iterations of odd(3*n+1) with output reduced
\\    mod m  and formatted in vectors of length l (default 10)
\\    when 0 < l <100 (otherwise no output) returning the number of
\\    iterations to 1, assuming this number is finite, as conjectured.
\\    Note: the iteration process itself is NOT reduced mod m.
{n32m(m,n,l,t,u,v,vm,w,x)=
if(l<0,l=100,);
if(l==0,l=10,);
v=vector(l,z,0);
u=1;
v[u]=compo(mod(n,m),2);
x=0;
t=n;
while(u >= l,\
  if(l < 100,print(v),);\
  v=vector(l,z,0);\
  u=0);
while(t > 1,\
  u=u+1;\
  x=x+1;\
  t=modd(3*t+1);\
  v[u]=compo(mod(t,m),2);\
  while(u >= l,\
    if(l < 100,print(v),);\
    v=vector(l,z,0);\
    u=0);\
  );
if((l < 100) && (u > 0), w=vector(u,j,v[j]);print(w),);
x}
\\
\\  n32p2(a,b): compute the no. of iterations of modd(3*x+1) to 1
\\     starting with x = 2^j-1, as j varies from a to b, and return a
\\     vector consisting of the j value and the no. of iterations that
\\     is largest among these values of j.
{n32p2(a,b,m,k,u)=m=0;
k=0;
for(j=a,b,u=n32((2^j)-1,-1);if(u>m,m=u;k=j,));
v=[k,m];v}
\\
\\ f32p2(k,l): compute the no. of iterations of modd(3*x+1) to 1
\\    for x = 2^j-1, with j = 2*k-1 and j = 2*k, formatting the n32
\\    output in vectors of length l.
{f32p2(k,l)=print("---  ",k,"  ---");\
for(j=(2*k)-1,2*k,print(n32((2^j)-1,l)))}
\\
\\ rmf(p,x): remove all factors of p from x (return x or 0 if p <=1)
{rmf(p,x,t)=
if(p>1,t=x;while(mod(t,p)==0,t=t/p);t,0)}
\\
\\ fit(a,b,p,x): return  a*x+b  with the highest possible power of
\\    p  factored out
{fit(a,b,p,x)=
if(p!=0,
rmf(p,(a*x)+b)
,0)
}
\\
\\ isinvec(v,n): return position of  n  in vector  v,  otherwise  0.
{isinvec(v,n,k,j,s)=
k=0;
j=0;
s=matsize(v)[2];
while(j<s,j=j+1;if(n==v[j],k=j;j=s,));
k}
\\
\\ augvec(v,n): return vector of size one more than the size of  v
\\    with last component  n.
{augvec(v,n,s,w)=
s=matsize(v)[2];
w=vector(s+1,j,0);
for(j=1,s,w[j]=v[j]);
w[s+1]=n;
w}
\\
\\ nit(a,b,p,x,verb):
\\             iterate z -> (a*z+b)/(highest power of p)  from  x;
\\    stop either when  z < x  or a cycle is detected.
\\    Return last iterate.
\\    Use arg verb (optional) to receive iteration information every
\\       verb iterations.
\\\    Note:  k = number of iterations; one less than size of the
\\\           vector  m  after application of augvec
{nit(a,b,p,x,verb,k,t,l,m,cyc,j)=
m=vector(1,z,x);
k=1;
t=fit(a,b,p,x);
cyc=isinvec(m,t);
while((t>=x) && (cyc==0),\
 m=augvec(m,t);
 if(verb>2,\
  if(mod(k,verb)==0,\
   print("x = ",x," iterations ",k," size ",size(t));,),);\
 k=k+1;\
 t=fit(a,b,p,t);\
 cyc=isinvec(m,t);\
);\
if(cyc!=0,\
 print("Cycle: ",cyc-1,"  ",k);\
 j=0;\
 k=k+1;\
 v=vector(k-cyc,z,0);\
 while(j<(k-cyc),j=j+1;l=cyc+(j-1);v[j]=m[l]);\
 print(v);\
,\
 print("Smaller in ",k," iterations"));\
t
}
\\
\\ rmv(m,x): return the "part" of x that is coprime to m
{rmv(m,x,t,u,v,f,s)=
if(m>1,t=x;\
 u=gcd(m,x);\
 f=factor(u);\
 s=matsize(f)[1];\
 for(j=1,s,\
  t=rmf(f[j,1],t)\
  );\
  t\
,0)}
\\
\\ fiv(a,b,p,x): return (a*x+b) made coprime to  p
{fiv(a,b,p,x)=
if(p!=0,rmv(p,(a*x)+b),0)}
\\\
\\\ nivo(a,b,p,x): iterate z -> [(a*z+b) made coprime to  p]  from  x;
\\\    stop either when  z < x  or a cycle is detected.
\\\    This version does not detect cycles correctly.  See niv.
{nivo(a,b,p,x,k,t,l,m,mold,j)=
l=0;\
mold=x;\
m=x;\
k=1;\
t=fiv(a,b,p,x);\
if(t>m,m=t;l=k,);\
while((t>=x) && (t!=mold),\
 if(mold<m,mold=m,);\
 k=k+1;\
 if(mod(k,1000)==0,print("Iterations: ",k);print(log(t)),);\
 t=fiv(a,b,p,t);\
 if(t>m,m=t;l=k,);\
);\
if(t==mold,\
 print("Repetition: ",l,"  ",k,"  ",t);\
 j=0;\
 v=vector(k-l,y,0);\
 while(j<(k-l),j=j+1;t=fiv(a,b,p,t);v[j]=t);\
 print(v);\
,\
 print("Iterations: ",k));\
t
}
\\
\\ niv(a,b,p,x,verb): iterate z -> [ a*z+b made coprime to p] from x;
\\    stop either when  z < x  or a cycle is detected.
\\    Print cycle or no. of iterations to end.  Return last iterate.
\\    Use arg verb (optional) to receive iteration information every
\\       verb iterations.
\\\    Note:  k = number of iterations; one less than size of the
\\\           vector  m  after application of augvec
{niv(a,b,p,x,verb,k,t,l,m,cyc,j,v)=
m=vector(1,z,x);
k=1;
t=fiv(a,b,p,x);
cyc=isinvec(m,t);
while((t>=x) && (cyc==0),\
 m=augvec(m,t);\
 if(verb>2,\
  if(mod(k,verb)==0,\
   print("x = ",x," iterations ",k," size ",size(t));,),);\
 k=k+1;\
 t=fiv(a,b,p,t);\
 cyc=isinvec(m,t);\
);\
if(cyc!=0,\
 print("Cycle: ",cyc-1,"  ",k);\
 j=0;\
 k=k+1;\
 v=vector(k-cyc,z,0);\
 while(j<(k-cyc),j=j+1;l=cyc+(j-1);v[j]=m[l]);\
 print(v);\
,\
 print("Smaller in ",k," iterations"));\
t}
\\
\\ nivc(a,b,p,x,verb):
\\                 iterate z -> [ a*z+b made coprime to p]  from  x;
\\    Stop when a cycle is detected.
\\    Use arg verb (optional) to receive iteration information every
\\       verb iterations.
\\    Return vector with four components:
\\       (1) offset for cycle start,
\\       (2) offset for cycle end,
\\       (3) smallest entry in cycle,
\\       (4) the cycle (itself a vector).
\\\    Note:  k = number of iterations; one less than size of the
\\\           vector  m  after application of augvec
{nivc(a,b,p,x,verb,k,t,l,m,cyc,j,v,q)=
m=vector(1,z,x);
k=1;
t=fiv(a,b,p,x);
cyc=isinvec(m,t);
while(cyc==0,\
 m=augvec(m,t);\
 if(verb>2,\
  if(mod(k,verb)==0,\
   print("x = ",x," iterations ",k," size ",size(t));,),);\
 k=k+1;\
 t=fiv(a,b,p,t);\
 cyc=isinvec(m,t);\
);
j=0;
k=k+1;
v=vector(k-cyc,z,0);
while(j<(k-cyc),j=j+1;l=cyc+(j-1);v[j]=m[l]);
j=0;
q=v[1];
while(j<(k-cyc),j=j+1;if(v[j]<q,q=v[j],););
[cyc-1,k-1,q,v]}
\\
\\ nivf(a,b,p,x,width):
\\                   iterate z -> [ a*z+b made coprime to p]  from  x;
\\    stop when a cycle is detected.
\\    Print all iterates as generated.
\\       width (optional) -- number of iterates per line (default 10)
\\    Return vector with four components:
\\       (1) offset for cycle start,
\\       (2) offset for cycle end,
\\       (3) smallest entry in cycle,
\\       (4) the cycle (itself a vector).
\\\    Note:  k = number of iterations; one less than size of the
\\\           vector  m  after application of augvec
{nivf(a,b,p,x,wid,k,t,l,m,cyc,j,v,r,u,rp,q)=
if(wid<0,wid=100,);
if(wid==0,wid=10,);
r=vector(wid,z,0);
u=1;
r[u]=x;
m=vector(1,z,x);
while(u>=wid,\
 if(wid<100,print(r),);\
 r=vector(wid,z,0);\
 u=0;\
);
k=1;
u=u+1;
t=fiv(a,b,p,x);
r[u]=t;
while(u>=wid,\
 if(wid<100,print(r),);\
 r=vector(wid,z,0);\
 u=0;\
);
cyc=isinvec(m,t);
while(cyc==0,\
 m=augvec(m,t);\
 k=k+1;\
 u=u+1;\
 t=fiv(a,b,p,t);\
 r[u]=t;\
 while(u>=wid,\
  if(wid<100,print(r),);\
  r=vector(wid,z,0);\
  u=0;\
 );\
 cyc=isinvec(m,t);\
);
if((u>0) && (wid < 100), rp=vector(u,z,r[z]);print(rp),);
j=0;
k=k+1;
v=vector(k-cyc,z,0);
while(j<(k-cyc),j=j+1;l=cyc+(j-1);v[j]=m[l]);
j=0;
q=v[1];
while(j<(k-cyc),j=j+1;if(v[j]<q,q=v[j],););
[cyc-1,k-1,q,v]}
\\
\\ nivck(a,b,p,vec,ia,ib,verb):
\\                 iterate z -> [ a*z+b made coprime to p]  from  x
\\       for all  ia <= x <= ib  until a value in the vector  vec
\\       is found.  This is based on the premise that successive such
\\       iterations always lead to a cycle belonging to a finite set
\\       and that each such cycle has at least one representative
\\       member among the coordinates of  vec .  Note: it is not always
\\       the case that such iterations lead to a cycle.
\\    Stop when a cycle is detected.
\\    Use arg verb (optional) to receive iteration information every
\\       verb iterations.
\\    Return the maximum number of iterations required along the way.
\\\    Note:  k = number of iterations; one less than size of the
\\\           vector  m  after application of augvec
{nivck(a,b,p,vv,ia,ib,verb,x,k,t,m,l,cyc,s)=
m=0;
l=0;
s=0;
for(x=ia,ib,\
k=0;\
cyc=0;\
t=x;\
while(cyc==0,\
 k=k+1;\
 if(verb>2,\
  if(mod(k,verb)==0,\
   print("x = ",x," iterations ",k," size ",size(t));,),);\
 t=fiv(a,b,p,t);\
 cyc=isinvec(vv,t);\
);\
if(k>m,l=x;m=k;s=t,);\
);
[l,s,m]}
\\
\\ rfiv(a,b,p,x,r,w): iterate  [ z -> (az + b) made coprime to p ]
\\      starting with  x  for  r  iterations, printing  w  iterates
\\      per line, w = 0 or 1 (default 0 = no output); return the r-th
\\      iterate
\\
{rfiv(a,b,p,x,r,w,t)=t=x;
if(w>1,print(" z -> ",a,"z + ",b," made coprime to ",p);\
 print("   from ",x," for ",r," iterations"),);
if(w>0,for(j=1,r,t=fiv(a,b,p,t);print(t)),for(j=1,r,t=fiv(a,b,p,t)));
t}
\\
\\ nive(a,b,p,s,x,verb):
\\                 iterate z -> [ a*z+b made coprime to p]  from  x;
\\    Stop when the second occurrence of  s  is found.
\\    Use arg verb (optional) to receive iteration information every
\\       verb iterations.
\\    Return vector with five components:
\\       (1) offset for the first occurrence of  s,
\\       (2) offset for the second occurrence of  s,
\\       (3) smallest entry in cycle delimited by (1) and (2),
\\       (4) the cycle (itself a vector) beginning with  s,
\\       (5) the initial sequence of iterates up to the first
\\           occurrence of  s.
\\\    Note:  k = number of iterations; one less than size of the
\\\           vector  m  after application of augvec
{nive(a,b,p,s,x,verb,k,t,l1,l2,m1,m2,j,q)=
m1=vector(0,z,0);
k=1;
t=x;
while(t!=s,\
 if(verb>2,\
  if(mod(k,verb)==0,\
   print("x = ",x," iterations ",k," size ",size(t));,),);\
 m1=augvec(m1,t);\
 k=k+1;\
 t=fiv(a,b,p,t);\
);
l1=k-1;
m2=vector(1,z,s);
k=k+1;
t=fiv(a,b,p,t);
while(t!=s,\
 if(verb>2,\
  if(mod(k,verb)==0,\
   print("x = ",x," iterations ",k," size ",size(t));,),);\
 m2=augvec(m2,t);\
 k=k+1;\
 t=fiv(a,b,p,t);\
);
l2=k-1;
j=0;
l=l2-l1;
q=m2[1];
if(matsize(m2)[2]==l,\
while(j<l,j=j+1;if(m2[j]<q,q=m2[j],););[l1,l2,q,m2,m1]\
,\
print("Error in function nive"));}
\\
\\ cfmat(v,e) = the continued fraction (2 X 2) matrix of the vector  v
\\              with (constant) "numerators"  e  (default e = 1)
{cfmat(v,e,n,a,m,r,u)=
r=matsize(v)[1];
if(r!=1,print("First arg. must be a (row) vector."),\
if(e==0,e=1,);\
n=matsize(v)[2];\
m=[1,0;0,1];\
for(j=1,n,\
u=[v[j],e;1,0];\
m=m*u;\
);\
);
m}
\\
\\ matact(m,q) = the number obtained by letting the 2 X 2 invertible
\\              matrix  m  act projectively on the number q
{matact(m,q,t,u,v,w)=
u=0;
v=m[1,1]*q+m[1,2];
w=m[2,1]*q+m[2,2];
if(w!=0,u=v/w;u,[v,w])}
\\
\\ q2cf(q,e)= the continued fraction of the quadratic number  q  with
\\           constant numerator  e  (default 1)  returned as a vector
\\           with two components, each a vector, the first being the
\\           cycle and the second being the initial sequence.
\\
{q2cf(q,e,halt,pr,t,a,vq,va,k,l,v1,v2,m,j)=
if(e==0,e=1,);
if((e>0)&&(e!=1),print("q2cf: Warning -- e not 1");e=1,);
if((e<0)&&(e!=-1),print("q2cf: Warning -- e not -1");e=-1,);
vq=vector(0,z,0);
va=vector(0,z,0);
k=0;
t=q;
halt=0;
cyc=0;
pr=setprecision(-1);
setprecision(1000);
while((cyc==0)&&(halt==0),\
  k=k+1;\
  vq=augvec(vq,t);\
  if(e>0,a=floor(1.*t),a=-floor(-1.*t));\
  va=augvec(va,a);\
  t=t-a;\
  if(t==0,halt=1,t=e/t;cyc=isinvec(vq,t););\
);
if(halt==1,v1=[];v2=va,\
  k=k+1;\
  l=cyc;\
  v2=vector(l-1,z,va[z]);\
  m=k-l;\
  v1=vector(m,z,0);\
  j=0;\
  while(j<m,\
    j=j+1;\
    v1[j]=va[l+j-1];\
  );\
);
setprecision(pr+9);
[v1,v2]}
\\
\\ n32c(n) = the element of the quadratic field of discriminant  5
\\           corresponding to the init-orbit of  n  under the map
\\           x -> [ (3x + 1) made odd ]
\\    Note this quadratic number is always of the form (u + w)/v
\\    where  u  and  v  are integers and  w = quadgen(5).  If  m =
\\    [a,b;c,d] = cfmat(the orbit), then [u,v]~ = (-1)^k m [-c c+d]~,
\\    where  k  is the number of terms in the init-orbit.
{n32c(n,q,x,m,y)=
q=quadgen(5);
x=nive(3,1,2,1,n)[5];
m=cfmat(x);
y=matact(m,q);
y}
\\
\\ n32cv(n) = the element of the quadratic field of discriminant  5
\\            corresponding to the orbit of  n  under the map
\\            x -> [ (3x + 1) made odd ]
\\   n32cv differs from n32c in that it is quite verbose: it prints
\\   (1) the part of the orbit of  n  up to the cycle  [1]
\\           hypothesizing the eventual existence of this cycle
\\   (2) the unimodular integer matrix corresponding to this
\\           orbit regarded as the beginning of a continued fraction
\\   (3) the reciprocal of the "imaginary" part of the quadratic number
\\   (4) the number (3) times the "real" part of the quadratic number
{n32cv(n,q,x,m,y,u,v)=
q=quadgen(5);
x=nive(3,1,2,1,n)[5];
print(x);
m=cfmat(x);
print(m);
y=matact(m,q);
v=1/imag(y);
print(v);
u=v*real(y);
print(u);
y}
\\
\\ n32ci(n) = the vector of two integers  [u,v]  where  (u + w)/v  is
\\            the element of the quadratic field of discriminant  5
\\            corresponding to the orbit of  n  under the map
\\            x -> [ (3x + 1) made odd ] .    (See n32c.)
{n32ci(n,q,x,m,y,u,v,w,s)=
q=quadgen(5);
x=nive(3,1,2,1,n)[5];
m=cfmat(x);
y=matact(m,q);
v=1/imag(y);
u=v*real(y);
if(matsize(x)[2]%2==0,s=1,s=-1);
w=trans(m*[-m[2,1],m[2,1]+m[2,2]]~);
w=s*w;
if(w[1]!=u,print("u = ",u,"  w1 = ",w[1]),);
if(w[2]!=v,print("v = ",v,"  w2 = ",w[2]),);
w}
\\
\\ n32cb(n) = the triple consisting of the sign of the determinant and
\\            the bottom row of the unimodular integer matrix
\\            corresponding to the init-orbit of  n  under the map
\\            x -> [ (3x + 1) made odd ] .    (See n32c and cfmat.)
{n32cb(n,q,x,m,y,u,v,w,s)=
q=quadgen(5);
x=nive(3,1,2,1,n)[5];
m=cfmat(x);
if(matsize(x)[2]%2==0,s=1,s=-1);
w=[s,m[2,1],m[2,2]];
w}
\\\
\\\ cyc2qo(v) = the real quadratic irrational number whose continued
\\\          fraction is purely cyclic with cycle  v  (a vector of
\\\          positive integers).
{cyc2qo(v,pr,m,a,b,c,d,t,w,srd,x1,x2,n1,n2,x)=
pr=setprecision(-1);
setprecision(1000);
m=cfmat(v);
a=m[2,1];
b=m[2,2]-m[1,1];
c=-m[1,2];
d=b^2-4*(a*c);
w=quadgen(d);
srd=2*w-(d%2);
x1=(srd-b)/(2*a);
x2=(-srd-b)/(2*a);
n1=floor(1.*x1);
n2=floor(1.*x2);
x=0;
if(n1==v[1],x=x1,);
if(n2==v[1],x=x2,);
setprecision(pr+9);
x}
\\
\\ cyc2q(v,e) = the real quadratic irrational number whose continued
\\              fraction with constant numerator  e  (default 1) is
\\              purely cyclic with cycle  v  (a vector of positive
\\              integers).  The returned quadratic is expressed in
\\              terms of the quadgen for the number field
\\              discriminant.
{cyc2q(v,e,pr,m,a,b,c,d,lw,sw,srd,srsmd,smd,quo,x1,x2,n1,n2,x)=
if(e==0,e=1,);
if((e>0)&&(e!=1),print("cyc2q: Warning -- e not 1");e=1,);
if((e<0)&&(e!=-1),print("cyc2q: Warning -- e not -1");e=-1,);
pr=setprecision(-1);
setprecision(1000);
m=cfmat(v,e);
a=m[2,1];
b=m[2,2]-m[1,1];
c=-m[1,2];
d=b^2-4*(a*c);
if(d==0,print("cyc2q: Warning -- discriminant is 0."),);
if(d<0,print("cyc2q: Warning -- discriminant is negative."),);
lw=quadgen(d);
smd=discf(compo(lw,1));
quo=d/smd;
sw=quadgen(smd);
srsmd=2*sw-(smd%2);
srd=isqrt(quo)*srsmd;
x1=(srd-b)/(2*a);
x2=(-srd-b)/(2*a);
x=x1;
if(d>0,\
  x=0;\
  n1=e*floor(1.*e*x1);n2=e*floor(1.*e*x2);\
  if(n1==v[1],x=x1,);\
  if(n2==v[1],x=x2;print("cyc2q: second root"),);\
  if(x==0,\
    print("cyc2q: Warning -- cannot match floor of quadratic root.")\
  ,);\
,);
setprecision(pr+9);
x}
\\\
\\\ cf2qo(c,v) = the real quadratic irrational number having continued
\\\              fraction with init  v  and cycle  c.
{cf2qo(c,v,t,m,x)=
t=cyc2qo(c);
m=cfmat(v);
x=matact(m,t);
x}
\\
\\ cf2q([c,v],e) =
\\ cf2q(c,v,e)   = the real quadratic irrational number having continued
\\                 fraction with init  v  and cycle  c  and constant
\\                 "numerator"  e  (default 1).
{cf2q(c,v,e,t,m,x,cs,vt)=
vt=type(v);
if(vt!=17,\
  e=v;\
  v=c[2];\
  c=c[1];\
,);
if(e==0,e=1,);
cs=matsize(c)[2];
if(cs==0,\
  m=cfmat(v,e);\
  if(m[2,1]==0,x=[m[1,1],m[2,1]],x=m[1,1]/m[2,1]);\
,\
  t=cyc2q(c,e);\
  m=cfmat(v,e);\
  x=matact(m,t);\
);
x}
\\
\\ smquad(q) = the "small discriminant" representation of the quadratic
\\             number  q
\\
{smquad(q,verb,f,d,smd,e,le,w,lw,quo,sq)=
if(verb!=0,prqv(q);,);
f=compo(q,1);
if(verb!=0,print("   f: ",f);,);
d=disc(q);
if(verb!=0,print("   d: ",d);,);
smd=discf(f);
if(verb!=0,print(" smd: ",smd);,);
quo=isqrt(d/smd);
if(verb!=0,print(" quo: ",quo);,);
w=quadgen(smd);
if(verb!=0,print("   w: ",compo(w,1),"  ",[compo(w,2),compo(w,3)]);,);
e=smd%2;
if(verb!=0,print("   e: ",e);,);
le=d%2;
if(verb!=0,print("  le: ",le);,);
lw=(((2*w-e)*quo)+le)/2;
if(verb!=0,print("   W: ",compo(lw,1),"  ",[compo(lw,2),compo(lw,3)]);,);
sq=real(q)+(imag(q)*lw);
if(verb!=0,prqv(sq);,);
sq}
\\
\\ prqv(q):  print verbose information about a quadratic number
\\
{prqv(q,f,t,n,d)=
if(type(q)==8,\
  f=compo(q,1);\
  t=trace(q);\
  n=norm(q);\
  d=(t^2)-(4*n);\
  print("trace ",t," norm ",n," disc ",\
       d,"  ",q,"   ",compo(q,1)," (",disc(f),")");\
,\
  print("prqv: Arg is not a quadratic number.");
);
q}
\\
\\ cf2qv(c,v,e) = the real quadratic irrational number having continued
\\                fraction with init  v ,  cycle  c  and constant
\\                "numerator"  e  (default 1).
{cf2qv(c,v,e,x)=
if(e==0,e=1,);
x=cf2q(c,v,e);
prqv(x);
x}
\\ r2cf(q) = the vector representing the (finite) continued fraction
\\            expansion of the rational number  q  with constant
\\            "numerator"  e  (default 1)
{r2cf(x,e,v,a,b,w,q,r,t)=
if(e==0,e=1,);
t=type(x);
if(((t!=4)&&(t!=5)),\
  if(t==1,v=[x];,print("ratcf: Arg. is type ", t));\
,\
  v=vector(0,z,0);\
  a=compo(x,1);\
  b=compo(x,2);\
  while(b!=0,\
    w=divres(a,b)~;\
    q=w[1];\
    r=w[2];\
    v=augvec(v,q);\
    a=b;\
    b=r;);\
);
v}
\\
\\ vhead(v,n) = vector consisting of the first n coordinates of v
\\              If n >= size of  v  ,  then  vhead(v) = v .
{vhead(v,n,s,w)=
s=matsize(v)[2];
if(n>=s,w=v,w=vector(n,z,v[z]));
w}
\\
\\ vtail(v,n) = vector consisting of the last n coordinates of v
\\              If n >= size of  v  ,  then  vtail(v) = v .
{vtail(v,n,s,w)=
s=matsize(v)[2];
if(n>=s,w=v,w=vector(n,z,v[z+s-n]));
w}
\\
\\ nivef(a,b,p,s,x,width):
\\                   iterate z -> [ a*z+b made coprime to p]  from  x;
\\    Stop when the second occurrence of  s  is found.
\\    Print all iterates as generated.
\\       width (optional) -- number of iterates per line (default 10)
\\    Return vector with five components:
\\       (1) offset for the first occurrence of  s,
\\       (2) offset for the second occurrence of  s,
\\       (3) smallest entry in the cycle delimited by (1) and (2),
\\       (4) the cycle (itself a vector) beginning with  s,
\\       (5) the initial sequence of iterates up to the first  s.
\\\    Note:  k = number of iterations; one less than size of the
\\\           vector  m  after application of augvec
\\\           r = running vector for printing in width wid
\\\           t = current iterate
\\\           u = index for  r
\\\          l1 = location (k-value) of first  s
\\\          l2 = location (k-value) of second  s
\\\          m1 = init vector
\\\          m2 = cycle (vector)
\\\           q = min. entry in cycle
\\\           j = misc. loop index
\\\          rp = shortened  r  for last printing
{nivef(a,b,p,s,x,wid,j,k,l1,l2,m1,m2,q,r,rp,t,u,l)=
if(wid<0,wid=100,);
if(wid==0,wid=10,);
r=vector(wid,z,0);
t=x;
k=0;
u=1;
r[u]=t;
m1=vector(0,z,0);
while(u>=wid,\
 if(wid<100,print(r),);\
 r=vector(wid,z,0);\
 u=0;\
);
while(t!=s,\
 m1=augvec(m1,t);\
 k=k+1;\
 t=fiv(a,b,p,t);\
 u=u+1;\
 r[u]=t;\
 while(u>=wid,\
  if(wid<100,print(r),);\
  r=vector(wid,z,0);\
  u=0;\
 );\
);
l1=k;
m2=vector(1,z,s);
k=k+1;
t=fiv(a,b,p,t);
u=u+1;
r[u]=t;
while(u>=wid,\
 if(wid<100,print(r),);\
 r=vector(wid,z,0);\
 u=0;\
);
while(t!=s,\
 m2=augvec(m2,t);\
 k=k+1;\
 t=fiv(a,b,p,t);\
 u=u+1;\
 r[u]=t;\
 while(u>=wid,\
  if(wid<100,print(r),);\
  r=vector(wid,z,0);\
  u=0;\
 );\
);
l2=k;
if((u>0) && (wid < 100), rp=vector(u,z,r[z]);print(rp),);
j=0;
l=l2-l1;
q=m2[1];
if(matsize(m2)[2]==l,\
while(j<l,j=j+1;if(m2[j]<q,q=m2[j],););[l1,l2,q,m2,m1]\
,\
print("Error in function nivef"))}
\\
\\ nivqf(a,b,p,s,x,width):
\\                   iterate z -> [ a*z+b made coprime to p]  from  x;
\\    Stop when the second occurrence of  s  is found.
\\    Print all iterates as generated.
\\       width (optional) -- number of iterates per line (default 10)
\\    Return vector with three components:
\\       (1) offset for the first occurrence of  s,
\\       (2) offset for the second occurrence of  s,
\\       (3) smallest entry in the cycle delimited by (1) and (2),
\\\    Note:  k = number of iterations; one less than size of the
\\\           vector  m  after application of augvec
\\\           r = running vector for printing in width wid
\\\           t = current iterate
\\\           u = index for  r
\\\          l1 = location (k-value) of first  s
\\\          l2 = location (k-value) of second  s
\\\           q = min. entry in cycle
\\\           j = misc. loop index
\\\          rp = shortened  r  for last printing
{nivqf(a,b,p,s,x,wid,j,k,l1,l2,q,r,rp,t,u,l)=
if(wid<0,wid=100,);
if(wid==0,wid=10,);
r=vector(wid,z,0);
t=x;
k=0;
u=1;
r[u]=t;
while(u>=wid,\
 if(wid<100,print(r),);\
 r=vector(wid,z,0);\
 u=0;\
);
while(t!=s,\
 k=k+1;\
 t=fiv(a,b,p,t);\
 u=u+1;\
 r[u]=t;\
 while(u>=wid,\
  if(wid<100,print(r),);\
  r=vector(wid,z,0);\
  u=0;\
 );\
);
l1=k;
q=s;
k=k+1;
t=fiv(a,b,p,t);
if(t<q,q=t,);
u=u+1;
r[u]=t;
while(u>=wid,\
 if(wid<100,print(r),);\
 r=vector(wid,z,0);\
 u=0;\
);
while(t!=s,\
 k=k+1;\
 t=fiv(a,b,p,t);\
 if(t<q,q=t,);\
 u=u+1;\
 r[u]=t;\
 while(u>=wid,\
  if(wid<100,print(r),);\
  r=vector(wid,z,0);\
  u=0;\
 );\
);
l2=k;
if((u>0) && (wid < 100), rp=vector(u,z,r[z]);print(rp),);
[l1,l2,q]}
\\
\\ prodp(n) = the product of all primes less than or equal to  n
{prodp(n,t)=t=1;forprime(p=1,n,t=t*p);
  print("prodp(",n,") has size ",size(t));t;}

{fivnorm(a,b,p,x)=
if(p!=0,rmv(p,(a*x)+b),0)}
{fivquad(a,b,p,x)=
if(p!=0,rmv(p,a*(x^2)+b),0)}
{fivcube(a,b,p,x)=
if(p!=0,rmv(p,a*(x^3)+b),0)}
{fivseven(a,b,p,x)=
if(p!=0,rmv(p,a*(x^7)+b),0)}
{fivn(n,a,b,p,x)=
if(p!=0,rmv(p,a*(x^n)+b),0)}
{fivp(polyn,p,x)=0}

\\
\\ nivqs(a,b,p,s,x,width):
\\                   iterate z -> [ a*z+b made coprime to p]  from  x;
\\    Stop when the second occurrence of  s  is found.
\\    Print the SIZES of all iterates as generated.
\\       width (optional) -- number of iterates per line (default 10)
\\    Return vector with three components:
\\       (1) offset for the first occurrence of  s,
\\       (2) offset for the second occurrence of  s,
\\       (3) smallest entry in the cycle delimited by (1) and (2),
\\\    Note:  k = number of iterations; one less than size of the
\\\           vector  m  after application of augvec
\\\           r = running vector for printing in width wid
\\\           t = current iterate
\\\           u = index for  r
\\\          l1 = location (k-value) of first  s
\\\          l2 = location (k-value) of second  s
\\\           q = min. entry in cycle
\\\           j = misc. loop index
\\\          rp = shortened  r  for last printing
{nivqs(a,b,p,s,x,wid,j,k,l1,l2,q,r,rp,t,u,l)=
if(wid<0,wid=100,);
if(wid==0,wid=10,);
r=vector(wid,z,0);
t=x;
k=0;
u=1;
r[u]=size(t);
while(u>=wid,\
 if(wid<100,print(r),);\
 r=vector(wid,z,0);\
 u=0;\
);
while(t!=s,\
 k=k+1;\
 t=fiv(a,b,p,t);\
 u=u+1;\
 r[u]=size(t);\
 while(u>=wid,\
  if(wid<100,print(r),);\
  r=vector(wid,z,0);\
  u=0;\
 );\
);
l1=k;
q=s;
k=k+1;
t=fiv(a,b,p,t);
if(t<q,q=t,);
u=u+1;
r[u]=size(t);
while(u>=wid,\
 if(wid<100,print(r),);\
 r=vector(wid,z,0);\
 u=0;\
);
while(t!=s,\
 k=k+1;\
 t=fiv(a,b,p,t);\
 if(t<q,q=t,);\
 u=u+1;\
 r[u]=size(t);\
 while(u>=wid,\
  if(wid<100,print(r),);\
  r=vector(wid,z,0);\
  u=0;\
 );\
);
l2=k;
if((u>0) && (wid < 100), rp=vector(u,z,r[z]);print(rp),);
[l1,l2,q]}
\\
\\ nivfmod(a,b,p,m,x,width):
\\                   iterate z -> [ a*z+b made coprime to p]  from  x;
\\    stop when a cycle is detected.
\\    Print all iterates mod  m  as generated.
\\       width (optional) -- number of iterates per line (default 10)
\\    Return vector with four components:
\\       (1) offset for cycle start,
\\       (2) offset for cycle end,
\\       (3) smallest entry in cycle,
\\       (4) the cycle mod  m  (itself a vector).
\\\    Note:  k = number of iterations; one less than size of the
\\\           vector  m  after application of augvec
{nivfmod(a,b,p,m,x,wid,k,t,l,vm,cyc,j,v,r,u,rp,q,tm)=
if(wid<0,wid=100,);
if(wid==0,wid=10,);
r=vector(wid,z,0);
u=1;
t=x;
tm=t%m;
r[u]=tm;
vm=vector(1,z,t);
while(u>=wid,\
 if(wid<100,print(r),);\
 r=vector(wid,z,0);\
 u=0;\
);
k=1;
u=u+1;
t=fiv(a,b,p,t);
tm=t%m;
r[u]=tm;
while(u>=wid,\
 if(wid<100,print(r),);\
 r=vector(wid,z,0);\
 u=0;\
);
cyc=isinvec(vm,t);
while(cyc==0,\
 vm=augvec(vm,t);\
 k=k+1;\
 u=u+1;\
 t=fiv(a,b,p,t);\
 tm=t%m;\
 r[u]=tm;\
 while(u>=wid,\
  if(wid<100,print(r),);\
  r=vector(wid,z,0);\
  u=0;\
 );\
 cyc=isinvec(vm,t);\
);
if((u>0) && (wid < 100), rp=vector(u,z,r[z]);print(rp),);
j=0;
k=k+1;
v=vector(k-cyc,z,0);
while(j<(k-cyc),j=j+1;l=cyc+(j-1);v[j]=vm[l]);
j=0;
q=v[1];
while(j<(k-cyc),j=j+1;if(v[j]<q,q=v[j],););
[cyc-1,k-1,q,v%m]}
\\
\\ nivcv(a,b,p,v,w): print the terminating cycle (see nivc) when the
\\                   map  x --> [ (ax + b) made coprime to  p ] is
\\                   iterated beginning with the successive
\\                   coordinates of the vector  v .  A non-zero value
\\                   for  w  causes iteration counts to be printed
\\                   every  w  iterations.
{nivcv(a,b,p,v,w,m,t)=
m=matsize(v)[2];
for(j=1,m,t=v[j];print(t,"  ",nivc(a,b,p,t,w)))}
\\
\\ nivfm(a,b,p,x,num,width):
\\                   iterate z -> [ a*z+b made coprime to p]  from  x;
\\    stop when a cycle is detected or otherwise after  num  iterations.
\\    Print all iterates as generated.
\\      width (optional) (default 10) is the number of iterates per row
\\    Return vector with four components:
\\       (1) offset for cycle start,
\\       (2) offset for cycle end,
\\       (3) smallest entry in cycle,
\\       (4) the cycle (itself a vector).
\\    The returned values will be null if no cycle has been detected.
\\\    Note:  k = number of iterations; one less than size of the
\\\           vector  m  after application of augvec
{nivfm(a,b,p,x,num,wid,k,t,l,m,cyc,j,v,r,u,rp,q)=
if(wid<0,wid=100,);
if(wid==0,wid=10,);
r=vector(wid,z,0);
u=1;
r[u]=x;
m=vector(1,z,x);
while(u>=wid,\
 if(wid<100,print(r),);\
 r=vector(wid,z,0);\
 u=0;\
);
k=1;
u=u+1;
t=fiv(a,b,p,x);
r[u]=t;
while(u>=wid,\
 if(wid<100,print(r),);\
 r=vector(wid,z,0);\
 u=0;\
);
cyc=isinvec(m,t);
while(((cyc==0) && (k<num)),\
 m=augvec(m,t);\
 k=k+1;\
 u=u+1;\
 t=fiv(a,b,p,t);\
 r[u]=t;\
 while(u>=wid,\
  if(wid<100,print(r),);\
  r=vector(wid,z,0);\
  u=0;\
 );\
 cyc=isinvec(m,t);\
);
if((u>0) && (wid < 100), rp=vector(u,z,r[z]);print(rp),);
j=0;
k=k+1;
if(cyc!=0,\
 v=vector(k-cyc,z,0);\
 while(j<(k-cyc),j=j+1;l=cyc+(j-1);v[j]=m[l]);\
 j=0;\
 q=v[1];\
 while(j<(k-cyc),j=j+1;if(v[j]<q,q=v[j],););\
,\
q=0;\
v=vector(0,z,0);\
);
[cyc-1,k-1,q,v]
}
\\
\\ nivfmmod(a,b,p,n,x,num,width):
\\                   iterate z -> [ a*z+b made coprime to p]  from  x;
\\    stop when a cycle is detected or otherwise after  num  iterations.
\\    Print all iterates mod  n  as generated.
\\      width (optional) (default 10) is the number of iterates per row
\\    Return vector with four components:
\\       (1) offset for cycle start,
\\       (2) offset for cycle end,
\\       (3) smallest entry in cycle,
\\       (4) the cycle mod  n  (itself a vector).
\\    The returned values will be null if no cycle has been detected.
\\\    Note:  k = number of iterations; one less than size of the
\\\           vector  m  after application of augvec
{nivfmmod(a,b,p,n,x,num,wid,k,t,tn,l,m,cyc,j,v,r,u,rp,q)=
if(wid<0,wid=100,);
if(wid==0,wid=10,);
r=vector(wid,z,0);
u=1;
t=x;
tn=t%n;
r[u]=tn;
m=vector(1,z,x);
while(u>=wid,\
 if(wid<100,print(r),);\
 r=vector(wid,z,0);\
 u=0;\
);
k=1;
u=u+1;
t=fiv(a,b,p,t);
tn=t%n;
r[u]=tn;
while(u>=wid,\
 if(wid<100,print(r),);\
 r=vector(wid,z,0);\
 u=0;\
);
cyc=isinvec(m,t);
while(((cyc==0) && (k<num)),\
 m=augvec(m,t);\
 k=k+1;\
 u=u+1;\
 t=fiv(a,b,p,t);\
 tn=t%n;\
 r[u]=tn;\
 while(u>=wid,\
  if(wid<100,print(r),);\
  r=vector(wid,z,0);\
  u=0;\
 );\
 cyc=isinvec(m,t);\
);
if((u>0) && (wid < 100), rp=vector(u,z,r[z]);print(rp),);
j=0;
k=k+1;
if(cyc!=0,\
 v=vector(k-cyc,z,0);\
 while(j<(k-cyc),j=j+1;l=cyc+(j-1);v[j]=m[l]);\
 j=0;\
 q=v[1];\
 while(j<(k-cyc),j=j+1;if(v[j]<q,q=v[j],););\
,\
q=0;\
v=vector(0,z,0);\
);
[cyc-1,k-1,q,v%n]
}