testw.pas

delphi矩阵运算、回归分析等数学计算

   - 1.000065950300622136103491886720203687457,
   - 1.00006995095046930174807034225078340539,
   - 1.000073734868993836022551364404248563489,
   - 1.007391489031309264813153180819941418531,
   - 1.010463846806564696239430620915659099361,
   - 1.012825626038880105597738761474228363542,
   - 1.014819592594577564927398399257428492725,
   - 1.016578512742400177512255407989807698099,
   - 1.018170476517182636083407324035097024753,
   - 1.019635932177973215702948823070039007361,
   - 1.021001233780440980009663527189756124983,
   - 1.022284686760270309618528459224732558767,
   - 1.023499619082082348038906498637836105447,
   - 1.033342436522109918536891072083243897233,
   - 1.040939194524944844012076988438386306843,
   - 1.047373634492196231421755878017895964061,
   - 1.053065496629607111615572634090884988897,
   - 1.058230030703619902820337106917657189605,
   - 1.06299509412938704964298950857825124135,
   - 1.067443986111355120560366087010834293872,
   - 1.071634561663924136735470389832541413862,
   - 1.07560894118662498941494486924522316597,
   - 1.108081880631165502564629660944418191031,
   - 1.133487001006868638317076487349933855181,
   - 1.155245851821528613609784258821055266176,
   - 1.174682608817289477552149783867901714817,
   - 1.192475850408615960644596781702366026321,
   - 1.209028378276581220769059281765172085749,
   - 1.224602449817731587352403997390766826335,
   - 1.239380103200799714836392811991400433357,
   - 1.253493791367214516100457405907304877145);

{ Lower branch, X close to 0 }
const WM2 : array[1..68] of Extended =
  (- 96.67475603368003636615083422832414231073,
   - 95.97433737593292677679699834708774152264,
   - 95.56459382507349364043974513871359837914,
   - 95.27386489130628866261760496192716897551,
   - 95.0483515329550645558378163981346085731,
   - 94.86408948075132603599669916724611501067,
   - 94.70829516116125928735505687861553800851,
   - 94.57333777268864473984718625104978190798,
   - 94.45429521137271454134108055166168787618,
   - 94.34780665137385269060032461028588648619,
   - 73.37311031382297679706747875812087452918,
   - 72.67033891766978907253811160121558649554,
   - 72.25920015786413889986246462168541066725,
   - 71.9674726772681844325410195162521230103,
   - 71.74117979478106456261839111929684980709,
   - 71.55627755942675731851469544735603279342,
   - 71.39993966508440988906136771384270319452,
   - 71.26450969134836299738230265916604879247,
   - 71.14504894849287026061378869987876478469,
   - 71.03818524971357411174259539994036186653,
   - 49.96298427667447244531514297262540669957,
   - 49.25557728489066973476436802404294348267,
   - 48.84167348449764278180827692232244878061,
   - 48.54795966722336777861228992424053274064,
   - 48.32011181512544381639923088992262632623,
   - 48.13392971864323755677656581326964166718,
   - 47.97650308277095858785998266172487372203,
   - 47.84012504158555569017852884133421231303,
   - 47.71982419568730714141619502661365943996,
   - 47.61220592218922310708330388890925734186,
   - 26.29523881924692569411012882185491823773,
   - 25.57429135222126159950976461862116397347,
   - 25.15218334705420805339928870686463192335,
   - 24.85251554543232259250342343156454440592,
   - 24.61997095867949438248843689371454503831,
   - 24.42989922074834124324823589226341161866,
   - 24.26914663885402405126372567664280388966,
   - 24.12985944288624210229238972590881794092,
   - 24.00697058168597098928369714836882130512,
   - 23.8970195845316574350263109196222825525,
   - 14.16360081581018300910955630361089957762,
   - 13.41624453595298662833544556875899262976,
   - 12.97753279184081358418630625949303360266,
   - 12.66551396826200331850774017793451747947,
   - 12.42304039760186078066171146072124458063,
   - 12.22461776385387453853455424320739669321,
   - 12.05663003490708840623665404674007291018,
   - 11.91094134143842011964821167497982287763,
   - 11.78229922740701885487699061601349928173,
   - 11.66711453256635441837882744697047370583,
   - 10.90655739570090676132157335673785028979,
   - 10.45921112040100393534625826514848865968,
   - 10.14059243262036578763968437893562720385,
   - 9.892699522704254067620287857665824159861,
   - 9.689637966382397752838283301312347921626,
   - 9.517569762038614935107630230444563521109,
   - 9.368222172408836799233763466046500781388,
   - 9.236251966692597369166416348621131600216,
   - 9.11800647040274012125833718204681427427,
   - 8.335081377982507150789361715143483020265,
   - 7.872521380098708883395239767904984410792,
   - 7.541940416432904084217222998374802941589,
   - 7.283997135099081646930521042317118095276,
   - 7.072162048994701667487346245044653243434,
   - 6.892241486671583156187212318718730022068,
   - 6.735741661607793269808533725369490789074,
   - 6.597171733627119347342347717832724288261,
   - 6.472775124394004694741057892724488037104);

var
  I, NDT : Integer;
  W : Float;

  function NDigits(W_appr : Float; W_exact : Extended) : Integer;
  { Returns the number of correct digits found in W_appr }
  begin
    if W_appr = W_exact then
      NDigits := NDT
    else
      NDigits := Round(Log10(Abs(W_exact / (W_appr - W_exact))) + 0.5);
  end;

  procedure PrintTitle(Msg : String; Offset : Boolean);
  { Prints the title for a group of results }
  var
    I : Integer;
    S : String;
  begin
    WriteLn;
    WriteLn(Msg);
    for I := 1 to Length(Msg) do
      Write('-');
    WriteLn; WriteLn;
    if Offset then
      S := ' Offset X'
    else
      S := '    X    ';
    WriteLn(S, '          W(X) (appr)       W(X) (exact)      Digits correct');
  end;

  procedure PrintResult(X, W_appr : Float; W_exact : Extended);
  { Prints one line of results }
  begin
    WriteLn(X:17, ' ', W_appr:17, ' ', W_exact:17, ' ',
            NDigits(W_appr, W_exact):9);
  end;

begin
  NDT := Round(0.5 - Log10(MACHEP));

  WriteLn('Lambert''s W-function');
  PrintTitle('Upper branch results for X near -1/e', True);

  for I := 1 to 68 do
    begin
      { Check whether underflow occurs }
      if DX1[I] = 0.0 then
        W := LambertW(- Exp(- 1.0), True, False)
      else
        W := LambertW(DX1[I], True, True);

      PrintResult(DX1[I], W, WP1[I]);
    end;

  PrintTitle('Upper branch results for X near 0', False);

  for I := 1 to 20 do
    begin
      W := LambertW(X2[I], True, False);
      PrintResult(X2[I], W, WP2[I]);
    end;

  for I := 1 to 20 do
    begin
      W := LambertW(- X2[I], True, False);
      PrintResult(- X2[I], W, WP2[20 + I]);
    end;

  PrintTitle('Other upper branch results', False);

  for I := 1 to 10 do
    begin
      W := LambertW(X3[I], True, False);
      PrintResult(X3[I], W, WP3[I]);
    end;

  PrintTitle('Lower branch results for X near -1/e', True);

  for I := 1 to 68 do
    begin
      W := LambertW(DX1[I], False, True);
      PrintResult(DX1[I], W, WM1[I]);
    end;

  PrintTitle('Lower branch results for X near 0', False);

  for I := 1 to 68 do
    { Check for underflow }
    if DX1[I] >= 0.0 then
      begin
        W := LambertW(- DX1[I], False, False);
        PrintResult(- DX1[I], W, WM2[I]);
      end;
end.