# segmented: a package to manipulate piecewise defined functions
#             Copyright 1994 by Robert K. Wright
#             Mathematics Dept, University of Vermont
#             Burlington, Vermont, USA, 05401
#             email: wright@emba.uvm.edu
#
# 7/19/94
#MODIFIED:
#  10/15/94: convert/Htoseg to handle Heaviside(-x+a)
#  10/26/94: evalsg/korder to handle the symbol Pi automatically
#  11/10/94: fix sgint/dirint to treat Dirac(x) as even
#  11/21/94: make sgD catch jumps at ends of supports
#  12/03/94: evalsg/sgrefine returns list, not expr. seq.
#  02/16/95: when f-g refines to const, return const not error
#  01/12/95: type/segmented now returns false for [[], []] 
#  02/23/96  begin conversion of basic type to table, got through
#            type/segmented, modifications to evalsg
#  04/28/96  fixed plotsg to handle +- infinity as knots
#  04/28/96  fixed convert/ifproc to handle user knotord correctly
#  05/15/96  big revision of evalsg. Fixed bug in type/segmented,
#            fixed evalsg/findtype to handle tables better
#  05/16/96  included segment utility
#  05/22/96  fixed type/segmented: set, not list indices.
#            changed default for sgmapf: now subs piece ranges too.
#  06/5/96   fixed sgint/dirint to handle factors (expr+Dirac)
#  08/11/96  changed evalsg to evaluate non-indexed procedure in the
#            point-evaluation section (for R4: R0 has no problem) 
#  09/4/96   changes to sgint and sgint/dirint to handle Dirac factors
#            better: essentially assume Heaviside is right continous
#  09/16/96  fixed convert/Htoseg to handle table format
#  09/16/96  changed name of "boxcar" to "segment"
#  11/15/96  `evalsg/hastype` replaces hastype in `evalsg/findtype`
#            to beat hastype bug. segment calls error if endpts not ordered
#  02/12/97  tinker `sgint/dirint` to better handle floats in Dirac args
#  02/17/97  more changes to `sgint/dirint` to avoid simplify problems
#            in Maple release 4
#  03/31/97  convert/Htoseg now handles indep. variable correctly
#  05/03/97  sgD now separates out distribution results
#  05/05/97  plotsg handles floating point knots better
#  05/08/97  evalsg/korder revised, treats polyns in infinity as non
#            archimedian above all other orders
#  05/08/97  evalsg/sgrefine now converts all floats in knots to rational
#            so that e.g. 1.0 and 1 are seen as the same knot
#  05/10/97  small changes to accomodate R4 handling of polyns in
#            infinity. Also changed sgint handling of case where
#            integration range is not comparable with knots: use
#            fund. theorem of calculus in all cases.
#  05/29/97  trimsg now indep of evalsg, thus faster
#  05/29/97  type/segmented now excludes ranges like [a,a]
#  05/29/97  if evalsg(expr) is not segmented or procedure, then
#            return just expr (treat like constant) on evalsg(expr,x)
#  10/03/97  sgD now returns non-distribution part on original support
#  11/04/98  included conversions to and from new maple piecewise
#            expressions
#  11/04/98  use `plotsg/pweval/pdum` instead of _pdum in plot routines.
#            (no problems with _pdum, just poor practice to use it).
#  12/11/98  Simpler method to check for x0 comparable with knots in
#            evalsg: avoid possible traperror bug in R5.
#  05/15/99  modified `plotsg/pweval` so that functions to be plotted
#            are left continuous at right end of plotting interval.
#  12/31/99  changed from formal package layout to single script layout
#            to make it easier to load in some situations. Included
#            help reader: uses help but not ?
#  1/26/00   modified `plotsg/pweval` to trap infinity as upper limit.
evalsg := proc(expr,x0,dir)
# If expr has a subexpression of type segmented then evaluate expr to
# a single segmented expression. If x0 exists evaluate at x0. If expr
# has no segmented subexpression and x0 exists, evaluate all proecedure
# subexpressions at x0. If dir exists and is left or right, do left (right)
# evaluation at knots.
# Extensive modification, 5/13/96
local i,j,np,ex1,exr, pwset,scpset,expr1,sl,
 dummy,s,l,rflag,evalpiece,pw,x,expr0,X,tipeseg, evalsubs,X1j;
#
tipeseg := type(expr,segmented);
if nargs=1 and tipeseg then RETURN(expr) fi;
#
# 1st section is for possible 'point evaluation'
if nargs > 1 then
#
  if not tipeseg then
    exr := eval(evalsg(expr));
   else
    exr := eval(expr);
  fi;
  if not type(exr,segmented) then
     if type(exr,procedure) then
      RETURN(exr(x0));
     else 
       RETURN(expr);
     fi;
  else
# Below determine whether x0 is variable or compat. with knots by
#  comparing it with highest knots: knots are totally ordered so
#  this provides a complete check.
#   
   np:=nops(exr[2]);
   s := traperror(`evalsg/korder`(x0,exr[2][np][2]));
   sl := traperror(`evalsg/korder`(x0,exr[2][np][1]));
   if ({s,sl} minus {`true`,`false`}) = {} then
    np := nops(exr[2]);
    if (`evalsg/korder`(x0,op(1,exr[2][1]))
      or `evalsg/korder`(op(2,exr[2][np]),x0) ) then
      RETURN(0);
    fi;
#  below assumes knots in ascending korder
    if ( nargs < 3 or dir = 'right') then
    
     for i from 1 to np do
      if `evalsg/korder`(x0,op(2,exr[2][i])) then
        ex1:= eval(exr[1][i]);
        RETURN(ex1(x0));
      fi;
     od; # fall thru => = last knot
     if nargs < 3 then
      ex1:=eval(exr[1][np]);
      RETURN(ex1(x0)); #left eval for this knot unless
     else # nargs=3 in which case right was called
      RETURN(0);
     fi;
    else # 3rd arg <> right assumed = left
     for i from np by -1 to 1 do
      if `evalsg/korder`(op(1,exr[2][i]),x0) then
        ex1:= eval(exr[1][i]);
        RETURN(ex1(x0));
      fi;
     od;
     RETURN(0); #fall thru => = first knot, and here val=left limit
    fi;
   fi; # here if x0 not korderable:
   if type(x0,name) then 
    s:= map(unapply, exr[1](x0), x0);
    RETURN(table([s, exr[2]]));
   else
    RETURN(expr);
   fi;
  fi; # from else exr is type segmented
fi; # from if nargs > 1
# end 'point evaluate' section
#
# 'symbolic evaluate'
#
if not hastype(expr,{algebraic,segmented}) then RETURN(expr); fi;
#
if not (tipeseg or type(expr,algebraic)) then
  expr1 := expr;
  for i from 1 to nops(expr) do
    ex1 := op(i,expr);
    expr1 := subs(ex1=evalsg(ex1),expr1);
  od;
  RETURN(expr1);
fi;
#
#
pwset := [op(`evalsg/findtype`(expr, segmented))];
#
if pwset = [] then
  RETURN(expr);
else
#
 expr0 := subs({seq(pwset[i]=table([X,i]),i=1..nops(pwset))},expr);
 scpset := `evalsg/findtype`(expr0,procedure);
# if necessary, replace all seg types by representations on common mesh

 expr1 := expr;

 dummy := table([map(0,pwset[1][1]), pwset[1][2]]);
 rflag := `false`;
 for i from 2 to nops(pwset) do
  if not (pwset[i][2] = dummy[2]) then
    dummy :=  `evalsg/sgrefine`(dummy,pwset[i])[1];
    rflag := `true`;
  fi;
 od;
 if rflag = `true` then
   for i from 1 to nops(pwset) do
     expr1 := subs(pwset[i]=`evalsg/sgrefine`(pwset[i],dummy)[1],expr1);
   od;
   pwset := [op(`evalsg/findtype`(expr1,segmented))];
 fi;
fi;
#
#
np := nops(pwset[1][1]);
exr := [];
X1j := (X,x,j)->X=X[1][j](x);
for j from 1 to np do
 evalsubs:= {op(map(X1j,pwset,x,j))};
 evalpiece := subs(evalsubs,eval(expr1));
 evalsubs := map( (X,x)-> X=X(x),scpset,x);
 evalpiece := subs(evalsubs, eval(evalpiece));
 exr := [op(exr),unapply(eval(evalpiece), x)];
od;
RETURN(table([exr,pwset[1][2]]));
end:

`evalsg/korder` := proc(a,b)
local ab, al, bl, bma, big, lc;
if has({a,b},infinity) then
  al:= subs(infinity=big,a);
  bl:= subs(infinity=big,b);
  bma := collect(bl-al, big);
  lc := traperror(lcoeff(bma,big));
  if lc = lasterror then
   ERROR(` can't determine`,a,`<`,b,`: check or define your evalsg/knotord `);
  fi;
  RETURN(`evalsg/korder`(0,lc));
fi;
if assigned(`evalsg/knotord`) then
 ab := traperror(`evalsg/knotord`(a,b));
 RETURN(ab);
fi;
#Here if no infinities and no special routines to order knots:
# try to treat knots as real numbers:
al := evalf(a); bl := evalf(b);
ab := evalb(al<bl);
if (ab=`true` or ab = `false`) then
 RETURN(ab);
else
 RETURN(`FAIL`);
fi;
end:

`evalsg/sgrefine` := proc(pw1,pw2)
# obtain common refinement for pw1 and pw2 partitions, then define
# representations of pw1,pw2 on that partition.
# extensions are by 0.
  local i,npa,la,npb,lb,l,nl,nlm,pwar,pwbr,ra,ja,jb,j,pwa,pwb,v1,v2,
        alo,bhi,korder,x;
 if not type(pw1,segmented)
        or not type(pw2,segmented) then ERROR(`wrong types`); fi;
#
  pwa := eval(pw1); pwb := eval(pw2);
# First convert all floating pts in knots to rational so, eg. 1.0 = 1 true
  pwa := table([pwa[1],convert(pwa[2],rational)]);
  pwb := table([pwb[1],convert(pwb[2],rational)]);
# Now handle potential unequal intervals
  npa := nops(pwa[1]);
  npb := nops(pwb[1]);
  v1 := op(1,pwa[2][1]); v2 := op(1,pwb[2][1]);
  if `evalsg/korder`(v1,v2) then
   npb := npb+1;
   pwb := [[x->0,op(pwb[1])], [[v1,v2],op(pwb[2])]];
  fi;
  if `evalsg/korder`(v2,v1) then
   npa := npa+1;
   pwa := [[x->0,op(pwa[1])], [[v2,v1],op(pwa[2])]];
  fi;
  v1 := op(2,pwa[2][npa]); v2 := op(2,pwb[2][npb]);
  if `evalsg/korder`(v1,v2) then
   npa := npa+1;
   pwa := [[op(pwa[1]),x->0], [op(pwa[2]),[v1,v2]]];
  fi;
  if `evalsg/korder`(v2,v1) then
   npb := npb+1;
   pwb := [[op(pwb[1]),x->0], [op(pwb[2]),[v2,v1]]];
  fi;
 # find the common refinement
 la := {seq(op(1,pwa[2][i]),i=1..npa)} union {op(2,pwa[2][npa])};
 lb := {seq(op(1,pwb[2][i]),i=1..npb)} union {op(2,pwb[2][npb])};
 l := la union lb;
 l := [op(l)];
 l := sort(l,`evalsg/korder`);
 alo := l[1]; bhi := l[nops(l)];
 nl := nops(l); nlm := nl-1;
 # construct the representations of pw1,pw2
 ra := [seq([l[i],l[i+1]],i=1..nlm)];
 ja := 1; jb := 1;
 while not `evalsg/korder`(alo,op(2,pwa[2][ja])) do
  ja := ja+1;
 od;
 while not `evalsg/korder`(alo,op(2,pwb[2][jb])) do
  jb := jb+1;
 od;
 pwar := [seq(0,i=1..nlm)]; pwbr := [seq(0,i=1..nlm)];
 for i from 1 to nlm do
  if `evalsg/korder`(op(2,pwa[2][ja]),op(2,ra[i])) then ja := ja+1; fi;
  if `evalsg/korder`(op(2,pwb[2][jb]),op(2,ra[i])) then jb := jb+1; fi;
  pwar := subsop(i=pwa[1][ja],pwar);
  pwbr := subsop(i=pwb[1][jb],pwbr);
 od;
 RETURN( table([pwar,ra]), table([pwbr,ra]) );
end:
`evalsg/findtype` := proc(expr,tipe)
 # returns set whose members are all distinct instances in expr of the
 #  type tipe
 #  intended to be used with types which don't
 #  have subexpresssions of the same type 
 #  patched 2/23
 local rex, n, i, subex;
 if type(expr,tipe) then RETURN({expr}); fi;
 rex := {};
 if `evalsg/hastype`(expr,tipe) then
   n := nops(expr);
   for i from 1 to n do
    subex := op(i,expr);
    if subex <> NULL then
     if type(subex,tipe) then
       rex := rex union {op(i,expr)};
     else
       rex := rex union `evalsg/findtype`(op(i,expr),tipe);
     fi;
    fi;
   od;
 fi;
 RETURN(rex);
end:

`evalsg/hastype` := proc(ex,tipe)
 local ht;
 if type(ex,tipe) then
   RETURN(`true`);
 elif op(ex)=ex then
   RETURN(`false`);
 else
   if member(`true`,map(`evalsg/hastype`,{op(ex)},tipe)) then
    ht := `true`;
   else
    ht := `false`;
   fi;
   RETURN(ht);
 fi;
end:

`type/segmented` := proc(expr)
# test for the type segmented:
# expr must  be a 2-entry table, both entries being lists:
# The 1st list must be a list of procedures, numeric constants,
# powers of procedures, or functions
# The 2nd list must be a list of pair lists, providing ranges for
# the 1st list. The ranges must be contiguous and must be in
# ascending order.
local i, np, rl, rr, e1, e2;
 if not type(expr,table) then RETURN(`false`); fi;
 if not {indices(expr)} = {[1],[2]} then RETURN(`false`); fi;
 np := nops(expr[1]);
 if np = 0 then RETURN(`false`); fi;
 if np <> nops(expr[2]) then RETURN(`false`); fi;
  for i from 1 to np do
   e1 := expr[1][i];
   if ( type(e1, procedure)
     or type(e1,numeric)
     or (type(e1,`^`) and type(op(1,e1),procedure))
     or (type(e1,`*`) and type(op(2,e1),procedure))
          )
     or type(e1, function)
   then 
     `o.k.`; # dummy branch for o.k.
   else
     RETURN(`false`);
   fi;
   e2 := expr[2][i];
   if not (type(e2, list) and nops(e2)=2) then RETURN(`false`);fi;
   rl := op(1,e2);   
   if i>1 and rl <> rr then RETURN(`false`); fi;
   rr := op(2,e2);
   if traperror(`evalsg/korder`(rl,rr))=`true` then
    `o.k.`;
   else
    if printlevel > 1 then print(lasterror); fi;
    RETURN(`false`);
   fi;
  od;
  RETURN(`true`);
end:

sgint := proc(pw1,r)
local np,nr,ipwf,a0,a1,a,b,c,i,s,ints,t,x,defint,userx,pw1l,Fa,Fb;
# Create segmented expression whose pieces are antiderivs of
# pieces of pw1 such that the result is continuous if possible.
# Return this unless range is given: if range given use antideriv
# to produce definite integral.
#
defint := `false`; # definite integral flag;
userx := `false`; # whether to use user x in eval. integrand
if nargs=2 then
  if type(r,name) then x:= r; userx := `true`; fi;
  if ( type(r,`=`) and type(op(1,r),name) ) then
    x:= op(1,r); userx := `true`; fi;
  if nops(r)=2 then
    if type(r,list) then
       a:= r[1]; b:= r[2]; defint:= `true`;
    elif type(op(2,r),range) then
       a:=op(1,op(2,r)); b:=op(2,op(2,r));  defint:= `true`;
    else fi;
   fi;
fi; # if r is garbage ignore it

if userx = `true` then
  pw1l := evalsg(pw1,x);
else
  pw1l := evalsg(pw1);
fi;
if not type(pw1l, segmented) then 
 RETURN(`sgint/dirint`(pw1l,r,infinity)); 
fi; 
# So now pw1l is segmented:
a0 := op(1,pw1l[2][1]);
if has(a0,infinity) then
   a1 := op(2,pw1l[2][1]);
   if has(a1,infinity) then
     a0 := 0;
   else
     a0 := a1;
   fi;
fi;
np := nops(pw1l[1]); c := 0; ipwf := [];
for i from 1 to np do
 s := pw1l[1][i](t);
 if type(s,polynom(anything,t)) then
   ints := unapply(`pint/sgint`(s,t=a0..x)+c,x);
 else
   a1 := op(2,pw1l[2][i]);
   ints := unapply(`sgint/dirint`(s,t=a0..x,a1)+c,x);
 fi;
 a0 := op(2,pw1l[2][i]); c := ints(a0);
 ipwf := [op(ipwf),eval(ints)];
od;
if a0 = infinity then
 ipwf :=table([ipwf, pw1l[2]]);
else
 ipwf:= table([[op(ipwf),unapply(c,x)], [op(pw1l[2]), [a0,infinity]]]);
fi;
if defint=`true` then
 Fb := evalsg(ipwf,b,left);
 Fa := evalsg(ipwf,a,right);
 RETURN(eval(Fb) - eval(Fa));
else
 RETURN(eval(ipwf));
fi;
end:
`sgint/dirint` := proc(exin0, r,a1)
 local noint, ex, exin, i, exi, n, nd, ar, xsi, xsi0, u, ex0,
              ad, factor, x, a, b, c, Dir, Dirn;
 exin := expand(exin0);
 if type(exin,`+`) then
   ex := 0;
   for i from 1 to nops(exin) do
    ex := ex + `sgint/dirint`(op(i,exin),r,a1);
   od;
   RETURN(ex);
 fi;
 if not ( has(exin,Dirac) and (type(exin,`*`) or type(exin,function)) ) then
   RETURN(int(exin, r));
 fi;
 ex := copy(exin)*c; #c is dummy to make ex have at least 2 factors
 x := op(1,r);
 if type(r,`=`) then
  a := op(1,op(2,r)); b := op(2,op(2,r));
 else
  a:= -infinity; b:= infinity;
 fi;
 for i from 1 to nops(ex) do
   noint := 1;
   exi := op(i,ex);
   if traperror(op(0,exi)) = Dirac then
     n := nops(exi); ar := op(n,exi);
     if degree(ar,x) = 1 then
      if n=1 then nd := 0 else nd := op(1,exi); fi;
       if type(nd,integer) and nd >=0 then
        xsi := [solve(ar=u,x)];
        xsi := xsi[1];
        xsi0 := subs(u=0,xsi); ad := diff(xsi,u);
        if type(nd,even) then
          if `evalsg/korder`(ad,0) then
            ad := -ad;
          fi;
        fi;
        if nd = 0 then
         ex0 := subsop(i=Dir,ex)*ad;
        else
         ex0 := subsop(i=Dirn,ex)*ad;
        fi;
        ex0 := simplify(subs(x=xsi, ex0));
        if a = xsi0 then
         factor := 1;
        elif `evalsg/korder`(a,xsi0) and `evalsg/korder`(xsi0,a1) then
         factor := Heaviside(x-xsi0);
        else
         factor := 0;
        fi;
        if nd=0 then
          ex0 := subs(Dir=factor,ex0);
          ex0 := traperror(subs(u=0, ex0));
        else
          ex0 := subs(Dirn=factor,ex0);
          ex0 := (-1)^nd*traperror(subs(u=0, diff(ex0, u$nd)));
        fi;
        if ex0 = lasterror then
         break; # couldn't substitute
        else
         ex := ex0; noint := 0; break;
        fi;
       fi;
     fi;
   fi;
 od;
 ex := subs(c=1,ex);
 if noint = 1 then
  RETURN(int(ex,r));
 else
  if type(r,`=`) then ex := eval(subs(x=b,ex)); fi;
  RETURN(ex);
 fi;
 end:

#
`pint/sgint` := proc(f, r)
# given polyomial f in x where r is x=a..b,
# produce the integral of
# f from a to b: hopefully avoid bogging down the way direct
# integration can with big problems
  local t, n, i, fl, a, b, s, jf, ap;
  t := op(1,r);
  a := op(1,op(2,r)); b := op(2,op(2,r));
  fl := collect(f, t);
  n := degree(fl,t);
  s := 1; ap := 1;
  jf := subs(t=0,fl);
  for i from 1 to n do
   ap := a*ap;
   s := b*s+ap;
   jf := jf+coeff(fl,t^i)*s/(i+1);
  od;
  RETURN((b-a)*jf);
  end:
#
sgD := proc (fin, x)
# differentiate (segmented) procedures:
# eval and diff w.r.t local z unless x is present, then
# diff w.r.t. x
 local farg, f, np, i, j, df1, w, z, l, r, jd, zro, ddr, dnodr;
 farg := evalsg(fin);
 if not type(farg,segmented) then ERROR( `invalid arguments`); fi;
 if not assigned(Dirac) then readlib(Dirac); fi;
 if nargs > 1 then
   w:= x;
 else
   w:= z;
 fi;
 zro := table([[0],[[-infinity,infinity]]]);
 f := evalsg(farg+zro);

   np := nops(f[1]);
   df1:=[]; ddr:= []; jd := 0;
   for i from 1 to np do
     if i > 1 then
        l := op(2,f[2][i-1]);
        r := op(1,f[2][i]);
        jd := (f[1][i](r)-f[1][i-1](l));
        jd := eval(subs(Dirac=0,jd));
        jd := jd*Dirac(z-r);
     fi;
     df1:= [op(df1), unapply(diff(f[1][i](z),w),z)];
     ddr :=[op(ddr), unapply(jd,z)];
   od;
   dnodr := table([df1,f[2]]);
   np:= nops(farg[2]);
   dnodr := trimsg(dnodr, farg[2][1][1],farg[2][np][2]);
   if {op(ddr)} = {0} then
     ddr := 0;
   else
     ddr := table([ddr,f[2]]);
   fi;
   RETURN(eval(dnodr)+eval(ddr));
 end:
#
sgmapf := proc(s,pwf,extras)
local efnlist,x,i,n,pwfl;
# apply the function s to the evaluated pieces of pwf:
# return resulting segmented procedure
# to apply to fun variable, enter pwf as [sg, x] and refer to x in extras
# 
if not type(pwf,segmented) then
  if (nops(pwf)=2 and type(pwf[1],segmented)) then
    pwfl := pwf[1]; x := pwf[2];
  else
    ERROR(pwf, `is the wrong type`);
  fi;
else
  pwfl := pwf;
fi;
n := nops(pwfl[1]); efnlist := [];
for i from 1 to n do
  efnlist := [op(efnlist),pwfl[1][i](x)];
od;
if nargs < 3 then
  efnlist := map(s,efnlist);
 else
  efnlist := map(s,efnlist,args[3..nargs]);
fi;
efnlist := traperror(map(unapply,efnlist,x));
if efnlist = lasterror then
  ERROR(` is `,x,` the variable in the segmented expression ?`);
fi;
RETURN(table([efnlist, pwfl[2]]));
end:

sgsubs := proc(sarg, pwf, x0)
# return segmented expression whose pieces are the
#  result of substituting into those of pwf according to the
#  equation(s) in s. Default now is subs in ranges too.
 local x, t, s, i, n, d, f, pwfs, rangeflag, r;
 if not type(pwf,segmented) then ERROR( ` wrong argument type `); fi;
 if nargs = 2 then
   t:=x
  elif nargs = 3 then
    t:=x0;
  else ERROR(` wrong # of arguments `); 
 fi;
 if type(sarg,list) then 
   s := sarg[1];
 else
   s := sarg; 
 fi;
 if type(sarg,list) and nops(sarg)=2 and sarg[2]=`nr` then
  rangeflag := 0;
 else
  rangeflag := 1;
 fi;
 pwfs := [];
 n := nops(pwf[1]);
 for i from 1 to n do
  f := pwf[1][i](t);
  d := subs(s,f);
  pwfs := [op(pwfs), unapply(d,t)];
 od;
 if rangeflag = 1 then
  r := subs(s, pwf[2]);
 else
  r := pwf[2];
 fi;
 RETURN(table([pwfs,r]));
 end:

`type/mmsp` := proc(expr)
# test for spline expressed in formalism of Maple bspline routine
# expr must  be a list np>=2 long of 2-member lists
# expr[i] must be of form [a,b] where a is [x < k[i]], i=1,2,..,np-1
# [x >= k[np]]
# where x is fixed
# and k[1],...,k[np] are the knots.
# and where b is an expression.
local i, np, x;
 if not type(expr,listlist) then
  RETURN(`false`);
 fi;
 np := nops(expr);
 if np < 2 then RETURN(`false`); fi;
 x := op(1,expr[1][1]);
 for i from 1 to np-1 do
   if not (type(expr[i][1],`<`) and op(1,expr[i][1])=x) then
    RETURN(`false`);
   fi;
 od;
   if not type(expr[np],list) then RETURN(`false`); fi;
   if not (type(expr[np][1],`<=`) and op(2,expr[i][1])=x) then
    RETURN(`false`);
   fi;
  RETURN(`true`);
end:

`convert/segtommsp` := proc(sp,v)
# convert splines from segmented format to Maple bspline format
local np, x, a, i, msp;
if type(sp,segmented) then
 np := nops(sp[1]);
 if nargs = 2 then x := v; fi;
 a := op(1,sp[2][1]);
 msp := [[x < a, 0]];
 for i from 1 to np do
  a := op(2,sp[2][i]);
  msp := [op(msp), [x < a, sp[1][i](x)]];
 od;
 msp := [op(msp), [a <= x, 0]];
 RETURN(msp);
else
 ERROR(` arg must be a segmented expression `);
fi;
end:

`convert/mmsptoseg` := proc(sp)
# convert splines from Maple bspline format to segmented format
# and vice versa
local np, x, fnpieces, ranges, i;
if type(sp,mmsp) then
 np := nops(sp);
 x := op(1,sp[1][1]);
 fnpieces := []; ranges := [];
 for i from 1 to np do
    fnpieces :=[op(fnpieces), unapply(sp[i][2],x)];
    if i = 1 then
     ranges := [[-infinity,op(2,sp[1][1])]];
    elif i=np then
     ranges := [op(ranges),[op(2,sp[np-1][1]),infinity]];
    else
     ranges := [op(ranges), [op(2,sp[i-1][1]),op(2,sp[i][1])]];
    fi;
 od;
 RETURN(table([fnpieces, ranges]));
else
 ERROR(` arg must be Maple bspline format spline `);
fi;
end:

plotsg := proc(s, moreargs)
#Extend the semantics of plot to segmented type
local sl, x, pwset, ifset, f, r, a, b;
if nargs >= 2 and type(args[2], name=range) then
  r := op(2, args[2]);
  a := op(1,r); b := op(2,r);
  x := op(1, args[2]);
elif nargs >= 2 and type(args[2],range) then
  r:= args[2];
  a:= op(1,r); b:= op(2,r);
else
  a := -10; b := 10;
fi;
pwset := `evalsg/findtype`(s,segmented);
ifset := map((f,a,b) -> f=`plotsg/pweval`(f,a,b)(x), pwset,a,b);
sl := subs(ifset, s);
plot(sl, args[2..nargs]);
end:

`plotsg/pweval` := proc(f,a,b)
# fix zero right limit problem, 5/99: extend f to the right
# by its left limit value at end knot. The extension is on
# tiny support, so won't show in a plot.
 local pp, ip, fl, al, bl, bf, flb, n, newend;
 global `plotsg/pweval/pdum`;
 `plotsg/pweval/pdum` := '`plotsg/pweval/pdum`';
  fl:= map(evalf,f);
  al:= evalf(a);
  bl:= evalf(b);
  n:= nops(fl[2]);
  bf:= fl[2][n][2];
  if has(bf,infinity) then
    ERROR(`can't handle infinity in this context`);
  elif bf = 0 then
    newend := Float(1,-8);
  else
    newend := bf +abs(bf)*Float(1,-8);
  fi;
  flb:= evalf(fl[1][n](bf));
  fl:= table([[op(fl[1]),flb], [op(fl[2]),
   [bf, newend]]]);
  pp := convert(eval(fl), ifproc);
  ip := subs(`plotsg/pweval/pdum`=eval(pp),
    proc(x) if type(x,numeric) then
      `plotsg/pweval/pdum`(x) else 
       '`plotsg/pweval/pdum`(x)'; fi; end
             );
  RETURN(eval(ip));
 end:

trimsg := proc(f,a,b,t)
#Returns segmented object whose support is the intersection
# of [a,b] with the support of f
local n, i, ia, ib, leftrng, leftfn, rtrng, rtfn, newfns, newrngs;
if not type(f,segmented) then ERROR(`f must be segmented`); fi;
if not traperror(`evalsg/korder`(a,b)) then
 ERROR(`a,b incorrectly ordered`);
fi;
n:= nops(f[2]);
ia:= n;
for i from n by -1 to 1 while traperror(`evalsg/korder`(a,op(1,f[2][i]))) do
 ia := ia-1;
od;
ib:= 1;
for i from 1 to n  while traperror(`evalsg/korder`(op(2,f[2][i]),b)) do
 ib := ib+1;
od;
if ia=ib then
 RETURN(table([ [f[1][ia]], [[a,b]] ]));
fi;
if ia=0 and nargs=3 then
  leftrng:= [a,op(1,f[2][1])]; leftfn := 0;
elif (ia=0 and nargs=4) or not `evalsg/korder`(a,op(2,f[2][ia])) then
  leftrng:= NULL; leftfn:= NULL;
else
    leftrng:= [a,op(2,f[2][ia])]; leftfn := f[1][ia];
fi;
if ib = n+1 and nargs=3 then
  rtrng := [op(2,f[2][n]),b]; rtfn := 0;
elif (ib=n+1 and nargs=4) or not `evalsg/korder`(op(1,f[2][ib]),b) then
  rtrng:= NULL; rtfn:= NULL;
else
  rtrng := [op(1,f[2][ib]),b]; rtfn := f[1][ib];
fi;
newfns := [eval(leftfn),seq(f[1][i],i=ia+1..ib-1),eval(rtfn)];
newrngs := [leftrng,seq(f[2][i],i=ia+1..ib-1),rtrng];
if newfns = [] then
  RETURN(0);
else
 RETURN(table([newfns,newrngs]));
fi;
end:

`convert/ifproc` := proc(f, symknots)
#convert type segmented expression to Maple procedure with if/then structure
#assume usual < unless called with 2 args
local i, iffun, xx, t, fm, rel;
if not type(f, segmented) then ERROR(` wrong type of argument`); fi;
fm := convert(f, segtommsp, t);
if nargs = 2 then
 for i from 1 to nops(fm)-1 do
   rel := fm[i][1];
   fm := subs(rel='`evalsg/korder`'(t, op(2,rel)), fm);
 od;
 rel := fm[nops(fm)][1];
 fm := subs(rel='not `evalsg/korder`'(t, op(1,rel)), fm);
fi;
iffun := &if(seq(xx[i],i=1..2*nops(fm)));
for i from 1 to nops(fm) do
  iffun := subs(xx[2*i-1]=fm[i][1], xx[2*i]=fm[i][2], iffun);
od;
readlib(procmake);
iffun := &proc([t], [], [], iffun);
RETURN(procmake(iffun));
end:

`convert/segtoH` := proc(sp,v)
# convert splines from segmented format to expression in terms of Heaviside
# functions
local np, x, a, b, i, u, w, hpw;  
if type(sp,segmented) then
 np := nops(sp[1]); 
 if nargs = 2 then x := v; fi;
 u := sp[1][1](x); a:= op(1,sp[2][1]); hpw := u;
 if not has(a,infinity) then hpw := Heaviside(x-a)*u; fi;
 for i from 2 to np do
  w := sp[1][i](x);
  a := op(1,sp[2][i]);
  hpw := hpw + Heaviside(x-a)*(w-u);
  u := copy(w);
 od;
 a := op(2,sp[2][np]);
 if not has(a,infinity) then
   hpw := hpw - Heaviside(x-op(2,sp[2][np]))*u;
 fi;
 RETURN(hpw);  
else
 ERROR(` arg must be a segmented expression`);
fi;
end: 

`convert/Htoseg` := proc(exin,x)
#convert expression in terms of Heaviside functions to
#segmented type. Heaviside args must be linear in x.
local ex0,ex1,i,n;
ex0 := `convert/Htoseg/uneval`(exin,x);
#ex0 := evalsg(ex0,x);
if type(ex0,segmented) then
 n := nops(ex0[1]);
 if ex0[1][n]=0 then
  ex0 := table([[seq(ex0[1][i],i=1..n-1)],[seq(ex0[2][i],i=1..n-1)]]);
  n := n-1;
 fi;
 if ex0[1][1]=0 then
  ex0 := table([[seq(ex0[1][i],i=2..n)],[seq(ex0[2][i],i=2..n)]]);
 fi;
fi;
RETURN(eval(ex0));
end:
`convert/Htoseg/uneval` := proc(exin,x)
 local ex, i, exi,ar, dar, HF, xsi,u, ex0,c,Hspln;
 if type(exin,`+`) then
   ex := 0; 
   for i from 1 to nops(exin) do
    ex := evalsg(ex + `convert/Htoseg/uneval`(op(i,exin),x),x); 
   od; 
   RETURN(eval(ex)); 
 fi; 
 if not (type(exin,`*`) or type(exin,function)) then
   RETURN(exin); 
 fi; 
 ex0 := copy(exin)*c; #c is dummy to make ex have at least 2 factors
 for i from 1 to nops(ex0) do
   exi := op(i,ex0); 
    if type(exi,function) then
      if op(0,exi)=Heaviside then
        ar := op(1,exi); xsi := [solve(ar=0,x)]; 
        xsi := traperror(xsi[1]); 
        if xsi = lasterror then break; fi; #couldn't solve
        dar := subs(x=xsi, diff(ar,x));
        if `evalsg/korder`(0,dar) then
         HF := [0,1];
        elif `evalsg/korder`(dar,0) then
         HF := [1,0];
        else break; #because can't find sign of ar derivative
        fi;
        Hspln := table([HF,[[-infinity,xsi], [xsi,infinity]]]);
        ex0 := subs(exi=eval(Hspln),ex0);
      fi;
    fi;
 od;
 ex0 := evalsg(subs(c=1,ex0),x);
 RETURN(eval(ex0));
 end:

segment := proc(a,b,f,x0)
#Attemp user friendly entry of segmented equivalent of
#(Heavside(b)-Heaviside(a))*f, represented as segmented object
#See help document
local x, X, fl, nameset, procset, evalsubs, i;
if not traperror(`evalsg/korder`(a,b)) = `true` then
  ERROR(a,b,` aren't ordered `);
fi;
if nargs=4 then
 if type(x0,name) and not member(x0, {constants}) then
  x := x0;
 else
  ERROR(x0, `is not a variable but should be`);
 fi;
fi;
procset := `evalsg/findtype`(f,procedure);
if nargs = 3 then
 nameset := `evalsg/findtype`(f,name);
 nameset := nameset minus {constants};
 nameset := nameset minus procset;
 if nameset <> {} then
  x := nameset[1];
 fi;
fi;
if type(f, procedure) then 
 fl:= f(x);
else
 evalsubs:= map((X,x) -> X = X(x), procset,x);
 fl := subs(evalsubs, f);
fi;
RETURN(table([[unapply(fl,x)],[[a,b]]]));
end:
#
`convert/segtopw` := proc(f,x)
#convert segmented to piecewise function a la Mohrenschildt
local fl, fo, i, m, n;
if nargs = 2 then
 fl := convert(f,segtommsp,x);
else
 fl := convert(f,segtommsp);
fi;
m:= 1;
n:= nops(fl);
fo := piecewise(seq(op(fl[i]),i=m..n));
RETURN(fo);
end:
#
`convert/pwtoseg`:= proc(f,x)
#convert from Mohrenschild format to segmented
 local fl, fo, i, fa, n, a, b;
 n:= nops(f);
 if not type(n,even) then
   ERROR(` wrong format for f `);
 else
   fa := [op(f)];
 fi;
 fl := [seq([fa[2*i-1],fa[2*i]],i=1..n/2)];
 if type(fl, mmsp) then
   fo := convert(fl,mmsptoseg,x);
   n:= nops(fo[1]);
   if fo[1][1]=0 then
     a:= fo[2][1][2];  b:= fo[2][n][2];
     fo := trimsg(fo,a,b);
   fi;
   n:= nops(fo[1]);
   if fo[1][n]=0 then
     a:= fo[2][1][1]; b:= fo[2][n][1];
     fo := trimsg(fo,a,b);
   fi;
   RETURN(eval(fo));
 else
   ERROR(` wrong type for f `);
 fi;
end:
               
`help/text/evalsg` := TEXT(
`FUNCTION: evalsg - evaluate expression involving segmented expressions`,
`   `,
`CALLING SEQUENCE:`,
`   evalsg(<expression>,<evaluation variable>,<direction>)`,
`   `,
`PARAMETERS:`,
`  <expression> - an expression possibly involving segmented expressions`,
`  <evaluation variable> - optional: if present it is used in final`,
`   evaluation`,
`  <direction> - optional: if present should be left or right`,
`   `,
`SYNOPSIS:   `,
`- The function evalsg evaluates an expression involving objects of type`,
`  segmented, in one of three ways, depending on <evaluation variable>:`,
`   1) with no evaluation variable, the expression is evaluated to an`,
`      object of type segmented whose pieces are in terms of x distinct`,
`      from an x defined external to evalsg.`,
`   2) with an evaluation variable x0 which is a name and which also is`,
`      not comparable with the breakpoints, the expression is`,
`      evaluated to a segmented object in terms of x0. This enables`,
`      non-procedure subexpressions in x0 to be included in <expression>`,
`      and folded into the result.`,
`   3) with an evaluation variable x0 which is comparable with the`,
`      breakpoints, the expression is evaluated to a segmented object`,
`      which is then evaluated at x0. If x0 is outside the explicit`,
`      domain of the evaluated segmented ojbect then 0 is returned.`,
`      If direction variable is present and is = left, then if x0 is a`,
`      breakpoint the value is that of the left piece procedure, while`,
`      if direction is = right, the value is that of the right procedure.`,
`      If direction variable is omitted or is nonsense, right evaluation`,
`      is performed at breakpoints.`,
`   `,
`   `,
`- An argument which is an ordinary procedure is evaluated at <evaluation`,
`  variable> if it is present, otherwise passed unevaluated.`,
`   `,
`- Segmented expressions with different breakpoints are admissible; the`,
`  set of breakpoints for the result is union of the breakpoints for the`,
`  operands in the initial expression. When procedures must be extended`,
`  beyond their original explicit domains they are extended by zero.`,
`   `,
`- By default all breakpoints are assumed to be real with the usual real`,
`  order. Users may employ their own orders on breakpoints by`,
`  defining the procedure ``evalsg/knotord``. It should implement a total`,
`  order on the user's set of breakpoints. Thus, for example, symbolic`,
`  breakpoints can be used. In case a,b are not comparable in the user's`,
`  order,``evalsg/knotord``(a,b) should return FAIL, rather than a`,
`  deliberate or unanticipated ERROR.`,
`   `,
`   `,
`EXAMPLES:   `,
`   `,
`   `,
`f := table([[x->2*x,x->3-x],[[0,1], [1,3] ]]);`,
`                         table([`,
`                             1 = [x -> 2 x, x -> 3 - x]`,
`                             2 = [[0, 1], [1, 3]]`,
`                         ])`,
`   `,
`g :=table([[x->1,x->2-x],[[-1,1], [1,2] ]]);`,
`                           table([`,
`                               1 = [1, x -> 2 - x]`,
`                               2 = [[-1, 1], [1, 2]]`,
`                           ])`,
`   `,
`h:=evalsg(f+g);`,
`              table([`,
`                  1 = [1, x -> 2 x + 1, x -> 5 - 2 x, x -> 3 - x]`,
`                  2 = [[-1, 0], [0, 1], [1, 2], [2, 3]]`,
`              ])`,
`   `,
`evalsg(f+t*g);`,
`       table([`,
`           1 = [x -> t, x -> 2 x + t, x -> 3 - x + t (2 - x), x -> 3 - x]`,
`           2 = [[-1, 0], [0, 1], [1, 2], [2, 3]]`,
`       ])   `,
`   `,
`evalsg(f+t*g,t);`,
`         table([`,
`             1 = [t -> t, t -> 3 t, t -> 3 - t + t (2 - t), t -> 3 - t]`,
`             2 = [[-1, 0], [0, 1], [1, 2], [2, 3]]`,
`         ])`,
`   `,
`evalsg(f^2);`,
`                       table([`,
`                                        2              2`,
`                           1 = [x -> 4 x , (x -> 3 - x) ]`,
`                           2 = [[0, 1], [1, 3]]`,
`                       ])`,
`   `,
`evalsg(sqrt(%),1.5);`,
`                                      1.5`,
`   `,
`#SEE ALSO:  ``type/segmented``, unapply, sgD, sgint, sgsubs, sgmapf `,
`   `
):
`help/text/sgint` := TEXT(
`FUNCTION: sgint - Integrator for segmented procedures`,
`   `,
`CALLING SEQUENCE:`,
`   sgint(f)`,
`   sgint(f,x)`,
`   sgint(f, x=a..b)`,
`   `,
`PARAMETERS:`,
`   f - algebraic expression which may be of type segmented`,
`   x - variable of integration`,
`   a..b - range of integration`,
`   `,
`SYNOPSIS:   `,
`- sgint first internally replaces f by evalsg(f,x). If this result`,
`  is not of type segmented then sgint passes it to the Maple`,
`  procedure int, with other arguments if supplied. If the evalsg result`,
`  is of type segmented, then sgint creates a segmented expression`,
`  which is an antiderivative of f, with respect to x if supplied. infinity`,
`  is the rightmost breakpoint of the antiderivative. With the possible`,
`  exception of this point the breakpoints are the same as those of f.`,
`  If no integration range is given this antiderivative is returned. If`,
`  a range is supplied then the antiderivative is used in the usual way`,
`  to give a definite integral. The endpoints a and b in the range should`,
`  be comparable with the breakpoints.`,
`   `,
`   `,
`EXAMPLES:   `,
`   `,
`   `,
`   `,
`   `,
`   `,
`ABS := table([[x-> -x, x->x], [[-infinity,0], [0,infinity]]]):`,
`sgint(exp(-x)*ABS,x=-1..1);`,
`                                       2`,
`                                  - ------ + 2`,
`                                    exp(1)`,
`   `,
`sgint(ABS^2,x);`,
`                    table([`,
`                                       3            3`,
`                        1 = [x -> 1/3 x , x -> 1/3 x ]`,
`                        2 = [[- infinity, 0], [0, infinity]]`,
`                    ])`,
`   `,
`   `,
`SEE ALSO:  sgD,  evalsg,  int`,
`   `
):
`help/text/sgsubs` := TEXT(
`FUNCTION: sgsubs - substitute into segmented expressions`,
`   `,
`CALLING SEQUENCE:`,
`   sgsubs(s, f)`,
`   `,
`PARAMETERS:`,
`   s - equation or set of equations acceptible to subs, or a list [e,`,
`       expr] where e is an equation or set of equations and expr is any`,
`       expression at all`,
`   f - segmented expression`,
`   `,
`SYNOPSIS:   `,
`- The pieces of f are evaluated at a variable x local to sgsubs and thus`,
`  distinct from any variables at the level from which sgsubs is called.`,
`  Substitution is performed on the resulting expressions according to the`,
`  equations in s (e). The segmented expression is formed and returned whose`,
`  pieces evaluate to the results of the substitutions. Note that the`,
`  dependence on x cannot be "substituted out" because x is local to sgsubs.`,
`  Presenting s as a list [e,``nr``] is used to signal that the equations`,
`  in e should not be used for substitution in the ranges of f; otherwise`,
`  the ranges are affected by the substitutions.`,
`   `,
`EXAMPLES:   `,
`   `,
`f:= table([[x->2*x, x->2*(4-x)], [[0,2],[2,4]]]):`,
`sgsubs({2=3,-2=-3}, f);`,
`                        table([`,
`                            1 = [x -> 3 x, x -> 8 - 3 x]`,
`                            2 = [[0, 3], [3, 4]]`,
`                        ])`,
`   `,
`sgsubs([{2=3,-2=-3},'nr'],f);`,
`                        table([`,
`                            1 = [x -> 3 x, x -> 8 - 3 x]`,
`                            2 = [[0, 2], [2, 4]]`,
`                        ])`,
`   `,
`SEE ALSO:  subs, evalsg`,
`   `,
`   `,
`   `,
`   `,
`   `,
`   `,
`   `
):

`help/text/sgD` := TEXT(
`FUNCTION: sgD - differentiator for segmented expressions`,
`   `,
`CALLING SEQUENCE:`,
`   sgD(f)   `,
`   sgD(f,x)`,
`   `,
`PARAMETERS:`,
`   f - expression that evaluates to type segmented`,
`   x - argument of type name to differentiate with respect to`,
`   `,
`SYNOPSIS:   `,
`- The call sgD(f) first forms evalsg(f). If the result is of type`,
`  segmented, sgD evaluates each piece procedure at z local to sgD`,
`  differentiates with respect to z or x if present, and returns the`,
`  segmented expression whose pieces correspond to these derivatives.`,
`  sgD uses shifted Dirac functions to handle discontinuities`,
`  at range boundaries.`,
`   `,
`EXAMPLES:   `,
`f:=table([[x -> 1/3 + 1/3*x, x -> 7 - x], [[2, 5], [5, 6]]]):`,
`sgD(f);`,
`       table([`,
`               1 = [1/3,-1]`,
`               2 = [[2,5],[5,6]]`,
`                                ]) + table([`,
`               1 = [0, z -> Dirac(z-2), 0, -(z -> Dirac(z-6))]`,
`               2 = [[-infinity, 2], [2, 5], [5, 6], [6, infinity]]`,

`                                                                 ])`,
`sgD(%);`,
`       table([`,
`               1 = [0, z -> Dirac(1, z - 2), 0, -(z -> Dirac(1, z - 6))]`,
`               2 = [[-infinity, 2], [2, 5], [5, 6], [6, infinity]]`,
`                                                            ]) + table([`,
`               1 = [0, z -> 1/3 Dirac(z - 2), z -> - 4/3 Dirac(z - 5),`,
`                    z -> Dirac(z - 6)]`,
`               2 = [[-infinity, 2], [2, 5], [5, 6], [6, infinity]]`,
`                                                                    ])`,
`tf:=evalsg(t*f);`,
`table([`,
`       1 = [x -> t (1/3 + 1/3 x), x -> t (7 - x)]`,
`       2 = [[2, 5], [5, 6]]`,
`                                                ])`,
` sgD(tf,t);`,
` table([`,
`       1 = [z -> 1/3 + 1/3 z, z -> 7 - z]`,
`       2 = [[2, 5], [5, 6]]`,
`                                          ]) + table([`,
`       1 = [0, z -> t Dirac(z - 2), 0, z -> -t Dirac(z - 6)]`,
`       2 = [[-infinity, 2], [2, 5], [5, 6], [6, infinity]]`,
`                                                            ])`,
`   `,
`SEE ALSO:  diff, sgint, Dirac, evalsg`
):
`help/text/sgmapf` := TEXT(
`   `,
`FUNCTION: sgmapf - apply transformations to segmented expressions`,
`   `,
`CALLING SEQUENCE:`,
`   sgmapf(proc, f)`,
`   sgmapf(proc, f, extras)`,
`   `,
`PARAMETERS:`,
`   proc   - a procedure to be apsgied to evaluated pieces of f`,
`   f      - segmented expression`,
`   extras - additional arguments to proc, if necessary`,
`   `,
`SYNOPSIS:   `,
`- The pieces of f are evaluated at a variable x local to sgmapf and thus`,
`  distinct from any variables at the level from which sgmapf is called.`,
`  The procedure proc is applied to the evaluated pieces. The segmented`,
`  expression is then created whose pieces evaluate to the proc-transformed`,
`  results. This segmented expression is returned.`,
`   `,
`EXAMPLES:   `,
`   `,
`f := table([[x->sin(x+a), x->cos(a-x)], [[0,1], [1,2]] ]);`,
`                   table([`,
`                       1 = [x -> sin(x + a), x -> cos(a - x)]`,
`                       2 = [[0, 1], [1, 2]]`,
`                   ])`,
`   `,
`sgmapf(expand,f,trig);`,
`table([   `,
`    1 =   `,
`   `,
`       [x -> sin(x) cos(a) + cos(x) sin(a), x -> cos(a) cos(x) + sin(a) sin(x)]`,
`    2 = [[0, 1], [1, 2]]`,
`])   `,
`   `,
`f := table([[x->(x^2-1)/(x-1), x->(x^3-1)/(x-1)], [[0,1], [1,2]] ]):`,
`sgmapf(simplify,f);`,
`                     table([`,
`                                                2`,
`                         1 = [x -> x + 1, x -> x  + x + 1]`,
`                         2 = [[0, 1], [1, 2]]`,
`                     ])`,
`   `,
`SEE ALSO:  evalsg, map`,
`   `
):
`help/text/type/segmented` := TEXT(
`   `,
`FUNCTION: type/segmented - check for a segmented expression`,
`   `,
`CALLING SEQUENCE:`,
`   type(f, segmented)`,
`   `,
`PARAMETERS:`,
`   f - an expression`,
`   segmented - name of the type`,
`   `,
`SYNOPSIS:   `,
`- This function checks whether f is a segmented expression.`,
`  A segmented expression is represented as a two-entry table whose`,
`  entries are each lists. The entries of the first list are procedures`,
`  with one argument, constants, or integer powers of procedures with`,
`  one argument. The entries of the second list are pair lists where the`,
`  i-th gives the range for the i-th piece procedure. The ranges`,
`  should be such that if the i-th is [a,b] and the i+1-st is [c,d] then`,
`  a < b = c < d holds.`,
`   `,
`   `,
`EXAMPLES:   `,
`   `,
`   `,
`f := table([[x-> x-1, x->3-x], [[1,2], [2,3]]]):`,
`type(f,segmented);`,
`                                      true`,
`   `,
`f := table([[x->3-x, x-> x-1], [[2,3], [1,2]]]):`,
`type(f,segmented);`,
`                                     false`,
`   `,
`f := table([[x^2, x^3], [[-1,0], [0,1]]]):`,
`type(f,segmented);`,
`                                     false`,
`   `,
`f := table([[x->x^2, x->x^3], [[-1,0],[0,1]]]):`,
`type(f,segmented);`,
`                                      true`,
`   `,
`SEE ALSO:  evalsg, convert/mmsptoseg, convert/segtommsp, convert/ifproc`,
`   `
):
`help/text/plotsg` := TEXT(
`FUNCTION: plotsg - 2-dimensional plots for segmented expressions`,
`   `,
`CALLING SEQUENCE:`,
`   plotsg(f, h, v)`,
`   plotsg(f, h, v,...)`,
`   `,
`PARAMETERS:`,
`   f - segmented expression(s) to be plotted`,
`   h - horizontal range`,
`   v - vertical range (optional)`,
`   `,
`SYNOPSIS:   `,
`- This function extends the functionality and calling syntax of Maple's`,
`  2D plot routine to the class of segmented expressions. Maple's`,
`  help documentation for plot applies. If f is a segmented expression it`,
`  should be referred to as f, not f(x), likewise for any segmented`,
`  expressions in f if f is a set or list (see plot[function]). Ordinary`,
`  functions and expressions may be referred to just as with Maple's plot.`,
`   `,
`SEE ALSO:  plot, evalsg.`
):
`help/text/type/mmsp` := TEXT(
`   `,
`   `,
`FUNCTION: type/mmsp - check for a spline in old Maple spline format`,
`   `,
`CALLING SEQUENCE:`,
`   type(f, mmsp)`,
`   `,
`PARAMETERS:`,
`   f - an expression`,
`   mmsp - name of the type`,
`   `,
`SYNOPSIS:   `,
`- This function checks whether f is in the following format: f is a list`,
`  used in the Maple bspline routine. Such a structure is a list`,
`  of n 2-member lists, n at least 2. f[i] must be of form [a,b] where`,
`  a is [x < k[i]], i=1,2,...,n-1, a is [x >= k[n]] for i=n, and where each`,
`  b is an expression.`,
`   `,
`f := table([[x-> x-1, x->3-x], [[1,2], [2,3]]]):`,
`type(f,mmsp);`,
`                                     false`,
`   `,
`f := [[x<1,0], [x<2, x-1], [x<3, 3-x], [x>=3,0]]:`,
`type(f,mmsp);`,
`                                      true`,
`   `,
`SEE ALSO:  type/segmented, evalsg, piecewise`
):
`help/text/convert/pwtoseg` := TEXT(
`   `,
`FUNCTION: convert/pwtoseg - convert Maple piecewise expression to`,
` object of type segmented.`,
`   `,
`CALLING SEQUENCE:`,
`   convert( g, pwtoseg,x)`,
`   `,
`PARAMETERS:`,
`   g - Maple piecewise expression`,
`   x - g is assumed to be a function of x`,
`   `,
`SYNOPSIS:   `,
`- convert/pwtoseg converts g in Maple piecewise format`,
`  into an object  h such that type(h,segmented) = true. This`,
`  conversion is "user defined" and thus doesn't appear in Maple's`,
`  list of topics known to the basic convert routine.`,
`   `,
`EXAMPLE:    `,
`   `,
`   `,
`g := piecewise(x < 1, 0, x < 2, x - 1, x < 3, 3 - x, 3 <= x, 0):`,
`convert(g, pwtoseg);`,
`   `,
`          table([`,
`               1 = [0, x -> x - 1, x -> 3 - x, 0]`,
`               2 = [[- infinity, 1], [1, 2], [2, 3], [3, infinity]]`,
`           ])`,
`   `,
`SEE ALSO: type/segmented, convert, bspline`
):
`help/text/convert/segtopw` := TEXT(
`   `,
`FUNCTION: convert/segtopw - convert segmented object to piecewise function`,
`   `,
`CALLING SEQUENCE:`,
`   convert( f, segtopw, x)`,
`   `,
`PARAMETERS:`,
`   f - segmented object`,
`   x - variable of type name`,
`   `,
`SYNOPSIS:   `,
`- convert/segtopw converts a segmented expression f (structure for which`,
`  type(f, segmented) is true) into a piecewise expression.`,
`  This conversion is "user defined" and thus doesn't appear in`,
`  Maple's default list of topics known to the basic convert routine.`,
`   `,
`EXAMPLE:   `,
`   `,
`f := table([[x -> x - 1, x -> 3 - x], [[1, 2], [2, 3]]]):`,
`g := convert(f,segtopw);`,
`          piecewise(x<1,0,x<2,x-1,x<3,-x+3,3<=x,0)`,
`   `,
`SEE ALSO: type/segmented, type/mmsp, convert, bspline`,
`   `
):
`help/text/convert/ifproc` := TEXT(
`   `,
`FUNCTION: convert/ifproc - convert segmented expression to procedure`,
`using if clauses`,
`   `,
`CALLING SEQUENCE:`,
`   convert( f, ifproc, symb)`,
`   `,
`PARAMETERS:`,
`   f      - segmented expression`,
`   ifproc - the string ifproc itself`,
`   symb   - if present triggers use of ``evalsg/korder```,
`   `,
`SYNOPSIS:   `,
`- converts a segmented expression f to a single procedure in whose`,
`  definition each piece of f corresponds to a clause. The latter format`,
`  is acceptible to some Maple routines which don't handle segmented`,
`  structures directly. An example is numerical solution of differential`,
`  equations using dsolve when coefficients or even the equation itself`,
`  is of type segmented. For numerical solution dsolve tries to use evalhf`,
`  which can't handle the type segmented directly. It can handle the result`,
`  of applying convert/ifproc when the breakpoints are numerical.`,
`  When called with a 3rd argument (whose value is immaterial), the`,
`  conditions are in terms of ``evalsg/korder``, allowing user-defined`,
`  orders on knots to be applied in the converted procedure. Such results`,
`  likely won't be acceptable to evalhf.`,
`  This conversion is 'user defined' and thus doesn't appear in the`,
`  default Maple help document on convert.`,
`   `,
`EXAMPLE:   `,
`   `,
`f := table([[x->x+1,x->x-1],[[-1,0],[0,1]]]):`,
`fproc := convert(f, ifproc);`,
`   `,
`proc(t)   `,
`    if t < -1 then 0`,
`    elif t < 0 then t+1`,
`    elif t < 1 then t-1`,
`    elif 1 <= t then 0`,
`    fi   `,
`end   `,
`   `,
`gproc := convert(f, ifproc, symb);`,
`   `,
`proc(t)   `,
`    if ``evalsg/korder``(t,-1) then 0`,
`    elif ``evalsg/korder``(t,0) then t+1`,
`    elif ``evalsg/korder``(t,1) then t-1`,
`    elif not ``evalsg/korder``(t,1) then 0`,
`    fi   `,
`end   `,
`   `,
`evalhf(evalsg(f, 0.1));`,
`Error:   unable to evaluate expression to hardware floats`,
`evalhf(fproc(0.1));`,
`                               -.9000000000000000`,
`   `,
`SEE ALSO: type[segmented], convert, evalhf`
):
`help/text/convert/segtoH` := TEXT(
`FUNCTION: convert/segtoH - convert segmented procedure to Heaviside formula`,
`   `,
`CALLING SEQUENCE:`,
`   convert(f, segtoH, x)`,
`   `,
`PARAMETERS:`,
`   `,
`   f - expression of type segmented`,
`   segtoH - indicates conversion from segmented to Heaviside`,
`   x - variable name to use in Heaviside expression`,
`   `,
`SYNOPSIS:   `,
`-Converts a segmented expression to a linear combination of shifted`,
` Heaviside functions in terms of the variable x. Such expressions can`,
` be used with the laplace and fourier functions in Maple. Caution should`,
` be used since valid segmented procedures may have piece procedures which`,
` are not defined outside their ranges. The corresponding Heaviside`,
` expressions can ``blow up`` if inappropriate values are substituted for`,
` x in such cases.`,
`   `,
`EXAMPLES:   `,
`   `,
`f := table([[1,x->exp(-x)],[[-1,0],[0,1]]]):`,
`convert(f,segtoH,t);`,
`   Heaviside(t + 1) + Heaviside(t) (exp(- t) - 1) - Heaviside(t - 1) exp(- t)`,
`   `,
`f := table([[1,x->exp(-x)],[[-infinity,0],[0,infinity]]]):`,
`convert(f,segtoH,t);`,
`                        1 + Heaviside(t) (exp(- t) - 1)`,
`   `,
`SEE ALSO: convert/Htoseg, evalsg, Heaviside, laplace, fourier`
):
`help/text/convert/Htoseg` := TEXT(
`FUNCTION: convert/Htoseg - convert Heaviside formula to segmented procedure`,
`   `,
`CALLING SEQUENCE:`,
`   convert(f, Htoseg, x)`,
`   `,
`PARAMETERS:`,
`   `,
`   f - expression involving Heaviside functions`,
`   Htoseg - indicates conversion from Heaviside to segmented`,
`   x - variable used in Heaviside expression`,
`   `,
`SYNOPSIS:   `,
`-Converts a linear combination of shifted Heaviside functions in terms`,
` of the variable x to a segmented procedure; this procedure then may`,
` be manipulated with the functions in the segmentd package.`,
`   `,
`EXAMPLES:   `,
`   `,
`f := 1-Heaviside(t) + Heaviside(t)*exp(-t):`,
`convert(f,Htoseg,t);`,
`                    table([`,
`                        1 = [1, t -> exp(- t)]`,
`                        2 = [[- infinity, 0], [0, infinity]]`,
`                    ])`,
`   `,
`SEE ALSO: convert/segtoH, evalsg, Heaviside, int, sgint`
):
`help/text/segment` := TEXT(
`FUNCTION: segment - create single-piece piecewise function in the segmented`,
`                   format, defined by an expression on an interval.`,
`  `,
`CALLING SEQUENCE:`,
`   segment(a, b, <expression>, <variable>)`,
`  `,
`PARAMETERS:`,
`a, b -  the left and right endpoints of the support interval;`,
`<expression> - the expression defining the function. It may be an`,
`  algebraic expression in a single variable and/or single variable`,
`  procedures.`,
`<variable> - an optional argument. If present it should be of Maple type`,
` name, excluding system constants: segment will automatically reject system`,
` constants. If present it is used as the variable for the segmented`,
` object which is created. If absent, then segment attempts to identify the`,
` independent variable in <expression>. If no variable is evident, because,`,
` for instance, <expression> is a procedure or a constant, then segment will`,
` use a variable which will look like x, but will be different from and`,
` independent of any x outside of segment.`,
`  `,
`  `,
`EXAMPLES:`,
`  `,
`  `,
`> segment(0,1, x^2+1);`,
`                            table([ `,
`                                           2`,
`                                1 = [x -> x  + 1]`,
`                                2 = [[0, 1]]`,
`                            ])`,
`  `,
`> lt:= segment(0, infinity, t*exp(-x*t),t);`,
`                         table([`,
`                             1 = [t -> t exp(- x t)]`,
`                             2 = [[0, infinity]]`,
`                         ])`,
`  `,
`> evalsg(lt,2);`,
`                                 2 exp(- 2 x)`,
`  `,
`  `,
`  `,
`> p := x -> x^2:`,
`> segment(-1, 1, t*p);`,
`                             table([`,
`                                              3 `,
`                                 1 = [(t -> t) ]`,
`                                 2 = [[-1, 1]]`,
`                             ])`,
`  `,
`> segment(-1, 1, p, 2.4);`,
`Error: (in segment), 2.4, is not a variable but should be`,
`  `,
`SEE ALSO: evalsg, convert/Htoseg `
):
unprotect(help):
oldhelp:= eval(help):
help := proc(who)
  local newnames;

newnames:= [`convert/Htoseg`, `convert/Htoseg/uneval`, `convert/ifproc`,
 `convert/mmsptoseg`, `convert/pwtoseg`, `convert/segtoH`,
 `convert/segtommsp`, `convert/segtopw`, evalsg, `evalsg/findtype`,
 `evalsg/hastype`, `evalsg/korder`, `evalsg/sgrefine`, `pint/sgint`,
  plotsg, `plotsg/pweval`, segment, sgint, `sgint/dirint`, sgmapf,
  sgsubs, trimsg, `type/mmsp`, `type/segmented`];
if member(who,newnames) then
  Help(who);
else
  oldhelp(who);
fi;
end:

Help := proc (u) local zz; zz := eval(cat(`help/text/`,u)); INTERFACE_HELP('
display',('topic') = u,('text') = zz) end:

protect(help):