blob: 37893792ac6f077e69fb297eedbc9325416dde8d [file] [log] [blame]
/*
* *****************************************************************************
*
* Copyright 2018 Gavin D. Howard
*
* Permission to use, copy, modify, and/or distribute this software for any
* purpose with or without fee is hereby granted.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH
* REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY
* AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,
* INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM
* LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR
* OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
* PERFORMANCE OF THIS SOFTWARE.
*
* *****************************************************************************
*
* The bc math library.
*
*/
scale=20
define e(x){
auto b,s,n,r,d,i,p,f,v
b=ibase
ibase=A
if(x<0){
n=1
x=-x
}
s=scale
r=6+s+0.44*x
scale=scale(x)+1
while(x>1){
d+=1
x/=2
scale+=1
}
scale=r
r=x+1
p=x
f=v=1
for(i=2;v!=0;++i){
p*=x
f*=i
v=p/f
r+=v
}
while((d--)!=0)r*=r
scale=s
ibase=b
if(n!=0)return(1/r)
return(r/1)
}
define l(x){
auto b,s,r,p,a,q,i,v
b=ibase
ibase=A
if(x<=0){
r=(1-10^scale)/1
ibase=b
return(r)
}
s=scale
scale+=6
p=2
while(x>=2){
p*=2
x=sqrt(x)
}
while(x<=0.5){
p*=2
x=sqrt(x)
}
r=a=(x-1)/(x+1)
q=a*a
v=1
for(i=3;v!=0;i+=2){
a*=q
v=a/i
r+=v
}
r*=p
scale=s
ibase=b
return(r/1)
}
define s(x){
auto b,s,r,n,a,q,i
b=ibase
ibase=A
s=scale
scale=1.1*s+2
a=a(1)
if(x<0){
n=1
x=-x
}
scale=0
q=(x/a+2)/4
x=x-4*q*a
if(q%2!=0)x=-x
scale=s+2
r=a=x
q=-x*x
for(i=3;a!=0;i+=2){
a*=q/(i*(i-1))
r+=a
}
scale=s
ibase=b
if(n!=0)return(-r/1)
return(r/1)
}
define c(x){
auto b,s
b=ibase
ibase=A
s=scale
scale*=1.2
x=s(2*a(1)+x)
scale=s
ibase=b
return(x/1)
}
define a(x){
auto b,s,r,n,a,m,t,f,i,u
b=ibase
ibase=A
n=1
if(x<0){
n=-1
x=-x
}
if(x==1){
if(scale<=64){
return(.7853981633974483096156608458198757210492923498437764552437361480/n)
}
}
if(x==.2){
if(scale<=64){
return(.1973955598498807583700497651947902934475851037878521015176889402/n)
}
}
s=scale
if(x>.2){
scale+=5
a=a(.2)
}
scale=s+3
while(x>.2){
m+=1
x=(x-.2)/(1+.2*x)
}
r=u=x
f=-x*x
t=1
for(i=3;t!=0;i+=2){
u*=f
t=u/i
r+=t
}
scale=s
ibase=b
return((m*a+r)/n)
}
define j(n,x){
auto b,s,o,a,i,v,f
b=ibase
ibase=A
s=scale
scale=0
n/=1
if(n<0){
n=-n
if(n%2==1)o=1
}
a=1
for(i=2;i<=n;++i)a*=i
scale=1.5*s
a=(x^n)/2^n/a
r=v=1
f=-x*x/4
scale=scale+length(a)-scale(a)
for(i=1;v!=0;++i){
v=v*f/i/(n+i)
r+=v
}
scale=s
ibase=b
if(o!=0)return(-a*r/1)
return(a*r/1)
}