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 plz(x){ 278 if(leading_zero())print x 279 else{ 280 if(x>-1&&x<1&&x!=0){ 281 if(x<0)print"-" 282 print 0,abs(x) 283 } 284 else print x 285 } 286} 287define void plznl(x){ 288 plz(x) 289 print"\n" 290} 291define void pnlz(x){ 292 auto s,i 293 if(leading_zero()){ 294 if(x>-1&&x<1&&x!=0){ 295 s=scale(x) 296 if(x<0)print"-" 297 print"." 298 x=abs(x) 299 for(i=0;i<s;++i){ 300 x<<=1 301 print x$ 302 x-=x$ 303 } 304 return 305 } 306 } 307 print x 308} 309define void pnlznl(x){ 310 pnlz(x) 311 print"\n" 312} 313define void output_byte(x,i){ 314 auto j,p,y,b 315 j=ibase 316 ibase=A 317 s=scale 318 scale=0 319 x=abs(x)$ 320 b=x/(2^(i*8)) 321 b%=256 322 y=log(256,obase) 323 if(b>1)p=log(b,obase)+1 324 else p=b 325 for(i=y-p;i>0;--i)print 0 326 if(b)print b 327 scale=s 328 ibase=j 329} 330define void output_uint(x,n){ 331 auto i 332 for(i=n-1;i>=0;--i){ 333 output_byte(x,i) 334 if(i)print" " 335 else print"\n" 336 } 337} 338define void hex_uint(x,n){ 339 auto o 340 o=obase 341 obase=G 342 output_uint(x,n) 343 obase=o 344} 345define void binary_uint(x,n){ 346 auto o 347 o=obase 348 obase=2 349 output_uint(x,n) 350 obase=o 351} 352define void uintn(x,n){ 353 if(scale(x)){ 354 print"Error: ",x," is not an integer.\n" 355 return 356 } 357 if(x<0){ 358 print"Error: ",x," is negative.\n" 359 return 360 } 361 if(x>=2^(n*8)){ 362 print"Error: ",x," cannot fit into ",n," unsigned byte(s).\n" 363 return 364 } 365 binary_uint(x,n) 366 hex_uint(x,n) 367} 368define void intn(x,n){ 369 auto t 370 if(scale(x)){ 371 print"Error: ",x," is not an integer.\n" 372 return 373 } 374 t=2^(n*8-1) 375 if(abs(x)>=t&&(x>0||x!=-t)){ 376 print "Error: ",x," cannot fit into ",n," signed byte(s).\n" 377 return 378 } 379 x=s2un(x,n) 380 binary_uint(x,n) 381 hex_uint(x,n) 382} 383define void uint8(x){uintn(x,1)} 384define void int8(x){intn(x,1)} 385define void uint16(x){uintn(x,2)} 386define void int16(x){intn(x,2)} 387define void uint32(x){uintn(x,4)} 388define void int32(x){intn(x,4)} 389define void uint64(x){uintn(x,8)} 390define void int64(x){intn(x,8)} 391define void uint(x){uintn(x,ubytes(x))} 392define void int(x){intn(x,sbytes(x))} 393define bunrev(t){ 394 auto a,s,m[] 395 s=scale 396 scale=0 397 t=abs(t)$ 398 while(t!=1){ 399 t=divmod(t,2,m[]) 400 a*=2 401 a+=m[0] 402 } 403 scale=s 404 return a 405} 406define band(a,b){ 407 auto s,t,m[],n[] 408 a=abs(a)$ 409 b=abs(b)$ 410 if(b>a){ 411 t=b 412 b=a 413 a=t 414 } 415 s=scale 416 scale=0 417 t=1 418 while(b){ 419 a=divmod(a,2,m[]) 420 b=divmod(b,2,n[]) 421 t*=2 422 t+=(m[0]&&n[0]) 423 } 424 scale=s 425 return bunrev(t) 426} 427define bor(a,b){ 428 auto s,t,m[],n[] 429 a=abs(a)$ 430 b=abs(b)$ 431 if(b>a){ 432 t=b 433 b=a 434 a=t 435 } 436 s=scale 437 scale=0 438 t=1 439 while(b){ 440 a=divmod(a,2,m[]) 441 b=divmod(b,2,n[]) 442 t*=2 443 t+=(m[0]||n[0]) 444 } 445 while(a){ 446 a=divmod(a,2,m[]) 447 t*=2 448 t+=m[0] 449 } 450 scale=s 451 return bunrev(t) 452} 453define bxor(a,b){ 454 auto s,t,m[],n[] 455 a=abs(a)$ 456 b=abs(b)$ 457 if(b>a){ 458 t=b 459 b=a 460 a=t 461 } 462 s=scale 463 scale=0 464 t=1 465 while(b){ 466 a=divmod(a,2,m[]) 467 b=divmod(b,2,n[]) 468 t*=2 469 t+=(m[0]+n[0]==1) 470 } 471 while(a){ 472 a=divmod(a,2,m[]) 473 t*=2 474 t+=m[0] 475 } 476 scale=s 477 return bunrev(t) 478} 479define bshl(a,b){return abs(a)$*2^abs(b)$} 480define bshr(a,b){return (abs(a)$/2^abs(b)$)$} 481define bnotn(x,n){ 482 auto s,t,m[] 483 s=scale 484 scale=0 485 t=2^(abs(n)$*8) 486 x=abs(x)$%t+t 487 t=1 488 while(x!=1){ 489 x=divmod(x,2,m[]) 490 t*=2 491 t+=!m[0] 492 } 493 scale=s 494 return bunrev(t) 495} 496define bnot8(x){return bnotn(x,1)} 497define bnot16(x){return bnotn(x,2)} 498define bnot32(x){return bnotn(x,4)} 499define bnot64(x){return bnotn(x,8)} 500define bnot(x){return bnotn(x,ubytes(x))} 501define brevn(x,n){ 502 auto s,t,m[] 503 s=scale 504 scale=0 505 t=2^(abs(n)$*8) 506 x=abs(x)$%t+t 507 scale=s 508 return bunrev(x) 509} 510define brev8(x){return brevn(x,1)} 511define brev16(x){return brevn(x,2)} 512define brev32(x){return brevn(x,4)} 513define brev64(x){return brevn(x,8)} 514define brev(x){return brevn(x,ubytes(x))} 515define broln(x,p,n){ 516 auto s,t,m[] 517 s=scale 518 scale=0 519 n=abs(n)$*8 520 p=abs(p)$%n 521 t=2^n 522 x=abs(x)$%t 523 if(!p)return x 524 x=divmod(x,2^(n-p),m[]) 525 x+=m[0]*2^p%t 526 scale=s 527 return x 528} 529define brol8(x,p){return broln(x,p,1)} 530define brol16(x,p){return broln(x,p,2)} 531define brol32(x,p){return broln(x,p,4)} 532define brol64(x,p){return broln(x,p,8)} 533define brol(x,p){return broln(x,p,ubytes(x))} 534define brorn(x,p,n){ 535 auto s,t,m[] 536 s=scale 537 scale=0 538 n=abs(n)$*8 539 p=abs(p)$%n 540 t=2^n 541 x=abs(x)$%t 542 if(!p)return x 543 x=divmod(x,2^p,m[]) 544 x+=m[0]*2^(n-p)%t 545 scale=s 546 return x 547} 548define bror8(x,p){return brorn(x,p,1)} 549define bror16(x,p){return brorn(x,p,2)} 550define bror32(x,p){return brorn(x,p,4)} 551define bror64(x,p){return brorn(x,p,8)} 552define brol(x,p){return brorn(x,p,ubytes(x))} 553define bmodn(x,n){ 554 auto s 555 s=scale 556 scale=0 557 x=abs(x)$%2^(abs(n)$*8) 558 scale=s 559 return x 560} 561define bmod8(x){return bmodn(x,1)} 562define bmod16(x){return bmodn(x,2)} 563define bmod32(x){return bmodn(x,4)} 564define bmod64(x){return bmodn(x,8)} 565