A signature inequality? Given two real symmetric matrices $A$ and $B$ of common square size $n$ with no strictly negative eigenvalues, can the symmetric matrix $AB+BA$ have strictly more than $n/2$ eigenvalues which are strictly negative?
The answer to this question is yes, thanks to Junkie. My random examples did not hit a counterexample since I did them in dimension considerably greater than $3$ (typically $8$
or something similar) where it seems quite hard to find counterexamples by taking generic matrices.
This suggests however a series of new questions (which I can unfortunately no longer accept since I
gave already credit to Junkie for a correct answer):
What is the maximal number of strictly negative eigenvalues of $AB+BA$ if $A$ and $B$ are 
definite positive symmetric matrices of common size $n\times n$?
This number is at least roughly $3n/4$ by Junkie's examples (put them along the diagonal). Can it be considerably higher?
I have for example currently no example with $4$ strictly negative eigenvalues for $n=5$.
($3$ strictly negative eigenvalues in dimension $n=5$ are easy to achieve by combining
Junkie's example with an example in dimension $2$ yielding signature $(1,1)$.)
It seems that there is always at least one non-negative eigenvalue (this is obvious if 
$A$ and $B$ have only positive coefficients by Perron-Frobenius and it is probably not very hard in the general case).
 A: Example found randomly:
EDIT: Make one with positive coefficients:
$A=\pmatrix{1&2&3\cr2&5&6\cr3&6&10}$ and
$B=\pmatrix{1&1&2\cr1&2&6\cr2&6&21}$.
EDIT: Here's how, with Magma
I get about a 25% probability with this Magma code:
R := RealField(30);
function FindCounterExample()
  S := RandomSLnZ(3,5,5); A := S*Transpose(S);                                       
  S := RandomSLnZ(3,5,5); B := S*Transpose(S);                                       
  ROOTS := Roots(CharacteristicPolynomial(A*B+B*A),R);                               
  ROOTS := [r[1] : r in ROOTS | r[1] ge 0];                                            
  if #ROOTS eq 1 then A; B; end if;                                                
  return #ROOTS;
  end function;                                                     

I get about a 25% probability of 1 positive eigenvalue, 75% of 2, and 0.15% of 3.
OUTPUT := [FindCounterExample() : i in [1..100000]];                                    
SequenceToMultiset(OUTPUT); // {* 1^^25563, 2^^74296, 3^^141 *}

EDIT: I think this can be described as saying that there is about a 75% chance of the determinant of $AB+BA$ being negative, and when it in the 25% positive case, the chance is not too great that all the eigenvalues are positive. It can also depends on what RandomSLnZ is doing. The split might only be close to 75-25 and not exact.
EDIT: Yes when I did it with RandomSLnZ(3,3,3) I get a split of about $114+844+42$, so the 75-25 is meaningless.
A: This review seems to imply that any symmetric real matrix $C$ with positive trace is the Jordan product $(AB+BA)/2$ of two positive definite real matrices $A,B$. If so, then the maximum number of negative eigenvalues of $(AB+BA)/2$ for $n\times n$ symmetric positive definite $A,B$ is $n-1$ (it cannot be $n$ because of the positive definiteness (hence positive trace) of $A^{1/2}BA^{1/2}$, which is conjugate to $AB$).
A: Overnight search (about 25 million random 5x5) found:
[ 835  791 -119   -1  981]
[ 791  755 -113    0  931]
[-119 -113   17    0 -140]
[  -1    0    0    1    2]
[ 981  931 -140    2 1166]

[   5   76   -4    2  -14]
[  76 2849    0   75 -531]
[  -4    0   17    0    0]
[   2   75    0    2  -14]
[ -14 -531    0  -14   99]

A: Here is Magma code that will construct a 10x10 example, with one eigenvalue about +10 and nine of -1. Note that Magma has little numerical linear algebra, so I do it myself.
As can be seen in the final $P$ and $Q$, the numbers become of size about $10^{40}$ here. You can round to the nearest integer after multiplying by $10^{100}$ and get an integral answer of course. It seems much harder to require that $P$ and $Q$ have positive coefficients.
RF:=RealField(100);                                                             
function ThetaRoot(M)
 h11,h12,h21,h22:=Explode(Eltseq(M));                                           
 t2:=h11^2; t1:=h22^2; t12:=(h12+h21)^2-2*h11*h22;                              
 a:=t1+t2+t12; b:=-2*t1-t12; c:=t1;                                             
 R:=Roots(Polynomial([c,b,a]),RF); if #R eq 0 then R:=[<-b/2/a,2>]; end if;     
 th:=Arccos(Sqrt(R[1][1])); return th; end function;                            

function UnitarilySimilar(M)
 t:=Trace(M); d:=Degree(Parent(M)); S:=SymmetricGroup(d);                       
 M:=M-t/d*Parent(M)!1;                                                          
 if d eq 1 or M[1][1] eq 0 then return Parent(M)!1; end if;                     
 D:=Diagonal(M);                                                                
 POS:=[i : i in [1..d] | D[i] gt 0]; NEG:=[i : i in [1..d] | D[i] lt 0];        
 p:=POS[1]; n:=NEG[1]; u:=S!1;                                                  
 if p ne 1 then u*:=S!(p,1); end if; if n ne 2 then u*:=S!(n,2); end if;        
 if p eq 2 and n eq 1 then u:=S!(1,2); end if;                                  
 T1:=PermutationMatrix(BaseRing(M),u);                                          
 M1:=T1*M*Transpose(T1);                                                        
 if Abs(M1[1][1]) gt Abs(M1[2][2]) then
  Tt:=Parent(M)!1; Tt[1][1]:=0; Tt[1][2]:=-1; Tt[2][1]:=1; Tt[2][2]:=0;         
  T1:=Tt*T1; M1:=T1*M*Transpose(T1); end if;                                    
 th:=ThetaRoot(Submatrix(M1,[1,2],[1,2]));                                      
 T2:=Parent(M)!1;                                                               
 T2[p][p]:=Cos(th); T2[p][n]:=-Sin(th);                                         
 T2[n][p]:=-T2[p][n]; T2[n][n]:=T2[p][p];                                       
 M2:=T2*M1*Transpose(T2);                                                       
 if Abs(M2[1][1]) gt 10^(-25) then T2[p][n]:=Sin(th); T2[n][p]:=-T2[p][n];      
  M2:=T2*M1*Transpose(T2); end if;                                              
 U:=UnitarilySimilar(Submatrix(M2,2,2,d-1,d-1));                                
 U:=DirectSum(<DiagonalMatrix([BaseRing(U)!1]),U>);                             
 return U*T2*T1; end function;                                                  

M:=DiagonalMatrix([RF!10,-1,-1,-1,-1,-1,-1,-1,-1,-1]); d:=Degree(Parent(M));    
U:=UnitarilySimilar(M); MT:=U*M*Transpose(U); I:=Parent(MT)!1;                  
for a in [1..d] do for b in [a+1..d] do // Ballantine gives "too diagonal"
 I[a][b]:=Random([-2^25..2^25])/2^32; end for; end for; // perturb it           
PERTURB:=I*MT*Transpose(I); B:=PERTURB; TRANS:=I*U;                             
Diagonal(PERTURB);                                                              
for a in [1..d] do for b in [a+1..d] do B[a][b]:=0; end for; end for;           
for a in [1..d] do B[a][a]:=B[a][a]/2; end for;                                 

function EpsilonKernel(M) d:=Degree(Parent(M));                                 
 S:=Vector([BaseRing(M)!0 : i in [1..d]]);                                      
 for e in [1..d] do T:=M;                                                       
  for f in Reverse([1..e-1]) do S[f]:=T[e][f]/T[f][f];                          
   T[e]-:=T[f]*S[f]; end for;                                                   
  if Abs(T[e][e]) lt 10^(-50) then S[e]:=-1; return S/Sqrt(Norm(S)); end if;    
  end for; end function;                                                        

function MagmaBrainDeadEigenvectors(M)
 R:=[r[1] : r in Roots(CharacteristicPolynomial(M))];                           
 return DiagonalMatrix(R),Matrix([EpsilonKernel(M-r*Parent(M)!1) : r in R]);    
end function;                                                                   

D,V:=MagmaBrainDeadEigenvectors(B);                                             

H:=Transpose(V)*D*V; P:=Transpose(B)*H; Q:=H^(-1);                              
P:=(P+Transpose(P))/2; Q:=(Q+Transpose(Q))/2; // cheat for symmetry             

Roots(CharacteristicPolynomial(P*Q+Q*P));

EDIT: In larger dimensions, the $Q=H^{-1}$ becomes time-consuming, perhaps for stability.
