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