function siam4(ndig)

mp.Digits(ndig);

sol = mp(strcat(['-3.30686864747523728007611377089851565716648236147628821750129308',...
'54963091998378882950358254880752834991861926308911364307171326954',...
'79183570527535481485323493119666055496491770987805562739806177408',...
'72901417026427239222708620733463947755147440345559705900071013093',...
'04845497900412573978770067606624031615977246112178932595362940539',...
'37533901074129148867476633116962279417362942396292658321079163381',...
'01316851063267855401990040152471171435568823442063752168543176730',...
'94579094257688823510971485048530835209863263399835579478491546301',...
'04547099679467311070004309817910593652674429408839748633882190364',...
'16544051720947665168360175353514147576012479164362997035635341769',...
'38957353505986354736421228225286888787496955098278981436285213904',...
'94115341245395385269727851257267362553188546604443742861873679538',...
'62838938586899796304777948631853716589905225050450767234347594968',...
'87804395245097967371042619199815779909795089672978863171733928079',...
'94861717751825429547302477283761015961319823829914040695514974564',...
'40511842541142406112230332341888455648573656622252085632075182302',...
'66083427054258008450368335076975967282094653382389414146131200851',...
'06299528299825023519889089340183007458722421466269942167133669278',...
'83341450327272827867357800483744332635448485548065564089408002517',...
'35139101672911979974958138274962496020119299267212446172560896708',...
'65468288638366468135129976048162880156657338839310067728263645720',...
'82538754562810788610398782908027199153649182587953477795007729568',...
'75072228092246347097796171462871788753134572655654328154498571431',...
'99587242798869150258121097516385474973713888845369120442844272495',...
'99180491213466992683866664276305657610970050905190579769822840120',...
'56339111743849434420963474067840643770565653300873268879367833545',...
'42409733851368687361556421467677329544042429770058745783723398865']));



 
[X,Y] = meshgrid(linspace(-4,4,1001),linspace(-4,4,1001));
Z     = problem4f([X(:) Y(:)] );
 
[I,J] = min(Z);
 
U0    = [X(J); Y(J)];
U0    = mp(U0);
 
options.Method = 'lbfgs';%'lbfgs';
options.progTol = mp(10^(-2*ndig-10));
options.optTol = mp(10^(-ndig-10));
options.MaxIter = 10000;
opts = [];
[U,fx] = minFunc(@problem4f,U0,options,opts);
 
fprintf(['Optimal value: %.',num2str(ndig),'f\n'],fx);
fprintf(['Known sol    : %.',num2str(ndig),'f\n'],sol);
fprintf('Error to known solution: %.1e\n',abs(sol-fx)); 

res= fx;

function [v,g] = problem4f(vec,opts);
 
if prod(size(vec))==2
  x = vec(1);
  y = vec(2);
else
  x = vec(:,1);
  y = vec(:,2);
end
 
v = exp(sin(50*x))+sin(60*exp(y))+sin(70*sin(x))+sin(sin(80*y))-sin(10*(x+y))+(x.^2+y.^2)/4;
 
 
 
if nargout>1
g = zeros(size(vec));
g(1) = exp(sin(50*x))*50.*cos(50*x)+cos(70*sin(x)).*70.*cos(x)-cos(10*(x+y))*10+x/2;
 
g(2) = cos(60*exp(y))*60.*exp(y)+cos(sin(80*y)).*cos(80*y)*80-cos(10*(x+y))*10+y/2;
 
g = g(:);
end
