# source:git/Singular/LIB/hyperel.lib@1a3911

spielwiese
Last change on this file since 1a3911 was 1a3911, checked in by Frank Seelisch <seelisch@…>, 14 years ago
removed some docu errors prior to release 3-1-0 git-svn-id: file:///usr/local/Singular/svn/trunk@11626 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100644`
File size: 24.4 KB
Line
1///////////////////////////////////////////////////////////////////////////////
2version="\$Id: hyperel.lib,v 1.2 2009-04-06 12:39:02 seelisch Exp \$";
3category="Teaching";
4info="
5LIBRARY:    hyperel.lib
6AUTHOR:     Markus Hochstetter, markushochstetter@gmx.de
7
8NOTE: The library provides procedures for computing with divisors in the
9      jacobian of hyperelliptic curves. In addition procedures are available
10      for computing the rational representation of divisors and vice versa.
11      The library is intended to be used for teaching and demonstrating
12      purposes but not for efficient computations.
13
14
15
16
17PROCEDURES:
18  ishyper(h,f)           test, if y^2+h(x)y=f(x) is hyperelliptic
19  isoncurve(P,h,f)       test, if point P is on C: y^2+h(x)y=f(x)
20  chinrestp(b,moduli)    compute polynom x, s.t. x=b[i] mod moduli[i]
21  norm(a,b,h,f)          norm of a(x)-b(x)y in IF[C]
22  multi(a,b,c,d,h,f)     (a(x)-b(x)y)*(c(x)-d(x)y)  in IF[C]
23  ratrep (P,h,f)         returns polynomials a,b, s.t. div(a,b)=P
24  divisor(a,b,h,f,[])    computes divisor of a(x)-b(x)y
25  gcddivisor(p,q)        gcd of the divisors p and q
26  semidiv(D,h,f)         semireduced divisor of the pair of polys D[1], D[2]
28  cantorred(D,h,f)       returns reduced divisor which is equivalent to D
29  double(D,h,f)          computes 2*D on y^2+h(x)y=f(x)
30  cantormult(m,D,h,f)    computes m*D on y^2+h(x)y=f(x)
31
32              [parameters in square brackets are optional]
33";
34///////////////////////////////////////////////////////////////////////////////
35
36
37//=============== Test, if a given curve is hyperelliptic =====================
38
39proc ishyper(poly h, poly f)
40"USAGE:   ishyper(h,f); h,f=poly
41RETURN:  1 if y^2+h(x)y=f(x) is hyperelliptic, 0 otherwise
42NOTE:    Tests, if y^2+h(x)y=f(x) is a hyperelliptic curve.
43         Curve is defined over basering. Additionally shows error-messages.
44EXAMPLE: example ishyper; shows an example
45"
46{
47  // constructing a copy of the basering (only variable x),
48  // with variables x,y.
49  def R=basering;
50  list l= ringlist(R);
51  list ll=l[2];
52  ll="x","y";
53  l[2]=ll;
54  intvec v= l[3][1][2];
55  v=v,1;
56  l[3][1][2]=v;
57  def s=ring(l);
58  setring s;
59
60  // test, if y^2 + hy - f is hyperelliptic.
61  int i=1;
62  poly h=imap(R,h);
63  poly f=imap(R,f);
64  poly F=y2 + h*y - f;
65  ideal I=F, diff(F,x) , diff(F,y);
66  ideal J=std(I);
67  if ( J != 1 )
68  {
69      i=0;
70      "The curve is singular!";
71  }
72  if ( deg(f) mod 2 != 1 )
73  {
74     i=0;
75      "The polynomial ",f," has even degree!";
76  }
77  if ( leadcoef(f) != 1 )
78  {
79     i=0;
80     "The polynomial ",f," is not monic!";
81  }
82  if ( 2*deg(h) > deg(f)-1 )
83  {
84    i=0;
85     "The polynomial ",h," has degree ",deg(h),"!";
86  }
87  setring(R);
88  return(i);
89}
90example
91{ "EXAMPLE:"; echo = 2;
92   ring R=7,x,dp;
93   // hyperelliptic curve y^2 + h*y = f
94   poly h=x;
95   poly f=x5+5x4+6x2+x+3;
96   ishyper(h,f);
97}
98
99
100//================= Test, if a given ponit is on the curve ====================
101
102proc isoncurve(list P, poly h, poly f)
103"USAGE:   isoncurve(P,h,f); h,f=poly; P=list
104RETURN:  1 or 0 (if P is on curve or not)
105NOTE:    Tests, if P=(P[1],P[2]) is on the hyperelliptic curve y^2+h(x)y=f(x).
106         Curve is defined over basering.
107EXAMPLE: example isoncurve; shows an example
108"
109{
110   if ( P[2]^2 + subst(h,var(1),P[1])*P[2] - subst(f,var(1),P[1]) == 0 )
111   {
112      return(1);
113   }
114   return(0);
115}
116example
117{ "EXAMPLE:"; echo = 2;
118   ring R=7,x,dp;
119   // hyperelliptic curve y^2 + h*y = f
120   poly h=x;
121   poly f=x5+5x4+6x2+x+3;
122   list P=2,3;
123   isoncurve(P,h,f);
124}
125
126
127//====================== Remainder of a polynomial division ===================
128
129proc divrem(f,g)
130"USAGE:   divrem(f,g);  f,g poly
131RETURN:  remainder of the division f/g
132NOTE:    Computes R, s.t. f=a*g + R, and deg(R) < deg(g)
133EXAMPLE: example divrem; shows an example
134"
135{
136   return(reduce(f,std(g)));
137}
138example
139{ "EXAMPLE:"; echo = 2;
140   ring R=0,x,dp;
141   divrem(x2+1,x2);
142}
143
144
145//================ chinese remainder theorem for polynomials ==================
146
147proc chinrestp(list b,list moduli)
148"USAGE:   chinrestp(b,moduli); moduli, b, moduli=list of polynomials
149RETURN:  poly x, s.t. x= b[i] mod moduli[i]
150NOTE:    chinese remainder theorem for polynomials
151EXAMPLE: example chinrestp; shows an example
152"
153{
154   int i;
155   int n=size(moduli);
156   poly M=1;
157   for(i=1;i<=n;i++)
158   {
159      M=M*moduli[i];
160   }
161   list m;
162   for(i=1;i<=n;i++)
163   {
164      m[i]=M/moduli[i];
165   }
166   list y;
167   for(i=1;i<=n;i++)
168   {
169      y[i]= extgcd(moduli[i],m[i])[3];
170   }
171   poly B=0;
172   for(i=1;i<=n;i++)
173   {
174      B=B+y[i]*m[i]*b[i];
175   }
176   B=divrem(B,M);
177   return(B);
178}
179example
180{ "EXAMPLE:"; echo = 2;
181   ring R=7,x,dp;
182   list b=3x-4, -3x2+1, 1, 4;
183   list moduli=(x-2)^2, (x-5)^3, x-1, x-6;
184   chinrestp(b,moduli);
185}
186
187
188//========================= norm of a polynomial ===============================
189
190proc norm(poly a, poly b, poly h, poly f)
191"USAGE:   norm(a,b,h,f);
192RETURN:  norm of a(x)-b(x)y in IF[C]
193NOTE:    The norm is a polynomial in just one variable.
194         Curve C: y^2+h(x)y=f(x) is defined over basering.
195EXAMPLE: example norm; shows an example
196"
197{
198   poly n=a^2+a*b*h-b^2*f;
199   return(n);
200}
201example
202{ "EXAMPLE:"; echo = 2;
203   ring R=7,x,dp;
204   // hyperelliptic curve y^2 + h*y = f
205   poly h=x;
206   poly f=x5+5x4+6x2+x+3;
207   poly a=x2+1;
208   poly b=x;
209   norm(a,b,h,f);
210}
211
212
213//========== multiplikation of polynomials in the coordinate ring =============
214
215proc multi(poly a, poly b, poly c, poly d, poly h, poly f)
216"USAGE:   multi(a,b,c,d,h,f);
217RETURN:  list L with L[1]-L[2]y=(a(x)-b(x)y)*(c(x)-d(x)y) in IF[C]
218NOTE:    Curve C: y^2+h(x)y=f(x) is defined over basering.
219EXAMPLE: example multi; shows an example
220"
221{
222   poly A=a*c + b*d*f;
223   poly B=b*c +a*d + b*h*d;
224   return (list(A,B));
225}
226example
227{ "EXAMPLE:"; echo = 2;
228   ring R=7,x,dp;
229   poly h=x;
230   poly f=x5+5x4+6x2+x+3;
231   // hyperelliptic curve y^2 + h*y = f
232   poly a=x2+1;
233   poly b=x;
234   poly c=5;
235   poly d=-x;
236   multi(a,b,c,d,h,f);
237}
238
239
240//================== polynomial expansion around a point ========================
241
242proc darst(list P,int k, poly h, poly f)
243"USAGE:   darst(P,k,h,f);
244RETURN:  list c of length k
245NOTE:    expansion around point P in IF[C], s.t.
246         y=c[1]+c[2]*(x-P[1]) +...+c[k]*(x-P[1])^k-1 + rest.
247         Curve C:y^2+h(x)y=f(x) is defined over basering.
248EXAMPLE: example darst; shows an example
249"
250{
251
252   if ( P[2] == -P[2]- subst(h,var(1),P[1]))
253   {
254    ERROR("no special points allowed");
255   }
256   list c;
257   list r;
258   list n;
259   poly N;
260   c[1]=P[2];
261   r[1]=list(0,-1,1,0);
262   poly r1,r2,r3,r4;
263   // rational function are represented as (r1 - r2*y) / (r3 - r4*y)
264
265   for (int i=1; i<k ; i++)
266   {
267      r1=r[i][1]-c[i]*r[i][3];
268      r2=r[i][2]-c[i]*r[i][4];
269      r3=r[i][3];
270      r4=r[i][4];
271      n=multi(r3,r4,r1+r2*h,-r2,h,f);
272      N=r1*r1 + r1*r2*h-r2*r2*f;
273      r[i+1]=list(N/(var(1)-P[1]),0,n[1],n[2]);
274      while ((divrem(r[i+1][1],var(1)-P[1]) ==0) and (divrem(r[i+1][2],var(1)-P[1]) ==0) and (divrem(r[i+1][3],var(1)-P[1]) ==0) and (divrem(r[i+1][4],var(1)-P[1]) ==0))
275      {
276          // reducing the rationl function
277          //(r[i+1][1] - r[i+1][2]*y)/(r[i+1][3] - r[i+1][4]) , otherwise there
278          // could be a pole, s.t. conditions are not fulfilled.
279          r[i+1][1]=(r[i+1][1]) / (var(1)-P[1]);
280          r[i+1][2]=(r[i+1][2]) / (var(1)-P[1]);
281          r[i+1][3]=(r[i+1][3]) / (var(1)-P[1]);
282          r[i+1][4]=(r[i+1][4]) / (var(1)-P[1]);
283      }
284      c[i+1]=(subst(r[i+1][1],var(1),P[1]) - subst(r[i+1][2],var(1),P[1])*P[2]) / (subst(r[i+1][3],var(1),P[1]) - subst(r[i+1][4],var(1),P[1])*P[2]);
285   }
286   return(c);
287}
288example
289{ "EXAMPLE:"; echo = 2;
290   ring R=7,x,dp;
291   // hyperelliptic curve y^2 + h*y = f
292   poly h=x;
293   poly f=x5+5x4+6x2+x+3;
294   list P=5,3;
295   darst(P,3,h,f);
296}
297
298
299//================ rational representation of a divisor =======================
300
301proc ratrep1 (list P, poly h, poly f)
302"USAGE:   ratrep1(P,k,h,f);
303RETURN:  list (a,b)
304NOTE:    Important: P has to be semireduced!
305         Computes rational representation of the divisor
306         P[1][3]*(P[1][1], P[1][2]) +...+ P[sizeof(P)][3]*
307         *(P[sizeof(P)][1], P[sizeof(P)][2]) - (*)infty=div(a,b)
308         Divisor P has to be semireduced.
309         Curve C:y^2+h(x)y=f(x) is defined over basering.
311EXAMPLE: example ratrep1; shows an example
312"
313{
314   poly a=1;
315   list b;
316   list m;
317   list koef;
318   int k;
319
320   // Determination of the polynomial b[i] for each point using procedure darst
321   for (int i=1 ; i<= size(P); i++)
322   {
323      a=a*(var(1)-P[i][1])^(P[i][3]);    // computing polynomial a
324      m[i]=(var(1)-P[i][1])^(P[i][3]);
325      b[i]=P[i][2];
326      k=1;
327
328      while (divrem(b[i]*b[i] + b[i] *h - f,(var(1)-P[i][1])^(P[i][3])) != 0)
329      {
330          k=k+1;   // b[i]=P[i][2];
331          koef=darst(list (P[i][1],P[i][2]), k, h,f);
332          // could be improved, if one doesn't compute list coef completely new
333          // every time
334          b[i]=b[i]+ koef[k]*(var(1)-P[i][1])^(k-1);
335      }
336   }
337   // Return polynomial a and b. Polynomial b is solution of the congruencies
338   // b[i] mod m[i] .
339   return(list(a,chinrestp(b,m)));
340
341}
342example
343{ "EXAMPLE:"; echo = 2;
344   ring R=7,x,dp;
345   // hyperelliptic curve y^2 + h*y = f
346   poly h=x;
347   poly f=x5+5x4+6x2+x+3;
348   //divisor P
349   list P=list(-1,-3,1),list(1,1,1);
350   ratrep1(P,h,f);
351}
352
353
354//================ rational representation of a divisor =======================
355
356proc ratrep (list P, poly h, poly f)
357"USAGE:   ratrep(P,k,h,f);
358RETURN:  list (a,b)
359NOTE:    Importatnt: P has to be semireduced!
360         Computes rational representation of the divisor
361         P[1][3]*(P[1][1], P[1][2]) +...+ P[sizeof(P)][3]*
362         *(P[sizeof(P)][1], P[sizeof(P)][2]) - (*)infty=div(a,b)
363         Divisor P has to be semireduced.
364         Curve C:y^2+h(x)y=f(x) is defined over basering.
365         Works faster than ratrep1.
367EXAMPLE: example ratrep; shows an example
368"
369{
370   poly a=1;
371   list b;
372   list m;
373   list koef;
374   int k;
375
376   poly c;
377   list r;
378   list n;
379   poly Norm;
380   poly r1,r2,r3,r4;
381
382   // Determination of the polynomial b[i] for each point using procedure darst
383   for (int i=1 ; i<= size(P); i++)
384   {
385      a=a*(var(1)-P[i][1])^(P[i][3]);     // computing polynomial a
386      m[i]=(var(1)-P[i][1])^(P[i][3]);
387      b[i]=P[i][2];
388      k=1;
389      c=P[i][2];
390      r=0,-1,1,0;
391      while (divrem(b[i]*b[i] + b[i] *h - f,(var(1)-P[i][1])^(P[i][3])) != 0)
392      {
393          k=k+1;
394          // here, the procedure darst was integrateg. In every pass a new
395          // coefficient c[i] is determined.
396          r1=r[1]-c*r[3];
397          r2=r[2]-c*r[4];
398          r3=r[3];
399          r4=r[4];
400          n=multi(r3,r4,r1+r2*h,-r2,h,f);
401          Norm=r1*r1 + r1*r2*h-r2*r2*f;
402          r=list(Norm/(var(1)-P[i][1]),0,n[1],n[2]);
403          while ((divrem(r[1],var(1)-P[i][1]) ==0) and (divrem(r[2],var(1)-P[i][1]) ==0) and (divrem(r[3],var(1)-P[i][1]) ==0) and (divrem(r[4],var(1)-P[i][1]) ==0))
404          {
405
406           // reducing the rationl function
407           // (r[1]-r[2]y)/(r[3]-r[4]y) , otherwise there
408           // could be a pole, s.t. conditions are not fulfilled.
409              r[1]=(r[1]) / (var(1)-P[i][1]);
410              r[2]=(r[2]) / (var(1)-P[i][1]);
411              r[3]=(r[3]) / (var(1)-P[i][1]);
412              r[4]=(r[4]) / (var(1)-P[i][1]);
413          }
414          c=(subst(r[1],var(1),P[i][1]) - subst(r[2],var(1),P[i][1])*P[i][2]) / (subst(r[3],var(1),P[i][1]) - subst(r[4],var(1),P[i][1])*P[i][2]);
415          b[i]=b[i]+ c*(var(1)-P[i][1])^(k-1);
416      }
417   }
418   // return polynomial a and b. Polynomial b is solution of the congruencies
419   // b[i] mod m[i] .
420   return(list(a,chinrestp(b,m)));
421}
422example
423{ "EXAMPLE:"; echo = 2;
424   ring R=7,x,dp;
425   // hyperelliptic curve y^2 + h*y = f
426   poly h=x;
427   poly f=x5+5x4+6x2+x+3;
428   //Divisor P
429   list P=list(-1,-3,1),list(1,1,1);
430   ratrep(P,h,f);
431}
432
433
434//============== Order of a zero in a polynomial ==============================
435
436proc ordnung(poly x0 , poly g)
437"USAGE:   ordnung(x0,g);
438RETURN:  int i
439NOTE:    i is maximal, s.t. (x-x0)^i divides g
440EXAMPLE: example ordnung; shows an example
441"
442{
443   poly gg=g;
444   int i;
445   while ( divrem(gg,var(1)-x0) ==0 )
446   {
447      i=i+1;
448      gg=gg/(var(1)-x0);
449   }
450   return(i);
451}
452example
453{ "EXAMPLE:"; echo = 2;
454   ring R=0,x,dp;
455   poly g=(x-5)^7*(x-3)^2;
456   number x0=5;
457   ordnung(x0,g);
458}
459
460
461//================== divisor of a polynomial function =========================
462
463proc divisor(poly a, poly b, poly h, poly f, list #)
464"USAGE:   divisor(a,b,h,f); optional: divisor(a,b,h,f,s); s=0,1
465RETURN:  list P
466NOTE:    P[1][3]*(P[1][1], P[1][2]) +...+ P[size(P)][3]*
467         *(P[size(P)][1], P[size(P)][2]) - (*)infty=div(a(x)-b(x)y)
468         if there is an optional parameter s!=0, then divisor additonally
469         returns a parameter, which says, whether irreducible polynomials
470         occured during computations or not. Otherwise only warnings are
471         displayed on the monitor. For s=0 nothing happens.
472         Curve C: y^2+h(x)y=f(x) is defined over basering.
473EXAMPLE: example divisor; shows an example
474"
475{
476   list p;
477   int j;
478   poly x0;
479   list y;
480   list fa=factorize(gcd(a,b));  // wanted: common roots of a and b
481   poly Norm=norm(a,b,h,f);
482
483   int s;
484   int irred=0;
485   if (size(#)>0)
486   {
487      s=#[1];
488   }
489   else
490   {
491      s=0;
492   }
493
494   for (int i=2; i<=size(fa[1]) ; i++)
495   {
496      // searching roots by finding polynomials of degree 1
497      if ( deg(fa[1][i]) !=1 )
498      {
499         if (s==0)
500         {
501         "WARNIG: ", fa[1][i], "is irreducible over this field !";
502         }
503         else
504         {
505            irred=1;
506         }
507      }
508      else
509      {
510         x0=var(1) - fa[1][i];
511         // finding the y-coordinates; max. 2
512         y= factorize(var(1)^2 + var(1)*subst(h,var(1),x0) - subst(f,var(1),x0));
513         if ( deg(y[1][2]) == 1)
514         // if root belongs to point on curve, then...
515         {
516            // compute order of a-b*y in the founded point
517            j=j+1;
518            p[j]=list(x0,var(1)-y[1][2],fa[2][i]);
519            if ( y[2][2]== 1)          // ordinary point
520            {
521               j=j+1;
522               p[j]=list(x0 , var(1)-y[1][3] , fa[2][i] );
523               if (a/(var(1)-x0)^(fa[2][i]) - b/(var(1)-x0)^(fa[2][i]) * p[j][2] ==0 )
524               {
525                  p[j][3]= p[j][3] + ordnung(x0,norm(a/(var(1)-x0)^(fa[2][i]) , b/(var(1)-x0)^(fa[2][i]),h,f));
526               }
527               if (a/(var(1)-x0)^(fa[2][i]) - b/(var(1)-x0)^(fa[2][i]) * p[j-1][2] ==0 )
528               {
529               p[j-1][3]=p[j-1][3] + ordnung(x0,norm(a/(var(1)-x0)^(fa[2][i]) , b/(var(1)-x0)^(fa[2][i]),h,f));
530               }
531            }
532            else        // special point
533            {
534               p[j][3]=p[j][3] *2  ;
535               if (a/(var(1)-x0)^(fa[2][i]) - b/(var(1)-x0)^(fa[2][i]) * p[j][2] ==0 )
536               {
537                  p[j][3]=p[j][3] + ordnung(x0,norm(a/(var(1)-x0)^(fa[2][i]) , b/(var(1)-x0)^(fa[2][i]),h,f));
538               }
539            }
540
541         }
542         // Norm of a-b*y is reduced by common root of a and b
543         // (is worked off)
544         Norm = Norm/((var(1)-x0)^(ordnung(x0,Norm)));
545      }
546   }
547
548   // some points are still missing; points for which a and b have no common
549   // roots, but norm(a-b*Y)=0 .
550   fa=factorize(Norm);
551   for ( i=2 ; i<=size(fa[1]) ; i++)
552   {
553      if ( deg(fa[1][i]) !=1)
554      {
555         if (s==0)
556         {
557         "WARNING: ", fa[1][i], "is irreducible over this field !";
558         }
559         else
560         {
561            irred=1;
562         }
563      }
564      else
565      {
566         x0=var(1) - fa[1][i];
567         y= factorize(var(1)^2 + var(1)*subst(h,var(1),x0) - subst(f,var(1),x0));
568         if ( deg(y[1][2]) == 1)
569         // if root belongs to point on curve, then...
570         {
571            if (subst(a,var(1),x0)- subst(b,var(1),x0)* (var(1)-y[1][2]) ==0)
572            {
573               p[size(p)+1]=list(x0,var(1)-y[1][2], ordnung(x0,Norm,h,f));
574            }
575            if ( y[2][2]== 1)     // ordinary point
576            {
577               if (subst(a,var(1),x0)- subst(b,var(1),x0)* (var(1)-y[1][3]) ==0)
578               {
579               p[size(p)+1]=list(x0 , var(1)-y[1][3] , ordnung(x0,Norm,h,f));
580               }
581            }
582         }
583      }
584   }
585   if (s==0)
586   {
587      return(p);
588   }
589   return(p,irred);
590}
591example
592{ "EXAMPLE:"; echo = 2;
593   ring R=7,x,dp;
594   // hyperelliptic curve y^2 + h*y = f
595   poly h=x;
596   poly f=x5+5x4+6x2+x+3;
597   poly a=(x-1)^2*(x-6);
598   poly b=0;
599   divisor(a,b,h,f,1);
600}
601
602
603//===================== gcd of two divisors ===================================
604
605proc gcddivisor(list p, list q)
606"USAGE:   gcddivisor(p,q);
607RETURN:  list P
608NOTE:    gcd of two divisors
609EXAMPLE: example gcddivisor; shows an example
610"
611{
612   list e;
613   int i,j;
614   for (i=1 ; i<= size(p) ; i++)
615   {
616      for (j=1 ; j<= size(q) ; j++)
617      {
618         if ( p[i][1] == q[j][1] and p[i][2] == q[j][2])
619         {
620            if ( p[i][3] <= q[j][3] )
621            {
622               e[size(e)+1]= list (p[i][1] , p[i][2] , p[i][3]);
623            }
624            else
625            {
626               e[size(e)+1]= list (q[j][1] , q[j][2] , q[j][3]);
627            }
628         }
629      }
630   }
631   return(e);
632}
633example
634{ "EXAMPLE:"; echo = 2;
635   ring R=7,x,dp;
636   // hyperelliptic curve y^2 + h*y = f
637   poly h=x;
638   poly f=x5+5x4+6x2+x+3;
639   // two divisors
640   list p=list(-1,-3,1),list(1,1,2);
641   list q=list(1,1,1),list(2,2,1);
642   gcddivisor(p,q);
643}
644
645
646//========== semireduced divisor from rational representation=================
647
648proc semidiv(list D,poly h, poly f)
649"USAGE:   semidiv(D,h,f);
650RETURN:  list P
651NOTE:    important: Divisor D has to be semireduced!
652         Computes semireduced divisor P[1][3]*(P[1][1], P[1][2]) +...+ P[size(P)][3]*
653         *(P[size(P)][1], P[size(P)][2]) - (*)infty=div(D[1],D[2])@*
654         Curve C:y^2+h(x)y=f(x) is defined over basering.
655EXAMPLE: example semidiv; shows an example
656"
657{
658   if ( deg(D[2]) >= deg(D[1]) or divrem(D[2]^2+D[2]*h-f,D[1]) != 0 )
659   {
660      ERROR("Pair of polynomials doesn't belong to semireduced divisor!");
661   }
662   list D1,D2;
663   int s1,s2;
664   D1,s1=divisor(D[1],0,h,f,1);
665   D2,s2=divisor(D[2],1,h,f,1);
666
667   // Only if irreducible polynomials occured in D1 !and! D2 a warning
668   // is necessary.
669   if (s1==1 and s2==1)
670   {
671      "Attention:
673   }
674   return(gcddivisor(D1,D2));
675}
676example
677{ "EXAMPLE:"; echo = 2;
678   ring R=7,x,dp;
679   // hyperelliptic curve y^2 + h*y = f
680   poly h=x;
681   poly f=x5+5x4+6x2+x+3;
682   // Divisor
683   list D=x2-1,2x-1;
684   semidiv(D,h,f);
685}
686
687
688//=============== Cantor's algorithm - composition ============================
689
690proc cantoradd(list D, list Q, poly h, poly f)
692RETURN:  list P
693NOTE:    Cantor's Algorithm - composition
694         important: D and Q have to be semireduced!
695         Computes semireduced divisor div(P[1],P[2])= div(D[1],D[2]) + div(Q[1],Q[2])
696         The divisors are defined over the basering.
697         Curve C: y^2+h(x)y=f(x) is defined over the basering.
698EXAMPLE: example cantoradd; shows an example
699"
700{
701   poly a;
702   poly b;
703   list e=extgcd(D[1],Q[1]);
704   if ( e[1]==1 )
705   {
706      a=D[1]*Q[1];
707      b=divrem( e[2]*D[1]*Q[2] + e[3]*Q[1]*D[2] ,a);
708      return(list(a,b));
709   }
710   list c=extgcd(e[1],D[2]+Q[2]+h);
711   poly s1=e[2]*c[2];
712   poly s2=c[2]*e[3];
713   poly s3=c[3];
714   a=D[1]*Q[1]/c[1]^2;
715   b=divrem((s1*D[1]*Q[2] + s2*Q[1]*D[2] + s3*(D[2]*Q[2] + f))/c[1],a);
716   return(list(a,b));
717}
718example
719{ "EXAMPLE:"; echo = 2;
720   ring R=7,x,dp;
721   // hyperelliptic curve y^2 + h*y = f
722   poly h=x;
723   poly f=x5+5x4+6x2+x+3;
724   // two divisors in rational representation
725   list D=x2-1,2x-1;
726   list Q=x2-3x+2,-3x+1;
728}
729
730
731//==================== Cantor's algorithm - reduction =========================
732
733proc cantorred(list D,poly h,poly f)
734"USAGE:   cantorred(D,h,f);
735RETURN:  list N
736NOTE:    Cantor's algorithm - reduction.
737         important: Divisor D has to be semireduced!
738         Computes reduced divisor div(N[1],N[2])= div(D[1],D[2]).@*
739         The divisors are defined over the basering.
740         Curve C: y^2+h(x)y=f(x) is defined over the basering.
741EXAMPLE: example cantorred; shows an example
742"
743{
744   list N=D;
745   while ( 2*deg(N[1]) > deg(f)-1 )
746   {
747      N[1]=(f - N[2]*h - N[2]^2)/N[1];
748      N[2]=divrem(-h-N[2],N[1]);
749   }
751   return(N);
752}
753example
754{ "EXAMPLE:"; echo = 2;
755   ring R=7,x,dp;
756   // hyperelliptic curve y^2 + h*y = f
757   poly h=x;
758   poly f=x5+5x4+6x2+x+3;
759   // semireduced divisor
760   list D=2x4+3x3-3x-2, -x3-2x2+3x+1;
761   cantorred(D,h,f);
762}
763
764
765//================= doubling a semireduced divisor ============================
766
767proc double(list D, poly h, poly f)
768"USAGE:   double(D,h,f);
769RETURN:  list Q=2*D
770NOTE:    important: Divisor D has to be semireduced!
771         Special case of Cantor's algorithm.
772         Computes reduced divisor div(Q[1],Q[2])= 2*div(D[1],D[2]).@*
773         The divisors are defined over the basering.
774         Curve C:y^2+h(x)y=f(x) is defined over the basering.
775EXAMPLE: example double; shows an example
776"
777{
778   list c=extgcd(D[1], 2*D[2] + h);
779   poly a=D[1]*D[1]/c[1]^2;
780   poly b=divrem((c[2]*D[1]*D[2] + c[3]*(D[2]*D[2] + f))/c[1],a);
781   return(cantorred(list(a,b),h,f));
782}
783example
784{ "EXAMPLE:"; echo = 2;
785   ring R=7,x,dp;
786   // hyperelliptic curve y^2 + h*y = f
787   poly h=x;
788   poly f=x5+5x4+6x2+x+3;
789   // reduced divisor
790   list D=x2-1,2x-1;
791   double(D,h,f);
792}
793
794
795//================ multiples of a semireduced divisor =========================
796
797proc cantormult(int m, list D, poly h, poly f)
798"USAGE:   cantormult(m,D,h,f);
799RETURN:  list res=m*D
800NOTE:    important: Divisor D has to be semireduced!
801         Uses repeated doublings for a faster computation
802         of the reduced divisor m*D.
803         Attention: Factor m=int, this means bounded.
804         For m<0 the inverse of m*D is returned.
805         The divisors are defined over the basering.
806         Curve C: y^2+h(x)y=f(x) is defined over the basering.
807EXAMPLE: example cantormult; shows an example
808"
809{
810   list res=1,0;
811   list bas=D;
812   int exp=m;
813   if (exp==0) { return(list(1,0)); }
814   if (exp==1) { return(D); }
815   if (exp==-1) { return(list(D[1],-D[2]-h)) ; }
816   if ( exp < 0)
817   {
818      exp=-exp;
819   }
820   while ( exp > 0 )
821   {
822     if ( (exp mod 2) !=0 )
823     {
825       exp=exp-1;
826     }
827     bas=double(bas,h,f);
828     exp=exp/2;
829   }
830   if ( m < 0 )
831   {
832      res[2]=-res[2]-h;
833   }
834   return(res);
835}
836example
837{ "EXAMPLE:"; echo = 2;
838   ring R=7,x,dp;
839   // hyperelliptic curve y^2 + h*y = f
840   poly h=x;
841   poly f=x5+5x4+6x2+x+3;
842   // reduced divisor
843   list D=x2-1,2x-1;
844   cantormult(34,D,h,f);
845}
846
847
848/*
849//=============================================================================
850//   In the following you find a large example, which demostrates the use of
851//   the most important procedures.
852//=============================================================================
853//---- field with 2^5=32 elements ----
854ring r=(2,a),x,dp;
855minpoly=a5+a2+1;
856
857//---- hyperelliptic curve y^2 + hy = f ----
858poly h=x2+x;
859poly f=x5+x3+1;
860
861//---- two divisors ----
862list l1=list(a30,0,1),list(0,1,1);
863list l2=list(a30,a16,1),list(1,1,1);
864
865//---- their rational representation ----
866list D1=ratrep(l1,h,f); D1;
867//[1]:
868//   x2+(a4+a)*x
869//[2]:
870//   (a)*x+1
871
872list D2=ratrep(l2,h,f); D2;
873//[1]:
874//   x2+(a4+a+1)*x+(a4+a)
875//[2]:
876//   (a3+a2+a+1)*x+(a3+a2+a)
877
878//---- back to the point-based-representation ----
879semidiv(D1,h,f);
880//[1]:
881//   [1]:
882//      (a4+a)
883//   [2]:
884//      0
885//   [3]:
886//      1
887//[2]:
888//   [1]:
889//      0
890//   [2]:
891//      1
892//   [3]:
893//      1
894
895semidiv(D2,h,f);
896//[1]:
897//   [1]:
898//      (a4+a)
899//   [2]:
900//      (a4+a3+a+1)
901//   [3]:
902//      1
903//[2]:
904//   [1]:
905//      1
906//   [2]:
907//      1
908//   [3]:
909//      1
910
911//---- adding D1 and D2 ----
913//[1]:
914//   x2+x
915//[2]:
916//   1
917
918//---- D1+D2 in point-based-representation ----
919semidiv(D12,h,f);
920//[1]:
921//   [1]:
922//      1
923//   [2]:
924//      1
925//   [3]:
926//      1
927//[2]:
928//   [1]:
929//      0
930//   [2]:
931//      1
932//   [3]:
933//      1
934
935//---- D1 + D1   (2 possible ways, same result) ----
937double(D1,h,f);
938//[1]:
939//   x2+(a3+1)
940//[2]:
941//   (a4+a3+a+1)*x+(a4+a3+a2+a+1)
942
943//---- order of D1 in the jacobian over the basering ----
944int i=1;
945list E=D1;
946while (E[1] != 1 or E[2] != 0 )
947{
949   i=i+1;
950}
951i;   // 482
952
953//---- proof with multiplikation validates the result ----
954cantormult(i,D1,h,f);
955//[1]:
956//   1
957//[2]:
958//   0
959
960//---- computing the inverse of D1 ----
961list d1= cantormult(-1,D1,h,f);  d1;
962//[1]:
963//   x2+(a4+a)*x
964//[2]:
965//   x2+(a+1)*x+1
966
967//---- proof validates the result ----