xref: /freebsd/contrib/bc/gen/lib.bc (revision d11f81afd5a4a71d5f725950b0592ca212084780)
1/*
2 * *****************************************************************************
3 *
4 * SPDX-License-Identifier: BSD-2-Clause
5 *
6 * Copyright (c) 2018-2021 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
36scale=2*A
37define 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}
69define 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}
98define 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}
121define 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}
132define 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}
175define 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}
202