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