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 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 log(x,b){ 97 auto p,s 98 s=scale 99 if(scale<K)scale=K 100 if(scale(x)>scale)scale=scale(x) 101 scale*=2 102 p=l(x)/l(b) 103 scale=s 104 return p@s 105} 106define l2(x){return log(x,2)} 107define l10(x){return log(x,A)} 108define root(x,n){ 109 auto s,m,r,q,p 110 if(n<0)sqrt(n) 111 n=n$ 112 if(n==0)x/n 113 if(x==0||n==1)return x 114 if(n==2)return sqrt(x) 115 s=scale 116 scale=0 117 if(x<0&&n%2==0)sqrt(x) 118 scale=s+2 119 m=(x<0) 120 x=abs(x) 121 p=n-1 122 q=10^ceil((length(x$)/n)$,0) 123 while(r!=q){ 124 r=q 125 q=(p*r+x/r^p)/n 126 } 127 if(m)r=-r 128 scale=s 129 return r@s 130} 131define cbrt(x){return root(x,3)} 132define pi(s){ 133 auto t,v 134 if(s==0)return 3 135 s=abs(s)$ 136 t=scale 137 scale=s+1 138 v=4*a(1) 139 scale=t 140 return v@s 141} 142define t(x){ 143 auto s,c,l 144 l=scale 145 scale+=2 146 s=s(x) 147 c=c(x) 148 scale=l 149 return s/c 150} 151define a2(y,x){ 152 auto a,p 153 if(!x&&!y)y/x 154 if(x<=0){ 155 p=pi(scale+2) 156 if(y<0)p=-p 157 } 158 if(x==0)a=p/2 159 else{ 160 scale+=2 161 a=a(y/x)+p 162 scale-=2 163 } 164 return a@scale 165} 166define sin(x){return s(x)} 167define cos(x){return c(x)} 168define atan(x){return a(x)} 169define tan(x){return t(x)} 170define atan2(y,x){return a2(y,x)} 171define r2d(x){ 172 auto r,i,s 173 s=scale 174 scale+=5 175 i=ibase 176 ibase=A 177 r=x*180/pi(scale) 178 ibase=i 179 scale=s 180 return r@s 181} 182define d2r(x){ 183 auto r,i,s 184 s=scale 185 scale+=5 186 i=ibase 187 ibase=A 188 r=x*pi(scale)/180 189 ibase=i 190 scale=s 191 return r@s 192} 193define frand(p){ 194 p=abs(p)$ 195 return irand(10^p)>>p 196} 197define ifrand(i,p){return irand(abs(i)$)+frand(p)} 198define srand(x){ 199 if(irand(2))return -x 200 return x 201} 202define brand(){return irand(2)} 203define void output(x,b){ 204 auto c 205 c=obase 206 obase=b 207 x 208 obase=c 209} 210define void hex(x){output(x,G)} 211define void binary(x){output(x,2)} 212define ubytes(x){ 213 auto p,b,i 214 b=ibase 215 ibase=A 216 x=abs(x)$ 217 i=2^8 218 for(p=1;i-1<x;p*=2){i*=i} 219 ibase=b 220 return p 221} 222define sbytes(x){ 223 auto p,b,n,z 224 z=(x<0) 225 x=abs(x) 226 x=x$ 227 n=ubytes(x) 228 b=ibase 229 ibase=A 230 p=2^(n*8-1) 231 if(x>p||(!z&&x==p))n*=2 232 ibase=b 233 return n 234} 235define void output_byte(x,i){ 236 auto j,p,y,b 237 j=ibase 238 ibase=A 239 s=scale 240 scale=0 241 x=abs(x)$ 242 b=x/(2^(i*8)) 243 b%=2^8 244 y=log(256,obase) 245 if(b>1)p=log(b,obase)+1 246 else p=b 247 for(i=y-p;i>0;--i)print 0 248 if(b)print b 249 scale=s 250 ibase=j 251} 252define void output_uint(x,n){ 253 auto i,b 254 b=ibase 255 ibase=A 256 for(i=n-1;i>=0;--i){ 257 output_byte(x,i) 258 if(i)print" " 259 else print"\n" 260 } 261 ibase=b 262} 263define void hex_uint(x,n){ 264 auto o 265 o=obase 266 obase=G 267 output_uint(x,n) 268 obase=o 269} 270define void binary_uint(x,n){ 271 auto o 272 o=obase 273 obase=2 274 output_uint(x,n) 275 obase=o 276} 277define void uintn(x,n){ 278 if(scale(x)){ 279 print"Error: ",x," is not an integer.\n" 280 return 281 } 282 if(x<0){ 283 print"Error: ",x," is negative.\n" 284 return 285 } 286 if(x>=2^(n*8)){ 287 print"Error: ",x," cannot fit into ",n," unsigned byte(s).\n" 288 return 289 } 290 binary_uint(x,n) 291 hex_uint(x,n) 292} 293define void intn(x,n){ 294 auto t 295 if(scale(x)){ 296 print"Error: ",x," is not an integer.\n" 297 return 298 } 299 t=2^(n*8-1) 300 if(abs(x)>=t&&(x>0||x!=-t)){ 301 print "Error: ",x," cannot fit into ",n," signed byte(s).\n" 302 return 303 } 304 if(x<0)x=2^(n*8)-(-x) 305 binary_uint(x,n) 306 hex_uint(x,n) 307} 308define void uint8(x){uintn(x,1)} 309define void int8(x){intn(x,1)} 310define void uint16(x){uintn(x,2)} 311define void int16(x){intn(x,2)} 312define void uint32(x){uintn(x,4)} 313define void int32(x){intn(x,4)} 314define void uint64(x){uintn(x,8)} 315define void int64(x){intn(x,8)} 316define void uint(x){uintn(x,ubytes(x))} 317define void int(x){intn(x,sbytes(x))} 318