95% 的人解不出这道题《史上最贱的数学题》能否用数学软件计算出来？

       lI=longInteger(); a=lI.new('154476802108746166441951315019919837485664325669565431700026634898253202035277999'); b=lI.new('36875131794129999827197811565225474825492979968971970996283137471637224634055579'); c=lI.new('4373612677928697257861252602371390152816537558161613618621437993378423467772036'); %a*(a+b)*(a+c)+b*(b+a)*(b+c)+c*(c+a)*(c+b)=4*(a+b)*(a+c)*(b+c) left=lI.add3(lI.mult3(a,lI.add(a,b),lI.add(a,c)),lI.mult3(b,lI.add(b,a),lI.add(b,c)),lI.mult3(c,lI.add(c,a),lI.add(c,b))); right=lI.mult4(lI.new('4'),lI.add(a,b),lI.add(a,c),lI.add(b,c)); lI.minus(left,right)     

       clear lI=longInteger(); lF=longFrac(lI); lM=longMatrix(lF); lP=longPoly(lF); a=lP.new({[1 1],[1 1]},1,1); b=lP.new({[1 1],[1 1]},2,1); c=lP.new({[1 1],[1 1]},3,1); x1=lP.mult3(a,lP.add(a,b),lP.add(a,c)); x2=lP.mult3(b,lP.add(b,a),lP.add(b,c)); x3=lP.mult3(c,lP.add(c,a),lP.add(c,b)); right=lP.mult4(lP.new({lI.new('4'),lI.new('1')},0,0),lP.add(a,b),lP.add(b,c),lP.add(c,a)); left=lP.add3(x1,x2,x3); equ0=lP.minus(left,right); equ0_string=lP.display(equ0)     

       P={[-1 1],[1 1];[1 1],[1 1];0,[1 1]}; lineZ=lP.new({0,[1 1]},0,0); for i=1:3     diff1=lP.diff(equ0,i);     diff1=lP.subtitConst(diff1,[1 2 3],P);     lineZ=lP.add(lineZ,lP.mult(diff1,lP.minus(lP.new({[1 1],[1 1]},i,1),lP.new(P(i,:),0,0)))); end lineZ_string=lP.display(lineZ)     

       Q={[1 1],[1 1];[1 1],[1 1];[1 2 1],[1 1]}; P3={[1 1],[1 1];0,[1,1];0,[1 1]}; M1=lM.zeros(3,3); M1.values([1 4 7],:)=Q; M1.values([2 5 8],:)=P; M1.values([3 6 9],:)=P3; M1inv=lM.inv(M1); d2abc=lP.add3(lP.new(M1inv.values(1,:),1,1),lP.new(M1inv.values(2,:),2,1),lP.new(M1inv.values(3,:),3,1)); e2abc=lP.add3(lP.new(M1inv.values(4,:),1,1),lP.new(M1inv.values(5,:),2,1),lP.new(M1inv.values(6,:),3,1)); f2abc=lP.add3(lP.new(M1inv.values(7,:),1,1),lP.new(M1inv.values(8,:),2,1),lP.new(M1inv.values(9,:),3,1)); a2def=lP.add3(lP.new(M1.values(1,:),4,1),lP.new(M1.values(2,:),5,1),lP.new(M1.values(3,:),6,1)); b2def=lP.add3(lP.new(M1.values(4,:),4,1),lP.new(M1.values(5,:),5,1),lP.new(M1.values(6,:),6,1)); c2def=lP.add3(lP.new(M1.values(7,:),4,1),lP.new(M1.values(8,:),5,1),lP.new(M1.values(9,:),6,1)); equ1=equ0; equ1=lP.subtitPoly(equ1,1,a2def); equ1=lP.subtitPoly(equ1,2,b2def); equ1=lP.subtitPoly(equ1,3,c2def); idx1=sum(equ1.power(:,4:6)==[3 0 0],2)==3; %set coeff of d^3 to 1 frac1=equ1.coeff(idx1,:); equ1=lP.mult(equ1,lP.new(frac1(1,[2 1]),0,0)); equ1_string=lP.display(equ1)     

       idx1=sum(equ1.power(:,4:6)==[0 2 1],2)==3; %make coeff of e^2*f to 1 f2g=lP.new(equ1.coeff(idx1,[2 1]),7,1); f2g=lP.neg(f2g); g2f=[]; g2f.coeff=lF.validate(f2g.coeff(1,[2 1])); g2f.power=[0 0 0 0 0 1 0]; g2abc=lP.subtitPoly(g2f,6,f2abc); equ2=lP.subtitPoly(equ1,6,f2g); equ2_string=lP.display(equ2)  idx1=sum(equ2.power(:,[5 7])==[1 2],2)==2; %take out e*g term frac1=lF.mult(equ2.coeff(idx1,:),{[1 1],[1 2]}); e2h=[]; e2h.coeff=[{[1 1],[1 1]};frac1]; e2h.power=[0 0 0 0 0 0 0 1;0 0 0 0 0 0 1 0]; h2e=[]; h2e.coeff=[{[1 1],[1 1]};lF.neg(frac1)]; h2e.power=[0 0 0 0 1 0 0 0;0 0 0 0 0 0 1 0]; equ3=lP.subtitPoly(equ2,5,e2h); equ3_string=lP.display(equ3) h2abc=lP.subtitPoly(h2e,5,e2abc); h2abc=lP.subtitPoly(h2abc,7,g2abc); d2abc_string=lP.display(d2abc) g2abc_string=lP.display(g2abc) h2abc_string=lP.display(h2abc)     

       equ4=equ3; equ4=lP.subtitPoly(equ4,4,d2abc); equ4=lP.subtitPoly(equ4,7,g2abc); equ4=lP.subtitPoly(equ4,8,h2abc); equ4_string=lP.display(equ4)     

       %save to longInteger.m function lI=longInteger() % longInteger x=x1*(x2+x3*10+x4*100+x5*1000+...) % x=0 means 0     function y=neg(x1)         y=x1;         y(1)=-y(1);     end     function y=validate(x1)         flag1=1;         y=x1;         y(1)=sign(y(1));         n1=numel(y);         while flag1             if n1==1                 y=0;                 flag1=0;             elseif y(n1)==0                 n1=n1-1;             elseif y(n1)<0                 y=-y;             elseif y(n1)>9                 y(n1+1)=floor(y(n1)/10);                 y(n1)=mod(y(n1),10);                 n1=n1+1;             else                 idx1=find(y(2:n1)>9 | y(2:n1)<0)+1;                 inc1=floor(y(idx1)/10);                 y(idx1+1)=y(idx1+1)+inc1;                 y(idx1)=y(idx1)-inc1*10;                 if isempty(idx1)                     flag1=0;                 end             end         end         y=y(1:n1);     end     function y=new(str1)         sign1=1;         if str1(1)=='+'             str1=str1(2:end);         elseif str1(1)=='-'             str1=str1(2:end);             sign1=-1;         end         n1=numel(str1)+1;         y=zeros(1,n1);         y(1)=sign1;         for i=2:n1             y(i)=str1(n1-i+1)-'0';         end         y=validate(y);             end     function y=display(x)         if sum(abs(x))==0             y='0';         else             y=[' ','0'+x(end:-1:2)];             if x(1)==1                 y(1)='+';             else                 y(1)='-';             end            end     end     function y=add(x1,x2)%x=x1*(x2+x3*10+x4*100+x5*1000+...)         n1=numel(x1);         n2=numel(x2);         if n1>=n2             y=x1;             y(2:n2)=y(2:n2)+sign(x1(1)*x2(1))*x2(2:n2);         else             y=x2;             y(2:n1)=y(2:n1)+sign(x1(1)*x2(1))*x1(2:n1);         end         y=validate(y);     end     function y=add3(x1,x2,x3)         y=add(add(x1,x2),x3);     end     function y=add4(x1,x2,x3,x4)         y=add(add3(x1,x2,x3),x4);     end     function y=add5(x1,x2,x3,x4,x5)         y=add(add4(x1,x2,x3,x4),x5);     end     function y=minus(x1,x2)         y=add(x1,neg(x2));     end     function y=mult(x1,x2)         n1=numel(x1);         n2=numel(x2);         y=zeros(1,n1+n2);         y(1)=sign(x1(1)*x2(1));         for i=2:n1             y(i+(2:n2)-2)=y(i+(2:n2)-2)+x1(i)*x2(2:n2);         end         y=validate(y);     end     function y=mult3(x1,x2,x3)         y=mult(mult(x1,x2),x3);     end     function y=mult4(x1,x2,x3,x4)         y=mult(mult3(x1,x2,x3),x4);     end     function y=mult5(x1,x2,x3,x4,x5)         y=mult(mult4(x1,x2,x3,x4),x5);     end     function y=compare(x1,x2) %x1>x2->1, x1<x2->-1, x1==x2->0         x11=validate(x1);         x22=validate(x2);         if x11(1)>x22(1)             y=1;         elseif x11(1)<x22(1)             y=-1;         else             n1=numel(x11);             n2=numel(x22);             y=0;             if n1>n2                 y=1;             elseif n1<n2                 y=-1;             else                 while n1>1 && y==0                     if x11(n1)>x22(n1)                         y=1;                     elseif x11(n1)<x22(n1)                         y=-1;                     else                         n1=n1-1;                     end                 end             end             y=y*x11(1);         end     end     function [q,r]=div(x1,x2)         n1=numel(x1);         n2=numel(x2);         if x2==0             q=NaN;             r=NaN;         else             q=zeros(1,n1);             q(1)=sign(x1(1)*x2(1));             r=x1;             n3=n1;             while n3>1                 if compare([1 r(n3:n1)],[1 x2(2:n2)])>=0                     q(n3)=q(n3)+1;                     r(n3+(0:n2-2))=r(n3+(0:n2-2))-x2(2:n2);                 else                     n3=n3-1;                 end             end             q=validate(q);             r=validate(r);         end     end     function y=gcd(x1,x2)         m=x1;         n=x2;         m(1)=1;         n(1)=1;         r=[1 1];         while numel(r)>1 %r~=0             [~,r]=div(m,n);             m=n;             n=r;         end         y=m;     end lI=[]; lI.validate=@validate; lI.new=@new; lI.display=@display; lI.add=@add; lI.add3=@add3; lI.add4=@add4; lI.add5=@add5; lI.minus=@minus; lI.neg=@neg; lI.mult=@mult; lI.mult3=@mult3; lI.mult4=@mult4; lI.mult5=@mult5; lI.compare=@compare; lI.div=@div; lI.gcd=@gcd; end     

       %save to longFrac.m function lF=longFrac(lI)     function y=validate(F1)         y=F1;         if y{2}(1)==0             y{1}=NaN;         else             if y{2}(1)==-1                 y{1}=lI.neg(y{1});                 y{2}=lI.neg(y{2});             end             I1=lI.gcd(y{1},y{2});             y{1}=lI.div(y{1},I1);             y{2}=lI.div(y{2},I1);         end      end     function y=new(I1,I2)         y=cell(1,2);         y{1}=I1;         y{2}=I2;         y=validate(y);     end     function y=neg(F1)         y=F1;         y{1}=lI.neg(y{1});     end     function y=add(F1,F2)         y=cell(1,2);         y{1}=lI.add(lI.mult(F1{1},F2{2}),lI.mult(F1{2},F2{1}));         y{2}=lI.mult(F1{2},F2{2});         y=validate(y);     end     function y=add3(F1,F2,F3)         y=add(add(F1,F2),F3);     end     function y=add4(F1,F2,F3,F4)         y=add(add3(F1,F2,F3),F4);     end     function y=add5(F1,F2,F3,F4,F5)         y=add(add4(F1,F2,F3,F4),F5);     end     function y=minus(F1,F2)         y=add(F1,neg(F2));     end     function y=mult(F1,F2)         y=cell(1,2);         y{1}=lI.mult(F1{1},F2{1});         y{2}=lI.mult(F1{2},F2{2});         y=validate(y);     end     function y=mult3(F1,F2,F3)         y=mult(mult(F1,F2),F3);     end     function y=mult4(F1,F2,F3,F4)         y=mult(mult3(F1,F2,F3),F4);     end     function y=mult5(F1,F2,F3,F4,F5)         y=mult(mult4(F1,F2,F3,F4),F5);     end     function y=div(F1,F2)         y=cell(1,2);         y{1}=lI.mult(F1{1},F2{2});         y{2}=lI.mult(F1{2},F2{1});         y=validate(y);     end     function y=display(F1)         if numel(F1{2})==2 && F1{2}(2)==1             y=lI.display(F1{1});         else             y=lI.display(F1{2});             y=[lI.display(F1{1}) '/' y(2:end)];         end     end lF=[]; lF.validate=@validate; lF.new=@new; lF.neg=@neg; lF.add=@add; lF.add3=@add3; lF.add4=@add4; lF.add5=@add5; lF.minus=@minus; lF.mult=@mult; lF.mult3=@mult3; lF.mult4=@mult4; lF.mult5=@mult5; lF.div=@div; lF.display=@display; end     

       %save to longMatrix.m function lM=longMatrix(lF)     function y=zeros(m,n)         y=[];         y.m=m;         y.n=n;         y.values=cell(m*n,2);         for i=1:m*n             y.values{i,1}=0;             y.values{i,2}=[1 1];         end             end     function y=ones(m,n)         y=[];         y.m=m;         y.n=n;         y.values=cell(m*n,2);         for i=1:m*n             y.values{i,1}=[1 1];             y.values{i,2}=[1 1];         end             end     function y=neg(M1)         y=M1;         for i=1:y.m*y.n             y.values(i,:)=lF.neg(y.values(i,:));         end     end     function y=validate(M1)         y=M1;         for i=1:M1.m*M1.n             y.values(i,:)=lF.validate(M1.values(i,:));         end     end     function y=inv(M1)         if M1.m~=M1.n             y=NaN;         else             m=M1.m;             values=cell(m*m*2,2);             values(1:m*m,:)=M1.values;             for i=1:m*m                 values{m*m+i,1}=0;                 values{m*m+i,2}=[1 1];             end             for i=1:m                 values{m*m+(i-1)*m+i,1}=[1 1];             end             for i=1:m                 flag1=0;                 idx1=i;                 while flag1==0                     if numel(values{(i-1)*m+idx1,1})>1                         flag1=1;                     elseif idx1==m                         flag1=2;                     else                         idx1=idx1+1;                     end                 end                 if flag1==1                     values([(0:m-1)*m+i (0:m-1)*m+idx1],:)=values([(0:m-1)*m+idx1 (0:m-1)*m+i],:);                     for k=i+1:m*2                         values((k-1)*m+i,:)=lF.div(values((k-1)*m+i,:),values((i-1)*m+i,:));                     end                     values{(i-1)*m+i,1}=[1 1];                     values{(i-1)*m+i,2}=[1 1];                     for j=[1:i-1 i+1:m]                         for k=i+1:m*2                             values((k-1)*m+j,:)=lF.minus(values((k-1)*m+j,:),lF.mult(values((i-1)*m+j,:),values((k-1)*m+i,:)));                         end                         values{(i-1)*m+j,1}=0;                         values{(i-1)*m+j,2}=[1 1];                     end                 else                     break                 end             end             y=[];             y.m=m;             y.n=m;             y.values=values(m*m+(1:m*m),:);             y=validate(y);         end      end  lM=[]; lM.zeros=@zeros; lM.ones=@ones; lM.neg=@neg; lM.inv=@inv; end     

       function lP=longPoly(lF)     function y=validate(x)         maxP=max(max(x.power));         [m,n]=size(x.power);         weight1=zeros(m,1);         for i=1:n             weight1=weight1+x.power(:,i)*(maxP+1)^i;         end         [~,idx1]=sort(weight1);         y=[];         y.coeff=cell(m,2);         y.power=zeros(m,n);         i=1;         power0=zeros(1,n)-1;         coeff0={0,[1 1]};         for j=1:m             if sum(power0==x.power(idx1(j),:))==n                 coeff0=lF.add(coeff0,x.coeff(idx1(j),:));                 y.coeff(i,:)=coeff0;             else                 if numel(coeff0{1})>1                     i=i+1;                 end                 power0=x.power(idx1(j),:);                 coeff0=x.coeff(idx1(j),:);                 y.power(i,:)=power0;                 y.coeff(i,:)=coeff0;             end         end         if i==1 && numel(y.coeff{1})==1             y.power=0;             y.coeff={0,[1 1]};         elseif numel(coeff0{1})==1             y.coeff=y.coeff(1:i-1,:);             y.power=y.power(1:i-1,:);         else             y.coeff=y.coeff(1:i,:);             y.power=y.power(1:i,:);         end     end     function y=new(coeff,nVar,power)         y=[];         y.coeff=coeff;         y.power=zeros(1,max(nVar,1));         y.power(1,max(nVar,1))=power;     end     function y=display(x)         if sum(sum(abs(x.power)))==0             y='0';         else             [m,n]=size(x.power);             maxL1=0;             maxL2=0;             coeff=cell(m,1);             power=cell(m,1);             for i=1:m                 coeff{i}=lF.display(x.coeff(i,:));                 if numel(coeff{i})>maxL1                     maxL1=numel(coeff{i});                 end                 power{i}='';                 flag1=0;                 for j=1:n                     if x.power(i,j)>1                         if flag1                             power{i}=[power{i} '*'];                         end                                                     flag1=1;                         power{i}=[power{i} 'a'+j-1 '^' num2str(x.power(i,j))];                     elseif x.power(i,j)>0                         if flag1                             power{i}=[power{i} '*'];                         end                                                     flag1=1;                         power{i}=[power{i} 'a'+j-1];                                           end                 end                 if numel(power{i})>maxL2                     maxL2=numel(power{i});                 end             end             y=zeros(m,maxL1+maxL2)+' ';             for i=1:m                 if numel(power{i})==0                     y(i,1:numel(coeff{i}))=coeff{i};                 elseif numel(coeff{i})==2 && coeff{i}(2)=='1'                     y(i,1)=coeff{i}(1);                     y(i,2:numel(power{i})+1)=power{i};                 else                     y(i,1)=coeff{i}(1);                     y(i,2:numel(power{i})+1)=power{i};                     y(i,numel(power{i})+2)='*';                     y(i,numel(power{i})+1+(2:numel(coeff{i})))=coeff{i}(2:end);                 end             end          end         y=char(y);     end     function y=neg(x)         y=x;         [m,~]=size(x.power);         for i=1:m             y.coeff(i,:)=lF.neg(y.coeff(i,:));         end     end     function y=add(x1,x2)         y=[];         y.coeff=[x1.coeff;x2.coeff];         n1=size(x1.power,2);         n2=size(x2.power,2);         if n1>=n2             y.power=[x1.power;x2.power x2.power(:,1)*zeros(1,n1-n2)];         else             y.power=[x1.power x1.power(:,1)*zeros(1,n2-n1);x2.power];         end         y=validate(y);     end     function y=add3(x1,x2,x3)         y=add(add(x1,x2),x3);     end     function y=add4(x1,x2,x3,x4)         y=add(add3(x1,x2,x3),x4);     end     function y=add5(x1,x2,x3,x5)         y=add(add4(x1,x2,x3,x4),x5);     end     function y=minus(x1,x2)         y=add(x1,neg(x2));     end     function y=mult(x1,x2)         y=[];         [m1,n1]=size(x1.power);         [m2,n2]=size(x2.power);         m=m1*m2;         n=max(n1,n2);         y.coeff=cell(m,2);         y.power=zeros(m,n);         for i=1:m1             for j=1:m2                 idx1=(i-1)*m2+j;                 y.coeff(idx1,:)=lF.mult(x1.coeff(i,:),x2.coeff(j,:));                 y.power(idx1,1:n1)=x1.power(i,:);                 y.power(idx1,1:n2)=y.power(idx1,1:n2)+x2.power(j,:);             end         end         y=validate(y);     end     function y=mult3(x1,x2,x3)         y=mult(mult(x1,x2),x3);     end     function y=mult4(x1,x2,x3,x4)         y=mult(mult3(x1,x2,x3),x4);     end     function y=mult5(x1,x2,x3,x4,x5)         y=mult(mult4(x1,x2,x3,x4),x5);     end     function y=subset(x1,idx1)         y=[];         y.coeff=x1.coeff(idx1,:);         y.power=x1.power(idx1,:);     end     function y=diff(x1,n)         y=x1;         [m,~]=size(x1.power);         for i=1:m             if y.power(i,n)>0                 y.coeff(i,:)=lF.mult({[1 y.power(i,n)],[1 1]},y.coeff(i,:));                 y.power(i,n)=y.power(i,n)-1;             else                 y.coeff(i,:)={0,[1 1]};                 y.power(i,:)=0;             end         end          y=validate(y);     end     function y=subtitConst(x1,varNo,values)         y=x1;         [m,~]=size(x1.power);         for i=1:m             for j=1:numel(varNo)                 for k=1:y.power(i,varNo(j))                     y.coeff(i,:)=lF.mult(y.coeff(i,:),values(j,:));                 end                 y.power(i,varNo(j))=0;             end         end         y=validate(y);     end     function y=subtitPoly(x1,varNo,poly1)         flag1=1;         if size(poly1.power,2)>=varNo             if sum(poly1.power(:,varNo)~=0)>1                 flag1=0;                 y=NaN;             end         end         if flag1             idx1=find(x1.power(:,varNo)==0);             y=subset(x1,idx1);             idx1=find(x1.power(:,varNo)~=0);             for i=1:numel(idx1)                 x2=subset(x1,idx1(i));                 x2.power(1,varNo)=0;                 for j=1:x1.power(idx1(i),varNo)                     x2=mult(x2,poly1);                 end                 y=add(y,x2);             end         end         y=validate(y);     end     function [q,r]=div(x1,x2,varOrder) %x2 highest order only constant coefficient          [m1,n1]=size(x1.power);         [m2,n2]=size(x2.power);         n=max(n1,n2);          if ~exist('varOrder','var')             varOrder=(n:-1:1).^2;%:-1:1;         else             varOrder=varOrder.^2;         end         r=x1;         x22=x2;         if n1<n             r.power=[r.power zeros(m1,n-n1)];         end         if n2<n             x22.power=[x22.power zeros(m2,n-n2)];         end         weight2=sum(x22.power,2)*(n+1)^3+x22.power*varOrder';         [~,idx2]=sort(weight2,'descend');         x22.power=x22.power(idx2,:);         x22.coeff=x22.coeff(idx2,:);         q=lP.new({0,[1 1]},0,0); %y=new(coeff,nVar,power)         if m2==1 && numel(x22.coeff{1,1})==1             q=NaN;             r=NaN;         else             flag1=1;              while flag1                 power=r.power-ones(m1,1)*x22.power(1,:);                 weight1=sum(power,2)*(n+1)^3+power*varOrder';                 [~,idx1]=max(weight1-(sum(power>=0,2)<n)*max(abs(weight1))*1000);                 power=power(idx1,:);                 if sum(power>=0)==n                      coeff=lF.div(r.coeff(idx1,:),x22.coeff(1,:));                     poly1=[];                     poly1.coeff=coeff;                     poly1.power=power;                     q=lP.add(q,poly1);                     poly2=lP.mult(x22,poly1);                     r=lP.minus(r,poly2);                     [m1,~]=size(r.power);                 else                     flag1=0;                 end             end             q=validate(q);             r=validate(r);         end     end lP.new=@new; lP.display=@display; lP.validate=@validate; lP.neg=@neg; lP.add=@add; lP.add3=@add3; lP.add4=@add4; lP.add5=@add5; lP.minus=@minus; lP.mult=@mult; lP.mult3=@mult3; lP.mult4=@mult4; lP.mult5=@mult5; lP.subset=@subset; lP.diff=@diff; lP.subtitConst=@subtitConst; lP.subtitPoly=@subtitPoly; lP.div=@div; end