(*
 * SharpSolver - Progetto di Programmazione e Calcolo a.a. 2018-19
 * Impl.fsi: implementazioni degli studenti
 * (C) 2018 Alvise Spano' @ Universita' Ca' Foscari di Venezia
 *)
module SharpSolver.Impl

open Absyn
open Prelude
open System
open System.Runtime.Remoting.Metadata.W3cXsd2001


let rationalize (x : float) : rational = //funzione che prende in input un float e restituisce la sua rappresentazione in tipo rational
    if float(int(x))=x then rational( int(x) ,1) // guardo se il numero non ha parte decimale e ritorno la parte intera / 1
    else   let x_splt = x.ToString().Split ',' // spezza la parte intera e la parte decimale
           rational (int(""+x_splt.[0]+x_splt.[1]),int(10. ** float(x_splt.[1].Length))) //converto in stringa la parte intera e la concateno alla parte decimale(anch'essa convertita in intero). Compongo il denominatore elevando la base 10 al numero di cifre della parte decimale
            // e poi ritorno la coppia numeratore fratto denominatore del tipo rational

let monomial_degree (m : monomial) : int = // mi restituisce il grado di un singolo monomio
     match m with // confronto m (monomio)
        Monomial(x, y) -> y //siccome y mi rappresenta il grado del monomio, restituisco y alla chiamata di monomial_degree.

let monomial_negate (m : monomial) : monomial = //Nego il monomio.
    match m with // confronto m (monomio)
        Monomial(x ,y) -> Monomial(-x, y) // Siccome x rappresenta il coefficiente del monomio, restituisco -x alla chiamata di monomial_negate
let rec polynomial_degree (p : polynomial) : int = //funzione ricorsiva che restituisce il grado di un polinomio sottoforma di intero.
        match p with
        |Polynomial([])->raise (NotImplementedException ()) // se il poliniomio dato in input è vuoto allora gestisco eccezione.
        |Polynomial([Monomial(x,y)])-> int(y) // se ho un singolo monomio allora il grado del polinomio sarà equivalente al grado del monomio: in questo caso y.
        |Polynomial(Monomial(x1,y1)::Monomial(x2,y2)::xs) -> if y2>y1 then polynomial_degree(Polynomial(Monomial(x2,y2)::xs)) else polynomial_degree(Polynomial(Monomial(x1,y1)::xs)) //Se ho un polinomio con due monomi come primi due elementi e una coda con altri monomi, allora confronto se il grado del secondo monomio è più grande del grado del primo monomio. Se sì allora assumo come massimo grado il secondo monomio e lo faccio diventare primo elemento concatenandolo alla coda. In caso contrario pongo il primo monomio come massimo grado il primo monomio e lo concateno alla coda facendolo restare il primo elemento. Così via fino alla fine del polinomio.

let polynomial_negate (p : polynomial) : polynomial = //nega un qualsiasi polinomio.
    let rec aux p= // utilizzo un'ausiliaria per comporre una lista di monomial_list
        match p with
            Polynomial([Monomial(x,y)])-> [monomial_negate (Monomial(x,y))] // caso base con un singolo monomio. nego il monomio usand la funzione monomial_negate precedente.
            |Polynomial(Monomial(x,y)::xs)->[monomial_negate (Monomial(x,y))]@(aux (Polynomial (xs)))  // nel caso in cui ho un polinomio con un primo monomio come primo elemento e una coda con altri monomi nego il primo monomio e lo tratto come fosse una lista concatenando la coda.
    in Polynomial (aux p) // richiamando poi il tipo Polynomial nella funzione ausiliaria

let normalized_polynomial_degree (np : normalized_polynomial) : int = //funzione che restituisce il grado di un polinomio normalizzato sottoforma di intero.
    let norm_to_array np =
        match np with
            NormalizedPolynomial (np) -> np // estraggo tramite pattern matching il vettore di rational
    let d =Array.toList (norm_to_array(np)) // tramite la funzione di libreria Array.toList trasformo l'array di rational in una lista di rational
    let rec rational_to_monomiallist ls c = // funzione che trasforma una lista di rational in una monomial list
        match ls with
            [x]-> [Monomial(x,c)] // se ho solo un elemento della lista ritorno un elemento monomial con coppia di rational e grado c passato a parametro
            |x::xs-> (Monomial(x,c)) :: rational_to_monomiallist (xs) (c+1) //se ho una lista rational di più elementi converto tutta lista rational in una lista di monomi con coppia di rational e grado c passato a parametro
    let rec mono_to_poli l= // funzione che trasforma una lista di  monomial list in un polinomio
        match l with
        [x]-> Polynomial([x]) // se la lista di monomi è di un solo elemento ritorno il polinomio con un solo elemento
        |x::xs ->Polynomial(x::xs) // se la lista di monomi ha più di un elemento ritorno il polinomio con una lista di più monomi
    in polynomial_degree (mono_to_poli(rational_to_monomiallist d 0)) // vado a richiamare la funzione polynomial_degree che ritorna il grado del polinomio passando a parametro il polinomio ottenuto tramite le due funzioni di appoggio

let normalize (p : polynomial) : normalized_polynomial = // funzione che permette di normalizzare un polinomio esempio: 3x+5x=0 -> 8x=0
    let rec polinomio_to_coppia lst= // trasformo il polinomio in una lista di coppie con rational e il grado del monomio
        match lst with
            Polynomial[Monomial(x,y)]-> [x,y]
            |Polynomial(Monomial(x,y)::xs)-> (x,y):: polinomio_to_coppia(Polynomial(xs))
    let array_poli = List.toArray(polinomio_to_coppia(p)) // trasformo la lista di coppie in un array
    let dim =polynomial_degree(p)+1 // dichiaro la dimensione del nuovo array che andrò ad inizializzare
    let arreo = Array.create(dim) (rational(0,1)) // creo il nuovo array di dimensione dim e inizializzo con rational a 0
    for i=0 to dim do // scorro il nuovo array dichiarato
            for j=0 to array_poli.Length-1 do // scorro l'array di coppie
                    let (a,b) =array_poli.[j] // dichiaro una coppia e li assegno l'elemento in indice j
                    if(b=i) then arreo.[i]<-arreo.[i]+a //se il grado è uguale alla posizione nel vettore creato inizialmente , sommo il valore già presente nella posizione i e sommo il coefficiente della coppia


    NormalizedPolynomial(arreo) // ritorno il vettore di rational


let derive (p : polynomial) : polynomial = //calcolo la derivata del polinomio dato in input
    let rec aux p = // funzione ausiliaria che permette il calcolo della derivata
        match p with
        Polynomial([Monomial(x,y)]) -> ([Monomial(((rational y)*x), y-1)]) //se ho un polinomio che ha un singolo monomio applico la definizione di derviata. Moltiplico l'esponente (facendo un cast di rational su y) per il coefficente noto che è x. Dopodiché sottraggo di 1 l'esponente.
        |Polynomial(Monomial(x,y)::xs)-> ((Monomial(((rational y)*x), y-1))::(aux(Polynomial(xs)))) //se ho una lista di monomi quindi un polinomio formato da un primo elemento (un monomio) e una coda che può contenere una lista di monomi, applico la definizione di derivata al primo monomio della lista e poi concateno la coda in modo da scorrere tutta la lista e applicare la  funzione a tutti i monomi.
    in Polynomial( aux p) //Ritorno la derivata del polinomio

let rec reduce (e : expr) : polynomial = //funzione che permette di distinguere un polinomio da una derivata
    match e with
    Poly e -> e // se ho un polinomio restituisco un polinomio
    |Derive e -> derive(reduce(e)) // se ho tipo derive applico la funzione derive in ricorsione e di conseguenza la funzione ritornerà il polinomio derivato

let solve0 (np : normalized_polynomial) : bool = // funzione che permette di risolvere le equazioni di grado 0
    let extract p =
        match np with
            NormalizedPolynomial(p) -> p // estraggo tramite pattern matching il vettore di rational
    let array = extract np // applico la funzione extract al polinomio normalizzato in input
    let a= array.[0] // accedo alla posizione 0 del vettore considerando 0 come il grado dell'equazione e quindi assegno ad a il coefficente del grado 0
    if(a=rational.Zero) then true else false //se il coefficiente è uguale a 0 ritorno true altrimenti false

let solve1 (np : normalized_polynomial) : rational = // funzione che permette di risolvere le equazioni di grado 1
    let extract p =
        match np with
        NormalizedPolynomial(p) -> p // estraggo tramite pattern matching il vettore di rational
    let array = extract np // applico la funzione extract al polinomio normalizzato in input
    let zero = array.[0] // accedo alla posizione 0 del vettore considerando 0 come il grado dell'equazione e quindi assegno ad a il coefficente del grado 0
    let one = array.[1] // accedo alla posizione 1 del vettore considerando 1 come il grado dell'equazione e quindi assegno ad a il coefficente del grado 1
    rational(((- zero.N)*(one.D)),((zero.D)*(one.N))) //moltiplico coefficiente di grado 1 e grado 0 ,alla fine  ritorno l'unica soluzione dell'equazione di grado 1 di tipo rational


let solve2 (np : normalized_polynomial) : (float * float option) option = // funzione che permette di risolvere le equazioni di grado 2
    let extract p =
        match np with
        NormalizedPolynomial(p) -> p // estraggo tramite pattern matching il vettore di rational
    let array = extract np // applico la funzione extract al polinomio normalizzato in input
    let zero = float(float(array.[0].N)/float(array.[0].D)) // accedo alla posizione 0 del vettore considerando 0 come il grado dell'equazione e quindi assegno a zero il coefficente del grado 0 , divido numeratore per il denominatore in modo da ottenere un numero float
    let one = float(float(array.[1].N)/float(array.[1].D)) // accedo alla posizione 1 del vettore considerando 1 come il grado dell'equazione e quindi assegno a one il coefficente del grado 1 ,  divido numeratore per il denominatore in modo da ottenere un numero float
    let two = float(float(array.[2].N)/float(array.[2].D)) // accedo alla posizione 2 del vettore considerando 2 come il grado dell'equazione e quindi assegno a two il coefficente del grado 2 ,  divido numeratore per il denominatore in modo da ottenere un numero float
    // il delta in quanto in matematica viene calcolato in questo modo delta = b^2 -4ac , ho deciso il procedimento in sottoprocedimenti sottostanti
    let delta1 = one*one // calcola il b^2 del delta
    let delta2 = 4.0 *zero*two // calcolo il 4ac del delta
    let delta = delta1 - delta2 // infine ottengono il delta effettuando b^2-4ac
    if(delta>0.) then // se il delta è maggiore di 0 andrò a ritornare 2 soluzioni
            let num1 = -one + sqrt(delta) // calcolo il -b+ la radice quadrata del delta (1a radice dell'equazione)
            let den1 = 2. *two //calcolo il denominatore 2*a
            let x1 =  num1/ den1 // divido il numeratore con il denominatore quindi (-b + radice del delta) /2a ed ottengo la prima radice dell'equazione
            let num2 = -one - sqrt(delta) // calcolo il -b-  la radice quadrata del delta (2a radice dell'equazione)
            let den2 = 2. *two //calcolo il denominatore 2*a
            let x2 = float(num2/den2) // divido il numeratore con il denominatore quindi (-b - radice del delta) /2a ed ottengo la seconda radice dell'equazione
            Some(x1,Some x2) //ritorno le due soluzioni
    elif(delta=0.) then // se il delta è uguale a 0 andrò a ritornare 2 soluzioni uguali.
            let num1 = -one // calcolo il -b+ la radice quadrata del delta (1a radice dell'equazione) , il delta è uguale a 0 e quindi non vado ne a sommare e ne a togliere la radice quadrata del delta
            let den1 = 2. *two //calcolo il denominatore 2*a
            let x1 =  float(num1/ den1) // divido il numeratore con il denominatore (-b/2a) ed ottengo una radice dell'equazione
            Some(x1,None) // ritorno l'unica soluzione
    else None // se il delta è minore di 0 non ho soluzioni

let solve3 (np: normalized_polynomial) : (float * float option * float option) option =
    let extract p =
        match np with
        NormalizedPolynomial(p) -> p
    let array = extract np // applico la funzione extract al polinomio normalizzato in input
    let d = float(float(array.[0].N)/float(array.[0].D)) // accedo alla posizione 0 del vettore considerando 0 come il grado dell'equazione e quindi assegno a zero il coefficente del grado 0 , divido numeratore per il denominatore in modo da ottenere un numero float
    let c =  float( float(array.[1].N)/ float(array.[1].D)) // accedo alla posizione 1 del vettore considerando 1 come il grado dell'equazione e quindi assegno a one il coefficente del grado 1 ,  divido numeratore per il denominatore in modo da ottenere un numero float
    let b =  float( float(array.[2].N)/ float(array.[2].D)) // accedo alla posizione 2 del vettore considerando 2 come il grado dell'equazione e quindi assegno a two il coefficente del grado 2 ,  divido numeratore per il denominatore in modo da ottenere un numero float
    let a =  float(float(array.[3].N)/ float(array.[3].D)) // accedo alla posizione 3 del vettore considerando 3 come grado dell'equazione e quindi assegno a three il coefficiente del grado 3, dividendo numeratore per denoinatore in modo da ottenere un numero float.
    let f = (((3.*c)/a)- ((b*b)/(a*a)))/3.
    let g = (((2.*(b*b*b))/(a*a*a)) - ((9.*b*c)/(a*a))+((27.*d)/a))/27.
    let h = ((g*g)/4.) + ((f*f*f)/27.0)
    if (f+g+h)=0. then   // all roots are real and equal
       let x1 = Math.Pow(float(d/a),0.33333333333333333333333333333333)*(-1.)
       let x2 = Math.Pow(float(d/a),0.33333333333333333333333333333333)*(-1.)
       let x3 = Math.Pow(float(d/a),0.33333333333333333333333333333333)*(-1.)
       Some(x1,Some x2,Some x3)
    elif(h<=0.) then       // all 3 roots are real
       let i = Math.Sqrt(((g*g)/4.) - h)
       let j = Math.Pow(i,0.33333333333333333333333333333333)
       let k = Math.Acos(-(g/(2.*i)))
       let l = j * -1.
       let m = Math.Cos(k/3.)
       let n = Math.Sqrt(3.) *Math.Sin(k/3.)
       let p = (b/(3.*a))* -1.
       let x1 = (2.*j)*Math.Cos(k/3.)-(b/(3.*a))
       let x2 = l * (m+n) + p
       let x3 = l * (m-n) + p
       Some(x1,Some x2,Some x3)
    elif(h>0.) then
       let r = (-1.*(g/2.0)) + Math.Sqrt(h)
       if r<0. then let s = Math.Pow(-r,0.33333333333333333333333333333333)
                    let t =  (-1.*(g/2.0)) - Math.Sqrt(h)
                    if t<0. then let u = Math.Pow(-t,0.33333333333333333333333333333333)
                                 let x1 = (-s-u) - (b/(3.0*a))
                                 Some(x1,None,None)
                    else
                        let u = Math.Pow(t,0.33333333333333333333333333333333)
                        let x1 = (-s+u) - (b/(3.0*a))
                        Some(x1,None,None)
       else let s = Math.Pow(r,0.33333333333333333333333333333333)
            let t =  (-1.*(g/2.0)) - Math.Sqrt(h)
            if t<0. then let u = Math.Pow(-t,0.33333333333333333333333333333333)
                         let x1 = (s-u) - (b/(3.0*a))
                         Some(x1,None,None)
                    else
                        let u = Math.Pow(t,0.33333333333333333333333333333333)
                        let x1 = (s+u) - (b/(3.0*a))
                        Some(x1,None,None)

    else None