source: git/Singular/LIB/findifs.lib

spielwiese
Last change on this file was 4689d0, checked in by Hans Schoenemann <hannes@…>, 3 years ago
rename poly.lib to polylib.lib
  • Property mode set to 100644
File size: 18.1 KB
Line 
1//////////////////////////////////////////////////////////////////////////////
2version="version findifs.lib 4.1.2.0 Feb_2019 "; // $Id$
3category="System and Control Theory";
4info="
5LIBRARY: findifs.lib  Tools for the finite difference schemes
6AUTHORS: Viktor Levandovskyy,     levandov@math.rwth-aachen.de
7
8OVERVIEW:
9        We provide the presentation of difference operators in a polynomial,
10        semi-factorized and a nodal form. Running @code{findifs_example();}
11        will demonstrate, how we generate finite difference schemes of linear PDEs
12        from given approximations.
13
14Theory: The method we use have been developed by V. Levandovskyy and Bernd Martin. The
15computation of a finite difference scheme of a given single linear partial
16differential equation with constant coefficients with a given approximation
17rules boils down to the computation of a Groebner basis of a submodule of
18a free module with respect to the ordering, eliminating module components.
19
20
21Support: SpezialForschungsBereich F1301 of the Austrian FWF
22
23PROCEDURES:
24findifs_example();  contains a guided explanation of our approach
25decoef(P,n);  decompose polynomial P into summands with respect to the number n
26difpoly2tex(S,P[,Q]); present the difference scheme in the nodal form
27
28
29exp2pt(P[,L]); convert a polynomial M into the TeX format, in nodal form
30texcoef(n);  converts the number n into TeX
31npar(n); search for 'n' among the parameters and returns its number
32magnitude(P); compute the square of the magnitude of a complex expression
33replace(s,what,with);  replace in s all the substrings with a given string
34xchange(w,a,b); exchange two substrings in a given string
35
36SEE ALSO: latex_lib, finitediff_lib
37";
38
39
40LIB "latex.lib";
41LIB "polylib.lib";
42
43proc tst_findif()
44{
45  example decoef;
46  example difpoly2tex;
47  example exp2pt;
48  example texcoef;
49  example npar;
50  example magnitude;
51  example replace;
52  example xchange;
53}
54
55// static procs:
56// par2tex(s);     convert special characters to TeX in s
57// mon2pt(P[,L]); convert a monomial M into the TeX format, in nodal form
58
59
60// 1. GLOBAL ASSUME: in the ring we have first Tx, then Tt: [FIXED, not needed anymore]!
61// 2. map vars other than Tx,Tt to parameters instead or just ignore them [?]
62// 3. clear the things with brackets
63// 4. todo: content resp lcmZ, gcdZ
64
65proc xchange(string where, string a, string b)
66"USAGE:  xchange(w,a,b);  w,a,b strings
67RETURN:  string
68PURPOSE:  exchanges substring 'a' with a substring 'b' in the string w
69NOTE:
70EXAMPLE: example xchange; shows examples
71"{
72  // replaces a<->b in where
73  // assume they are of the same size [? seems to work]
74  string s = "H";
75  string t;
76  t = replace(where,a,s);
77  t = replace(t,b,a);
78  t = replace(t,s,b);
79  return(t);
80}
81example
82{
83 " EXAMPLE:"; echo=2;
84 ring r = (0,dt,dh,A),Tt,dp;
85 poly p = (Tt*dt+dh+1)^2+2*A;
86 string s = texpoly("",p);
87 s;
88 string t = xchange(s,"dh","dt");
89 t;
90}
91
92static proc par2tex(string s)
93"USAGE:  par2tex(s);  s a string
94RETURN:  string
95PURPOSE: converts special characters to TeX in s
96NOTE: the convention is the following:
97      'Tx' goes to 'T_x',
98      'dx' to '\\tri x' (the same for dt, dy, dz),
99      'theta', 'ro', 'A', 'V' are converted to greek letters.
100EXAMPLE: example par2tex; shows examples
101"{
102  // can be done with the help of latex_lib
103
104// s is a tex string with a poly
105// replace theta with \theta
106// A with \lambda
107// dt with \tri t
108// dh with \tri h
109// Tx with T_x, Ty with T_y
110// Tt with T_t
111// V with \nu
112// ro with \rho
113// dx with \tri x
114// dy with \tri y
115  string t = s;
116  t = replace(t,"Tt","T_t");
117  t = replace(t,"Tx","T_x");
118  t = replace(t,"Ty","T_y");
119  t = replace(t,"dt","\\tri t");
120  t = replace(t,"dh","\\tri h");
121  t = replace(t,"dx","\\tri x");
122  t = replace(t,"dy","\\tri y");
123  t = replace(t,"theta","\\theta");
124  t = replace(t,"A","\\lambda");
125  t = replace(t,"V","\\nu");
126  t = replace(t,"ro","\\rho");
127  return(t);
128}
129example
130{
131 " EXAMPLE:"; echo=2;
132 ring r = (0,dt,theta,A),Tt,dp;
133 poly p = (Tt*dt+theta+1)^2+2*A;
134 string s = texfactorize("",p);
135 s;
136 par2tex(s);
137 string T = texfactorize("",p*(-theta*A));
138 par2tex(T);
139}
140
141proc replace(string s, string what, string with)
142"USAGE:  replace(s,what,with);  s,what,with strings
143RETURN:  string
144PURPOSE: replaces in 's' all the substrings 'what' with substring 'with'
145NOTE:
146EXAMPLE: example replace; shows examples
147"{
148  // clear: replace in s, "what" with "with"
149  int ss = size(s);
150  int cn = find(s,what);
151  if ( (cn==0) || (cn>ss))
152  {
153    return(s);
154  }
155  int gn = 0; // global counter
156  int sw = size(what);
157  int swith = size(with);
158  string out="";
159  string tmp;
160  gn  = 0;
161  while(cn!=0)
162  {
163//    "cn:";    cn;
164//    "gn";    gn;
165    tmp = "";
166    if (cn>gn)
167    {
168      tmp = s[gn..cn-1];
169    }
170//    "tmp:";tmp;
171//    out = out+tmp+" "+with;
172    out = out+tmp+with;
173//    "out:";out;
174    gn  = cn + sw;
175    if (gn>ss)
176    {
177      // ( (gn>ss) || ((sw>1) && (gn >= ss)) )
178      // no need to append smth
179      return(out);
180    }
181//     if (gn == ss)
182//     {
183
184//     }
185    cn  = find(s,what,gn);
186  }
187  // and now, append the rest of s
188  //  out = out + " "+ s[gn..ss];
189  out = out + s[gn..ss];
190  return(out);
191}
192example
193{
194 " EXAMPLE:"; echo=2;
195 ring r = (0,dt,theta),Tt,dp;
196 poly p = (Tt*dt+theta+1)^2+2;
197 string s = texfactorize("",p);
198 s;
199 s = replace(s,"Tt","T_t"); s;
200 s = replace(s,"dt","\\tri t"); s;
201 s = replace(s,"theta","\\theta"); s;
202}
203
204proc exp2pt(poly P, list #)
205"USAGE:  exp2pt(P[,L]);  P poly, L an optional list of strings
206RETURN:  string
207PURPOSE: convert a polynomial M into the TeX format, in nodal form
208ASSUME: coefficients must not be fractional
209NOTE: an optional list L contains a string, which will replace the default
210value 'u' for the discretized function
211EXAMPLE: example exp2pt; shows examples
212"{
213// given poly in vars [now Tx,Tt are fixed],
214// create Tex expression for points of lattice
215// coeffs must not be fractional
216  string varnm = "u";
217  if (size(#) > 0)
218  {
219    if (typeof(#[1])=="string")
220    {
221      varnm = string(#[1]);
222    }
223  }
224  //  varnm;
225  string rz,mz;
226  while (P!=0)
227  {
228    mz = mon2pt(P,varnm);
229    if (mz[1]=="-")
230    {
231      rz = rz+mz;
232    }
233    else
234    {
235      rz = rz + "+" + mz;
236    }
237    P = P-lead(P);
238  }
239  rz  = rz[2..size(rz)];
240  return(rz);
241}
242example
243{
244  " EXAMPLE:"; echo=2;
245  ring r = (0,dh,dt),(Tx,Tt),dp;
246  poly M = (4*dh*Tx^2+1)*(Tt-1)^2;
247  print(exp2pt(M));
248  print(exp2pt(M,"F"));
249}
250
251static proc mon2pt(poly M, string V)
252"USAGE:  mon2pt(M,V);  M poly, V a string
253RETURN:  string
254PURPOSE: convert a monomial M into the TeX format, nodal form
255EXAMPLE: example mon2pt; shows examples
256"{
257  // searches for Tx, then Tt
258  // monomial to the lattice point conversion
259  // c*X^a*Y^b --> c*U^{n+a}_{j+b}
260  number cM = leadcoef(M);
261  intvec e = leadexp(M);
262  //  int a = e[2]; // convention: first Tx, then Tt
263  //  int b = e[1];
264  int i;
265  int a , b, c = 0,0,0;
266  int ia,ib,ic = 0,0,0;
267  int nv = nvars(basering);
268  string s;
269  for (i=1; i<=nv ; i++)
270  {
271    s = string(var(i));
272    if (s=="Tt") { a = e[i]; ia = i;}
273    if (s=="Tx") { b = e[i]; ib = i;}
274    if (s=="Ty") { c = e[i]; ic = i;}
275  }
276  //  if (ia==0) {"Error:Tt not found!"; return("");}
277  //  if (ib==0) {"Error:Tx not found!"; return("");}
278  //  if (ic==0)   {"Error:Ty not found!"; return("");}
279  //  string tc = texobj("",c); // why not texpoly?
280  string tc = texcoef(cM);
281  string rs;
282  if (cM==-1)
283  {
284    rs = "-";
285  }
286  if (cM^2 != 1)
287  {
288// we don't need 1 or -1 as coeffs
289//    rs = clTex(tc)+" ";
290//    rs = par2tex(rmDol(tc))+" ";
291    rs = par2tex(tc)+" ";
292  }
293// a = 0 or b = 0
294  rs = rs + V +"^{n";
295  if (a!=0)
296  {
297    rs = rs +"+"+string(a);
298  }
299  rs = rs +"}_{j";
300  if (b!=0)
301  {
302    rs = rs +"+"+string(b);
303  }
304  if (c!=0)
305  {
306    rs = rs + ",k+";
307    rs = rs + string(c);
308  }
309  rs = rs +"}";
310  return(rs);
311}
312example
313{
314 "EXAMPLE:"; echo=2;
315  ring r = (0,dh,dt),(Tx,Tt),dp;
316  poly M = (4*dh^2-dt)*Tx^3*Tt;
317  print(mon2pt(M,"u"));
318  poly N = ((dh-dt)/(dh+dt))*Tx^2*Tt^2;
319  print(mon2pt(N,"f"));
320  ring r2 = (0,dh,dt),(Tx,Ty,Tt),dp;
321  poly M  = (4*dh^2-dt)*Tx^3*Ty^2*Tt;
322  print(mon2pt(M,"u"));
323}
324
325proc npar(number n)
326"USAGE:  npar(n);  n a number
327RETURN:  int
328PURPOSE: searches for 'n' among the parameters and returns its number
329EXAMPLE: example npar; shows examples
330"{
331  // searches for n amongst parameters
332  // and returns its number
333  int i,j=0,0;
334  list L = ringlist(basering);
335  list M = L[1][2]; // pars
336  string sn = string(n);
337  sn = sn[2..size(sn)-1];
338  for (i=1; i<=size(M);i++)
339  {
340    if (M[i] == sn)
341    {
342      j = i;
343    }
344  }
345  if (j==0)
346  {
347    "Incorrect parameter";
348  }
349  return(j);
350}
351example
352{
353 "EXAMPLE:"; echo=2;
354  ring r = (0,dh,dt,theta,A),t,dp;
355  npar(dh);
356  number T = theta;
357  npar(T);
358  npar(dh^2);
359}
360
361proc decoef(poly P, number n)
362"USAGE:  decoef(P,n);  P a poly, n a number
363RETURN:  ideal
364PURPOSE:  decompose polynomial P into summands with respect
365to the  presence of the number n in the coefficients
366NOTE:    n is usually a parameter with no power
367EXAMPLE: example decoef; shows examples
368"{
369  // decomposes polynomial into summands
370  // w.r.t. the presence of a number n in coeffs
371  // returns ideal
372  def br = basering;
373  int i,j=0,0;
374  int pos = npar(n);
375  if ((pos==0) || (P==0))
376  {
377    return(0);
378  }
379  pos = pos + nvars(basering);
380// map all pars except to vars, provided no things are in denominator
381  number con = content(P);
382  con = numerator(con);
383  P = cleardenom(P); //destroys content!
384  P = con*P; // restore the numerator part of the content
385  list M = ringlist(basering);
386  list L = M[1..4];
387  list Pars = L[1][2];
388  list Vars = L[2] + Pars;
389  L[1] = L[1][1]; // characteristic
390  L[2] = Vars;
391  // for non-comm things: don't need nc but graded algebra
392  //  list templ;
393  //  L[5] = templ;
394  //  L[6] = templ;
395  def @R = ring(L);
396  setring @R;
397  poly P = imap(br,P);
398  poly P0 = subst(P,var(pos),0);
399  poly P1 = P - P0;
400  ideal I = P0,P1;
401  setring br;
402  ideal I = imap(@R,I);
403  kill @R;
404  // check: P0+P1==P
405  poly Q = I[1]+I[2];
406  if (P!=Q)
407  {
408    "Warning: problem in decoef";
409  }
410  return(I);
411// subtract the pure part from orig and check if n is remained there
412}
413example
414{
415 " EXAMPLE:"; echo=2;
416  ring r = (0,dh,dt),(Tx,Tt),dp;
417  poly P = (4*dh^2-dt)*Tx^3*Tt + dt*dh*Tt^2 + dh*Tt;
418  decoef(P,dt);
419  decoef(P,dh);
420}
421
422proc texcoef(number n)
423"USAGE:  texcoef(n);  n a number
424RETURN:  string
425PURPOSE:  converts the number n into TeX format
426NOTE: if n is a polynomial, texcoef adds extra brackets and performs some space substitutions
427EXAMPLE: example texcoef; shows examples
428"{
429  // makes tex from n
430  // and uses substitutions
431  // if n is a polynomial, adds brackets
432  number D  = denominator(n);
433  int DenIsOne = 0;
434  if ( D==number(1) )
435  {
436    DenIsOne = 1;
437  }
438  string sd = texpoly("",D);
439  sd = rmDol(sd);
440  sd = par2tex(sd);
441  number N  = numerator(n);
442  string sn = texpoly("",N);
443  sn = rmDol(sn);
444  sn = par2tex(sn);
445  string sout="";
446  int i;
447  int NisPoly = 0;
448  if (DenIsOne)
449  {
450    sout = sn;
451    for(i=1; i<=size(sout); i++)
452    {
453      if ( (sout[i]=="+") || (sout[i]=="-") )
454      {
455        NisPoly = 1;
456      }
457    }
458    if (NisPoly)
459    {
460      sout = "("+sout+")";
461    }
462  }
463  else
464  {
465    sout = "\\frac{"+sn+"}{"+sd+"}";
466  }
467  return(sout);
468}
469example
470{
471 " EXAMPLE:"; echo=2;
472  ring r = (0,dh,dt),(Tx,Tt),dp;
473  number n1,n2,n3 = dt/(4*dh^2-dt),(dt+dh)^2, 1/dh;
474  n1; texcoef(n1);
475  n2; texcoef(n2);
476  n3; texcoef(n3);
477}
478
479static proc rmDol(string s)
480{
481  // removes $s and _no_ (s on appearance
482  int i = size(s);
483  if (s[1] == "$") { s = s[2..i]; i--;}
484  if (s[1] == "(") { s = s[2..i]; i--;}
485  if (s[i] == "$") { s = s[1..i-1]; i--;}
486  if (s[i] == ")") { s = s[1..i-1];}
487  return(s);
488}
489
490proc difpoly2tex(ideal S, list P, list #)
491"USAGE:  difpoly2tex(S,P[,Q]);  S an ideal, P and optional Q are lists
492RETURN:  string
493PURPOSE: present the difference scheme in the nodal form
494ASSUME: ideal S is the result of @code{decoef} procedure
495NOTE: a list P may be empty or may contain parameters, which will not
496appear in denominators
497@* an optional list Q represents the part of the scheme, depending
498on other function, than the major part
499EXAMPLE: example difpoly2tex; shows examples
500"
501{
502  // S = sum s_i = orig diff poly or
503  // the result of decoef
504  // P = list of pars (numbers) not to be divided with, may be empty
505  // # is an optional list of polys, repr. the part dep. on "f", not on "u"
506  //  S = simplify(S,2); // destroys the leadcoef
507  // rescan S and remove 0s from it
508  int i;
509  ideal T;
510  int ss = ncols(S);
511  int j=1;
512  for(i=1; i<=ss; i++)
513  {
514    if (S[i]!=0)
515    {
516      T[j]=S[i];
517      j++;
518    }
519  }
520  S  = T;
521  ss = j-1;
522  int GotF = 1;
523  list F;
524  if (size(#)>0)
525  {
526    F = #;
527    if ( (size(F)==1) && (F[1]==0) )
528    {
529      GotF = 0;
530    }
531  }
532  else
533  {
534    GotF = 0;
535  }
536  int sf = size(F);
537
538  ideal SC;
539  int  sp = size(P);
540  intvec np;
541  int GotP = 1;
542  if (sp==0)
543  {
544    GotP = 0;
545  }
546  if (sp==1)
547  {
548    if (P[1]==0)
549    {
550      GotP = 0;
551    }
552  }
553  if (GotP)
554  {
555    for (i=1; i<=sp; i++)
556    {
557      np[i] = npar(P[i])+ nvars(basering);
558    }
559  }
560  for (i=1; i<=ss; i++)
561  {
562    SC[i] = leadcoef(S[i]);
563  }
564  if (GotF)
565  {
566    for (i=1; i<=sf; i++)
567    {
568      SC[ss+i] = leadcoef(F[i]);
569    }
570  }
571  def br = basering;
572// map all pars except to vars, provided no things are in denominator
573  list M = ringlist(basering);
574  list L = M[1..4]; // erase nc part
575  list Pars = L[1][2];
576  list Vars = L[2] + Pars;
577  L[1] = L[1][1]; // characteristic
578  L[2] = Vars;
579
580  def @R = ring(L);
581  setring @R;
582  ideal SC = imap(br,SC);
583  if (GotP)
584  {
585    for (i=1; i<=sp; i++)
586    {
587      SC = subst(SC,var(np[i]),1);
588    }
589  }
590  poly q=1;
591  q = lcm(q,SC);
592  setring br;
593  poly  q = imap(@R,q);
594  number lq = leadcoef(q);
595  //  lq;
596  number tmp;
597  string sout="";
598  string vname = "u";
599  for (i=1; i<=ss; i++)
600  {
601    tmp  = leadcoef(S[i]);
602    S[i] = S[i]/tmp;
603    tmp = tmp/lq;
604    sout = sout +"+ "+texcoef(tmp)+"\\cdot ("+exp2pt(S[i])+")";
605  }
606  if (GotF)
607  {
608    vname = "p"; //"f";
609    for (i=1; i<=sf; i++)
610    {
611      tmp  = leadcoef(F[i]);
612      F[i] = F[i]/tmp;
613      tmp = tmp/lq;
614      sout = sout +"+ "+texcoef(tmp)+"\\cdot ("+exp2pt(F[i],vname)+")";
615    }
616  }
617  sout = sout[3..size(sout)]; //rm first +
618  return(sout);
619}
620example
621{
622  "EXAMPLE:"; echo=2;
623  ring r = (0,dh,dt,V),(Tx,Tt),dp;
624  poly M = (4*dh*Tx+dt)^2*(Tt-1) + V*Tt*Tx;
625  ideal I = decoef(M,dt);
626  list L; L[1] = V;
627  difpoly2tex(I,L);
628  poly G = V*dh^2*(Tt-Tx)^2;
629  difpoly2tex(I,L,G);
630}
631
632
633proc magnitude(poly P)
634"USAGE:  magnitude(P);  P a poly
635RETURN:  poly
636PURPOSE: compute the square of the magnitude of a complex expression
637ASSUME: i is the variable of a basering
638EXAMPLE: example magnitude; shows examples
639"
640{
641  // check whether i is present among the vars
642  list L = ringlist(basering)[2]; // vars
643  int j; int cnt = 0;
644  for(j=size(L);j>0;j--)
645  {
646    if (L[j] == "i")
647    {
648      cnt = 1; break;
649    }
650  }
651  if (!cnt)
652  {
653    ERROR("a variable called i is expected in basering");
654  }
655  // i is present, check that i^2+1=0;
656//   if (NF(i^2+1,std(0)) != 0)
657//   {
658//     "Warning: i^2+1=0 does not hold. Reduce the output manually";
659//   }
660  poly re = subst(P,i,0);
661  poly im = (P - re)/i;
662  return(re^2+im^2);
663}
664example
665{
666  "EXAMPLE:"; echo=2;
667  ring r = (0,d),(g,i,sin,cos),dp;
668  poly P = d*i*sin - g*cos +d^2*i;
669  NF( magnitude(P), std(i^2+1) );
670}
671
672
673static proc clTex(string s)
674// removes beginning and ending $'s
675{
676  string t;
677  if (size(s)>2)
678  {
679    // why -3?
680     t = s[2..(size(s)-3)];
681  }
682  return(t);
683}
684
685static proc simfrac(poly up, poly down)
686{
687  // simplifies a fraction up/down
688  // into the form up/down = RT[1] + RT[2]/down
689  list LL = division(up,down);
690  list RT;
691  RT[1] =  LL[1][1,1]; // integer part
692  RT[2] = L[2][1]; // new numerator
693  return(RT);
694}
695
696proc findifs_example()
697"USAGE:  findifs_example();
698RETURN:  nothing (demo)
699PURPOSE:  demonstration of our approach and this library
700EXAMPLE: example findifs_example; shows examples
701"
702{
703
704  "* Equation: u_tt - A^2 u_xx -B^2 u_yy = 0; A,B are constants";
705  "* we employ three central differences";
706  "* the vector we act on is (u_xx, u_yy, u_tt, u)^T";
707  "* Set up the ring: ";
708  "ring r = (0,A,B,dt,dx,dy),(Tx,Ty,Tt),(c,dp);";
709  ring r = (0,A,B,dt,dx,dy),(Tx,Ty,Tt),(c,dp);
710  "* Set up the matrix with equation and approximations: ";
711  "matrix M[4][4] =";
712    "    // direct equation:";
713  "    -A^2, -B^2, 1, 0,";
714  "    // central difference u_tt";
715  "    0, 0,  -dt^2*Tt, (Tt-1)^2,";
716  "    // central difference u_xx";
717  "    -dx^2*Tx, 0, 0, (Tx-1)^2,";
718  "    // central difference u_yy";
719  "    0, -dy^2*Ty, 0, (Ty-1)^2;";
720  matrix M[4][4] =
721    // direct equation:
722    -A^2, -B^2, 1, 0,
723    // central difference u_tt
724    0, 0,  -dt^2*Tt, (Tt-1)^2,
725    // central difference u_xx
726    -dx^2*Tx, 0, 0, (Tx-1)^2,
727    // central difference u_yy
728    0, -dy^2*Ty, 0, (Ty-1)^2;
729//=========================================
730// CHECK THE CORRECTNESS OF EQUATIONS AS INPUT:
731ring rX = (0,A,B,dt,dx,dy,Tx,Ty,Tt),(Uxx, Uyy,Utt, U),(c,Dp);
732matrix M = imap(r,M);
733vector X = [Uxx, Uyy, Utt, U];
734 "* Print the differential form of equations: ";
735print(M*X);
736// END CHECK
737//=========================================
738 setring r;
739 "* Perform the elimination of module components:";
740 " module R = transpose(M);";
741 " module S = std(R);";
742 module R = transpose(M);
743 module S = std(R);
744 " * See the result of Groebner bases: generators are columns";
745 " print(S);";
746 print(S);
747 " * So, only the first column has its nonzero element in the last component";
748 " * Hence, this polynomial is the scheme";
749 " poly   p = S[4,1];" ;
750 poly   p = S[4,1]; // by elimination of module components
751 " print(p); ";
752 print(p);
753 list L; L[1]=A;L[2] = B;
754 ideal I = decoef(p,dt); // make splitting w.r.t. the appearance of dt
755 "* Create the nodal of the scheme in TeX format: ";
756 " ideal I = decoef(p,dt);";
757 " difpoly2tex(I,L);";
758 difpoly2tex(I,L); // the nodal form of the scheme in TeX
759 "* Preparations for the semi-factorized form: ";
760 poly pi1 = subst(I[2],B,0);
761 poly pi2 = I[2] - pi1;
762 " poly pi1 = subst(I[2],B,0);";
763 " poly pi2 = I[2] - pi1;";
764 "* Show the semi-factorized form of the scheme: 1st summand";
765 " factorize(I[1]); ";
766 factorize(I[1]); // semi-factorized form of the scheme: 1st summand
767 "* Show the semi-factorized form of the scheme: 2nd summand";
768 " factorize(pi1);";
769 factorize(pi1); // semi-factorized form of the scheme: 2nd summand
770 "* Show the semi-factorized form of the scheme: 3rd summand";
771 " factorize(pi1);";
772 factorize(pi2); // semi-factorized form of the scheme: 3rd summand
773}
774example
775{
776  "EXAMPLE:"; echo=1;
777 findifs_example();
778}
Note: See TracBrowser for help on using the repository browser.