xref: /freebsd/contrib/bc/gen/lib.bc (revision 9bc300465e48e19d794d88d0c158a2adb92c7197)
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
36define 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}
68define 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}
97define 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}
120define 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}
131define 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}
174define 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}
201