home *** CD-ROM | disk | FTP | other *** search
/ PC Pro 1998 April / DPPCPRO0498.ISO / April / MathCad / SETUP / DATA.Z / DIGAMMA.C < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-31  |  10.4 KB  |  629 lines
  1. /* $Id: digamma.c,v 1.1 1994/08/23 13:34:54 lorczak Exp $   */
  2. /*                 $Revision: 1 $        */
  3. #include "mcadincl.h"
  4.  
  5. #include "math.h"
  6.  
  7.  
  8.  
  9. #define PI 3.141592653589793
  10.  
  11.  
  12.  
  13. #define  DOMAIN_ERROR           1
  14.  
  15. #define  NUMBER_OF_ERRORS       1   
  16.  
  17.  
  18.  
  19.     // Table of error messages.
  20.  
  21.     // 
  22.  
  23.     // There is one message to signal the case of an input
  24.  
  25.     // for which digamma is undefined.  Note there is no interrupt
  26.  
  27.     // message, since this function simply does a little 
  28.  
  29.     // arithmetic and so will be be too fast to interrupt.
  30.  
  31.  
  32.  
  33.     char * myErrorMessageTable[NUMBER_OF_ERRORS] =  
  34.  
  35.     {    "domain error"
  36.  
  37.     };
  38.  
  39.     
  40.  
  41.     
  42.  
  43.     // This code executes the digamma function computation.
  44.  
  45.  
  46.  
  47.     LRESULT  DigammaFunction(  COMPLEXSCALAR * const p, 
  48.  
  49.         const COMPLEXSCALAR * const z)
  50.  
  51.     {
  52.  
  53.         double Ber[16];
  54.  
  55.         double xr,xi,c1,c2,norm,xm12r,xm12i,tr,ti,sumr,sumi,trtemp;
  56.  
  57.         double sr,cr,shi,chi,fac;
  58.  
  59.  
  60.  
  61.         unsigned int i, Situation;
  62.  
  63.  
  64.  
  65.         // check for integer inputs less than or equal to 0
  66.  
  67.  
  68.  
  69.         if ( (z->imag == 0.0) && (z->real <= 0.0) && (floor(z->real) == ceil(z->real)))
  70.  
  71.                         return MAKELRESULT(DOMAIN_ERROR, 1 );
  72.  
  73.  
  74.  
  75.  
  76.  
  77.         // Create a vector containing the values of the necessary 
  78.  
  79.         // even-order Bernoulli numbers.
  80.  
  81.  
  82.  
  83.         Ber[1] = 0.166666666666667;
  84.  
  85.         Ber[2] = -0.033333333333333;
  86.  
  87.         Ber[3] = 0.023809523809524;
  88.  
  89.         Ber[4] = -0.033333333333334;
  90.  
  91.         Ber[5] = 0.075757575757579;
  92.  
  93.         Ber[6] = -0.253113553113593;
  94.  
  95.         Ber[7] = 1.166666666667406;
  96.  
  97.         Ber[8] = -7.092156862763144;
  98.  
  99.         Ber[9] = 54.97117794542229;
  100.  
  101.         Ber[10] = -529.1242424458143;
  102.  
  103.         Ber[11] = 6.192123189415645e+3;
  104.  
  105.         Ber[12] = -8.658025317003414e+4;
  106.  
  107.         Ber[13] = 1.425517170386453e+6;
  108.  
  109.         Ber[14] = -2.729823135274775e+7;
  110.  
  111.         Ber[15] = 6.015808990172174e+8;
  112.  
  113.  
  114.  
  115.         // Given the argument to digamma determine which situation
  116.  
  117.         // is present (large real value, small real value, etc.).
  118.  
  119.  
  120.  
  121.               Situation = 3;
  122.  
  123.         if (z->real >= 7.0)
  124.  
  125.               Situation = 0;
  126.  
  127.         if (z->real <= -6.0)
  128.  
  129.               Situation = 1;
  130.  
  131.         if  ((z->real >=  0.0) && (z->real <7.0))
  132.  
  133.               Situation = 2;
  134.  
  135.  
  136.  
  137.         // If you examine the Mathcad implementation of the digamma algorithm
  138.  
  139.         // you'll see that for a given input z, the basic large value model
  140.  
  141.         // is subsequently applied to either
  142.  
  143.         //          x = z
  144.  
  145.         //          x = 1 - z
  146.  
  147.         //          x =  z + 7
  148.  
  149.         //          x = 8 - z
  150.  
  151.         // These are the 4 situations tested for above.  The corresponding value of
  152.  
  153.         // x is computed below.  
  154.  
  155.  
  156.  
  157.  
  158.  
  159.         switch(Situation)
  160.  
  161.          {
  162.  
  163.              case 0:
  164.  
  165.                       xr = z->real;
  166.  
  167.                       xi = z->imag;
  168.  
  169.                       break;
  170.  
  171.  
  172.  
  173.              case 1:
  174.  
  175.                       xr = 1.0 - z->real;
  176.  
  177.                       xi =  - z->imag;
  178.  
  179.                       break;
  180.  
  181.  
  182.  
  183.              case 2:
  184.  
  185.                       xr = 7.0 + z->real;
  186.  
  187.                       xi = z->imag;
  188.  
  189.                       break;
  190.  
  191.  
  192.  
  193.              case 3:
  194.  
  195.                       xr = 8.0 - z->real;
  196.  
  197.                       xi = -z->imag;
  198.  
  199.                       break;
  200.  
  201.          }
  202.  
  203.  
  204.  
  205.         // The value of x will now have real part xr >= 7.0 so
  206.  
  207.         // compute digamma(x) using the large x model.
  208.  
  209.  
  210.  
  211.         // Start with the summation involving the Bernoulli numbers.
  212.  
  213.         // The statements below compute the first term of the sum
  214.  
  215.         //              Ber[1]/(2*(x-1)^2)
  216.  
  217.         // Note we have to allow for complex arithmetic.
  218.  
  219.  
  220.  
  221.         c1 =  2.0*xi*(xr-1.0);
  222.  
  223.         c2 =  xr*xr - 2.0*xr - xi*xi +1.0;
  224.  
  225.         norm = c1*c1 + c2*c2;
  226.  
  227.  
  228.  
  229.         xm12r = c2/norm;
  230.  
  231.         xm12i = -c1/norm;
  232.  
  233.  
  234.  
  235.         tr = xm12r;
  236.  
  237.         ti = xm12i;
  238.  
  239.         sumr = (Ber[1]/2.0) * tr;
  240.  
  241.         sumi = (Ber[1]/2.0) * ti;
  242.  
  243.  
  244.  
  245.         // Compute the rest of the summation piece.
  246.  
  247.  
  248.  
  249.         for (i=2;i<16;i++)
  250.  
  251.             {
  252.  
  253.                   trtemp = xm12r*tr - xm12i*ti;
  254.  
  255.                   ti = xm12r*ti + xm12i*tr;
  256.  
  257.                   tr = trtemp;
  258.  
  259.  
  260.  
  261.                   sumr = sumr + (Ber[i]/(2.0 * (double) i))*tr;
  262.  
  263.                   sumi = sumi + (Ber[i]/(2.0 * (double) i))*ti;
  264.  
  265.  
  266.  
  267.              }
  268.  
  269.  
  270.  
  271.  
  272.  
  273.         // Compute 1/(2*(x-1)).
  274.  
  275.  
  276.  
  277.         tr = (xr - 1.0) /(2.0 * ((xr-1.0) * (xr-1.0) + xi*xi)) - sumr;
  278.  
  279.         ti = - xi /(2.0 * ((xr-1.0) * (xr-1.0) + xi*xi)) - sumi;
  280.  
  281.  
  282.  
  283.         // Finally, add the possibly complex natural log(x - 1).
  284.  
  285.  
  286.  
  287.         p->real = tr + 0.5*log((xr-1.0) * (xr-1.0) + xi*xi) ;
  288.  
  289.         p->imag = ti + atan2(xi,xr-1.0);
  290.  
  291.  
  292.  
  293.  
  294.  
  295.         // The structure p now contains the value Psi(x).
  296.  
  297.         // Depending on the situation, p must be adjusted by either
  298.  
  299.         //               Pi/tan(Pi*z)
  300.  
  301.         //  or           sum(1/(z+i),i=1..6)
  302.  
  303.         //  or both to produce Psi(z).
  304.  
  305.         //  The appropriate adjustment is done below.   
  306.  
  307.  
  308.  
  309.         switch(Situation)
  310.  
  311.          {
  312.  
  313.              case 0:          // Psi(x) = Psi(z)
  314.  
  315.                       break;
  316.  
  317.  
  318.  
  319.              case 1:         // adjust by Pi(tan(Pi*z)
  320.  
  321.  
  322.  
  323.                       // compute possibly complex Pi/tan(Pi*z)
  324.  
  325.  
  326.  
  327.                       sr = sin(PI*z->real);
  328.  
  329.                       cr = cos(PI*z->real);
  330.  
  331.                       shi = sinh(PI*z->imag);
  332.  
  333.                       chi = cosh(PI*z->imag);
  334.  
  335.  
  336.  
  337.                       fac = PI*(cr*cr + shi*shi)/(cr*cr*sr*sr + chi*chi*shi*shi);
  338.  
  339.                       tr = sr*cr*fac;
  340.  
  341.                       ti = -shi*chi*fac;
  342.  
  343.  
  344.  
  345.                       p->real -= tr;
  346.  
  347.                       p->imag -= ti;
  348.  
  349.  
  350.  
  351.                       break;
  352.  
  353.  
  354.  
  355.              case 2:          // adjust by summation
  356.  
  357.  
  358.  
  359.                       tr = 0.0;
  360.  
  361.                       ti = 0.0;
  362.  
  363.                       for (i=0;i<7;i++)
  364.  
  365.                            { 
  366.  
  367.                               tr = tr + (z->real + i)/ ((z->real + i)* (z->real + i) + (z->imag)*(z->imag));
  368.  
  369.                               ti = ti + 1.0/ ((z->real + i)* (z->real + i) + (z->imag)*(z->imag));
  370.  
  371.                             }
  372.  
  373.  
  374.  
  375.                       ti = -(z->imag)*ti;
  376.  
  377.  
  378.  
  379.                       p->real -= tr;
  380.  
  381.                       p->imag -= ti;
  382.  
  383.  
  384.  
  385.                       break;
  386.  
  387.  
  388.  
  389.              case 3:         // adjust by summation and value of Pi/tan(Pi*z)
  390.  
  391.  
  392.  
  393.                       sr = sin(PI*z->real);
  394.  
  395.                       cr = cos(PI*z->real);
  396.  
  397.                       shi = sinh(PI*z->imag);
  398.  
  399.                       chi = cosh(PI*z->imag);
  400.  
  401.  
  402.  
  403.                       fac = PI*(cr*cr + shi*shi)/(cr*cr*sr*sr + chi*chi*shi*shi);
  404.  
  405.                       tr = sr*cr*fac;
  406.  
  407.                       ti = -shi*chi*fac;
  408.  
  409.  
  410.  
  411.                       p->real -= tr;
  412.  
  413.                       p->imag -= ti;
  414.  
  415.  
  416.  
  417.  
  418.  
  419.                       tr = 0.0;
  420.  
  421.                       ti = 0.0;
  422.  
  423.                       for (i=1;i<8;i++)
  424.  
  425.                            { 
  426.  
  427.                               tr = tr + (i - z->real)/ ((i - z->real )* (i - z->real) + (z->imag)*(z->imag));
  428.  
  429.                               ti = ti + 1.0/ ((i - z->real )* (i - z->real)  + (z->imag)*(z->imag));
  430.  
  431.                             }
  432.  
  433.  
  434.  
  435.                       ti = (z->imag)*ti;
  436.  
  437.  
  438.  
  439.                       p->real -= tr;
  440.  
  441.                       p->imag -= ti;
  442.  
  443.  
  444.  
  445.  
  446.  
  447.                       break;
  448.  
  449.          }
  450.  
  451.  
  452.  
  453.         // All done.
  454.  
  455.  
  456.  
  457.         return 0;
  458.  
  459.     }
  460.  
  461.  
  462.  
  463.     // Fill out a FUNCTIONINFO structure with
  464.  
  465.     // the information needed for registering the
  466.  
  467.     // function with Mathcad.
  468.  
  469.     FUNCTIONINFO Psi = 
  470.  
  471.     {
  472.  
  473.     // First, the name by which Mathcad will recognize the function.
  474.  
  475.     "Psi",        
  476.  
  477.     
  478.  
  479.     // Second, adescription of "Psi" parameters to be used
  480.  
  481.     // by the Insert Function dialog box.
  482.  
  483.     "z",   
  484.  
  485.     
  486.  
  487.     // Then a description of the function for the Insert Function dialog box.       
  488.  
  489.     "Digamma function for complex z",    
  490.  
  491.     
  492.  
  493.     // Now a pointer to the executable code,
  494.  
  495.     // i.e., code that should be executed
  496.  
  497.     // when a user types in "Psi(z)=".
  498.  
  499.     (LPCFUNCTION)DigammaFunction,  
  500.  
  501.         
  502.  
  503.     // Psi(z) returns a complex array.
  504.  
  505.     COMPLEX_SCALAR,
  506.  
  507.         
  508.  
  509.     // Psi(z) takes one argument.
  510.  
  511.     1,   
  512.  
  513.     
  514.  
  515.     // a complex scalar.
  516.  
  517.     COMPLEX_SCALAR 
  518.  
  519.     };
  520.  
  521.     
  522.  
  523.  
  524.  
  525.     // DLL entry point code.
  526.  
  527.     // The _CRT_INIT function is needed
  528.  
  529.     // if you are using Microsoft's 32-bit compiler
  530.  
  531.  
  532.  
  533.     BOOL WINAPI _CRT_INIT(HINSTANCE hinstDLL, DWORD dwReason, LPVOID lpReserved);
  534.  
  535.  
  536.  
  537.     BOOL WINAPI  DllEntryPoint (HINSTANCE hDLL, DWORD dwReason, LPVOID lpReserved)
  538.  
  539.     {
  540.  
  541.         switch (dwReason)
  542.  
  543.         {
  544.  
  545.                 case DLL_PROCESS_ATTACH:
  546.  
  547.  
  548.  
  549.                 //
  550.  
  551.                 // DLL is attaching to the address space of 
  552.  
  553.                 // the current process.
  554.  
  555.                 //
  556.  
  557.                     if (!_CRT_INIT(hDLL, dwReason, lpReserved)) 
  558.  
  559.                         return FALSE;
  560.  
  561.                 
  562.  
  563.                 
  564.  
  565.                     // Register the error message table.
  566.  
  567.                     // Note that if your function never returns
  568.  
  569.                     // an error, you do not need to 
  570.  
  571.                     // register an error message table.
  572.  
  573.                     if ( CreateUserErrorMessageTable(
  574.  
  575.                             hDLL, NUMBER_OF_ERRORS, myErrorMessageTable ) )
  576.  
  577.                         // And if the errors register OK
  578.  
  579.                         // go ahead and
  580.  
  581.                         // register user function.
  582.  
  583.                         CreateUserFunction( hDLL, &Psi );
  584.  
  585.                     break;
  586.  
  587.  
  588.  
  589.  
  590.  
  591.                 case DLL_THREAD_ATTACH:
  592.  
  593.                 case DLL_THREAD_DETACH:
  594.  
  595.                 case DLL_PROCESS_DETACH:
  596.  
  597.  
  598.  
  599.                     if (!_CRT_INIT(hDLL, dwReason, lpReserved)) 
  600.  
  601.                         return FALSE;
  602.  
  603.                     break;
  604.  
  605.                     
  606.  
  607.         }
  608.  
  609.         return TRUE;
  610.  
  611.     }
  612.  
  613.  
  614.  
  615. #undef INTERRUPTED    
  616.  
  617. #undef INSUFFICIENT_MEMORY
  618.  
  619. #undef MUST_BE_REAL   
  620.  
  621. #undef NUMBER_OF_ERRORS     
  622.  
  623.  
  624.  
  625.     
  626.  
  627.  
  628.  
  629.