function siam1(numd)

if nargin<1
  numd=20;
end

numd = numd+16;

sol = mp(strcat(['0.323367431677778761399370087952170446651046625725469661681036443',...
'43179033721067289443193037046410245138028052733121715128433295510',...
'03343518853814345933849959100496203187230243755305850326867826084',...
'81414362100828454472232543319742328956220838303605180613817022887',...
'59277962238894141948088977862776567200916363254809300626123688307',...
'09515617259626227109122796181848626492520354425805878303512878922',...
'71350617697250436342732427806269819545993593793494870337700399388',...
'78958444941938532524701051263491411668692808429670505450842283874',...
'06907420047518022017247554109230583360447095045817596946489788630',...
'16776842609173614450859675049401159693181697158702857210715805288',...
'89123103966740195477112287059791090468098755993831428538993148977',...
'23348536769454368361557868526285246345537463616787033605541989910',...
'53033730029055952049091013813195212753251996770518810078365111887',...
'72835687475009609176911865529334596642093059266995096076632168382',...
'51160999002569152976749531082493203607852959693352367456372093026',...
'71623952826823560013238254339251717524303362820709002066320838622',...
'10170860746188444225429402690637206001424734653924223735409965579',...
'36053703204816654826851632409400237764054524939112312211532967241',...
'96140044673395516138166098905000786282187631131821194253064572759',...
'61754712960298839710460678160366323608154407900171959697298305937',...
'83549197633852783311222104880117569307839989168003909313155386494']));

mp.Digits(numd);
%DD = integral(@func,mp(1e-19),mp(1),'AbsTol',mp(10^(-numd)),'RelTol',mp(10^(-numd)));
DD = quadgk(@func,mp(0),mp(1),'AbsTol',mp(10^(-numd)),'RelTol',mp(10^(-numd)),'MaxIntervalCount',10000) 
fprintf('Computed result correct (hopefully) up to %d digits:\n',numd-16);
fprintf(['%.',num2str(numd-16),'f\n'],DD);
fprintf('Error to known result: %.1e\n',abs(DD-sol));


function y = func(t)

z = t+1i*t.*(1-t);
y = real(z.^(1i./z-1).*(1+1i*(1-2*t)));

