gen/lib.bc

   1 /*
   2  * *****************************************************************************
   3  *
   4  * SPDX-License-Identifier: BSD-2-Clause
   5  *
   6  * Copyright (c) 2018-2020 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 scale=20
  37 define e(x){
  38         auto b,s,n,r,d,i,p,f,v
  39         b=ibase
  40         ibase=A
  41         if(x<0){
  42                 n=1
  43                 x=-x
  44         }
  45         s=scale
  46         r=6+s+.44*x
  47         scale=scale(x)+1
  48         while(x>1){
  49                 d+=1
  50                 x/=2
  51                 scale+=1
  52         }
  53         scale=r
  54         r=x+1
  55         p=x
  56         f=v=1
  57         for(i=2;v;++i){
  58                 p*=x
  59                 f*=i
  60                 v=p/f
  61                 r+=v
  62         }
  63         while(d--)r*=r
  64         scale=s
  65         ibase=b
  66         if(n)return(1/r)
  67         return(r/1)
  68 }
  69 define l(x){
  70         auto b,s,r,p,a,q,i,v
  71         if(x<=0)return((1-A^scale)/1)
  72         b=ibase
  73         ibase=A
  74         s=scale
  75         scale+=6
  76         p=2
  77         while(x>=2){
  78                 p*=2
  79                 x=sqrt(x)
  80         }
  81         while(x<=.5){
  82                 p*=2
  83                 x=sqrt(x)
  84         }
  85         r=a=(x-1)/(x+1)
  86         q=a*a
  87         v=1
  88         for(i=3;v;i+=2){
  89                 a*=q
  90                 v=a/i
  91                 r+=v
  92         }
  93         r*=p
  94         scale=s
  95         ibase=b
  96         return(r/1)
  97 }
  98 define s(x){
  99         auto b,s,r,a,q,i
 100         if(x<0)return(-s(-x))
 101         b=ibase
 102         ibase=A
 103         s=scale
 104         scale=1.1*s+2
 105         a=a(1)
 106         scale=0
 107         q=(x/a+2)/4
 108         x-=4*q*a
 109         if(q%2)x=-x
 110         scale=s+2
 111         r=a=x
 112         q=-x*x
 113         for(i=3;a;i+=2){
 114                 a*=q/(i*(i-1))
 115                 r+=a
 116         }
 117         scale=s
 118         ibase=b
 119         return(r/1)
 120 }
 121 define c(x){
 122         auto b,s
 123         b=ibase
 124         ibase=A
 125         s=scale
 126         scale*=1.2
 127         x=s(2*a(1)+x)
 128         scale=s
 129         ibase=b
 130         return(x/1)
 131 }
 132 define a(x){
 133         auto b,s,r,n,a,m,t,f,i,u
 134         b=ibase
 135         ibase=A
 136         n=1
 137         if(x<0){
 138                 n=-1
 139                 x=-x
 140         }
 141         if(scale<65){
 142                 if(x==1){
 143                         r=.7853981633974483096156608458198757210492923498437764552437361480/n
 144                         ibase=b
 145                         return(r)
 146                 }
 147                 if(x==.2){
 148                         r=.1973955598498807583700497651947902934475851037878521015176889402/n
 149                         ibase=b
 150                         return(r)
 151                 }
 152         }
 153         s=scale
 154         if(x>.2){
 155                 scale+=5
 156                 a=a(.2)
 157         }
 158         scale=s+3
 159         while(x>.2){
 160                 m+=1
 161                 x=(x-.2)/(1+.2*x)
 162         }
 163         r=u=x
 164         f=-x*x
 165         t=1
 166         for(i=3;t;i+=2){
 167                 u*=f
 168                 t=u/i
 169                 r+=t
 170         }
 171         scale=s
 172         ibase=b
 173         return((m*a+r)/n)
 174 }
 175 define j(n,x){
 176         auto b,s,o,a,i,r,v,f
 177         b=ibase
 178         ibase=A
 179         s=scale
 180         scale=0
 181         n/=1
 182         if(n<0){
 183                 n=-n
 184                 o=n%2
 185         }
 186         a=1
 187         for(i=2;i<=n;++i)a*=i
 188         scale=1.5*s
 189         a=(x^n)/2^n/a
 190         r=v=1
 191         f=-x*x/4
 192         scale+=length(a)-scale(a)
 193         for(i=1;v;++i){
 194                 v=v*f/i/(n+i)
 195                 r+=v
 196         }
 197         scale=s
 198         ibase=b
 199         if(o)a=-a
 200         return(a*r/1)
 201 }