###############################################################################
#                                                                             #
#                 Studio di funzioni razionali fratte                         #
#                 -----------------------------------                         #
#                                                                             #
# Autore: Claudio Marsan                                                      #
# Ultima modifica: 16 marzo 2003                                              #
# Versione: Maple V Release 6.02 for Windows 2000                             #
#                                                                             #
###############################################################################

studio_funzione_razionale := proc(funz_rat :: ratpoly, x, x1, x2, y1, y2)
  local load_lib, f, num, den, poli, num_poli, zeri, pn, num_pn, i, grv, t,
	as_or, grh, m_obl, q_obl, gro, f1, f2, aux_crit, crit, num_crit,
	flessi_or, k, cresc, num_cresc, aux_conc, conc, num_conc, conc2,
	num_conc2, grf;

  # Controllo input
  if x1 >= x2 then 
    ERROR(`Intervallo delle ascisse non valido!`);
  end if;
  if y1 >= y2 then 
    ERROR(`Intervallo delle ordinate non valido!`);
  end if;

  # Carica librerie ausiliarie
  if load_lib <> 1 then
    with(plots):
    load_lib := 1;
  end if;

  # Semplificazione della funzione data
  f := unapply(simplify(funz_rat(x)), x);

  # Numeratore e denominatore di f(x)
  num := numer(f(x));
  den := denom(f(x));
  lprint(`------------------------------------------------------------`);

  # Zeri del denominatore di f(x) (poli)
  poli := {fsolve(den = 0, x, fulldigits)};
  num_poli := nops(poli);
  if num_poli = 0 then
    lprint(`Il dominio della funzione e' l'insieme dei numeri reali`);
  else
    lprint(`Punti estranei al dominio della funzione: `);
    print(poli);
  end if;
  lprint(`------------------------------------------------------------`);

  # Zeri del numeratore di f(x)
  zeri := {fsolve(num = 0, x, fulldigits)};
  if nops(zeri) = 0 then
    lprint(`La funzione non ha zeri reali`);
  else
    lprint(`Zeri della funzione:`);
    print(zeri);
  end if;
  lprint(`------------------------------------------------------------`);

  # Ordinata all'origine
  if {0.0} intersect poli = {0.0} then
    lprint(`Il grafico della funzione non interseca l'asse y`);
  else
    lprint(`Ordinata all'origine:`);
    print(f(0));
  end if;
  lprint(`------------------------------------------------------------`);

  # Pari o dispari
  if eval(f(x)-f(-x)) = 0 then
    lprint(`La funzione e' pari`);
  else
    lprint(`La funzione non e' pari`);
  end if;
  if eval(f(x)+f(-x)) = 0 then
    lprint(`La funzione e' dispari`);
  else
    lprint(`La funzione non e' dispari`);
  end if;
  lprint(`------------------------------------------------------------`);

  # Segno della funzione negli intervalli determinati dai punti notevoli
  pn := sort(convert(poli union zeri, list));
  num_pn := nops(pn);
  if num_pn = 0 then
    if f(0) > 0 then
      lprint(`La funzione e' sempre positiva`);
    else
      lprint(`La funzione e' sempre negativa`);
    end if;
  else
    if f(pn[1]-1) < 0 then
      lprint(`La funzione e' negativa nell'intervallo aperto:`);
      print(-infinity..pn[1]);
    else
      lprint(`La funzione e' positiva nell'intervallo aperto:`);
      print(-infinity..pn[1]);
    end if;
    if num_pn <> 1 then
      for i from 1 to num_pn-1 do
        if f(0.5*(pn[i] + pn[i+1])) < 0 then
          lprint(`La funzione e' negativa nell'intervallo aperto:`);
	  print(pn[i]..pn[i+1]);
        else
          lprint(`La funzione e' positiva nell'intervallo aperto:`);
	  print(pn[i]..pn[i+1]);
	end if;
      end do;
    end if;
      if f(pn[num_pn]+1) < 0 then
        lprint(`La funzione e' negativa nell'intervallo aperto:`);
        print(pn[num_pn]..+infinity);
      else
        lprint(`La funzione e' positiva nell'intervallo aperto:`);
        print(pn[num_pn]..+infinity);
      end if;
  end if;
  lprint(`------------------------------------------------------------`);

  # Asintoti verticali
  if num_poli = 0 then
    lprint(`Non ci sono asintoti verticali`);
    grv := plot(0, x1..x2, color=BLACK);
  else
    lprint(`Asintoti verticali: `);
    for i from 1 to num_poli do
      print(`x` = poli[i]);
      lprint();
    end do;
    grv := seq(plot([poli[i], t, t=y1..y2], color=BLUE), i=1..num_poli);
  end if;
  lprint(`------------------------------------------------------------`);

  # Asintoti orizzontali
  as_or := limit(f(x), x = infinity);
  if {as_or} intersect {infinity, -infinity, undefined} = {} and type(as_or, realcons) then
    lprint(`Asintoto orizzontale:`);
    print(`y` = as_or);
    grh := plot(as_or, x1..x2, color=BLUE);
  else
    lprint(`Non ci sono asintoti orizzontali`);
    grh := plot(0, x1..x2, color=BLACK);
  end if;
  lprint(`------------------------------------------------------------`);

  # Asintoti obliqui
  m_obl := limit(f(x)/x, x = infinity);
  if {m_obl} intersect {0, infinity, -infinity, undefined} = {} and type(m_obl, realcons) then
    q_obl := limit(f(x) - m_obl*x, x = infinity);
    lprint(`Asintoto obliquo:`);
    print(`y` = m_obl*x + q_obl);
    gro := plot(m_obl*x + q_obl, x=x1..x2, color=BLUE);
  else
    lprint(`Non ci sono asintoti obliqui`);
    gro := plot(0, x1..x2, color=BLACK);
  end if;
  lprint(`------------------------------------------------------------`);

  # Calcolo della prima derivata di f(x);
  f1 := D(f);
  lprint(`La derivata prima della funzione:`);
  print(factor(f1(x)));
  lprint(`------------------------------------------------------------`);

  # Calcolo della seconda derivata della funzione
  f2 := D(f1);
  lprint(`La derivata seconda della funzione:`);
  print(factor(f2(x)));
  lprint();
  lprint(`------------------------------------------------------------`);

  # Punti critici
  aux_crit := {fsolve(numer(factor(f1(x))) = 0, x, fulldigits)};
  crit := sort(convert(aux_crit, list));
  num_crit := nops(crit);
  flessi_or := {};
  if num_crit = 0 then
    lprint(`Non ci sono punti critici`);
  else
    lprint(`Punti critici:`);
    print(seq([crit[i], f(crit[i])], i = 1..num_crit));
    lprint(`------------------------------------------------------------`);
    for i from 1 to num_crit do
      k := 2;
      while subs(x = crit[i], diff(f(x), x$k)) = 0 do 
         k := k + 1; 
      end do;
      if type(k, odd) then
        flessi_or := flessi_or union {crit[i]};
        lprint(`Flesso orizzontale in:`);
        print([crit[i], f(crit[i])]);
      else
        if subs(x = crit[i], diff(f(x), x$k)) > 0 then
          lprint(`Minimo locale in:`);
	  print([crit[i], f(crit[i])]);
	else
	  lprint(`Massimo locale in:`);
	  print([crit[i], f(crit[i])]);
	end if;
      end if;
    end do;
  end if;
  lprint(`------------------------------------------------------------`);

  # Intervalli in cui la funzione e' crescente o decrescente
  cresc := aux_crit union poli;
  cresc := sort(convert(cresc, list));
  num_cresc := nops(cresc);
  if num_cresc = 0 then
    if f1(0) > 0 then
      lprint(`La funzione e' sempre crescente`);
    else
      lprint(`La funzione e' sempre decrescente`);
    end if;
  else
    if f1(cresc[1]-1) < 0 then
      lprint(`La funzione e' decrescente nell'intervallo aperto:`);
      print(-infinity..cresc[1]);
    else
      lprint(`La funzione e' crescente nell'intervallo aperto:`);
      print(-infinity..cresc[1]);
    end if;
    if num_cresc <> 1	then
      for i from 1 to num_cresc-1 do
        if f1(0.5*(cresc[i] + cresc[i+1])) < 0 then
          lprint(`La funzione e' decrescente nell'intervallo aperto:`);
          print(cresc[i]..cresc[i+1]);
        else
          lprint(`La funzione e' crescente nell'intervallo aperto:`);
          print(cresc[i]..cresc[i+1]);
        end if;
      end do;
    end if;
    if f1(cresc[num_cresc]+1) < 0 then
      lprint(`La funzione e' decrescente nell'intervallo aperto:`);
      print(cresc[num_cresc]..infinity);
    else
      lprint(`La funzione e' crescente nell'intervallo aperto:`);
      print(cresc[num_cresc]..infinity);
    end if;
  end if;
  lprint(`------------------------------------------------------------`);

  # Flessi obliqui
  aux_conc := {fsolve(numer(factor(f2(x))) = 0, x, fulldigits)};
  aux_conc := aux_conc minus flessi_or;
  conc := sort(convert(aux_conc, list));
  num_conc := nops(conc);
  if num_conc = 0 then
    lprint(`Non ci sono flessi obliqui`);
  else
    lprint(`Flessi obliqui:`);
    print(seq([conc[i], f(conc[i])], i = 1..num_conc));
  end if;
  lprint(`------------------------------------------------------------`);

  # Intervalli nei quali la funzione e' convessa/concava
  conc2 := poli union aux_conc union flessi_or;
  conc2 := sort(convert(conc2, list));
  num_conc2 := nops(conc2);
  if num_conc2 = 0 then
    if f2(0) > 0 then
      lprint(`La funzione e' convessa`);
    else
      lprint(`La funzione e' concava`);
    end if;
  else
    if f2(conc2[1]-1) > 0 then
      lprint(`La funzione e' convessa nell'intervallo aperto:`);
      print(-infinity..conc2[1]);
    else
      lprint(`La funzione e' concava nell'intervallo aperto:`);
      print(-infinity..conc2[1]);
    end if;
    if num_conc2 <> 1	then
      for i from 1 to num_conc2 - 1 do
        if f2(0.5*(conc2[i] + conc2[i+1])) > 0 then
          lprint(`La funzione e' convessa nell'intervallo aperto:`);
          print(conc2[i]..conc2[i+1]);
        else
	  lprint(`La funzione e' concava nell'intervallo aperto:`);
	  print(conc2[i]..conc2[i+1]);
	end if;
      end do;
    end if;
    if f2(conc2[num_conc2]+1) > 0 then
      lprint(`La funzione e' convessa nell'intervallo aperto:`);
      print(conc2[num_conc2]..infinity);
    else
      lprint(`La funzione e' concava nell'intervallo aperto:`);
      print(conc2[num_conc2]..infinity);
    end if;
  end if;
  lprint(`------------------------------------------------------------`);

  # Visualizzazione grafica
  grf := plot(f(x), x=x1..x2, y=y1..y2, discont=true, color=RED);
  display([grh, gro, grf, grv]);

end: