share/contrib/sarag/topology.mac

   1 /* ------------------------------------------------------------------- */
   2 /* SARAG - Topology                                                    */
   3 /* by Fabrizio Caruso   modified by Alexandre Le Meur and Marie-Françoise Roy                 */
   4
   5 SARAG_LINE_COLOR : 1;
   6 SARAG_POINT_SIZE : 3;
   7 SARAG_POINT_STYLE_GENERIC : 3;
   8 SARAG_POINT_STYLE_CRITICAL : 2;
   9 SARAG_POINT_STYLE_NONCRITICAL : 1;
  10
  11 signAtRoot(pol,x,polyInX,rootInX) :=
  12   signsAtRoots([rootInX],polyInX,pol,x)[1][2];
  13
  14 bivariateSignChanges(seq,x,polyInX,rootInX) :=
  15   signChanges(map(lambda([arg],signAtRoot(arg,x,polyInX,rootInX)),seq));
  16
  17
  18
  19 upperBoundWithSignAt(polyInX,interval,x,sgnAtRoot) :=
  20   block([signFound,bound,lhsBound,rhsBound],
  21    signFound : false,
  22    while(not(signFound)) do
  23      (
  24      bound : upperBoundAt(polyInX,interval,x),
  25      lhsBound : bound[1],
  26      rhsBound : bound[2],
  27      if sgn(lhsBound)=sgnAtRoot and sgn(rhsBound)=sgnAtRoot then
  28        signFound : true
  29      else
  30        interval : refineRoot(interval,polyInX,x)
  31      ), /* end while */
  32    return([interval,[lhsBound,rhsBound]])
  33    ); /* end proc */
  34
  35
  36
  37 bivariateSSubResSignChanges(seq,x,polyInX,rootInX) :=
  38   genPermMVar(map(lambda([arg],signAtRoot(arg,x,polyInX,rootInX)),seq));
  39
  40 bivariateSSubResTarskiQuery(q,p,var,x,polyInX,rootInX) :=
  41   block([degQ, degP, sSRRes, sSR, seq, coeffSeq, j, b1, b0, Rbar],
  42
  43   p : expandIf(p),
  44   q : expandIf(q),
  45   degQ : degree(q,var),
  46   degP : degree(p,var),
  47
  48   if degQ = 0 then
  49     (
  50     coeffSeq : sSubResCoeff(p,diff(p,var),var),
  51     return(signAtRoot(leadCoeff(q,var),x,polyInX,rootInX)*
  52            bivariateSSubResSignChanges(coeffSeq,x,polyInX,rootInX))
  53     ) /* end if */
  54   else
  55     if degQ = 1 then
  56       (
  57       b0 : coeff(q,var,0),
  58       b1 : coeff(q,var,1),
  59
  60       Rbar : expand(diff(p,var)*q-degP*b1*p),
  61
  62       coeffSeq : sSubResCoeff(p,Rbar,var), /* SR(p b1 P - P' Q) */
  63
  64       return(bivariateSSubResSignChanges(coeffSeq,x,polyInX,rootInX))
  65       ) /* end if */
  66     else
  67       (
  68
  69       coeffSeq : sSubResCoeff(expand(-diff(p,var)*q,p,var)), /* SR(-P'Q,P) */
  70       if oddp(degQ-1) then
  71          return(bivariateSSubResSignChanges(coeffSeq,x,polyInX,rootInX) + signAtRoot(leadCoeff(q,var),x,polyInX,rootInX))
  72       else
  73          return(bivariateSSubResSignChanges(coeffSeq,x,polyInX,rootInX))
  74       ) /* end else */
  75   ); /* end proc */
  76
  77
  78 multiplicity(sol,pol,var) :=
  79   block([count],
  80   count : 0,
  81   while(expand(factor(subst(sol,var,pol)))=0) do
  82     (
  83     pol : diff(pol,var),
  84     count : count + 1
  85     ), /* end while */
  86   return(count)
  87   ); /* end proc */
  88
  89 bivariateMultiplicity(sol,pol,var,x,polyInX,rootInX) :=
  90   block([count],
  91   count : 0,
  92   while(signAtRoot(expand(factor(subst(sol,var,pol))),x,polyInX,rootInX)=0) do
  93     (
  94     pol : diff(pol,var),
  95     count : count + 1
  96     ), /* end while */
  97   return(count)
  98   ); /* end proc */
  99
 100
 101
 102 archimedeanTopology(pol,isolAlg,x,y) :=
 103   block([sRes,sSR,sSRCoeff,discr,isolDiscr,i,j,jList,epsList,flag,rootSign,degP,
 104          numOfRoots,rootsAtDiscr,numAbove,rootsAtDiscrAboveCr,taQ,
 105          numBetween,rootsBetween,res,
 106          coSys,generic,critY,multi,aboveCr],
 107
 108   pol : expand(pol),
 109
 110   if constantp(pol) then
 111     return([]),
 112
 113
 114   coSys : 0,
 115   while(not(degree(leadCoeff(pol,y),x)=0)) do
 116     (
 117     coSys : coSys + 1,
 118     if WARNINGS then
 119       (
 120       print("topology) WARNING: the function is not quasi-monomial in ", y),
 121       print("topology) WARNING: new curve is : ", expand(subst(x+y,x,pol)))
 122       ), /* end if */
 123     pol : expand(subst(x+y,x,pol))
 124     ), /* end while */
 125   generic : false,
 126
 127
 128
 129   while (not(generic)) do
 130     (
 131
 132     sSR : sSubResPol(pol,diff(pol,y),y),
 133     sSRCoeff : sSubResCoeff(pol,diff(pol,y),y),
 134     discr : last(sSRCoeff),
 135
 136
 137     isolDiscr : isListSort(isolAlg(discr,x)),
 138
 139     jList : [],
 140     epsList : [],
 141     rootsAtDiscr : [],
 142     rootsAtDiscrAboveCr : [],
 143     aboveCr : [],
 144     rootsBetween : [],
 145     degP : degree(expand(pol),y),
 146
 147
 148     generic : true, /* It "hopes" the curve is in generic position but it will check during computation */
 149     for i : 1 thru length(isolDiscr) do /* should insert check on generic position */
 150       (
 151       if generic then
 152         (
 153
 154         if length(isolDiscr[i]) = 1 then
 155           (
 156           j : 0,
 157           flag : false,
 158           while(not(flag)) do
 159             (
 160             rootSign : sgn(subst(isolDiscr[i][1],x,sSRCoeff[degP+1-j])),
 161
 162             if not(rootSign=0) then
 163               flag : true
 164             else
 165               j : j + 1
 166             ), /* end while */
 167
 168           jList : endcons(j,jList),
 169           epsList : endcons(rootSign,epsList),
 170
 171           if not(ASSUME_GENERIC_POSITION) then
 172             (
 173             critY : -1/j * coeff(subst(isolDiscr[i][1],x,sSR[degP+1-j]),y,j-1)/
 174                            subst(isolDiscr[i][1],x,sSRCoeff[degP+1-j]),
 175
 176
 177             multi : multiplicity(critY,subst(isolDiscr[i][1],x,sSR[degP+1-j]),y),
 178             if multi<j then
 179                (
 180                generic:false,
 181
 182                if WARNINGS then
 183                  (
 184                  print("topology) WARNING : Non-generic position --- multiplicity : ",
 185                         multi),
 186                  print("topology) WARNING: new curve is : ", expand(subst(x+y,x,pol)))
 187                  ), /* end if */
 188                pol : expand(subst(x+y,x,pol)),
 189                coSys : coSys + 1,
 190                while(not(degree(leadCoeff(pol,y),x)=0)) do
 191                  (
 192                  coSys : coSys + 1,
 193                  if WARNINGS then
 194                    (
 195                    print("topology) WARNING: the function is not quasi-monomial in ", y),
 196                    print("topology) WARNING: new curve is : ", expand(subst(x+y,x,pol)))
 197                    ), /* end if */
 198                  pol : expand(subst(x+y,x,pol))
 199                  ) /* end while */
 200
 201                ) /* end if */
 202             ), /* end if */
 203
 204           if generic then
 205             (
 206             numOfRoots : genPermMVar(subst(isolDiscr[i][1],x,sSRCoeff)),
 207
 208             rootsAtDiscr : endcons(numOfRoots,rootsAtDiscr),
 209             taQ : sSubResTarskiQuery(
 210                                expand(subst(isolDiscr[i][1],x,
 211                                  rootSign*(j*sSRCoeff[degP+1-j]*y+coeff(sSR[degP+1-j],y,j-1)))),
 212                                expand(subst(isolDiscr[i][1],x,pol)),y),
 213
 214             numAbove : (numOfRoots+taQ-1)/2,
 215
 216             rootsAtDiscrAboveCr : endcons(numAbove,rootsAtDiscrAboveCr)
 217             ) /* end if */
 218
 219           ) /* end if */
 220         else /* case when the root of the dicriminant is described by an open interval */
 221           (
 222
 223           j : 0,
 224           flag : false,
 225           while(not(flag)and j <= degP+1) do
 226             (
 227
 228             rootSign : signAtRoot(sSRCoeff[degP+1-j],
 229                                          x,discr,isolDiscr[i]),
 230
 231             if not(rootSign=0) then
 232               flag:true
 233             else
 234               j : j + 1
 235             ), /* end while */
 236
 237           jList : endcons(j,jList),
 238           epsList : endcons(rootSign,epsList),
 239
 240           if not(ASSUME_GENERIC_POSITION) then
 241             (
 242             critY : -1/j * coeff(sSR[degP+1-j],y,j-1)/sSRCoeff[degP+1-j],
 243
 244
 245
 246             multi : bivariateMultiplicity(critY,sSR[degP+1-j],y,x,discr,isolDiscr[i]),
 247
 248             if multi<j then
 249                (
 250                generic:false,
 251
 252                if WARNINGS then
 253                  (
 254                  print("topology) WARNING : Non-generic position --- multiplicity : ",
 255                         multi),
 256
 257                  print("topology) WARNING: new curve is : ", expand(subst(x+y,x,pol)))
 258                  ), /* end if */
 259                pol : expand(subst(x+y,x,pol)),
 260                coSys : coSys + 1,
 261                while(not(degree(leadCoeff(pol,y),x)=0)) do
 262                  (
 263                  coSys : coSys + 1,
 264                  if WARNINGS then
 265                    (
 266                    print("topology) WARNING: the function is not quasi-monomial in ", y),
 267                    print("topology) WARNING: new curve is : ", expand(subst(x+y,x,pol)))
 268                    ), /* end if */
 269                  pol : expand(subst(x+y,x,pol))
 270                  ) /* end while */
 271                ) /* end if */
 272             ), /* end if */
 273
 274           if generic then
 275             (
 276             numOfRoots : bivariateSSubResSignChanges(sSRCoeff,x,discr,isolDiscr[i]),
 277             rootsAtDiscr : endcons(numOfRoots,rootsAtDiscr),
 278
 279             taQ : bivariateSSubResTarskiQuery(
 280                                       rootSign*(j*sSRCoeff[degP+1-j]*y+coeff(sSR[degP+1-j],y,j-1)),
 281                                       pol,y,x,expand(squareFree(discr,x)),isolDiscr[i]),
 282
 283             numAbove : (numOfRoots+taQ-1)/2,
 284
 285             rootsAtDiscrAboveCr : endcons(numAbove,rootsAtDiscrAboveCr)
 286             ) /* end if */
 287           ) /* end else */
 288         ) /* end if */
 289       ) /* end for */
 290   ), /* end while */
 291
 292 if length(isolDiscr)>0 then
 293   (
 294
 295   rootsBetween : [genPermMVar(subst(isolDiscr[1][1]-1,x,sSRCoeff))],
 296
 297   for i : 1 thru length(isolDiscr)-1 do
 298     (
 299
 300
 301     rootsBetween : endcons(genPermMVar(
 302                               subst(intermidiatePoint(isolDiscr[i],isolDiscr[i+1],discr,x),
 303                                     x,sSRCoeff)),
 304                           rootsBetween)
 305
 306     ), /* end for */
 307
 308   rootsBetween : endcons(genPermMVar(subst(last(last(isolDiscr))+1,x,sSRCoeff)),rootsBetween),
 309   i : 1,
 310   res : [rootsBetween[1] ],
 311   while i<=length(isolDiscr) do
 312     (
 313     res : append(res,[[rootsAtDiscr[i],rootsAtDiscr[i]-rootsAtDiscrAboveCr[i]],rootsBetween[i+1]]),
 314     i : i + 1
 315     ), /* end while */
 316   if DRAW_TOPOLOGY then
 317      drawTopology(res),
 318   return([coSys,res])
 319   ) /* end if */
 320 else
 321   (
 322   res : [sSubResTarskiQuery(1,subst(0,x,pol),y)],
 323   if DRAW_TOPOLOGY then
 324      drawTopology(res),
 325   return([coSys,res])
 326   )
 327 ); /* end proc */
 328
 329
 330
 331 /* ----------------------------------------------------------- */
 332 /* Plotting of the topological graph of a planar curve */
 333 /* Very preliminary */
 334
 335
 336 getPointsOnProj(num,projX,height):=makelist([projX,i*height/(num+1)],i,1,num);
 337
 338
 339 drawPoints(ptList) :=
 340   plot2d([discrete,float(ptList)],
 341 ['style,[points, SARAG_POINT_STYLE_GENERIC]],
 342 [gnuplot_preamble, "set grid; set xrange[-25:25]; set yrange[-10:10]; unset key"]);
 343
 344
 345
 346
 347 drawPointsOnProj(num,projX,height) :=
 348   drawPoints(getPointsOnProj(num,projX,height));
 349
 350
 351
 352 getPointsOnCrit(num,pos,projX,height):=
 353   block([i,crit,nonCritBelow,nonCritAbove],
 354   nonCritBelow : [],
 355   nonCritAbove : [],
 356   for i : 1 thru num do
 357     if i=pos then
 358       crit : [projX,i*height/(num+1)]
 359     else
 360       if i<pos then
 361         nonCritBelow : append(nonCritBelow,[[projX,i*height/(num+1)]])
 362       else
 363         nonCritAbove : append(nonCritAbove,[[projX,i*height/(num+1)]]),
 364   return([crit,nonCritBelow,nonCritAbove])
 365   ); /* end proc */
 366
 367
 368 getTopologyPoints(tpg) :=
 369   block([i,j,nonCritPts,critPts,ptsOnCritProj,previous,jump,newNonCritPts,lastCritPts,lastNonCritPts,
 370          correct,newPtsAtInd,ptsAtInf,lineSegsAtInf, lineSegs, newPtsAtInf],
 371    nonCritPts : [],
 372    lastCritPts : [0,[],[]],
 373    lastNonCritPts : [],
 374    critPts : [],
 375    ptsOnCritProj : [0,[],[]],
 376    lineSegs : [],
 377
 378
 379 /* It checks the form of the input */
 380    i : 1,
 381    correct : true,
 382    while correct and i<= length(tpg) do
 383      (
 384
 385      if evenp(i) then
 386        (
 387        if not(listp(tpg[i])) then
 388           (
 389           correct : false
 390           ) /* end if */
 391        else
 392          if not(length(tpg[i])=2) or
 393             not(numberp(first(tpg[i]))) or
 394             not(numberp(second(tpg[i]))) then
 395            (
 396
 397            correct : false
 398            ) /* end if */
 399        ) /* end if */
 400      else
 401        (
 402        if not(numberp(tpg[i])) then
 403          correct : false
 404
 405        ), /* end else */
 406      i : i + 1
 407      ), /* end of check of the input format */
 408
 409    if not(correct) then
 410      (
 411      return(false)
 412      ) /* end if */
 413    else
 414    (
 415    ptsAtInf : [],
 416    lineSegsAtInf : [],
 417    if PLOT_AT_INFINITY then
 418     (
 419     newPtsAtInf : getPointsOnProj(tpg[1],0,PLOT_HEIGHT),
 420     ptsAtInf : append(ptsAtInf,newPtsAtInf),
 421     newNonCritPts : getPointsOnProj(tpg[1],1*PLOT_SPACING,PLOT_HEIGHT),
 422     for j : 1 thru tpg[1] do
 423       lineSegsAtInf : append(lineSegsAtInf,[[discrete,
 424                                    float([newPtsAtInf[j],
 425                                           newNonCritPts[j]    ])]])
 426
 427
 428     ), /* end if */
 429    for i : 1 thru length(tpg) do
 430     (
 431     if oddp(i) then /* non critical points */
 432       (
 433
 434       newNonCritPts : getPointsOnProj(tpg[i],i*PLOT_SPACING,PLOT_HEIGHT),
 435       nonCritPts : append(nonCritPts,newNonCritPts),
 436
 437       if not(ptsOnCritProj = [0,[],[]]) then
 438         (
 439
 440         jump : length(newNonCritPts)-length(ptsOnCritProj[2])-length(ptsOnCritProj[3]),
 441
 442
 443         for j : 1 thru length(ptsOnCritProj[2]) do
 444           lineSegs : append(lineSegs,[[discrete,
 445                                        float([ptsOnCritProj[2][j],
 446                                               newNonCritPts[j]])]]),
 447
 448         for j : 1 thru length(ptsOnCritProj[3]) do
 449           lineSegs : append(lineSegs,[[discrete,
 450                                        float([ptsOnCritProj[3][j],
 451                                               newNonCritPts[j+jump+length(ptsOnCritProj[2])]])]]),
 452
 453         for j : 1 thru jump do
 454           lineSegs : append(lineSegs,[[discrete,
 455                                        float([ptsOnCritProj[1],newNonCritPts[j+length(ptsOnCritProj[2])]])]])
 456
 457         ) /* end if */
 458
 459
 460       ) /* end if */
 461     else /* critical points */
 462       (
 463
 464
 465       ptsOnCritProj : getPointsOnCrit(tpg[i][1],tpg[i][2],i*PLOT_SPACING,PLOT_HEIGHT),
 466       jump : length(newNonCritPts)-(tpg[i][1]-1),
 467
 468
 469
 470       for j : 1 thru length(ptsOnCritProj[3]) do
 471           (
 472           lineSegs : append(lineSegs,[[discrete,
 473                                        float([newNonCritPts[j+jump+length(ptsOnCritProj[2])],
 474                                               ptsOnCritProj[3][j]])]])
 475
 476           ), /* end for */
 477
 478       for j : 1 thru length(ptsOnCritProj[2]) do
 479           lineSegs : append(lineSegs, [[discrete,
 480                                        float([newNonCritPts[j],ptsOnCritProj[2][j]])]]),
 481
 482       for j : 1 thru jump do
 483           lineSegs : append(lineSegs,[[discrete,
 484                                       float([newNonCritPts[j+length(ptsOnCritProj[2])],
 485                                              ptsOnCritProj[1]])]]),
 486
 487
 488       critPts : append(critPts,[ptsOnCritProj[1]]),
 489       nonCritPts : append(nonCritPts,ptsOnCritProj[2],ptsOnCritProj[3])
 490       ) /* end else */
 491     ), /* end for */
 492     if PLOT_AT_INFINITY then
 493     (
 494     newPtsAtInf : getPointsOnProj(tpg[length(tpg)],(length(tpg)+1)*PLOT_SPACING,PLOT_HEIGHT),
 495     ptsAtInf : append(ptsAtInf,newPtsAtInf),
 496
 497     newNonCritPts : getPointsOnProj(tpg[length(tpg)],length(tpg)*PLOT_SPACING,PLOT_HEIGHT),
 498     for j : 1 thru tpg[length(tpg)] do
 499       lineSegsAtInf : append(lineSegsAtInf,[[discrete,
 500                                              float([newPtsAtInf[j],
 501                                              newNonCritPts[j]    ])]])
 502
 503     )/* end if */
 504
 505   ), /* end while */
 506
 507    return([nonCritPts,critPts,lineSegs,ptsAtInf,lineSegsAtInf])
 508    ); /* end proc */
 509
 510 getCritPoint(critLine,projX,height) :=
 511   [projX,critLine[2]*height/critLine[1]];
 512
 513 drawTopology(tpg) :=
 514   block([nonCritPoints,critPoints,firstAttemptToGetPoints,secondAttemptToGetPoints, topologyPoints,style,preamble,i,
 515          leftMost,rightMost,top,bottom,ptsAtInf,fixedTpg, lineSegs],
 516
 517     firstAttemptToGetPoints: getTopologyPoints(tpg),
 518     if not firstAttemptToGetPoints=false then
 519        (
 520        topologyPoints : firstAttemptToGetPoints,
 521        fixedTpg : tpg
 522        )
 523     else
 524        (
 525        secondAttemptToGetPoints: getTopologyPoints(second(tpg)),
 526        if not secondAttemptToGetPoints=false then
 527           (
 528           topologyPoints : secondAttemptToGetPoints,
 529           fixedTpg : second(tpg)
 530           )
 531        ), /* end else */
 532     if topologyPoints=false then
 533       (
 534       print("Wrong input format!"),
 535       return(false)
 536       ) /* end if */
 537     else
 538      (
 539      nonCritPoints : topologyPoints[1],
 540      critPoints : topologyPoints[2],
 541      lineSegs : topologyPoints[3],
 542
 543      lineSegs : append(lineSegs,topologyPoints[5]),
 544
 545      ptsAtInf : topologyPoints[4],
 546
 547      leftMost : -2,
 548      rightMost : PLOT_SPACING*length(fixedTpg)+3,
 549      bottom : -1,
 550      top : PLOT_HEIGHT+1,
 551      preamble : concat(PLOT_STYLE, "set xrange[-1:",rightMost,
 552                        "]; set yrange[-1:",PLOT_HEIGHT,"];"),
 553
 554      style : [],
 555      for i : 1 thru length(lineSegs) do
 556         style : cons([lines, 1, SARAG_LINE_COLOR],style),
 557
 558      if not(PS_OUTPUT) then
 559        if critPoints=[] then
 560          plot2d(append([[discrete,float(nonCritPoints)]],
 561                        lineSegs),
 562                   cons('style,cons([points, SARAG_POINT_SIZE, SARAG_POINT_STYLE_NONCRITICAL],style)),
 563                   [gnuplot_preamble, preamble])
 564        else
 565          if nonCritPoints = [] then
 566             plot2d(append([[discrete,float(critPoints)]],
 567                           lineSegs),
 568                    cons('style,cons([points, SARAG_POINT_SIZE, SARAG_POINT_STYLE_CRITICAL],style)),
 569                    [gnuplot_preamble,preamble])
 570          else
 571          plot2d(append([[discrete,float(nonCritPoints)],
 572                         [discrete,float(critPoints)]],
 573                        lineSegs),
 574                   cons('style,append([[points, SARAG_POINT_SIZE, SARAG_POINT_STYLE_NONCRITICAL],
 575                                       [points, SARAG_POINT_SIZE, SARAG_POINT_STYLE_CRITICAL]],style)),
 576                   [gnuplot_preamble,preamble])
 577      else
 578        if critPoints=[] then
 579          plot2d(append([[discrete,float(nonCritPoints)]],
 580                        lineSegs),
 581                   cons('style,cons([points, SARAG_POINT_SIZE, SARAG_POINT_STYLE_NONCRITICAL],style)),
 582                   [gnuplot_preamble, preamble],
 583                   [gnuplot_term, ps], [gnuplot_out_file, PS_OUTPUT_FILE_NAME])
 584        else
 585          if nonCritPoints = [] then
 586             plot2d(append([[discrete,float(critPoints)]],
 587                           lineSegs),
 588                    cons('style,cons([points, SARAG_POINT_SIZE, SARAG_POINT_STYLE_CRITICAL],style)),
 589                    [gnuplot_preamble,preamble],
 590                    [gnuplot_term, ps], [gnuplot_out_file, PS_OUTPUT_FILE_NAME])
 591          else
 592          plot2d(append([[discrete,float(nonCritPoints)],
 593                         [discrete,float(critPoints)]],
 594                        lineSegs),
 595                   cons('style,append([[points, SARAG_POINT_SIZE, SARAG_POINT_STYLE_NONCRITICAL],
 596                                       [points, SARAG_POINT_SIZE, SARAG_POINT_STYLE_CRITICAL]],style)),
 597                   [gnuplot_preamble,preamble],
 598                   [gnuplot_term, ps], [gnuplot_out_file, PS_OUTPUT_FILE_NAME])
 599      ), /* end else */
 600
 601 return(true)
 602 ); /* end proc */
 603