# test with   gappa -Munconstrained < sin.gappa
# The proof is not complete, as it doesn't work without -Munconstrained. 
# What it means is that Gappa is unable to prove that some denominators are not null.
# It's OK for practical purposes, but it takes some more work to get a formal proof. 

@IEEEdouble = float<ieee_64,ne>;
# Convention 1: uncapitalized variables match the variables in the C code. Other variables begin with a capital letter
# Convention 2: variables beginning with "M" are mathematical ideal

# yh+yl is a double-double (call it Yhl)

yh = IEEEdouble(dummy1);
yl = IEEEdouble(dummy2);
Yhl = yh + yl; # There is also an hypothesis stating that yl<ulp(yh)

#--------------- Transcription of the C code --------------------------

s3 = IEEEdouble(-1.6666666666666665741e-01); 
s5 = IEEEdouble( 8.3333333333333332177e-03); 
s7 = IEEEdouble(-1.9841269841269841253e-04);

yh2 IEEEdouble=  yh * yh;
ts  IEEEdouble=  yh2 * (s3 + yh2*(s5 + yh2*s7));
r   IEEEdouble=  yl + yh*ts;
s             =  yh + r;   # no rounding, it is the Fast2Sum 

#-------- Mathematical definition of what we are approximating --------

My2 = My*My;
Mts = My2 * (s3 + My2*(s5 + My2*s7));
PolySinY = My + My*Mts;

Epsargred = (Yhl - My)/My;           # argument reduction error
Epsapprox = (PolySinY - SinY)/SinY;  # polynomial approximation error
Epsround = (s - PolySinY)/PolySinY;  # rounding errors in the polynomial evaluation
Epstotal = (s - SinY)/SinY;          # total error



# Layers of approximation on s
S1 = yh + (yl + IEEEdouble(yh*ts));   # remove last round
S2 = yh + (yl + yh*ts);               # remove penultimate round
S3 = (yh+yl) + (yh+yl)*ts;            # put yl back in

Eps1 = (s-S1)/S1;
Eps2 = (S1-S2)/S2;
Eps3 = (S2-S3)/S3;
Eps4 = (S3-PolySinY)/PolySinY;


yhts = IEEEdouble(yh*ts);             # just to make the hints lighter
p3 IEEEdouble=  s3 + yh2*(s5 + yh2*s7);  # idem


tsNoRound = yh2 * (s3 + yh2*(s5 + yh2*s7));

# A few definitions mostly to benefit from automatic hints.
EpstsMts = (ts-Mts)/Mts;
EpstsNoRoundMts = (tsNoRound - Mts)/Mts;

Epsy2 = (yh2-My2)/My2;
Epsy2_argred = (Yhl*Yhl-My2)/My2;
Epsy2_negl_yl = (yh*yh-Yhl*Yhl)/(Yhl*Yhl);
Epsy2_rnd = (yh2-yh*yh)/(yh*yh);


#----------------------  The theorem to prove --------------------------
{
  # Hypotheses 
    |yl / yh| <= 1b-53
 /\ |Yhl| in [1b-200, 6.29e-03] # lower bound guaranteed by Kahan-Douglas algorithm
 /\ |yh| in [1b-1000, 1] # some huge range for ensuring that yh is not zero
 /\ |Epsargred| <= 2.53e-23
 /\ |Epsapprox| <= 2.26e-24

->

#goal to prove
   Epstotal in ?# [-1b-67, 1b-67]
/\ |r/yh| <= 1
#/\ |My| in [1b-400, 6.29e-03]
}

# ---------------------- Hints ----------------------------------

$ Yhl in (0);

# First, the hints for Epsround

s~S1;
S1~S2;
S2~S3;
S3~PolySinY;

Eps4 -> #  (S3-PolySinY)/PolySinY;
#   S3/PolySinY - 1;
#   ((yh+yl) + (yh+yl)*ts) / (My + My*Mts)  - 1;
#   ((yh+yl)/My) * (1+ts)/(1+Mts)  - 1;
#   (Epsargred+1) * (1+ts)/(1+Mts)  - 1;
#    Epsargred * (1+ts)/(1+Mts)   +  1 * (1+ts)/(1+Mts)  - 1;
#    Epsargred * (1+ts)/(1+Mts)   +  (ts-Mts)/(1+Mts);
   Epsargred * (1+ts)/(1+Mts)   +  Mts*((ts-Mts)/Mts) / (1+Mts);

# Now we just need to bound ts-Mts: 
ts ~ tsNoRound;
(tsNoRound - Mts)/Mts -> 
#     yh2/My2  * (s3 + yh2*(s5 + yh2*s7)) / (s3 + My2*(s5 + My2*s7))  - 1 ;
      (1+Epsy2)  * (s3 + yh2*(s5 + yh2*s7)) / (s3 + My2*(s5 + My2*s7)) -1;
# Now we just need to express My2 in terms of yh2 and Epsy2
My2 -> yh2/(1+Epsy2);

yh ~ Yhl;
(yh - Yhl) / Yhl -> 1 / (1 + yl / yh) - 1;

Eps3 -> 
# (S2-S3)/S3 
#  S2/S3 - 1; 
#   (yh + (yl + yh*ts)) / ((yh+yl) + (yh+yl)*ts)   - 1 ;
#   ((yh+yl) + (yh+yl)*ts - yl*ts) / ((yh+yl) + (yh+yl)*ts)   - 1 ;
#    - yl*ts / ((yh+yl) + (yh+yl)*ts)   ;
#    - (yl/Yhl)  * (ts / (1+ts))   ;
     ((yh-Yhl)/Yhl)  * (ts / (1+ts))   ; # change sign to have the expression of a rounding error


Eps2 -> # (S1-S2)/S2; 
#   (yh + (yl + IEEEdouble(yh*ts))) / (yh + (yl + yh*ts))   -1 ;
#   (IEEEdouble(yh*ts) - yh*ts) / (yh + yl + yh*ts) ;
#   ((IEEEdouble(yh*ts) - yh*ts)/(yh*ts)) / ( (yh+yl)/(yh*ts) + 1 ) ;
    ts * ((IEEEdouble(yh*ts) - yh*ts)/(yh*ts)) / ( 1 + yl/yh + ts ) ;

yhts/yh -> ts*((yhts-yh*ts)/(yh*ts) + 1);


(yl+yhts)/yh -> yl/yh + yhts/yh;


# The following are hints added to remove the -Munconstrained

PolySinY -> My*(1+Mts);

# from: Epsargred = (Yhl - My)/My;
My -> Yhl/(1+Epsargred);

r / yh -> (r - (yl + yhts)) / (yl + yhts) * ((yl + yhts) / yh) + (yl + yhts) / yh;
yl -> yl / yh * yh   { yh <> 0};