gen/lib.bc

   1 /*
   2  * *****************************************************************************
   3  *
   4  * SPDX-License-Identifier: BSD-2-Clause
   5  *
   6  * Copyright (c) 2018-2024 Gavin D. Howard and contributors.
   7  *
   8  * Redistribution and use in source and binary forms, with or without
   9  * modification, are permitted provided that the following conditions are met:
  10  *
  11  * * Redistributions of source code must retain the above copyright notice, this
  12  *   list of conditions and the following disclaimer.
  13  *
  14  * * Redistributions in binary form must reproduce the above copyright notice,
  15  *   this list of conditions and the following disclaimer in the documentation
  16  *   and/or other materials provided with the distribution.
  17  *
  18  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  19  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  20  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  21  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
  22  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  23  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  24  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  25  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  26  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  27  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  28  * POSSIBILITY OF SUCH DAMAGE.
  29  *
  30  * *****************************************************************************
  31  *
  32  * The bc math library.
  33  *
  34  */
  35
  36 define e(x){
  37         auto b,s,n,r,d,i,p,f,v
  38         b=ibase
  39         ibase=A
  40         if(x<0){
  41                 n=1
  42                 x=-x
  43         }
  44         s=scale
  45         r=6+s+.44*x
  46         scale=scale(x)+1
  47         while(x>1){
  48                 d+=1
  49                 x/=2
  50                 scale+=1
  51         }
  52         scale=r
  53         r=x+1
  54         p=x
  55         f=v=1
  56         for(i=2;v;++i){
  57                 p*=x
  58                 f*=i
  59                 v=p/f
  60                 r+=v
  61         }
  62         while(d--)r*=r
  63         scale=s
  64         ibase=b
  65         if(n)return(1/r)
  66         return(r/1)
  67 }
  68 define l(x){
  69         auto b,s,r,p,a,q,i,v
  70         if(x<=0)return((1-A^scale)/1)
  71         b=ibase
  72         ibase=A
  73         s=scale
  74         scale+=6
  75         p=2
  76         while(x>=2){
  77                 p*=2
  78                 x=sqrt(x)
  79         }
  80         while(x<=.5){
  81                 p*=2
  82                 x=sqrt(x)
  83         }
  84         r=a=(x-1)/(x+1)
  85         q=a*a
  86         v=1
  87         for(i=3;v;i+=2){
  88                 a*=q
  89                 v=a/i
  90                 r+=v
  91         }
  92         r*=p
  93         scale=s
  94         ibase=b
  95         return(r/1)
  96 }
  97 define s(x){
  98         auto b,s,r,a,q,i
  99         if(x<0)return(-s(-x))
 100         b=ibase
 101         ibase=A
 102         s=scale
 103         scale=1.1*s+2
 104         a=a(1)
 105         scale=0
 106         q=(x/a+2)/4
 107         x-=4*q*a
 108         if(q%2)x=-x
 109         scale=s+2
 110         r=a=x
 111         q=-x*x
 112         for(i=3;a;i+=2){
 113                 a*=q/(i*(i-1))
 114                 r+=a
 115         }
 116         scale=s
 117         ibase=b
 118         return(r/1)
 119 }
 120 define c(x){
 121         auto b,s
 122         b=ibase
 123         ibase=A
 124         s=scale
 125         scale*=1.2
 126         x=s(2*a(1)+x)
 127         scale=s
 128         ibase=b
 129         return(x/1)
 130 }
 131 define a(x){
 132         auto b,s,r,n,a,m,t,f,i,u
 133         b=ibase
 134         ibase=A
 135         n=1
 136         if(x<0){
 137                 n=-1
 138                 x=-x
 139         }
 140         if(scale<65){
 141                 if(x==1){
 142                         r=.7853981633974483096156608458198757210492923498437764552437361480/n
 143                         ibase=b
 144                         return(r)
 145                 }
 146                 if(x==.2){
 147                         r=.1973955598498807583700497651947902934475851037878521015176889402/n
 148                         ibase=b
 149                         return(r)
 150                 }
 151         }
 152         s=scale
 153         if(x>.2){
 154                 scale+=5
 155                 a=a(.2)
 156         }
 157         scale=s+3
 158         while(x>.2){
 159                 m+=1
 160                 x=(x-.2)/(1+.2*x)
 161         }
 162         r=u=x
 163         f=-x*x
 164         t=1
 165         for(i=3;t;i+=2){
 166                 u*=f
 167                 t=u/i
 168                 r+=t
 169         }
 170         scale=s
 171         ibase=b
 172         return((m*a+r)/n)
 173 }
 174 define j(n,x){
 175         auto b,s,o,a,i,r,v,f
 176         b=ibase
 177         ibase=A
 178         s=scale
 179         scale=0
 180         n/=1
 181         if(n<0){
 182                 n=-n
 183                 o=n%2
 184         }
 185         a=1
 186         for(i=2;i<=n;++i)a*=i
 187         scale=1.5*s
 188         a=(x^n)/2^n/a
 189         r=v=1
 190         f=-x*x/4
 191         scale+=length(a)-scale(a)
 192         for(i=1;v;++i){
 193                 v=v*f/i/(n+i)
 194                 r+=v
 195         }
 196         scale=s
 197         ibase=b
 198         if(o)a=-a
 199         return(a*r/1)
 200 }