xref: /freebsd/contrib/bc/gen/lib2.bc (revision 59c8e88e72633afbc47a4ace0d2170d00d51f7dc)
1/*
2 * *****************************************************************************
3 *
4 * SPDX-License-Identifier: BSD-2-Clause
5 *
6 * Copyright (c) 2018-2023 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 second bc math library.
33 *
34 */
35
36define p(x,y){
37	auto a,i,s,z
38	if(y==0)return 1@scale
39	if(x==0){
40		if(y>0)return 0
41		return 1/0
42	}
43	a=y$
44	if(y==a)return(x^a)@scale
45	z=0
46	if(x<1){
47		y=-y
48		a=-a
49		z=x
50		x=1/x
51	}
52	if(y<0){
53		return e(y*l(x))
54	}
55	i=x^a
56	s=scale
57	scale+=length(i)+5
58	if(z){
59		x=1/z
60		i=x^a
61	}
62	i*=e((y-a)*l(x))
63	scale=s
64	return i@scale
65}
66define r(x,p){
67	auto t,n
68	if(x==0)return x
69	p=abs(p)$
70	n=(x<0)
71	x=abs(x)
72	t=x@p
73	if(p<scale(x)&&x-t>=5>>p+1)t+=1>>p
74	if(n)t=-t
75	return t
76}
77define ceil(x,p){
78	auto t,n
79	if(x==0)return x
80	p=abs(p)$
81	n=(x<0)
82	x=abs(x)
83	t=(x+((x@p<x)>>p))@p
84	if(n)t=-t
85	return t
86}
87define f(n){
88	auto r
89	n=abs(n)$
90	for(r=1;n>1;--n)r*=n
91	return r
92}
93define max(a,b){
94	if(a>b)return a
95	return b
96}
97define min(a,b){
98	if(a<b)return a
99	return b
100}
101define perm(n,k){
102	auto f,g,s
103	if(k>n)return 0
104	n=abs(n)$
105	k=abs(k)$
106	f=f(n)
107	g=f(n-k)
108	s=scale
109	scale=0
110	f/=g
111	scale=s
112	return f
113}
114define comb(n,r){
115	auto s,f,g,h
116	if(r>n)return 0
117	n=abs(n)$
118	r=abs(r)$
119	s=scale
120	scale=0
121	f=f(n)
122	h=f(r)
123	g=f(n-r)
124	f/=h*g
125	scale=s
126	return f
127}
128define fib(n){
129	auto i,t,p,r
130	if(!n)return 0
131	n=abs(n)$
132	t=1
133	for (i=1;i<n;++i){
134		r=p
135		p=t
136		t+=r
137	}
138	return t
139}
140define log(x,b){
141	auto p,s
142	s=scale
143	if(scale<K)scale=K
144	if(scale(x)>scale)scale=scale(x)
145	scale*=2
146	p=l(x)/l(b)
147	scale=s
148	return p@s
149}
150define l2(x){return log(x,2)}
151define l10(x){return log(x,A)}
152define root(x,n){
153	auto s,t,m,r,q,p
154	if(n<0)sqrt(n)
155	n=n$
156	if(n==0)x/n
157	if(x==0||n==1)return x
158	if(n==2)return sqrt(x)
159	s=scale
160	scale=0
161	if(x<0&&n%2==0){
162		scale=s
163		sqrt(x)
164	}
165	scale=s+scale(x)+5
166	t=s+5
167	m=(x<0)
168	x=abs(x)
169	p=n-1
170	q=A^ceil((length(x$)/n)$,0)
171	while(r@t!=q@t){
172		r=q
173		q=(p*r+x/r^p)/n
174	}
175	if(m)r=-r
176	scale=s
177	return r@s
178}
179define cbrt(x){return root(x,3)}
180define gcd(a,b){
181	auto g,s
182	if(!b)return a
183	s=scale
184	scale=0
185	a=abs(a)$
186	b=abs(b)$
187	if(a<b){
188		g=a
189		a=b
190		b=g
191	}
192	while(b){
193		g=a%b
194		a=b
195		b=g
196	}
197	scale=s
198	return a
199}
200define lcm(a,b){
201	auto r,s
202	if(!a&&!b)return 0
203	s=scale
204	scale=0
205	a=abs(a)$
206	b=abs(b)$
207	r=a*b/gcd(a,b)
208	scale=s
209	return r
210}
211define pi(s){
212	auto t,v
213	if(s==0)return 3
214	s=abs(s)$
215	t=scale
216	scale=s+1
217	v=4*a(1)
218	scale=t
219	return v@s
220}
221define t(x){
222	auto s,c
223	l=scale
224	scale+=2
225	s=s(x)
226	c=c(x)
227	scale-=2
228	return s/c
229}
230define a2(y,x){
231	auto a,p
232	if(!x&&!y)y/x
233	if(x<=0){
234		p=pi(scale+2)
235		if(y<0)p=-p
236	}
237	if(x==0)a=p/2
238	else{
239		scale+=2
240		a=a(y/x)+p
241		scale-=2
242	}
243	return a@scale
244}
245define sin(x){return s(x)}
246define cos(x){return c(x)}
247define atan(x){return a(x)}
248define tan(x){return t(x)}
249define atan2(y,x){return a2(y,x)}
250define r2d(x){
251	auto r,i,s
252	s=scale
253	scale+=5
254	i=ibase
255	ibase=A
256	r=x*180/pi(scale)
257	ibase=i
258	scale=s
259	return r@s
260}
261define d2r(x){
262	auto r,i,s
263	s=scale
264	scale+=5
265	i=ibase
266	ibase=A
267	r=x*pi(scale)/180
268	ibase=i
269	scale=s
270	return r@s
271}
272define frand(p){
273	p=abs(p)$
274	return irand(A^p)>>p
275}
276define ifrand(i,p){return irand(abs(i)$)+frand(p)}
277define i2rand(a,b){
278	auto n,x
279	a=a$
280	b=b$
281	if(a==b)return a
282	n=min(a,b)
283	x=max(a,b)
284	return irand(x-n+1)+n
285}
286define srand(x){
287	if(irand(2))return -x
288	return x
289}
290define brand(){return irand(2)}
291define void output(x,b){
292	auto c
293	c=obase
294	obase=b
295	x
296	obase=c
297}
298define void hex(x){output(x,G)}
299define void binary(x){output(x,2)}
300define ubytes(x){
301	auto p,i
302	x=abs(x)$
303	i=2^8
304	for(p=1;i-1<x;p*=2){i*=i}
305	return p
306}
307define sbytes(x){
308	auto p,n,z
309	z=(x<0)
310	x=abs(x)$
311	n=ubytes(x)
312	p=2^(n*8-1)
313	if(x>p||(!z&&x==p))n*=2
314	return n
315}
316define s2un(x,n){
317	auto t,u,s
318	x=x$
319	if(x<0){
320		x=abs(x)
321		s=scale
322		scale=0
323		t=n*8
324		u=2^(t-1)
325		if(x==u)return x
326		else if(x>u)x%=u
327		scale=s
328		return 2^(t)-x
329	}
330	return x
331}
332define s2u(x){return s2un(x,sbytes(x))}
333define void plz(x){
334	if(leading_zero())print x
335	else{
336		if(x>-1&&x<1&&x!=0){
337			if(x<0)print"-"
338			print 0,abs(x)
339		}
340		else print x
341	}
342}
343define void plznl(x){
344	plz(x)
345	print"\n"
346}
347define void pnlz(x){
348	auto s,i
349	if(leading_zero()){
350		if(x>-1&&x<1&&x!=0){
351			s=scale(x)
352			if(x<0)print"-"
353			print"."
354			x=abs(x)
355			for(i=0;i<s;++i){
356				x<<=1
357				print x$
358				x-=x$
359			}
360			return
361		}
362	}
363	print x
364}
365define void pnlznl(x){
366	pnlz(x)
367	print"\n"
368}
369define void output_byte(x,i){
370	auto j,p,y,b,s
371	s=scale
372	scale=0
373	x=abs(x)$
374	b=x/(2^(i*8))
375	j=2^8
376	b%=j
377	y=log(j,obase)
378	if(b>1)p=log(b,obase)+1
379	else p=b
380	for(i=y-p;i>0;--i)print 0
381	if(b)print b
382	scale=s
383}
384define void output_uint(x,n){
385	auto i
386	for(i=n-1;i>=0;--i){
387		output_byte(x,i)
388		if(i)print" "
389		else print"\n"
390	}
391}
392define void hex_uint(x,n){
393	auto o
394	o=obase
395	obase=G
396	output_uint(x,n)
397	obase=o
398}
399define void binary_uint(x,n){
400	auto o
401	o=obase
402	obase=2
403	output_uint(x,n)
404	obase=o
405}
406define void uintn(x,n){
407	if(scale(x)){
408		print"Error: ",x," is not an integer.\n"
409		return
410	}
411	if(x<0){
412		print"Error: ",x," is negative.\n"
413		return
414	}
415	if(x>=2^(n*8)){
416		print"Error: ",x," cannot fit into ",n," unsigned byte(s).\n"
417		return
418	}
419	binary_uint(x,n)
420	hex_uint(x,n)
421}
422define void intn(x,n){
423	auto t
424	if(scale(x)){
425		print"Error: ",x," is not an integer.\n"
426		return
427	}
428	t=2^(n*8-1)
429	if(abs(x)>=t&&(x>0||x!=-t)){
430		print "Error: ",x," cannot fit into ",n," signed byte(s).\n"
431		return
432	}
433	x=s2un(x,n)
434	binary_uint(x,n)
435	hex_uint(x,n)
436}
437define void uint8(x){uintn(x,1)}
438define void int8(x){intn(x,1)}
439define void uint16(x){uintn(x,2)}
440define void int16(x){intn(x,2)}
441define void uint32(x){uintn(x,4)}
442define void int32(x){intn(x,4)}
443define void uint64(x){uintn(x,8)}
444define void int64(x){intn(x,8)}
445define void uint(x){uintn(x,ubytes(x))}
446define void int(x){intn(x,sbytes(x))}
447define bunrev(t){
448	auto a,s,m[]
449	s=scale
450	scale=0
451	t=abs(t)$
452	while(t!=1){
453		t=divmod(t,2,m[])
454		a*=2
455		a+=m[0]
456	}
457	scale=s
458	return a
459}
460define band(a,b){
461	auto s,t,m[],n[]
462	a=abs(a)$
463	b=abs(b)$
464	if(b>a){
465		t=b
466		b=a
467		a=t
468	}
469	s=scale
470	scale=0
471	t=1
472	while(b){
473		a=divmod(a,2,m[])
474		b=divmod(b,2,n[])
475		t*=2
476		t+=(m[0]&&n[0])
477	}
478	scale=s
479	return bunrev(t)
480}
481define bor(a,b){
482	auto s,t,m[],n[]
483	a=abs(a)$
484	b=abs(b)$
485	if(b>a){
486		t=b
487		b=a
488		a=t
489	}
490	s=scale
491	scale=0
492	t=1
493	while(b){
494		a=divmod(a,2,m[])
495		b=divmod(b,2,n[])
496		t*=2
497		t+=(m[0]||n[0])
498	}
499	while(a){
500		a=divmod(a,2,m[])
501		t*=2
502		t+=m[0]
503	}
504	scale=s
505	return bunrev(t)
506}
507define bxor(a,b){
508	auto s,t,m[],n[]
509	a=abs(a)$
510	b=abs(b)$
511	if(b>a){
512		t=b
513		b=a
514		a=t
515	}
516	s=scale
517	scale=0
518	t=1
519	while(b){
520		a=divmod(a,2,m[])
521		b=divmod(b,2,n[])
522		t*=2
523		t+=(m[0]+n[0]==1)
524	}
525	while(a){
526		a=divmod(a,2,m[])
527		t*=2
528		t+=m[0]
529	}
530	scale=s
531	return bunrev(t)
532}
533define bshl(a,b){return abs(a)$*2^abs(b)$}
534define bshr(a,b){return(abs(a)$/2^abs(b)$)$}
535define bnotn(x,n){
536	auto s,t,m[]
537	s=scale
538	scale=0
539	t=2^(abs(n)$*8)
540	x=abs(x)$%t+t
541	t=1
542	while(x!=1){
543		x=divmod(x,2,m[])
544		t*=2
545		t+=!m[0]
546	}
547	scale=s
548	return bunrev(t)
549}
550define bnot8(x){return bnotn(x,1)}
551define bnot16(x){return bnotn(x,2)}
552define bnot32(x){return bnotn(x,4)}
553define bnot64(x){return bnotn(x,8)}
554define bnot(x){return bnotn(x,ubytes(x))}
555define brevn(x,n){
556	auto s,t,m[]
557	s=scale
558	scale=0
559	t=2^(abs(n)$*8)
560	x=abs(x)$%t+t
561	scale=s
562	return bunrev(x)
563}
564define brev8(x){return brevn(x,1)}
565define brev16(x){return brevn(x,2)}
566define brev32(x){return brevn(x,4)}
567define brev64(x){return brevn(x,8)}
568define brev(x){return brevn(x,ubytes(x))}
569define broln(x,p,n){
570	auto s,t,m[]
571	s=scale
572	scale=0
573	n=abs(n)$*8
574	p=abs(p)$%n
575	t=2^n
576	x=abs(x)$%t
577	if(!p)return x
578	x=divmod(x,2^(n-p),m[])
579	x+=m[0]*2^p%t
580	scale=s
581	return x
582}
583define brol8(x,p){return broln(x,p,1)}
584define brol16(x,p){return broln(x,p,2)}
585define brol32(x,p){return broln(x,p,4)}
586define brol64(x,p){return broln(x,p,8)}
587define brol(x,p){return broln(x,p,ubytes(x))}
588define brorn(x,p,n){
589	auto s,t,m[]
590	s=scale
591	scale=0
592	n=abs(n)$*8
593	p=abs(p)$%n
594	t=2^n
595	x=abs(x)$%t
596	if(!p)return x
597	x=divmod(x,2^p,m[])
598	x+=m[0]*2^(n-p)%t
599	scale=s
600	return x
601}
602define bror8(x,p){return brorn(x,p,1)}
603define bror16(x,p){return brorn(x,p,2)}
604define bror32(x,p){return brorn(x,p,4)}
605define bror64(x,p){return brorn(x,p,8)}
606define brol(x,p){return brorn(x,p,ubytes(x))}
607define bmodn(x,n){
608	auto s
609	s=scale
610	scale=0
611	x=abs(x)$%2^(abs(n)$*8)
612	scale=s
613	return x
614}
615define bmod8(x){return bmodn(x,1)}
616define bmod16(x){return bmodn(x,2)}
617define bmod32(x){return bmodn(x,4)}
618define bmod64(x){return bmodn(x,8)}
619