{isreduced(a,b,c) = 
   if(b^2-4*a*c>=0 || a<0, 
      error("reduce: (a,b,c) must be positive definite."));
   if(!(abs(b)<=a && a<=c), return(0));
   if(abs(b)==a || a==c, return(b>=0));
   return(1);
}

{reduce(f) = 
   local(D, k, t, a,b,c); 
   a=f[1]; b=f[2]; c=f[3];
   D=b^2-4*a*c;
   if(D>=0 || a<0, error("reduce: (a,b,c) must be positive definite."));
   while(!isreduced(a,b,c),       \\ ! means ``not''
      if( c<a || (c==a && b<0), 
         b = -b; t = a; a = c; c = t,
      \\ else
         if (abs(b)>a || -b==a || a==c, 
            k = floor((a-b)/(2*a)); 
            b = b+2*k*a;
            c = (b^2-D)/(4*a);
         )
      )
   );
   return([a,b,c])
}

{reducedforms(D)=
   local(bound, forms, b, r);
   if (D > 0 || D%4 == 2 || D%4==3, error("Invalid discriminant"));
   bound = floor(-D/3);
   forms = [];
   for(a = 1, bound,
      for(c = 1, bound,
         if(3*a*c<=-D && issquare(4*a*c+D),
            b = floor(sqrt(4*a*c+D));
            r = reduce([a,b,c]);
            print1([a,b,c], " ----> ", r);
            if (gcd(r[1],gcd(r[2],r[3])) == 1,
               forms = setunion(forms,[r]); print(""),
               \\ else
               print ("  \t(not primitive)")
            )
         )
      )
   );
   return(eval(forms));    \\ eval gets rid of the annoying quotes.
}

{composition(f1, f2)=
   local(a1,b1,c1,a2,b2,c2,D,s,n,bz0,bz1,u,v,w);
   a1=f1[1]; b1=f1[2]; c1=f1[3];
   a2=f2[1]; b2=f2[2]; c2=f2[3];   
   D = b1^2 - 4*a1*c1;
   if(b2^2 - 4*a2*c2 != D, error("Forms must have the same discriminant."));
   s = (b1+b2)/2;
   n = (b1-b2)/2;
   bz0 = bezout(a1,a2);
   bz1 = bezout(bz0[3],s);
   u = bz1[1]*bz0[1];
   v = bz1[1]*bz0[2];
   w = bz1[2];
   d = bz1[3];
   d0 = gcd(gcd(gcd(d,c1),c2),n);
   a3 = d0*a1*a2/d^2;
   b3 = b2+2*a2*(v*(s-b2)-w*c2)/d;
   c3 = (b3^2-D)/(4*a3);
   f3 = reduce([a3,b3,c3]);
   return(f3);
}

{orderform(a, iden)=
         local(product, order);
         order = 1;
         product = a;
         if (a == iden, return(order), order++);
         while(order > 0,
            if ((product = composition(a, product)) == iden, 
            return(order), order++);
         )
}