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