Changeset cfdc527 in git


Ignore:
Timestamp:
Sep 15, 2023, 12:17:32 PM (8 months ago)
Author:
slap <slaplagne@…>
Branches:
(u'spielwiese', 'd08f5f0bb3329b8ca19f23b74cb1473686415c3a')
Children:
05c66a9bb25a83fdfa3f3e8391519d12c0250e3b
Parents:
eb2ffedb2b98fc8bd9435f62972bddf37e7f82b3
git-author:
slap <slaplagne@gmail.com>2023-09-15 12:17:32+02:00
git-committer:
Hans Schoenemann <hannes@mathematik.uni-kl.de>2023-11-07 16:28:40+01:00
Message:
Only one expansion is returned

The puiseux command returns one representative of each class of Puiseux expansions. The procedures for computation of the integral basis are modified accordingly.
Location:
Singular/LIB
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/integralbasis.lib

    reb2ffe rcfdc527  
    12201220  dbprint(dbg, "--Building factors of different degrees...");
    12211221  list bF = buildFactors(classes);
     1222  //"Factors: "; bF;
    12221223
    12231224  debug_log_intbas(3, "------Time for the factors: ", timer - tt);
    … …  
    12341235
    12351236  list bestElem = vS[2];
     1237 
     1238  "classes";
     1239  ~;
    12361240
    12371241  // The highest exponent of the denominator
    … …  
    12661270  debug_log_intbas(3, "------Time Hensel Blocks: ", timer - tt);
    12671271  tt = timer;
    1268 
    12691272
    12701273////////////////////////////////////////////////////////////////////////
    … …  
    13991402  if((ifOut[4] == 1))  // Wrong number of factors, recompute orders
    14001403  {
     1404    "ifOut4 = 1";
    14011405    kill ordsFull;
    14021406    kill ordsBest;
    … …  
    15061510    //debug_log_intbas(3, "------Time for recomputation: ", timer - tt);
    15071511  }
    1508 
     1512 
    15091513  ideal tmp_var_mdm=std(var(1)^mdm);
    15101514  for(i = 1; i <= size(I2LiftedFull); i++)
    … …  
    15221526
    15231527  tt = timer;
     1528 
     1529  "Before locLoc";
     1530  ~;
    15241531
    15251532  if(optimize == 0)
    … …  
    19721979//         The third element is a list of the degrees of the corresponding
    19731980//         factors.
    1974 static proc buildFactors(list classes)
    1975 {
     1981
     1982// We adapt buildFactors using only one representative of the class.
     1983
     1984proc buildFactors(list classes)
     1985{
     1986 
     1987  "buildFactors";
     1988  "classes";
     1989  classes;
     1990 
     1991 
     1992  "here -1";
     1993
    19761994  int i, j;
    19771995  int d = -1;
    … …  
    19922010
    19932011  def R = basering;
    1994   list facs;
    1995   list ords;
    1996   list degs;
     2012  list facs = list();
     2013  list ords = list();
     2014  list degs = list();
    19972015  list gfChecks;
    19982016  int den;
    … …  
    20022020
    20032021  list polyGround;
    2004 
     2022 
    20052023  int isGround = 1;
    2006 
     2024   
     2025  "before for";
    20072026  for(i = 1; i <= size(classes); i++)
    20082027  {
     2028    "after for";
    20092029    facs[i] = list();
    20102030    ords[i] = list();
    … …  
    20152035    first = 1;
    20162036    stop1 = 0;
    2017 
    2018     // If there is only one expansion in the class, there is nothing to do
    2019     if((typeof(classes[i][1]) == "ring") or (size(classes[i]) > 1))
    2020     {
    2021       for(j = 1; j <= size(classes[i]); j++)
    2022       {
    2023         if(typeof(classes[i][j]) == "ring")
     2037   
     2038    // If the denominator is 1 there is nothing to do
     2039   
     2040    int algExt = 0;
     2041    if(typeof(classes[i][1]) == "ring")
     2042    {
     2043      algExt = 1;
     2044    }
     2045    if(typeof(classes[i][1]) <> "ring")
     2046    {
     2047      if(classes[i][1][2] > 1) {algExt = 1;}
     2048    }
     2049   
     2050    if(algExt == 1)
     2051    {
     2052      if(typeof(classes[i][1]) == "ring")
     2053      {
     2054        S = classes[i][1];
     2055        setring S;
     2056        if(typeof(PE[1][7]) != "none")
    20242057        {
    2025           S = classes[i][j];
    2026           setring S;
    2027           if(typeof(PE[1][7])!= "none")
    2028           {
    2029 
    2030             debug_log_intbas(6, "Computation of the characterisitic exponents from the integral basis. Check if correct or characteristic orders are needed.");
    2031 
    2032             //expsT = poly2intvec(PE[1][1]);
    2033             //"Check expsT 1: ", expsT;
    2034             expsT = list2intvec(PE[1][7]);
    2035             //"Check expsT 2: ", expsT;
    2036             //"Compare exps: ";
    2037             //poly2intvec(PE[1][1]);
    2038             //charExp(poly2intvec(PE[1][1]), PE[1][2]);
    2039             //expsT;
    2040           } else
    2041           {
    2042             expsT = deg(PE[1][1]);
    2043           }
    2044           if(j == 1)
    2045           {
    2046             den = PE[1][2];
    2047           }
    2048           setring R;
     2058          debug_log_intbas(6, "Computation of the characterisitic exponents from the integral basis. Check if correct or characteristic orders are needed.");
     2059          expsT = list2intvec(PE[1][7]);
    20492060        } else
    20502061        {
    2051           if(typeof(classes[i][j][7]) != "none")
    2052           {
    2053             expsT = list2intvec(classes[i][j][7]);
    2054           } else
    2055           {
    2056             expsT = deg(classes[i][j][1]);
    2057           }
    2058           if(j == 1)
    2059           {
    2060             den = classes[i][1][2];
    2061           }
     2062          expsT = deg(PE[1][1]);
    20622063        }
    2063         if(j > 1)
     2064        den = PE[1][2];
     2065        setring R;
     2066      } else
     2067      {
     2068        if(typeof(classes[i][1][7]) != "none")
    20642069        {
    2065           exps = appendIntvecs(exps, expsT);
     2070          expsT = list2intvec(classes[i][1][7]);
    20662071        } else
    20672072        {
    2068           exps = expsT;
     2073          expsT = deg(classes[i][1][1]);
    20692074        }
    2070       }
     2075        den = classes[i][1][2];
     2076      }
     2077      exps = expsT;
     2078     
    20712079
    20722080      for(k = 1; k <= size(exps); k++)
    … …  
    20782086        }
    20792087        t = timer;
    2080         for(j = 1; j <= size(classes[i]); j++)
     2088        if(typeof(classes[i][1]) == "ring")
    20812089        {
    2082           if(typeof(classes[i][j]) == "ring")
     2090          S = classes[i][1];
     2091          setring S;
     2092          // d = -1 for computing polynomial of full degree
     2093          if(d >= 0)
    20832094          {
    2084             if(!defined(dMP))
    2085             {
    2086               int dMP;
    2087             }
    2088             dMP = pardeg(minpoly);
    2089             S = classes[i][j];
    2090 
    2091             setring S;
    2092             // d = -1 for computing polynomial of full degree
    2093             if(d >= 0){
    2094               poly fF = jet(PE[1][1], d-1, vx);
    2095             } else
    2096             {
    2097               poly fF = PE[1][1];
    2098             }
    2099             if(dMP > 1)
    2100             {
    2101               poly mp = composePolys(minPolys);
    2102               poly mpX = subst(mp, var(2), var(1));
    2103               poly fFB = buildPolyFracNew(fF, mpX);
    2104               ideal rel = extendBack(fFB, erg[1], dMP);
    2105               matrix CC = coef(fFB, var(1)*var(2));
    2106             } else
    2107             {
    2108               poly fFB = buildPolyFrac(fF);
    2109             }
    2110             setring R;
    2111 
    2112             if(dMP > 1)
    2113             {
    2114               ideal rel = imap(S, rel);
    2115               matrix CC = imap(S, CC);
    2116               int kk;
    2117               fFrac[j] = 0;
    2118               for(kk = 1; kk <= ncols(CC); kk++)
    2119               {
    2120                 fFrac[j] = fFrac[j] + subst(rel[kk], var(1), par(1)) * CC[1, kk];
    2121               }
    2122             } else
    2123             {
    2124               fFrac[j] = imap(S, fFB);
    2125             }
    2126 
    2127             // Cleaning
    2128             setring S;
    2129             kill fF;
    2130             kill fFB;
    2131             setring R;
     2095            poly fF = jet(PE[1][1], d-1, vx);
    21322096          } else
    21332097          {
    2134             if(d>=0){
    2135               fFrac[j] = var(2) - jet(classes[i][j][1], d-1, vx);
    2136             } else
    2137             {
    2138               fFrac[j] = var(2) - classes[i][j][1];
    2139             }
    2140 
     2098            poly fF = PE[1][1];
    21412099          }
     2100          list pGL = buildPolyGroundXRoot(fF, den);
     2101          poly pG = pGL[1];
     2102          int gfCheck = pGL[2];
     2103          kill fF;
     2104          setring R;
     2105          poly pG = fetch(S, pG);
     2106          polyGround = list(pG, gfCheck);
     2107        } else
     2108        {
     2109          if(d>=0)
     2110          {
     2111            poly fF = jet(classes[i][1][1], d-1, vx);
     2112          } else
     2113          {
     2114            poly fF = classes[i][1][1];
     2115          }
     2116          polyGround = buildPolyGroundXRoot(fF, den);
    21422117        }
    21432118
    2144         polyGround = buildPolyGround(fFrac, den);
    2145 
    21462119        fGround = squarefree(polyGround[1]);
     2120
    21472121        if (polyGround[2] == 0)
    21482122        {
    … …  
    21502124          ind++;
    21512125          isGround = 0;
    2152         } else {
    2153 
     2126        } else
     2127        {
    21542128          debug_log_intbas(4, "------Time for this factor: ", timer - t);
    21552129          if(dbg - 1 > 0)
    … …  
    21662140            if(deg(fGround, vy) > degs[i][ind])
    21672141            {
    2168             ind++;
     2142              ind++;
    21692143            } else
    21702144            {
    2171             dbprint(dbg - 1, "Factor discarded");
     2145              dbprint(dbg - 1, "Factor discarded");
    21722146            }
    21732147          } else
    … …  
    21772151        }
    21782152        facs[i][ind] = fGround;
    2179         if(d>=0){
     2153        if(d>=0)
     2154        {
    21802155          ords[i][ind] = number(d) / number(den);
    2181         } else {
     2156        } else
     2157        {
    21822158          ords[i][ind] = number(deg(facs[i][ind],vx)) / number(den);
    21832159        }
    21842160        degs[i][ind] = deg(fGround, vy);
    21852161        gfChecks[i][ind] = polyGround[2];
    2186 
    21872162      }
    21882163    }
    … …  
    22392214      fTemp = subst(f, var(3), T(i));
    22402215      fNewTemp = fNewTemp * (var(2) - fTemp);
    2241       fNew = reduce(fNewTemp, I);
    2242     }
     2216    }
     2217    fNew = reduce(fNewTemp, I);
    22432218    setring R;
    22442219    poly fNew = imap(Q, fNew);
    … …  
    22542229
    22552230// Product of conjugate factors
    2256 // Input:  a polynomial f over an algebraic extension of the base ring.
    2257 // Output: polynomial = the product of (y-f_i) for all conjugate polynomials
    2258 //         f_i of f.
    2259 static proc buildPolyFrac(poly f)
    2260 {
     2231// Input:  a polynomial f representing a Puiseux expansions, where the denominator in the exponent is den.
     2232// Output: polynomial = the product of (y-f_i) for all conjugate polynomials, taking the conjugates of x^(1/den).
     2233proc buildPolyGroundXRoot(poly f, int den)
     2234{
     2235  "buildPolyGroundXRoot";
     2236  f, den;
    22612237  int i;
     2238  int gfCheck = 1;
    22622239  def R = basering;
    22632240
    … …  
    22652242  int degX = deg(f, dx);
    22662243
    2267   poly mp = minpoly;
    2268 
    2269 
    2270   mp = subst(minpoly, @a, var(1));
    2271   int d = pardeg(minpoly);
    2272   def Q = ring(0, (var(1), var(2), @a, T(1..d)), lp);
    2273   ring S = 0, (var(1), @a), dp;
    2274   poly f = imap(R, f);
    2275   intvec da = (0,1);
    2276 
    2277   if(deg(f, da) > 0)
    2278   {
    2279     poly mp = imap(R, mp);
    2280     mp = subst(mp,var(1),@a);
    2281     setring Q;
    2282     poly f = imap(S, f);
    2283     poly mp = imap(S, mp);
    2284     mp = mp / leadcoef(mp);
    2285     poly fNew = 1;
    2286     poly rels = 1;
    2287     //ideal I = var(1)^(degX + 1);  // We use this to reduce modulo x^maxDeg
    2288     ideal I;  // We use this to reduce modulo x^maxDeg
    2289     for(i = 1; i <= d; i++)
    2290     {
    2291       rels = rels * (var(1) - T(i));
    2292     }
    2293     matrix c1 = coeffs(rels, var(1));
    2294     matrix c2 = coeffs(mp, @a);
    2295     for(i = 1; i<=d; i++)
    2296     {
    2297       I[i] = c1[i,1] - c2[i,1];
    2298     }
    2299     I = groebner(I);
    2300 
    2301     int t = timer;
    2302     poly fTemp;
    2303     for(i = 1; i<=d; i++)
    2304     {
    2305       fTemp = subst(f, @a, T(i));
    2306       fNew = fNew * (var(2) - fTemp);
    2307       fNew = reduce(fNew, I);
    2308     }
    2309     debug_log_intbas(4,"--------Time for product reduction: ", timer - t);
    2310 
    2311     setring R;
    2312     poly fNew = imap(Q, fNew);
    2313 
    2314   } else
    2315   {
    2316     setring R;
    2317     poly fNew = (var(2) - f)^d;
    2318   }
    2319   return(fNew);
     2244  list l = ringlist(basering);
     2245  for(i = 1; i<=den; i++)
     2246  {
     2247    l[2][size(l[2])+1]="T(" + string(i) + ")";
     2248  }
     2249  def S = ring(l);
     2250  setring S;
     2251  list pols;
     2252  poly ff;
     2253  for(i = 1; i <= den; i++)
     2254  {
     2255    ff = fetch(R, f);
     2256    ff = subst(ff, var(1), var(1) * var(2+i));
     2257    pols[i] = ff;
     2258  }
     2259 
     2260  // Newton relations between the variables
     2261  poly rels = 1;
     2262  for(i = 1; i <= den; i++)
     2263  {
     2264    rels = rels * (var(1) - T(i));
     2265  }
     2266 
     2267  poly pExt = var(1)^den - 1; // A root of unity
     2268  matrix c1 = coeffs(rels, var(1));
     2269  matrix c2 = coeffs(pExt, var(1));
     2270 
     2271  ideal I;
     2272  for(i = 1; i<=den; i++)
     2273  {
     2274    I[i] = c1[i,1] - c2[i,1];
     2275  }
     2276  I = groebner(I);
     2277
     2278  poly fNewTemp = 1;
     2279  for(i = 1; i<=den; i++)
     2280  {
     2281    fNewTemp = fNewTemp * (var(2) - pols[i]);
     2282    fNewTemp = reduce(fNewTemp, I);
     2283  }
     2284  poly fNew = reduce(fNewTemp, I);
     2285 
     2286  // We replace x^den by x
     2287  ring P = (0, @a), (@x, @y, @X), dp;
     2288  poly fNew = fetch(S, fNew);
     2289  poly f2 = reduce(fNew, std(@X-@x^den));
     2290  if(deg(f2, (1,0,0)) > 0)
     2291  {
     2292    //"f2", f2;
     2293    ERROR("Polynomial is not over the ground field");
     2294    gfCheck = 0;
     2295  }
     2296  f2 = subst(f2, @X, @x);
     2297  setring R;
     2298  poly fNew = fetch(P, f2);
     2299 
     2300  "Output of polyGroundFrac", fNew;
     2301
     2302  return(list(fNew,gfCheck));
    23202303}
    23212304
    … …  
    23302313static proc buildPolyGround(list fFrac, int sD)
    23312314{
     2315  "buildPolyGround";
     2316  "fFrac, sD";
     2317  fFrac, sD;
     2318 
    23322319  int i;
    23332320  int gfCheck = 1;  // The polynomial computed is over the ground field.
    … …  
    23542341  return(list(fNew, gfCheck));
    23552342}
     2343
     2344
    23562345
    23572346////////////////////////////////////////////////////////////////////////
    … …  
    25572546static proc ibElement(matrix M, list MSelf, int d, intvec md)
    25582547{
     2548  "ibElement";
    25592549  // Best i
    25602550  int i, j, k;
    … …  
    26102600    }
    26112601  }
     2602 
     2603  "maxExpi";
     2604  ~;
    26122605  intvec elem = sums[maxExpi];
    26132606  return(list(maxExp, elem));
    … …  
    28042797//         Puiseux expansions (given also either as expansions over the base
    28052798//         rings or defined in extension rings).
    2806 static proc getClasses(list l)
     2799proc getClasses(list l)
    28072800{
    28082801  int i;
    … …  
    28942887  {
    28952888    degs[i] = 0;
    2896     for(j = 1; j <= size(classes[i]); j++)
    2897     {
    2898       if(typeof(classes[i][j]) == "ring")
    2899       {
    2900         def S = classes[i][j];
    2901         setring S;
    2902         degs[i] = degs[i] + (pardeg(minpoly) div degBase)*size(PE);
    2903         setring R;
    2904         kill S;
    2905       } else
    2906       {
    2907         degs[i] = degs[i] + 1;
    2908       }
     2889    if(typeof(classes[i][1]) == "ring")
     2890    {
     2891      def S = classes[i][1];
     2892      setring S;
     2893      // The number of conjugated expansions is the common demoninator of the exponents times the degree of the extension,
     2894      // except when the algebraic extension comes from taking root of unity.
     2895      // This degree should be checked if it is always correct computation, or if we should compute this number
     2896      // while computing the Puiseux expansions.
     2897      degs[i] = lcm((pardeg(minpoly) div degBase),PE[1][2]);
     2898      setring R;
     2899      kill S;
     2900    } else
     2901    {
     2902      degs[i] = classes[i][1][2];
    29092903    }
    29102904  }
    29112905  return(degs);
    29122906}
     2907
    29132908
    29142909////////////////////////////////////////////////////////////////////////
    … …  
    30503045 
    30513046  list I2Lifted = henselBlocks(f, degExpand, 1);
    3052 
    3053 
     3047 
    30543048  list updatedClasses = classes;
    30553049  int nClasses;
    … …  
    31113105
    31123106        newL = puiseux(I2Red, degExpand, 1);
    3113         dbprint(dbg, "Puiseux expansions finished");
    3114 
    31153107       
    31163108        classes2 = getClasses(newL);
    … …  
    31283120
    31293121
    3130           for(cl = 1; cl <= size(classes2[j]); cl++)
     3122          if(typeof(classes2[j][1]) == "ring")
    31313123          {
    3132 
    3133             if(typeof(classes2[j][cl]) == "ring")
    3134             {
    3135               def S = classes2[j][cl];
    3136               setring S;
    3137 
    3138               poly fu = buildPolyFrac(PE[1][1]);
    3139 
    3140               PEPE = PE[1][2];
    3141 
    3142               setring R;
    3143               fu = imap(S, fu);
    3144               kill S;
    3145               fuProd = fuProd * fu;
    3146             } else
    3147             {
    3148               fuProd = fuProd * (var(2) - classes2[j][cl][1]);
    3149             }
     3124            def S = classes2[j][1];
     3125            setring S;
     3126           
     3127            list bP = buildPolyGroundXRoot(PE[1][1], PE[1][2]);
     3128           
     3129            "bP ring", bP;
     3130
     3131
     3132            poly fuProd = bP[1];
     3133            int gfCheck = bP[2];
     3134            //
     3135            //PEPE = PE[1][2];
     3136
     3137            setring R;
     3138            poly fuProd = imap(S, fuProd);
     3139            kill S;
     3140          } else
     3141          {
     3142            list bP = buildPolyGroundXRoot(classes2[j][1][1], classes2[j][1][2]);
     3143            poly fuProd = bP[1];
     3144            int gfCheck = bP[2];
     3145           
     3146            "bP poly", bP;
     3147           
    31503148          }
    31513149
    3152           bPG = buildPolyGround(fuProd, PEPE);
    3153           dbprint(dbg, "buildPolyGround finished");
    3154 
     3150          bPG = list(fuProd, gfCheck);
     3151          //"bPG"; bPG;
    31553152
    31563153          if(bPG[2] == 0)
    … …  
    31663163
    31673164
    3168           I2LiftedNew[ind] = fu;
     3165          I2LiftedNew[ind] = bPG[1];
    31693166          gfCheckList[ind] = bPG[2];
    31703167
    … …  
    31803177  }
    31813178  list ll = list(I2Lifted, gfCheckList, gfCheck, wrongNumber, updatedClasses);
    3182 
    3183   dbprint(dbg, "END - irreducibleFactors");
    3184 
    31853179  return(ll);
    31863180}

  • Singular/LIB/puiseuxexpansions.lib

    reb2ffe rcfdc527  
    272272  if(atOrigin == 1)
    273273  {
    274     list p = puiseuxMain(f, maxDeg, 0, 1, 1);
     274    list p = puiseuxMainOneExpansion(f, maxDeg, 0, 1, 1);
     275    //list p = puiseuxMain(f, maxDeg, 0, 1, 1);
    275276  } else
    276277  {
    277     list p = puiseuxMain(f, maxDeg, 0, 0, 1);
     278    list p = puiseuxMainOneExpansion(f, maxDeg, 0, 0, 1);
     279    //list p = puiseuxMain(f, maxDeg, 0, 0, 1);
    278280  }
    279281
    … …  
    477479        @c = puiseuxStep(fTemp, slN, 1);
    478480        fc = factorize(@c);
    479 
     481       
    480482        for(j = 2; j <= size(fc[1]); j++)
    481483        {
    … …  
    493495          }
    494496
    495           if(fc[1][j] != var(2)){
     497          if(fc[1][j] != var(2))
     498          {
    496499            minP = fc[1][j];
    497500            if(deg(minP)==1)
    … …  
    499502              cofs = coeffs(minP, var(2));
    500503              c1 = number(-cofs[1,1]/cofs[2,1]);
    501               if(defined(erg)){
     504              if(defined(erg))
     505              {
    502506                sizeErg = size(erg);
    503               } else {
     507              } else
     508              {
    504509                sizeErg = 1;
    505510              }
    … …  
    533538                    {
    534539                      cs[size(cs) + 1] = list((csT[k][1] + c1)*(var(1)^(csT[k][2]*slN)), csT[k][2]*slD);
    535                     } else {
     540                    } else
     541                    {
    536542                      cs[size(cs) + 1] = list(csT[k][1] + c1, csT[k][2]*slD);
    537543                      cs[size(cs)][8] = list("Denominator", var(1) ^ (-slN*csT[k][2]));
    … …  
    541547                    // NEW CODING FOR CONJUGATED EXPANSIONS
    542548                    cs[size(cs)][6] = insert(csT[k][6], cod);
    543 
     549                   
    544550                    if(newExt == 1)
    545551                    {
    546                       if(typeof(csT[k][7])=="none"){
     552                      if(typeof(csT[k][7])=="none")
     553                      {
    547554                        cs[size(cs)][7] = list(0);
    548555                      } else
    … …  
    553560                    } else
    554561                    {
    555                       if(typeof(csT[k][7]) != "none"){
     562                      if(typeof(csT[k][7]) != "none")
     563                      {
    556564                        cs[size(cs)][7] = sum2All(csT[k][7], slN * csT[k][2]);
    557565                      }
    558566                    }
    559                   } else {
     567                  } else
     568                  {
    560569                    def RT = csT[k];
    561570                    setring RT;
    … …  
    784793  return(cs);
    785794}
     795
     796
     797
     798///////////////////////////////////////////////////////////////////////
     799
     800// Computes the Puiseux expansions starting with slope >= sN / sD;
     801
     802// Output:
     803// cs[1] = Puiseux expansion
     804// cs[2] = denominator of all the exponents
     805// cs[6] = code for identifying different classes
     806// cs[7] = exponents where new branches appear (this is used for computing
     807//         polynomials over the ground field).
     808
     809// It computes only one expansion in each class
     810static proc puiseuxMainOneExpansion(poly f, int maxDeg, int sN, int sD, int firstTime)
     811{
     812  //if(firstTime == 1)
     813  //{
     814  //  "Computing Puiseux expansions for f = ", f;
     815  //}
     816  //else
     817  //{
     818  //  "PuiseuxMain - not first time - f", f;
     819  //}
     820  //
     821  //~;
     822
     823  int dbg = 0;
     824
     825  def R = basering;
     826  intvec vy = (0,1);
     827  int d = deg(f, vy);
     828  int h,i,j, ii;
     829  int k = 1;
     830  int slN, slD;
     831  int g;
     832  int sizeErg;
     833  int stop;
     834  int newMaxDeg;
     835
     836  poly @c;
     837  int newExt;
     838  list fc;
     839  poly minP;
     840  list cs;    // coefficients of the p. exp
     841  list csT;
     842  list out;
     843  poly je;    // Minimal equation
     844  list fJe;   // Factorization of the initial equation
     845  matrix cofs;
     846  number c1;
     847  poly fNew1, fNew2;
     848  poly fTemp;
     849  int cod = 0;
     850  int mi;
     851
     852  // Case of Puiseux expansions with finite number of terms
     853  list dd = Integralbasis::divideBy(f, var(2));
     854  if(dd[2]>0)
     855  {
     856    for(i = 1; i <= dd[2]; i++)
     857    {
     858      cs[size(cs) + 1] = list(0, 1);
     859      cs[size(cs)][6] = list(cod);
     860      cod++;
     861    }
     862    f = dd[1];
     863  }
     864
     865
     866  list l = newtonpoly2(f);
     867  "newton poly", l;
     868
     869  //list l = newtonpoly(f);
     870
     871  //if((l[1][1] == 0) && (l[1][2] == 0))
     872  //{
     873  //  "ERROR: The polynomial must pass through the origin.";
     874  //  list cs = list();
     875  //} else
     876  //{
     877
     878
     879  //for(i = 1; i<size(l); i++)
     880  for(i = size(l)-1; i>=1; i--)
     881  {
     882    slN = l[i+1][1] - l[i][1];
     883    slD = l[i][2] - l[i+1][2];
     884
     885    // We always use positive denominator
     886    if(slD < 0){
     887      slN = (-1)*slN;
     888      slD = (-1)*slD;
     889    }
     890
     891    //if(firstTime == 1)
     892    //{
     893    //  "Computing for segment ", i, " with slope ", slN, " / ", slD;
     894    //} else
     895    //{
     896    //  "Further developing at segment ", i, " with slope ", slN, " / ", slD;
     897    //}
     898
     899    if(slN != 0)
     900    {
     901      g = gcd(slD, slN);
     902    } else {
     903      if(slD != 0)
     904      {
     905        g = slD;
     906      } else
     907      {
     908        g = 1;
     909      }
     910    }
     911
     912    slN = slN div g;
     913    slD = slD div g;
     914
     915
     916    // sD = sN = 0 indicates all slopes must be used.
     917
     918    //if(bigint(sD) * slN > bigint(sN) * slD)
     919    //if((bigint(sD) * slN > bigint(sN) * slD) or ((sD == 0) and (sN == 0)))
     920    if((slN > 0) or ((sD == 0) and (sN == 0)))
     921    {
     922      je = minEqNewton(f, slN, slD);
     923
     924      fJe = factorize(je);
     925      fJe = Integralbasis::sortFactors(fJe);
     926
     927      fNew1 = subst(f, var(1), var(1)^slD);
     928     
     929      "size(fJe[1])";
     930      size(fJe[1]);
     931     
     932      "fJe", fJe;
     933     
     934      for(ii = 2; ii <= size(fJe[1]); ii++)
     935      {
     936        if((ii == 2) || (nonRatCoeff(fJe[1][ii]) == 0))  // If the factor is over the ground field then it is a new conjugacy class.
     937        {
     938          cod++;
     939        }
     940        fTemp = subst(fJe[1][ii]^(fJe[2][ii]), var(1), var(1)^slD);
     941        // Checks if a new extension is needed.
     942        // If so, we have a new cutting point for building the factors.
     943        if((deg(fJe[1][ii], vy) > 1) or (size(fJe[1]) > 2))
     944        {
     945          newExt = 1;
     946        } else
     947        {
     948          newExt = 0;
     949        }
     950
     951        @c = puiseuxStep(fTemp, slN, 1);
     952        fc = factorize(@c);
     953       
     954        "here";
     955        fJe;
     956
     957          // Stopping criterium
     958          stop = 0;
     959          if(maxDeg <= 0)
     960          {
     961            if(fJe[2][ii] == 1)
     962            {
     963              stop = 1;
     964            }
     965          } else
     966          {
     967            if(slN >= slD * maxDeg)
     968            {
     969              stop = 1;
     970            }
     971          }
     972
     973          if(fJe[1][ii] != var(2))
     974          {
     975            //minP = fc[1][j];
     976            // The minpoly of Y^5-1 is Y-1 and the minpoly of Y^5-2 is Y-2.
     977            poly minPBefore = subst(fJe[1][ii],var(1),1);
     978            list minPFacs = factorize(minPBefore);
     979            minP = minPFacs[1][2];   // Check if this always the smallest degree
     980            "minPFacs", minPFacs;
     981            "minP = ", minP;
     982            if(deg(minP)==1)
     983            {
     984              cofs = coeffs(minP, var(2));
     985              c1 = number(-cofs[1,1]/cofs[2,1]);
     986              if(defined(erg))
     987              {
     988                sizeErg = size(erg);
     989              } else
     990              {
     991                sizeErg = 1;
     992              }
     993              if(stop == 0)
     994              {
     995                if(slN >= 0)
     996                {
     997                  "check this";
     998                  fNew2 = subst(fNew1, var(2), (c1+var(2))*(var(1)^slN));
     999                } else
     1000                {
     1001                  // Negative exponent case
     1002                  fNew2 = negExp(fNew1, c1, slN);
     1003                }
     1004                fNew2= sat(ideal(fNew2), var(1))[1];
     1005
     1006                if(maxDeg > 0)
     1007                {
     1008                  newMaxDeg = maxDeg * slD - slN;
     1009                } else {
     1010                  newMaxDeg = maxDeg;
     1011                }
     1012
     1013                csT = puiseuxMainOneExpansion(fNew2, newMaxDeg, slN, 1, 0);
     1014
     1015                for(k = 1; k<=size(csT); k++)
     1016                {
     1017                  if(typeof(csT[k]) != "ring")
     1018                  {
     1019                    // Case of a polynomial with the corresponding denominator;
     1020                    if(slN >= 0)
     1021                    {
     1022                      cs[size(cs) + 1] = list((csT[k][1] + c1)*(var(1)^(csT[k][2]*slN)), csT[k][2]*slD);
     1023                    } else
     1024                    {
     1025                      cs[size(cs) + 1] = list(csT[k][1] + c1, csT[k][2]*slD);
     1026                      cs[size(cs)][8] = list("Denominator", var(1) ^ (-slN*csT[k][2]));
     1027                    }
     1028
     1029
     1030                    // NEW CODING FOR CONJUGATED EXPANSIONS
     1031                    cs[size(cs)][6] = insert(csT[k][6], cod);
     1032                   
     1033                    if(newExt == 1)
     1034                    {
     1035                      if(typeof(csT[k][7])=="none"){
     1036                        cs[size(cs)][7] = list(0);
     1037                      } else
     1038                      {
     1039                        cs[size(cs)][7] = insert(csT[k][7], 0);
     1040                      }
     1041                      cs[size(cs)][7] = sum2All(cs[size(cs)][7], slN * csT[k][2]);
     1042                    } else
     1043                    {
     1044                      if(typeof(csT[k][7]) != "none"){
     1045                        cs[size(cs)][7] = sum2All(csT[k][7], slN * csT[k][2]);
     1046                      }
     1047                    }
     1048                  } else
     1049                  {
     1050                    def RT = csT[k];
     1051                    setring RT;
     1052                    number c1 = fetch(R, c1);
     1053                    // Change @a by par(1)?
     1054                    c1 = number(subst(c1, @a, number(erg[sizeErg])));
     1055
     1056                    for(h = 1; h <= size(PE); h++)
     1057                    {
     1058                      if(slN >= 0)
     1059                      {
     1060                        PE[h][1] = (PE[h][1] + c1) * (var(1)^(slN*PE[h][2]));
     1061                        if(newExt == 1)
     1062                        {
     1063                          PE[h][7] = insert(PE[h][7], 0);
     1064                        }
     1065                        PE[h][7] = sum2All(PE[h][7], slN * PE[h][2]);
     1066                        PE[h][2] = PE[h][2] * slD;
     1067
     1068                        // NEW CODING FOR CONJUGATED EXPANSIONS
     1069                        PE[h][6] = insert(PE[h][6], cod);
     1070                      } else {
     1071                        PE[h][1] = PE[h][1] + c1;
     1072                        if(newExt == 1)
     1073                        {
     1074                          PE[h][7] = insert(PE[h][7], 0);
     1075                        }
     1076
     1077                        // Check if this is correct
     1078                        // (before PE[h][8] was defined after PE[h][2] was modified)
     1079                        PE[h][8] = list("Denominator", var(1) ^ (-slN*PE[h][2]));
     1080                        PE[h][2] = PE[h][2] * slD;
     1081
     1082                        // NEW CODING FOR CONJUGATED EXPANSIONS
     1083                        PE[h][6] = insert(PE[h][6], cod);
     1084                      }
     1085                    }
     1086                    setring R;
     1087                    cs[size(cs) + 1] = RT;
     1088                    kill RT;
     1089                  }
     1090                }
     1091              } else {
     1092                cs[size(cs) + 1] = list(c1 * var(1)^slN, slD);
     1093                cs[size(cs)][6] = list(cod);
     1094                if(newExt == 1)
     1095                {
     1096                  cs[size(cs)][7] = list(slN);
     1097                }
     1098              }
     1099            } else {
     1100              if(npars(R) == 1)
     1101              {
     1102                if(!defined(erg))
     1103                {
     1104                  list erg;
     1105                  erg[1] = par(1);
     1106                }
     1107                if(!defined(minPolys))
     1108                {
     1109                  list minPolys;
     1110                  minPolys[1] = 1;
     1111                }
     1112
     1113                def S = Integralbasis::splitRingAt(minP, erg);
     1114                setring S;
     1115                sizeErg = size(erg);
     1116                poly newA = erg[size(erg)];  // The root founded. That is, newA is the root of minP in S.
     1117                poly fNew1 = fetch(R, fNew1);
     1118
     1119                if(!defined(minPolys))
     1120                {
     1121                  list minPolys = fetch(R, minPolys);
     1122                }
     1123
     1124                // We replace the old a by its image in the new ring.
     1125                minPolys[size(minPolys)+1] = fetch(R, minP);
     1126                for(mi = 1; mi <= size(minPolys); mi++)
     1127                {
     1128                  minPolys[mi] = subst(minPolys[mi], par(1), erg[size(erg)-1]);
     1129                }
     1130                fNew1 = subst(fNew1, par(1), erg[size(erg)-1]);
     1131
     1132              } else {
     1133                def S = Integralbasis::splitRingAt(minP);
     1134                setring S;
     1135                list minPolys;
     1136                minPolys[1] = 1;
     1137                minPolys[size(minPolys)+1] = fetch(R, minP);
     1138                poly fNew1 = fetch(R, fNew1);
     1139                poly newA = par(1);      // In this case, it is the new variable.
     1140                list erg = par(1);
     1141                export erg;
     1142                sizeErg = 1;
     1143              }
     1144              export(minPolys);
     1145
     1146              if(defined(PES))
     1147              {
     1148                PES = list();
     1149              } else
     1150              {
     1151                list PES;
     1152              }
     1153              if(stop == 0){
     1154                // We use the image of the root of minP in the new ring.
     1155                if(slN >= 0)
     1156                {
     1157                  poly fNew2 = subst(fNew1, var(2), (newA + var(2))*var(1)^slN);
     1158                } else
     1159                {
     1160                  poly fNew2 = negExp(fNew1, newA, slN);
     1161                }
     1162                fNew2= sat(ideal(fNew2), var(1))[1];
     1163                if(maxDeg > 0)
     1164                {
     1165                  newMaxDeg = maxDeg * slD - slN;
     1166                } else {
     1167                  newMaxDeg = maxDeg;
     1168                }
     1169                list csTS = puiseuxMainOneExpansion(fNew2, newMaxDeg, slN, 1, 0);
     1170                for(k = 1; k<=size(csTS); k++)
     1171                {
     1172                  if(typeof(csTS[k]) == "ring")
     1173                  {
     1174                    def RT = csTS[k];
     1175                    setring RT;
     1176                    for(h = 1; h<=size(PE); h++)
     1177                    {
     1178                      // The new coefficient is the image of a in the new ring.
     1179                      if(slN >= 0)
     1180                      {
     1181                        PE[h][1] = (PE[h][1] + erg[sizeErg]) * var(1)^(PE[h][2]*slN);
     1182                        if(newExt == 1)
     1183                        {
     1184                          PE[h][7] = insert(PE[h][7], 0);
     1185                        }
     1186                        PE[h][7] = sum2All(PE[h][7], slN * PE[h][2]);
     1187                        PE[h][2] = PE[h][2] * slD;
     1188
     1189                        // NEW CODING FOR CONJUGATED EXPANSIONS
     1190                        PE[h][6] = insert(PE[h][6], cod);
     1191                      } else
     1192                      {
     1193                        PE[h][1] = (PE[h][1] + erg[sizeErg]);
     1194                        if(newExt == 1)
     1195                        {
     1196                          PE[h][7] = insert(PE[h][7], 0);
     1197                        }
     1198                        PE[h][7] = sum2All(PE[h][7], slN * PE[h][2]);
     1199                        PE[h][2] = PE[h][2] * slD;
     1200
     1201                        // NEW CODING FOR CONJUGATED EXPANSIONS
     1202                        PE[h][6] = insert(PE[h][6], cod);
     1203
     1204                        PE[h][8] = list("Denominator", var(1) ^ (-slN*PE[h][2]));
     1205                      }
     1206                    }
     1207                    setring R;
     1208                    cs[size(cs) + 1] = RT;
     1209                    kill RT;
     1210                    setring S;
     1211                  } else {
     1212                    if(slN >= 0)
     1213                    {
     1214                      PES[size(PES) + 1] = list((csTS[k][1] + newA)*var(1)^(csTS[k][2]*slN), csTS[k][2]*slD);
     1215                    } else {
     1216                      PES[size(PES) + 1] = list(csTS[k][1] + newA, csTS[k][2]*slD);
     1217                      PES[size(PES)][8] = list("Denominator", var(1) ^ (-slN*csTS[k][2]));
     1218                    }
     1219
     1220                    // NEW CODING FOR CONJUGATED EXPANSIONS
     1221                    PES[size(PES)][6] = insert(csTS[k][6], cod);
     1222
     1223                    if(newExt == 1)
     1224                    {
     1225                      if(typeof(csTS[k][7]) == "none")
     1226                      {
     1227                        PES[size(PES)][7] = list(0);
     1228
     1229                      } else
     1230                      {
     1231                        PES[size(PES)][7] = insert(csTS[k][7], 0);
     1232                      }
     1233                      PES[size(PES)][7] = sum2All(PES[size(PES)][7], slN * csTS[k][2]);
     1234                    } else
     1235                    {
     1236                      PES[size(PES)][7] = csTS[k][7];
     1237                    }
     1238
     1239                  }
     1240                }
     1241                if(size(PES) > 0){
     1242                  list PE = PES;
     1243                  export PE;
     1244                  setring R;
     1245                  cs[size(cs) + 1] = S;
     1246                }
     1247              } else {
     1248                poly PE1 = newA*var(1)^slN;
     1249                list PE = list(list(PE1, slD));
     1250                PE[1][6] = list(cod);
     1251                if(newExt == 1)
     1252                {
     1253                  PE[1][7] = list(slN);
     1254                }
     1255                export PE;
     1256
     1257                setring R;
     1258                cs[size(cs) + 1] = S;
     1259              }
     1260              kill S;
     1261            }
     1262          }
     1263          setring R;
     1264       
     1265      }
     1266    }
     1267  }
     1268
     1269  if(size(cs) == 0){
     1270    cs = list(list(poly(0),1));
     1271    cs[1][6] = list(1);
     1272  }
     1273
     1274  return(cs);
     1275}
     1276
     1277
     1278
     1279
     1280
     1281
     1282
     1283
     1284
     1285
     1286
     1287
     1288
     1289
     1290
     1291
     1292
     1293
     1294
     1295
     1296
     1297
    7861298
    7871299
Note: See TracChangeset for help on using the changeset viewer.