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