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=A^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 gcd(a,b){ 133 auto g,s 134 if(!b)return a 135 s=scale 136 scale=0 137 a=abs(a)$ 138 b=abs(b)$ 139 if(a<b){ 140 g=a 141 a=b 142 b=g 143 } 144 while(b){ 145 g=a%b 146 a=b 147 b=g 148 } 149 scale=s 150 return a 151} 152define lcm(a,b){ 153 auto r,s 154 if(!a&&!b)return 0 155 s=scale 156 scale=0 157 a=abs(a)$ 158 b=abs(b)$ 159 r=a*b/gcd(a,b) 160 scale=s 161 return r 162} 163define pi(s){ 164 auto t,v 165 if(s==0)return 3 166 s=abs(s)$ 167 t=scale 168 scale=s+1 169 v=4*a(1) 170 scale=t 171 return v@s 172} 173define t(x){ 174 auto s,c 175 l=scale 176 scale+=2 177 s=s(x) 178 c=c(x) 179 scale-=2 180 return s/c 181} 182define a2(y,x){ 183 auto a,p 184 if(!x&&!y)y/x 185 if(x<=0){ 186 p=pi(scale+2) 187 if(y<0)p=-p 188 } 189 if(x==0)a=p/2 190 else{ 191 scale+=2 192 a=a(y/x)+p 193 scale-=2 194 } 195 return a@scale 196} 197define sin(x){return s(x)} 198define cos(x){return c(x)} 199define atan(x){return a(x)} 200define tan(x){return t(x)} 201define atan2(y,x){return a2(y,x)} 202define r2d(x){ 203 auto r,i,s 204 s=scale 205 scale+=5 206 i=ibase 207 ibase=A 208 r=x*180/pi(scale) 209 ibase=i 210 scale=s 211 return r@s 212} 213define d2r(x){ 214 auto r,i,s 215 s=scale 216 scale+=5 217 i=ibase 218 ibase=A 219 r=x*pi(scale)/180 220 ibase=i 221 scale=s 222 return r@s 223} 224define frand(p){ 225 p=abs(p)$ 226 return irand(A^p)>>p 227} 228define ifrand(i,p){return irand(abs(i)$)+frand(p)} 229define srand(x){ 230 if(irand(2))return -x 231 return x 232} 233define brand(){return irand(2)} 234define void output(x,b){ 235 auto c 236 c=obase 237 obase=b 238 x 239 obase=c 240} 241define void hex(x){output(x,G)} 242define void binary(x){output(x,2)} 243define ubytes(x){ 244 auto p,i 245 x=abs(x)$ 246 i=2^8 247 for(p=1;i-1<x;p*=2){i*=i} 248 return p 249} 250define sbytes(x){ 251 auto p,n,z 252 z=(x<0) 253 x=abs(x) 254 x=x$ 255 n=ubytes(x) 256 p=2^(n*8-1) 257 if(x>p||(!z&&x==p))n*=2 258 return n 259} 260define s2un(x,n){ 261 auto t,u,s 262 x=x$ 263 if(x<0){ 264 x=abs(x) 265 s=scale 266 scale=0 267 t=n*8 268 u=2^(t-1) 269 if(x==u)return x 270 else if(x>u)x%=u 271 scale=s 272 return 2^(t)-x 273 } 274 return x 275} 276define s2u(x){return s2un(x,sbytes(x))} 277define void output_byte(x,i){ 278 auto j,p,y,b 279 j=ibase 280 ibase=A 281 s=scale 282 scale=0 283 x=abs(x)$ 284 b=x/(2^(i*8)) 285 b%=256 286 y=log(256,obase) 287 if(b>1)p=log(b,obase)+1 288 else p=b 289 for(i=y-p;i>0;--i)print 0 290 if(b)print b 291 scale=s 292 ibase=j 293} 294define void output_uint(x,n){ 295 auto i 296 for(i=n-1;i>=0;--i){ 297 output_byte(x,i) 298 if(i)print" " 299 else print"\n" 300 } 301} 302define void hex_uint(x,n){ 303 auto o 304 o=obase 305 obase=G 306 output_uint(x,n) 307 obase=o 308} 309define void binary_uint(x,n){ 310 auto o 311 o=obase 312 obase=2 313 output_uint(x,n) 314 obase=o 315} 316define void uintn(x,n){ 317 if(scale(x)){ 318 print"Error: ",x," is not an integer.\n" 319 return 320 } 321 if(x<0){ 322 print"Error: ",x," is negative.\n" 323 return 324 } 325 if(x>=2^(n*8)){ 326 print"Error: ",x," cannot fit into ",n," unsigned byte(s).\n" 327 return 328 } 329 binary_uint(x,n) 330 hex_uint(x,n) 331} 332define void intn(x,n){ 333 auto t 334 if(scale(x)){ 335 print"Error: ",x," is not an integer.\n" 336 return 337 } 338 t=2^(n*8-1) 339 if(abs(x)>=t&&(x>0||x!=-t)){ 340 print "Error: ",x," cannot fit into ",n," signed byte(s).\n" 341 return 342 } 343 x=s2un(x,n) 344 binary_uint(x,n) 345 hex_uint(x,n) 346} 347define void uint8(x){uintn(x,1)} 348define void int8(x){intn(x,1)} 349define void uint16(x){uintn(x,2)} 350define void int16(x){intn(x,2)} 351define void uint32(x){uintn(x,4)} 352define void int32(x){intn(x,4)} 353define void uint64(x){uintn(x,8)} 354define void int64(x){intn(x,8)} 355define void uint(x){uintn(x,ubytes(x))} 356define void int(x){intn(x,sbytes(x))} 357define bunrev(t){ 358 auto a,s,m[] 359 s=scale 360 scale=0 361 t=abs(t)$ 362 while(t!=1){ 363 t=divmod(t,2,m[]) 364 a*=2 365 a+=m[0] 366 } 367 scale=s 368 return a 369} 370define band(a,b){ 371 auto s,t,m[],n[] 372 a=abs(a)$ 373 b=abs(b)$ 374 if(b>a){ 375 t=b 376 b=a 377 a=t 378 } 379 s=scale 380 scale=0 381 t=1 382 while(b){ 383 a=divmod(a,2,m[]) 384 b=divmod(b,2,n[]) 385 t*=2 386 t+=(m[0]&&n[0]) 387 } 388 scale=s 389 return bunrev(t) 390} 391define bor(a,b){ 392 auto s,t,m[],n[] 393 a=abs(a)$ 394 b=abs(b)$ 395 if(b>a){ 396 t=b 397 b=a 398 a=t 399 } 400 s=scale 401 scale=0 402 t=1 403 while(b){ 404 a=divmod(a,2,m[]) 405 b=divmod(b,2,n[]) 406 t*=2 407 t+=(m[0]||n[0]) 408 } 409 while(a){ 410 a=divmod(a,2,m[]) 411 t*=2 412 t+=m[0] 413 } 414 scale=s 415 return bunrev(t) 416} 417define bxor(a,b){ 418 auto s,t,m[],n[] 419 a=abs(a)$ 420 b=abs(b)$ 421 if(b>a){ 422 t=b 423 b=a 424 a=t 425 } 426 s=scale 427 scale=0 428 t=1 429 while(b){ 430 a=divmod(a,2,m[]) 431 b=divmod(b,2,n[]) 432 t*=2 433 t+=(m[0]+n[0]==1) 434 } 435 while(a){ 436 a=divmod(a,2,m[]) 437 t*=2 438 t+=m[0] 439 } 440 scale=s 441 return bunrev(t) 442} 443define bshl(a,b){return abs(a)$*2^abs(b)$} 444define bshr(a,b){return (abs(a)$/2^abs(b)$)$} 445define bnotn(x,n){ 446 auto s,t,m[] 447 s=scale 448 scale=0 449 t=2^(abs(n)$*8) 450 x=abs(x)$%t+t 451 t=1 452 while(x!=1){ 453 x=divmod(x,2,m[]) 454 t*=2 455 t+=!m[0] 456 } 457 scale=s 458 return bunrev(t) 459} 460define bnot8(x){return bnotn(x,1)} 461define bnot16(x){return bnotn(x,2)} 462define bnot32(x){return bnotn(x,4)} 463define bnot64(x){return bnotn(x,8)} 464define bnot(x){return bnotn(x,ubytes(x))} 465define brevn(x,n){ 466 auto s,t,m[] 467 s=scale 468 scale=0 469 t=2^(abs(n)$*8) 470 x=abs(x)$%t+t 471 scale=s 472 return bunrev(x) 473} 474define brev8(x){return brevn(x,1)} 475define brev16(x){return brevn(x,2)} 476define brev32(x){return brevn(x,4)} 477define brev64(x){return brevn(x,8)} 478define brev(x){return brevn(x,ubytes(x))} 479define broln(x,p,n){ 480 auto s,t,m[] 481 s=scale 482 scale=0 483 n=abs(n)$*8 484 p=abs(p)$%n 485 t=2^n 486 x=abs(x)$%t 487 if(!p)return x 488 x=divmod(x,2^(n-p),m[]) 489 x+=m[0]*2^p%t 490 scale=s 491 return x 492} 493define brol8(x,p){return broln(x,p,1)} 494define brol16(x,p){return broln(x,p,2)} 495define brol32(x,p){return broln(x,p,4)} 496define brol64(x,p){return broln(x,p,8)} 497define brol(x,p){return broln(x,p,ubytes(x))} 498define brorn(x,p,n){ 499 auto s,t,m[] 500 s=scale 501 scale=0 502 n=abs(n)$*8 503 p=abs(p)$%n 504 t=2^n 505 x=abs(x)$%t 506 if(!p)return x 507 x=divmod(x,2^p,m[]) 508 x+=m[0]*2^(n-p)%t 509 scale=s 510 return x 511} 512define bror8(x,p){return brorn(x,p,1)} 513define bror16(x,p){return brorn(x,p,2)} 514define bror32(x,p){return brorn(x,p,4)} 515define bror64(x,p){return brorn(x,p,8)} 516define brol(x,p){return brorn(x,p,ubytes(x))} 517define bmodn(x,n){ 518 auto s 519 s=scale 520 scale=0 521 x=abs(x)$%2^(abs(n)$*8) 522 scale=s 523 return x 524} 525define bmod8(x){return bmodn(x,1)} 526define bmod16(x){return bmodn(x,2)} 527define bmod32(x){return bmodn(x,4)} 528define bmod64(x){return bmodn(x,8)} 529