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