Print this page
11210 libm should be cstyle(1ONBLD) clean

Split Close
Expand all
Collapse all
          --- old/usr/src/lib/libm/common/Q/__lgammal.c
          +++ new/usr/src/lib/libm/common/Q/__lgammal.c
↓ open down ↓ 14 lines elided ↑ open up ↑
  15   15   * If applicable, add the following below this CDDL HEADER, with the
  16   16   * fields enclosed by brackets "[]" replaced with your own identifying
  17   17   * information: Portions Copyright [yyyy] [name of copyright owner]
  18   18   *
  19   19   * CDDL HEADER END
  20   20   */
  21   21  
  22   22  /*
  23   23   * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24   24   */
       25 +
  25   26  /*
  26   27   * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27   28   * Use is subject to license terms.
  28   29   */
  29   30  
  30   31  /*
  31   32   * long double __k_lgammal(long double x, int *signgamlp);
  32   33   * K.C. Ng, August, 1989.
  33   34   *
  34   35   * We choose [1.5,2.5] to be the primary interval. Our algorithms
  35   36   * are mainly derived from
  36   37   *
  37   38   *
  38   39   *                             zeta(2)-1    2    zeta(3)-1    3
  39   40   * lgamma(2+s) = s*(1-euler) + --------- * s  -  --------- * s  + ...
  40   41   *                                 2                 3
  41   42   *
  42   43   *
  43   44   * Note 1. Since gamma(1+s)=s*gamma(s), hence
  44      - *              lgamma(1+s) = log(s) + lgamma(s), or
  45      - *              lgamma(s) = lgamma(1+s) - log(s).
       45 + *              lgamma(1+s) = log(s) + lgamma(s), or
       46 + *              lgamma(s) = lgamma(1+s) - log(s).
  46   47   *         When s is really tiny (like roundoff), lgamma(1+s) ~ s(1-enler)
  47   48   *         Hence lgamma(s) ~ -log(s) for tiny s
  48   49   *
  49   50   */
  50   51  
  51   52  #include "libm.h"
  52   53  #include "longdouble.h"
  53   54  
  54   55  static long double neg(long double, int *);
  55   56  static long double poly(long double, const long double *, int);
  56   57  static long double polytail(long double);
  57   58  static long double primary(long double);
  58   59  
  59      -static const long double
  60      -c0 =     0.0L,
  61      -ch =     0.5L,
  62      -c1 =     1.0L,
  63      -c2 =     2.0L,
  64      -c3 =     3.0L,
  65      -c4 =     4.0L,
  66      -c5 =     5.0L,
  67      -c6 =     6.0L,
  68      -pi =     3.1415926535897932384626433832795028841971L,
  69      -tiny =   1.0e-40L;
       60 +static const long double c0 = 0.0L,
       61 +        ch = 0.5L,
       62 +        c1 = 1.0L,
       63 +        c2 = 2.0L,
       64 +        c3 = 3.0L,
       65 +        c4 = 4.0L,
       66 +        c5 = 5.0L,
       67 +        c6 = 6.0L,
       68 +        pi = 3.1415926535897932384626433832795028841971L,
       69 +        tiny = 1.0e-40L;
  70   70  
  71   71  long double
  72      -__k_lgammal(long double x, int *signgamlp) {
  73      -        long double t,y;
       72 +__k_lgammal(long double x, int *signgamlp)
       73 +{
       74 +        long double t, y;
  74   75          int i;
  75   76  
  76      -    /* purge off +-inf, NaN and negative arguments */
  77      -        if (!finitel(x)) return x*x;
       77 +        /* purge off +-inf, NaN and negative arguments */
       78 +        if (!finitel(x))
       79 +                return (x * x);
       80 +
  78   81          *signgamlp = 1;
  79      -        if (signbitl(x)) return (neg(x,signgamlp));
  80   82  
  81      -    /* for x < 8.0 */
  82      -        if (x<8.0L) {
  83      -            y = anintl(x);
  84      -            i = (int) y;
  85      -            switch(i) {
  86      -            case 0:
  87      -                if (x<1.0e-40L) return -logl(x); else
  88      -                return (primary(x)-log1pl(x))-logl(x);
  89      -            case 1:
  90      -                return primary(x-y)-logl(x);
  91      -            case 2:
  92      -                return primary(x-y);
  93      -            case 3:
  94      -                return primary(x-y)+logl(x-c1);
  95      -            case 4:
  96      -                return primary(x-y)+logl((x-c1)*(x-c2));
  97      -            case 5:
  98      -                return primary(x-y)+logl((x-c1)*(x-c2)*(x-c3));
  99      -            case 6:
 100      -                return primary(x-y)+logl((x-c1)*(x-c2)*(x-c3)*(x-c4));
 101      -            case 7:
 102      -                return primary(x-y)+logl((x-c1)*(x-c2)*(x-c3)*(x-c4)*(x-c5));
 103      -            case 8:
 104      -                return primary(x-y)+
 105      -                        logl((x-c1)*(x-c2)*(x-c3)*(x-c4)*(x-c5)*(x-c6));
 106      -            }
       83 +        if (signbitl(x))
       84 +                return (neg(x, signgamlp));
       85 +
       86 +        /* for x < 8.0 */
       87 +        if (x < 8.0L) {
       88 +                y = anintl(x);
       89 +                i = (int)y;
       90 +
       91 +                switch (i) {
       92 +                case 0:
       93 +
       94 +                        if (x < 1.0e-40L)
       95 +                                return (-logl(x));
       96 +                        else
       97 +                                return ((primary(x) - log1pl(x)) - logl(x));
       98 +
       99 +                case 1:
      100 +                        return (primary(x - y) - logl(x));
      101 +                case 2:
      102 +                        return (primary(x - y));
      103 +                case 3:
      104 +                        return (primary(x - y) + logl(x - c1));
      105 +                case 4:
      106 +                        return (primary(x - y) + logl((x - c1) * (x - c2)));
      107 +                case 5:
      108 +                        return (primary(x - y) + logl((x - c1) * (x - c2) * (x -
      109 +                            c3)));
      110 +                case 6:
      111 +                        return (primary(x - y) + logl((x - c1) * (x - c2) * (x -
      112 +                            c3) * (x - c4)));
      113 +                case 7:
      114 +                        return (primary(x - y) + logl((x - c1) * (x - c2) * (x -
      115 +                            c3) * (x - c4) * (x - c5)));
      116 +                case 8:
      117 +                        return (primary(x - y) + logl((x - c1) * (x - c2) * (x -
      118 +                            c3) * (x - c4) * (x - c5) * (x - c6)));
      119 +                }
 107  120          }
 108  121  
 109      -    /* 8.0 <= x < 1.0e40 */
      122 +        /* 8.0 <= x < 1.0e40 */
 110  123          if (x < 1.0e40L) {
 111      -            t = logl(x);
 112      -            return x*(t-c1)-(ch*t-polytail(c1/x));
      124 +                t = logl(x);
      125 +                return (x * (t - c1) - (ch * t - polytail(c1 / x)));
 113  126          }
 114  127  
 115      -    /* 1.0e40 <= x <= inf */
 116      -        return x*(logl(x)-c1);
      128 +        /* 1.0e40 <= x <= inf */
      129 +        return (x * (logl(x) - c1));
 117  130  }
 118  131  
 119      -static const long double an1[] = {              /* 20 terms */
 120      -  -0.0772156649015328606065120900824024309741L,
 121      -   3.224670334241132182362075833230130289059e-0001L,
 122      -  -6.735230105319809513324605383668929964120e-0002L,
 123      -   2.058080842778454787900092432928910226297e-0002L,
 124      -  -7.385551028673985266273054086081102125704e-0003L,
 125      -   2.890510330741523285758867304409628648727e-0003L,
 126      -  -1.192753911703260976581414338096267498555e-0003L,
 127      -   5.096695247430424562831956662855697824035e-0004L,
 128      -  -2.231547584535777978926798502084300123638e-0004L,
 129      -   9.945751278186384670278268034322157947635e-0005L,
 130      -  -4.492623673665547726647838474125147631082e-0005L,
 131      -   2.050721280617796810096993154281561168706e-0005L,
 132      -  -9.439487785617396552092393234044767313568e-0006L,
 133      -   4.374872903516051510689234173139793159340e-0006L,
 134      -  -2.039156676413643091040459825776029327487e-0006L,
 135      -   9.555777181318621470466563543806211523634e-0007L,
 136      -  -4.468344919709630637558538313482398989638e-0007L,
 137      -   2.216738086090045781773004477831059444178e-0007L,
 138      -  -7.472783403418388455860445842543843485916e-0008L,
 139      -   8.777317930927149922056782132706238921648e-0008L,
      132 +static const long double an1[] = {      /* 20 terms */
      133 +        -0.0772156649015328606065120900824024309741L,
      134 +        3.224670334241132182362075833230130289059e-0001L,
      135 +        -6.735230105319809513324605383668929964120e-0002L,
      136 +        2.058080842778454787900092432928910226297e-0002L,
      137 +        -7.385551028673985266273054086081102125704e-0003L,
      138 +        2.890510330741523285758867304409628648727e-0003L,
      139 +        -1.192753911703260976581414338096267498555e-0003L,
      140 +        5.096695247430424562831956662855697824035e-0004L,
      141 +        -2.231547584535777978926798502084300123638e-0004L,
      142 +        9.945751278186384670278268034322157947635e-0005L,
      143 +        -4.492623673665547726647838474125147631082e-0005L,
      144 +        2.050721280617796810096993154281561168706e-0005L,
      145 +        -9.439487785617396552092393234044767313568e-0006L,
      146 +        4.374872903516051510689234173139793159340e-0006L,
      147 +        -2.039156676413643091040459825776029327487e-0006L,
      148 +        9.555777181318621470466563543806211523634e-0007L,
      149 +        -4.468344919709630637558538313482398989638e-0007L,
      150 +        2.216738086090045781773004477831059444178e-0007L,
      151 +        -7.472783403418388455860445842543843485916e-0008L,
      152 +        8.777317930927149922056782132706238921648e-0008L,
 140  153  };
 141  154  
 142      -static const long double an2[] = {              /* 20 terms */
 143      -  -.0772156649015328606062692723698127607018L,
 144      -   3.224670334241132182635552349060279118047e-0001L,
 145      -  -6.735230105319809367555642883133994818325e-0002L,
 146      -   2.058080842778459676880822202762143671813e-0002L,
 147      -  -7.385551028672828216011343150077846918930e-0003L,
 148      -   2.890510330762060607399561536905727853178e-0003L,
 149      -  -1.192753911419623262328187532759756368041e-0003L,
 150      -   5.096695278636456678258091134532258618614e-0004L,
 151      -  -2.231547306817535743052975194022893369135e-0004L,
 152      -   9.945771461633313282744264853986643877087e-0005L,
 153      -  -4.492503279458972037926876061257489481619e-0005L,
 154      -   2.051311416812082875492678651369394595613e-0005L,
 155      -  -9.415778282365955203915850761537462941165e-0006L,
 156      -   4.452428829045147098722932981088650055919e-0006L,
 157      -  -1.835024727987632579886951760650722695781e-0006L,
 158      -   1.379783080658545009579060714946381462565e-0006L,
 159      -   2.282637532109775156769736768748402175238e-0007L,
 160      -   1.002577375515900191362119718128149880168e-0006L,
 161      -   5.177028794262638311939991106423220002463e-0007L,
 162      -   3.127947245174847104122426445937830555755e-0007L,
      155 +static const long double an2[] = {      /* 20 terms */
      156 +        -.0772156649015328606062692723698127607018L,
      157 +        3.224670334241132182635552349060279118047e-0001L,
      158 +        -6.735230105319809367555642883133994818325e-0002L,
      159 +        2.058080842778459676880822202762143671813e-0002L,
      160 +        -7.385551028672828216011343150077846918930e-0003L,
      161 +        2.890510330762060607399561536905727853178e-0003L,
      162 +        -1.192753911419623262328187532759756368041e-0003L,
      163 +        5.096695278636456678258091134532258618614e-0004L,
      164 +        -2.231547306817535743052975194022893369135e-0004L,
      165 +        9.945771461633313282744264853986643877087e-0005L,
      166 +        -4.492503279458972037926876061257489481619e-0005L,
      167 +        2.051311416812082875492678651369394595613e-0005L,
      168 +        -9.415778282365955203915850761537462941165e-0006L,
      169 +        4.452428829045147098722932981088650055919e-0006L,
      170 +        -1.835024727987632579886951760650722695781e-0006L,
      171 +        1.379783080658545009579060714946381462565e-0006L,
      172 +        2.282637532109775156769736768748402175238e-0007L,
      173 +        1.002577375515900191362119718128149880168e-0006L,
      174 +        5.177028794262638311939991106423220002463e-0007L,
      175 +        3.127947245174847104122426445937830555755e-0007L,
 163  176  };
 164  177  
 165      -static const long double an3[] = {              /* 20 terms */
 166      -  -.0772156649015328227870646417729220690875L,
 167      -   3.224670334241156699881788955959915250365e-0001L,
 168      -  -6.735230105312273571375431059744975563170e-0002L,
 169      -   2.058080842924464587662846071337083809005e-0002L,
 170      -  -7.385551008677271654723604653956131791619e-0003L,
 171      -   2.890510536479782086197110272583833176602e-0003L,
 172      -  -1.192752262076857692740571567808259138697e-0003L,
 173      -   5.096800771149805289371135155128380707889e-0004L,
 174      -  -2.231000836682831335505058492409860123647e-0004L,
 175      -   9.968912171073936803871803966360595275047e-0005L,
 176      -  -4.412020779327746243544387946167256187258e-0005L,
 177      -   2.281374113541454151067016632998630209049e-0005L,
 178      -  -4.028361291428629491824694655287954266830e-0006L,
 179      -   1.470694920619518924598956849226530750139e-0005L,
 180      -   1.381686137617987197975289545582377713772e-0005L,
 181      -   2.012493539265777728944759982054970441601e-0005L,
 182      -   1.723917864208965490251560644681933675799e-0005L,
 183      -   1.202954035243788300138608765425123713395e-0005L,
 184      -   5.079851887558623092776296577030850938146e-0006L,
 185      -   1.220657945824153751555138592006604026282e-0006L,
      178 +static const long double an3[] = {      /* 20 terms */
      179 +        -.0772156649015328227870646417729220690875L,
      180 +        3.224670334241156699881788955959915250365e-0001L,
      181 +        -6.735230105312273571375431059744975563170e-0002L,
      182 +        2.058080842924464587662846071337083809005e-0002L,
      183 +        -7.385551008677271654723604653956131791619e-0003L,
      184 +        2.890510536479782086197110272583833176602e-0003L,
      185 +        -1.192752262076857692740571567808259138697e-0003L,
      186 +        5.096800771149805289371135155128380707889e-0004L,
      187 +        -2.231000836682831335505058492409860123647e-0004L,
      188 +        9.968912171073936803871803966360595275047e-0005L,
      189 +        -4.412020779327746243544387946167256187258e-0005L,
      190 +        2.281374113541454151067016632998630209049e-0005L,
      191 +        -4.028361291428629491824694655287954266830e-0006L,
      192 +        1.470694920619518924598956849226530750139e-0005L,
      193 +        1.381686137617987197975289545582377713772e-0005L,
      194 +        2.012493539265777728944759982054970441601e-0005L,
      195 +        1.723917864208965490251560644681933675799e-0005L,
      196 +        1.202954035243788300138608765425123713395e-0005L,
      197 +        5.079851887558623092776296577030850938146e-0006L,
      198 +        1.220657945824153751555138592006604026282e-0006L,
 186  199  };
 187  200  
 188      -static const long double an4[] = {              /* 21 terms */
 189      -  -.0772156649015732285350261816697540392371L,
 190      -   3.224670334221752060691751340365212226097e-0001L,
 191      -  -6.735230109744009693977755991488196368279e-0002L,
 192      -   2.058080778913037626909954141611580783216e-0002L,
 193      -  -7.385557567931505621170483708950557506819e-0003L,
 194      -   2.890459838416254326340844289785254883436e-0003L,
 195      -  -1.193059036207136762877351596966718455737e-0003L,
 196      -   5.081914708100372836613371356529568937869e-0004L,
 197      -  -2.289855016133600313131553005982542045338e-0004L,
 198      -   8.053454537980585879620331053833498511491e-0005L,
 199      -  -9.574620532104845821243493405855672438998e-0005L,
 200      -  -9.269085628207107155601445001196317715686e-0005L,
 201      -  -2.183276779859490461716196344776208220180e-0004L,
 202      -  -3.134834305597571096452454999737269668868e-0004L,
 203      -  -3.973878894951937437018305986901392888619e-0004L,
 204      -  -3.953352414899222799161275564386488057119e-0004L,
 205      -  -3.136740932204038779362660900621212816511e-0004L,
 206      -  -1.884502253819634073946130825196078627664e-0004L,
 207      -  -8.192655799958926853585332542123631379301e-0005L,
 208      -  -2.292183750010571062891605074281744854436e-0005L,
 209      -  -3.223980628729716864927724265781406614294e-0006L,
      201 +static const long double an4[] = {      /* 21 terms */
      202 +        -.0772156649015732285350261816697540392371L,
      203 +        3.224670334221752060691751340365212226097e-0001L,
      204 +        -6.735230109744009693977755991488196368279e-0002L,
      205 +        2.058080778913037626909954141611580783216e-0002L,
      206 +        -7.385557567931505621170483708950557506819e-0003L,
      207 +        2.890459838416254326340844289785254883436e-0003L,
      208 +        -1.193059036207136762877351596966718455737e-0003L,
      209 +        5.081914708100372836613371356529568937869e-0004L,
      210 +        -2.289855016133600313131553005982542045338e-0004L,
      211 +        8.053454537980585879620331053833498511491e-0005L,
      212 +        -9.574620532104845821243493405855672438998e-0005L,
      213 +        -9.269085628207107155601445001196317715686e-0005L,
      214 +        -2.183276779859490461716196344776208220180e-0004L,
      215 +        -3.134834305597571096452454999737269668868e-0004L,
      216 +        -3.973878894951937437018305986901392888619e-0004L,
      217 +        -3.953352414899222799161275564386488057119e-0004L,
      218 +        -3.136740932204038779362660900621212816511e-0004L,
      219 +        -1.884502253819634073946130825196078627664e-0004L,
      220 +        -8.192655799958926853585332542123631379301e-0005L,
      221 +        -2.292183750010571062891605074281744854436e-0005L,
      222 +        -3.223980628729716864927724265781406614294e-0006L,
 210  223  };
 211  224  
 212      -static const long double ap1[] = {                      /* 19 terms */
 213      -  -0.0772156649015328606065120900824024296961L,
 214      -   3.224670334241132182362075833230047956465e-0001L,
 215      -  -6.735230105319809513324605382963943777301e-0002L,
 216      -   2.058080842778454787900092126606252375465e-0002L,
 217      -  -7.385551028673985266272518231365020063941e-0003L,
 218      -   2.890510330741523285681704570797770736423e-0003L,
 219      -  -1.192753911703260971285304221165990244515e-0003L,
 220      -   5.096695247430420878696018188830886972245e-0004L,
 221      -  -2.231547584535654004647639737841526025095e-0004L,
 222      -   9.945751278137201960636098805852315982919e-0005L,
 223      -  -4.492623672777606053587919463929044226280e-0005L,
 224      -   2.050721258703289487603702670753053765201e-0005L,
 225      -  -9.439485626565616989352750672499008021041e-0006L,
 226      -   4.374838162403994645138200419356844574219e-0006L,
 227      -  -2.038979492862555348577006944451002161496e-0006L,
 228      -   9.536763152382263548086981191378885102802e-0007L,
 229      -  -4.426111214332434049863595231916564014913e-0007L,
 230      -   1.911148847512947464234633846270287546882e-0007L,
 231      -  -5.788673944861923038157839080272303519671e-0008L,
      225 +static const long double ap1[] = {      /* 19 terms */
      226 +        -0.0772156649015328606065120900824024296961L,
      227 +        3.224670334241132182362075833230047956465e-0001L,
      228 +        -6.735230105319809513324605382963943777301e-0002L,
      229 +        2.058080842778454787900092126606252375465e-0002L,
      230 +        -7.385551028673985266272518231365020063941e-0003L,
      231 +        2.890510330741523285681704570797770736423e-0003L,
      232 +        -1.192753911703260971285304221165990244515e-0003L,
      233 +        5.096695247430420878696018188830886972245e-0004L,
      234 +        -2.231547584535654004647639737841526025095e-0004L,
      235 +        9.945751278137201960636098805852315982919e-0005L,
      236 +        -4.492623672777606053587919463929044226280e-0005L,
      237 +        2.050721258703289487603702670753053765201e-0005L,
      238 +        -9.439485626565616989352750672499008021041e-0006L,
      239 +        4.374838162403994645138200419356844574219e-0006L,
      240 +        -2.038979492862555348577006944451002161496e-0006L,
      241 +        9.536763152382263548086981191378885102802e-0007L,
      242 +        -4.426111214332434049863595231916564014913e-0007L,
      243 +        1.911148847512947464234633846270287546882e-0007L,
      244 +        -5.788673944861923038157839080272303519671e-0008L,
 232  245  };
 233  246  
 234      -static const long double ap2[] = {                      /* 19 terms */
 235      -  -0.077215664901532860606428624449354836087L,
 236      -   3.224670334241132182271948744265855440139e-0001L,
 237      -  -6.735230105319809467356126599005051676203e-0002L,
 238      -   2.058080842778453315716389815213496002588e-0002L,
 239      -  -7.385551028673653323064118422580096222959e-0003L,
 240      -   2.890510330735923572088003424849289006039e-0003L,
 241      -  -1.192753911629952368606185543945790688144e-0003L,
 242      -   5.096695239806718875364547587043220998766e-0004L,
 243      -  -2.231547520600616108991867127392089144886e-0004L,
 244      -   9.945746913898151120612322833059416008973e-0005L,
 245      -  -4.492599307461977003570224943054585729684e-0005L,
 246      -   2.050609891889165453592046505651759999090e-0005L,
 247      -  -9.435329866734193796540515247917165988579e-0006L,
 248      -   4.362267138522223236241016136585565144581e-0006L,
 249      -  -2.008556356653246579300491601497510230557e-0006L,
 250      -   8.961498103387207161105347118042844354395e-0007L,
 251      -  -3.614187228330216282235692806488341157741e-0007L,
 252      -   1.136978988247816860500420915014777753153e-0007L,
 253      -  -2.000532786387196664019286514899782691776e-0008L,
      247 +static const long double ap2[] = {      /* 19 terms */
      248 +        -0.077215664901532860606428624449354836087L,
      249 +        3.224670334241132182271948744265855440139e-0001L,
      250 +        -6.735230105319809467356126599005051676203e-0002L,
      251 +        2.058080842778453315716389815213496002588e-0002L,
      252 +        -7.385551028673653323064118422580096222959e-0003L,
      253 +        2.890510330735923572088003424849289006039e-0003L,
      254 +        -1.192753911629952368606185543945790688144e-0003L,
      255 +        5.096695239806718875364547587043220998766e-0004L,
      256 +        -2.231547520600616108991867127392089144886e-0004L,
      257 +        9.945746913898151120612322833059416008973e-0005L,
      258 +        -4.492599307461977003570224943054585729684e-0005L,
      259 +        2.050609891889165453592046505651759999090e-0005L,
      260 +        -9.435329866734193796540515247917165988579e-0006L,
      261 +        4.362267138522223236241016136585565144581e-0006L,
      262 +        -2.008556356653246579300491601497510230557e-0006L,
      263 +        8.961498103387207161105347118042844354395e-0007L,
      264 +        -3.614187228330216282235692806488341157741e-0007L,
      265 +        1.136978988247816860500420915014777753153e-0007L,
      266 +        -2.000532786387196664019286514899782691776e-0008L,
 254  267  };
 255  268  
 256      -static const long double ap3[] = {                      /* 19 terms */
 257      -  -0.077215664901532859888521470795348856446L,
 258      -   3.224670334241131733364048614484228443077e-0001L,
 259      -  -6.735230105319676541660495145259038151576e-0002L,
 260      -   2.058080842775975461837768839015444273830e-0002L,
 261      -  -7.385551028347615729728618066663566606906e-0003L,
 262      -   2.890510327517954083379032008643080256676e-0003L,
 263      -  -1.192753886919470728001821137439430882603e-0003L,
 264      -   5.096693728898932234814903769146577482912e-0004L,
 265      -  -2.231540055048827662528594010961874258037e-0004L,
 266      -   9.945446210018649311491619999438833843723e-0005L,
 267      -  -4.491608206598064519190236245753867697750e-0005L,
 268      -   2.047939071322271016498065052853746466669e-0005L,
 269      -  -9.376824046522786006677541036631536790762e-0006L,
 270      -   4.259329829498149111582277209189150127347e-0006L,
 271      -  -1.866064770421594266702176289764212873428e-0006L,
 272      -   7.462066721137579592928128104534957135669e-0007L,
 273      -  -2.483546217529077735074007138457678727371e-0007L,
 274      -   5.915166576378161473299324673649144297574e-0008L,
 275      -  -7.334139641706988966966252333759604701905e-0009L,
      269 +static const long double ap3[] = {      /* 19 terms */
      270 +        -0.077215664901532859888521470795348856446L,
      271 +        3.224670334241131733364048614484228443077e-0001L,
      272 +        -6.735230105319676541660495145259038151576e-0002L,
      273 +        2.058080842775975461837768839015444273830e-0002L,
      274 +        -7.385551028347615729728618066663566606906e-0003L,
      275 +        2.890510327517954083379032008643080256676e-0003L,
      276 +        -1.192753886919470728001821137439430882603e-0003L,
      277 +        5.096693728898932234814903769146577482912e-0004L,
      278 +        -2.231540055048827662528594010961874258037e-0004L,
      279 +        9.945446210018649311491619999438833843723e-0005L,
      280 +        -4.491608206598064519190236245753867697750e-0005L,
      281 +        2.047939071322271016498065052853746466669e-0005L,
      282 +        -9.376824046522786006677541036631536790762e-0006L,
      283 +        4.259329829498149111582277209189150127347e-0006L,
      284 +        -1.866064770421594266702176289764212873428e-0006L,
      285 +        7.462066721137579592928128104534957135669e-0007L,
      286 +        -2.483546217529077735074007138457678727371e-0007L,
      287 +        5.915166576378161473299324673649144297574e-0008L,
      288 +        -7.334139641706988966966252333759604701905e-0009L,
 276  289  };
 277  290  
 278      -static const long double ap4[] = {                      /* 19 terms */
 279      -  -0.0772156649015326785569313252637238673675L,
 280      -   3.224670334241051435008842685722468344822e-0001L,
 281      -  -6.735230105302832007479431772160948499254e-0002L,
 282      -   2.058080842553481183648529360967441889912e-0002L,
 283      -  -7.385551007602909242024706804659879199244e-0003L,
 284      -   2.890510182473907253939821312248303471206e-0003L,
 285      -  -1.192753098427856770847894497586825614450e-0003L,
 286      -   5.096659636418811568063339214203693550804e-0004L,
 287      -  -2.231421144004355691166194259675004483639e-0004L,
 288      -   9.942073842343832132754332881883387625136e-0005L,
 289      -  -4.483809261973204531263252655050701205397e-0005L,
 290      -   2.033260142610284888319116654931994447173e-0005L,
 291      -  -9.153539544026646699870528191410440585796e-0006L,
 292      -   3.988460469925482725894144688699584997971e-0006L,
 293      -  -1.609692980087029172567957221850825977621e-0006L,
 294      -   5.634916377249975825399706694496688803488e-0007L,
 295      -  -1.560065465929518563549083208482591437696e-0007L,
 296      -   2.961350193868935325526962209019387821584e-0008L,
 297      -  -2.834602215195368130104649234505033159842e-0009L,
      291 +static const long double ap4[] = {      /* 19 terms */
      292 +        -0.0772156649015326785569313252637238673675L,
      293 +        3.224670334241051435008842685722468344822e-0001L,
      294 +        -6.735230105302832007479431772160948499254e-0002L,
      295 +        2.058080842553481183648529360967441889912e-0002L,
      296 +        -7.385551007602909242024706804659879199244e-0003L,
      297 +        2.890510182473907253939821312248303471206e-0003L,
      298 +        -1.192753098427856770847894497586825614450e-0003L,
      299 +        5.096659636418811568063339214203693550804e-0004L,
      300 +        -2.231421144004355691166194259675004483639e-0004L,
      301 +        9.942073842343832132754332881883387625136e-0005L,
      302 +        -4.483809261973204531263252655050701205397e-0005L,
      303 +        2.033260142610284888319116654931994447173e-0005L,
      304 +        -9.153539544026646699870528191410440585796e-0006L,
      305 +        3.988460469925482725894144688699584997971e-0006L,
      306 +        -1.609692980087029172567957221850825977621e-0006L,
      307 +        5.634916377249975825399706694496688803488e-0007L,
      308 +        -1.560065465929518563549083208482591437696e-0007L,
      309 +        2.961350193868935325526962209019387821584e-0008L,
      310 +        -2.834602215195368130104649234505033159842e-0009L,
 298  311  };
 299  312  
 300  313  static long double
 301      -primary(long double s) {        /* assume |s|<=0.5 */
      314 +primary(long double s)                  /* assume |s|<=0.5 */
      315 +{
 302  316          int i;
 303  317  
 304      -        i = (int) (8.0L * (s + 0.5L));
 305      -        switch(i) {
 306      -        case 0: return ch*s+s*poly(s,an4,21);
 307      -        case 1: return ch*s+s*poly(s,an3,20);
 308      -        case 2: return ch*s+s*poly(s,an2,20);
 309      -        case 3: return ch*s+s*poly(s,an1,20);
 310      -        case 4: return ch*s+s*poly(s,ap1,19);
 311      -        case 5: return ch*s+s*poly(s,ap2,19);
 312      -        case 6: return ch*s+s*poly(s,ap3,19);
 313      -        case 7: return ch*s+s*poly(s,ap4,19);
      318 +        i = (int)(8.0L * (s + 0.5L));
      319 +
      320 +        switch (i) {
      321 +        case 0:
      322 +                return (ch * s + s * poly(s, an4, 21));
      323 +        case 1:
      324 +                return (ch * s + s * poly(s, an3, 20));
      325 +        case 2:
      326 +                return (ch * s + s * poly(s, an2, 20));
      327 +        case 3:
      328 +                return (ch * s + s * poly(s, an1, 20));
      329 +        case 4:
      330 +                return (ch * s + s * poly(s, ap1, 19));
      331 +        case 5:
      332 +                return (ch * s + s * poly(s, ap2, 19));
      333 +        case 6:
      334 +                return (ch * s + s * poly(s, ap3, 19));
      335 +        case 7:
      336 +                return (ch * s + s * poly(s, ap4, 19));
 314  337          }
      338 +
 315  339          /* NOTREACHED */
 316      -    return 0.0L;
      340 +        return (0.0L);
 317  341  }
 318  342  
 319  343  static long double
 320      -poly(long double s, const long double *p, int n) {
      344 +poly(long double s, const long double *p, int n)
      345 +{
 321  346          long double y;
 322  347          int i;
 323      -        y = p[n-1];
 324      -        for (i=n-2;i>=0;i--) y = p[i]+s*y;
 325      -        return y;
      348 +
      349 +        y = p[n - 1];
      350 +
      351 +        for (i = n - 2; i >= 0; i--)
      352 +                y = p[i] + s * y;
      353 +
      354 +        return (y);
 326  355  }
 327  356  
 328  357  static const long double pt[] = {
 329      -   9.189385332046727417803297364056176804663e-0001L,
 330      -   8.333333333333333333333333333331286969123e-0002L,
 331      -  -2.777777777777777777777777553194796036402e-0003L,
 332      -   7.936507936507936507927283071433584248176e-0004L,
 333      -  -5.952380952380952362351042163192634108297e-0004L,
 334      -   8.417508417508395661774286645578379460131e-0004L,
 335      -  -1.917526917525263651186066417934685675649e-0003L,
 336      -   6.410256409395203164659292973142293199083e-0003L,
 337      -  -2.955065327248303301763594514012418438188e-0002L,
 338      -   1.796442830099067542945998615411893822886e-0001L,
 339      -  -1.392413465829723742489974310411118662919e+0000L,
 340      -   1.339984238037267658352656597960492029261e+0001L,
 341      -  -1.564707657605373662425785904278645727813e+0002L,
 342      -   2.156323807499211356127813962223067079300e+0003L,
 343      -  -3.330486427626223184647299834137041307569e+0004L,
 344      -   5.235535072011889213611369254140123518699e+0005L,
 345      -  -7.258160984602220710491988573430212593080e+0006L,
 346      -   7.316526934569686459641438882340322673357e+0007L,
 347      -  -3.806450279064900548836571789284896711473e+0008L,
      358 +        9.189385332046727417803297364056176804663e-0001L,
      359 +        8.333333333333333333333333333331286969123e-0002L,
      360 +        -2.777777777777777777777777553194796036402e-0003L,
      361 +        7.936507936507936507927283071433584248176e-0004L,
      362 +        -5.952380952380952362351042163192634108297e-0004L,
      363 +        8.417508417508395661774286645578379460131e-0004L,
      364 +        -1.917526917525263651186066417934685675649e-0003L,
      365 +        6.410256409395203164659292973142293199083e-0003L,
      366 +        -2.955065327248303301763594514012418438188e-0002L,
      367 +        1.796442830099067542945998615411893822886e-0001L,
      368 +        -1.392413465829723742489974310411118662919e+0000L,
      369 +        1.339984238037267658352656597960492029261e+0001L,
      370 +        -1.564707657605373662425785904278645727813e+0002L,
      371 +        2.156323807499211356127813962223067079300e+0003L,
      372 +        -3.330486427626223184647299834137041307569e+0004L,
      373 +        5.235535072011889213611369254140123518699e+0005L,
      374 +        -7.258160984602220710491988573430212593080e+0006L,
      375 +        7.316526934569686459641438882340322673357e+0007L,
      376 +        -3.806450279064900548836571789284896711473e+0008L,
 348  377  };
 349  378  
 350  379  static long double
 351      -polytail(long double s) {
 352      -        long double t,z;
      380 +polytail(long double s)
      381 +{
      382 +        long double t, z;
 353  383          int i;
 354      -        z = s*s;
      384 +
      385 +        z = s * s;
 355  386          t = pt[18];
 356      -        for (i=17;i>=1;i--) t = pt[i]+z*t;
 357      -        return pt[0]+s*t;
      387 +
      388 +        for (i = 17; i >= 1; i--)
      389 +                t = pt[i] + z * t;
      390 +
      391 +        return (pt[0] + s * t);
 358  392  }
 359  393  
 360  394  static long double
 361      -neg(long double z, int *signgamlp) {
 362      -        long double t,p;
 363      -
 364      -     /*
 365      -      * written by K.C. Ng,  Feb 2, 1989.
 366      -      *
 367      -      * Since
 368      -      *         -z*G(-z)*G(z) = pi/sin(pi*z),
 369      -      * we have
 370      -      *         G(-z) = -pi/(sin(pi*z)*G(z)*z)
 371      -      *               =  pi/(sin(pi*(-z))*G(z)*z)
 372      -      * Algorithm
 373      -      *         z = |z|
 374      -      *         t = sinpi(z); ...note that when z>2**112, z is an int
 375      -      *         and hence t=0.
 376      -      *
 377      -      *         if (t == 0.0) return 1.0/0.0;
 378      -      *         if (t< 0.0) *signgamlp = -1; else t= -t;
 379      -      *         if (z<1.0e-40)  ...tiny z
 380      -      *             return -log(z);
 381      -      *         else
 382      -      *             return log(pi/(t*z))-lgamma(z);
 383      -      *
 384      -      */
      395 +neg(long double z, int *signgamlp)
      396 +{
      397 +        long double t, p;
      398 +
      399 +        /* BEGIN CSTYLED */
      400 +        /*
      401 +         * written by K.C. Ng,  Feb 2, 1989.
      402 +         *
      403 +         * Since
      404 +         *              -z*G(-z)*G(z) = pi/sin(pi*z),
      405 +         * we have
      406 +         *      G(-z) = -pi/(sin(pi*z)*G(z)*z)
      407 +         *            =  pi/(sin(pi*(-z))*G(z)*z)
      408 +         * Algorithm
      409 +         *              z = |z|
      410 +         *              t = sinpi(z); ...note that when z>2**112, z is an int
      411 +         *              and hence t=0.
      412 +         *
      413 +         *              if (t == 0.0) return 1.0/0.0;
      414 +         *              if (t< 0.0) *signgamlp = -1; else t= -t;
      415 +         *              if (z<1.0e-40)  ...tiny z
      416 +         *                  return -log(z);
      417 +         *              else
      418 +         *                  return log(pi/(t*z))-lgamma(z);
      419 +         *
      420 +         */
      421 +        /* END CSTYLED */
 385  422  
 386  423          t = sinpil(z);                  /* t := sin(pi*z) */
 387      -        if (t == c0)                    /* return   1.0/0.0 =  +INF */
 388      -            return c1/c0;
      424 +
      425 +        if (t == c0)                    /* return   1.0/0.0 =  +INF */
      426 +                return (c1 / c0);
 389  427  
 390  428          z = -z;
 391      -        if (z<=tiny)
 392      -            p = -logl(z);
      429 +
      430 +        if (z <= tiny)
      431 +                p = -logl(z);
 393  432          else
 394      -            p = logl(pi/(fabsl(t)*z))-__k_lgammal(z,signgamlp);
 395      -        if (t<c0) *signgamlp = -1;
 396      -        return p;
      433 +                p = logl(pi / (fabsl(t) * z)) - __k_lgammal(z, signgamlp);
      434 +
      435 +        if (t < c0)
      436 +                *signgamlp = -1;
      437 +
      438 +        return (p);
 397  439  }
    
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX