contrib/bc/gen/lib2.bc

   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
  36 define p(x,y){
  37         auto a
  38         a=y$
  39         if(y==a)return(x^a)@scale
  40         return e(y*l(x))
  41 }
  42 define 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 }
  53 define 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 }
  63 define f(n){
  64         auto r
  65         n=abs(n)$
  66         for(r=1;n>1;--n)r*=n
  67         return r
  68 }
  69 define 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 }
  82 define 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 }
  96 define 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 }
 108 define 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 }
 118 define l2(x){return log(x,2)}
 119 define l10(x){return log(x,A)}
 120 define 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 }
 147 define cbrt(x){return root(x,3)}
 148 define 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 }
 168 define 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 }
 179 define 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 }
 189 define 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 }
 198 define 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 }
 213 define sin(x){return s(x)}
 214 define cos(x){return c(x)}
 215 define atan(x){return a(x)}
 216 define tan(x){return t(x)}
 217 define atan2(y,x){return a2(y,x)}
 218 define 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 }
 229 define 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 }
 240 define frand(p){
 241         p=abs(p)$
 242         return irand(A^p)>>p
 243 }
 244 define ifrand(i,p){return irand(abs(i)$)+frand(p)}
 245 define srand(x){
 246         if(irand(2))return -x
 247         return x
 248 }
 249 define brand(){return irand(2)}
 250 define void output(x,b){
 251         auto c
 252         c=obase
 253         obase=b
 254         x
 255         obase=c
 256 }
 257 define void hex(x){output(x,G)}
 258 define void binary(x){output(x,2)}
 259 define 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 }
 266 define 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 }
 275 define 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 }
 291 define s2u(x){return s2un(x,sbytes(x))}
 292 define 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 }
 302 define void plznl(x){
 303         plz(x)
 304         print"\n"
 305 }
 306 define 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 }
 324 define void pnlznl(x){
 325         pnlz(x)
 326         print"\n"
 327 }
 328 define 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 }
 343 define 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 }
 351 define void hex_uint(x,n){
 352         auto o
 353         o=obase
 354         obase=G
 355         output_uint(x,n)
 356         obase=o
 357 }
 358 define void binary_uint(x,n){
 359         auto o
 360         o=obase
 361         obase=2
 362         output_uint(x,n)
 363         obase=o
 364 }
 365 define 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 }
 381 define 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 }
 396 define void uint8(x){uintn(x,1)}
 397 define void int8(x){intn(x,1)}
 398 define void uint16(x){uintn(x,2)}
 399 define void int16(x){intn(x,2)}
 400 define void uint32(x){uintn(x,4)}
 401 define void int32(x){intn(x,4)}
 402 define void uint64(x){uintn(x,8)}
 403 define void int64(x){intn(x,8)}
 404 define void uint(x){uintn(x,ubytes(x))}
 405 define void int(x){intn(x,sbytes(x))}
 406 define 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 }
 419 define 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 }
 440 define 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 }
 466 define 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 }
 492 define bshl(a,b){return abs(a)$*2^abs(b)$}
 493 define bshr(a,b){return(abs(a)$/2^abs(b)$)$}
 494 define 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 }
 509 define bnot8(x){return bnotn(x,1)}
 510 define bnot16(x){return bnotn(x,2)}
 511 define bnot32(x){return bnotn(x,4)}
 512 define bnot64(x){return bnotn(x,8)}
 513 define bnot(x){return bnotn(x,ubytes(x))}
 514 define 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 }
 523 define brev8(x){return brevn(x,1)}
 524 define brev16(x){return brevn(x,2)}
 525 define brev32(x){return brevn(x,4)}
 526 define brev64(x){return brevn(x,8)}
 527 define brev(x){return brevn(x,ubytes(x))}
 528 define 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 }
 542 define brol8(x,p){return broln(x,p,1)}
 543 define brol16(x,p){return broln(x,p,2)}
 544 define brol32(x,p){return broln(x,p,4)}
 545 define brol64(x,p){return broln(x,p,8)}
 546 define brol(x,p){return broln(x,p,ubytes(x))}
 547 define 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 }
 561 define bror8(x,p){return brorn(x,p,1)}
 562 define bror16(x,p){return brorn(x,p,2)}
 563 define bror32(x,p){return brorn(x,p,4)}
 564 define bror64(x,p){return brorn(x,p,8)}
 565 define brol(x,p){return brorn(x,p,ubytes(x))}
 566 define 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 }
 574 define bmod8(x){return bmodn(x,1)}
 575 define bmod16(x){return bmodn(x,2)}
 576 define bmod32(x){return bmodn(x,4)}
 577 define bmod64(x){return bmodn(x,8)}