New post
Question
· Jul 9, 2019

Calculate PI

#Caché

Hi,

I am looking for an Objectscript example to calculate PI (3,14) with at least 100 Decimals.

There are many examples on the internet (https://rosettacode.org/wiki/Pi), but I want to use it for testing purposes.

Regards,

Peter van Voorthuysen

Discussion (9)2
Log in or sign up to continue
Vitaliy Serdtsev · Jul 10, 2019

Translation of: BASIC256

calcPI(n=1000) public { 
  s len = 10*n\4,
    needdecimal = $$$YES,
    nines = 0,
    predigit = 0  ;# {First predigit is a 0}
   
  f j=1:1:len s a(j-1)=2 ;# {Start with 2s}
  
  f j=1:1:n {
    s q=0
    f i=len:-1:1 {
      ;#  {Work backwards}
      s x = 10*a(i-1) + (q*i),
        a(i-1)=x#(2*i-1),
        q=x\(2*i-1)
    }
    s a(0)=q#10,
      q=q\10
    i q=9 {
      s nines = nines + 1
    }elseif q=10 {
      s d = predigit+1 d outputd
      i nines>0 f k = 1:1:nines s d =  0 d outputd
      s predigit = 0, nines = 0
    }else{
      s d = predigit,predigit = q d outputd
      i nines {
        f k = 1:1:nines s d = 9 d outputd
        s nines = 0
      }
    }
  }
  w predigit
  q
   
outputd()
  if needdecimal {
     q:d=0
     w d_"."
     s needdecimal = $$$NO
  } else {
     w d
  }
}

The greater "n", the higher the accuracy.

If there is a lack of RAM, you can easily replace the local array "a" with globals ^||a or ^a.

2 0
Vitaliy Serdtsev · Jul 10, 2019

Translation of: LUA

calcPILua(n=1000) public { 
  s len = 10*n\3,
    nines = 0,
    predigit = 0
  
  f j=1:1:len s a(j)=2
  
  f j=1:1:n {
    s q=0
    f i=len:-1:1 {
      s x = 10*a(i) + (q*i),
        a(i)=x#(2*i-1),
        q=x\(2*i-1)
    }
    s a(1)=q#10,
      q=q\10
    i q=9 {
      s nines = nines + 1
    }elseif q=10 {
      w predigit+1
      f k = 1:1:nines w 0
      s predigit = 0, nines = 0
    }else{
      w predigit
      s predigit = q
      i nines {
        f k = 1:1:nines w 9
        s nines = 0
      }
    }
  }
  w predigit
}

The result of this example is exactly the same as the result of the program C# (tested at n=10000).

2 0
Jon Willeke · Jul 10, 2019

It occurs to me that the Visual Basic solution works with minimal modification as Caché Basic: weaken the DIM statements, change Debug.Print to Print, and replace the call to Format$().

Option Explicit

Const VECSIZE = 3350
Const BUFSIZE = 201
Dim buffer
Dim vect
Dim more, karray, num, k, l, n

  For n = 1 To VECSIZE
    vect(n) = 2
  Next n
  For n = 1 To BUFSIZE
    karray = 0
    For l = VECSIZE To 1 Step -1
      num = 100000 * vect(l) + karray * l
      karray = num \ (2 * l - 1)
      vect(l) = num - karray * (2 * l - 1)
    Next l
    k = karray \ 100000
    buffer(n) = more + k
    more = karray - k * 100000
  Next n
  Print buffer(1)
  Print "."
  l = 0
  For n = 2 To BUFSIZE
    Print Right("00000" & buffer(n), 5)
    l = l + 1
    If l = 10 Then
      l = 0
      Print 'line feed
    End If
  Next n
1 0
Stuart Salzer · Jul 15, 2019

I wrote this a while ago. Just ran on a MacBook Pro Mid 2015, 2.8 GHz Intel Core i7.

It computed the first 1000 digits in less than 1 minute, 2000 digits in 5 minute, 3000 digits, in 14 minute, 4000 digits in 32 minute, 5000 digits in 60 minute...

PI    ; SRS 2011-07-18
    ; Compute PI using Plouffe and Bellard algorithm, based upon "C" code
    ; found at <http://iweb.dl.sourceforge.net/project/projectpi/
    ; Digit%20Extraction%20Methods/Plouffe%20and%20Bellard%20v1/pi1_f.c>.
    KILL ^PI,^PRIME
    SET ^PI=-1
    FOR I=0:1:9 JOB C($INCREMENT(^PI))
    QUIT
C(p)    NEW $ETRAP SET $ETRAP="DO ^%ETN"
    SET n=p*9+1
    SET bign=+$TRANSLATE(n+20*$ZLN(10)/$ZLN(2),".","!")
    SET sum=0
    SET a=3 FOR  QUIT:a>(2*bign)  DO  SET a=$$NP(a)
    . SET vmax=+$PIECE($ZLN(2*bign)/$ZLN(a),".")
    . SET av=1 FOR i=1:1:vmax SET av=av*a
    . SET s=0,num=1,den=1,v=0,kq=1,kq2=1
    . FOR k=1:1:bign DO
    . . SET t=k
    . . IF kq'<a FOR  SET t=t\a,v=v-1 IF t#a SET kq=0 QUIT
    . . SET kq=kq+1,num=num*t#av
    . . SET t=2*k-1
    . . DO:kq2'<a 
    . . . IF kq2=a FOR  SET t=t\a,v=v+1 QUIT:t#a
    . . . SET kq2=kq2-a
    . . SET den=den*t#av,kq2=kq2+2
    . . DO:v>0
    . . . SET t=$$IM(den,av),t=t*num#av,t=t*k#av
    . . . FOR i=v+1:1:vmax SET t=t*a#av
    . . . SET s=s+t
    . . . SET:s'<av s=s-av
    . SET t=$$PM(10,n-1,av),s=s*t#av
    . SET sum=sum+(s/av),sum=+("."_$PIECE(sum,".",2))
    SET ^PI(p)=$EXTRACT($PIECE(sum,".",2)_"000000000",1,9)
    JOB C($INCREMENT(^PI))
    QUIT
NP(a)    NEW (a) FOR  SET r=$ORDER(^PRIME(a)) QUIT:r'=""  DO  QUIT:r'=""
    . LOCK ^PRIME 
    . SET r=$ORDER(^PRIME(a)) IF r'="" LOCK  QUIT
    . IF $DATA(^PRIME)=0 SET ^PRIME=3,^PRIME(2)="",^PRIME(3)="" LOCK  QUIT
    . FOR r=^PRIME:2 DO  IF pr SET ^PRIME(r)="" QUIT:r>a
    . . SET pr=1 FOR p=3:2:$ZSQR(r) IF r#p=0 SET pr=0 QUIT
    . SET ^PRIME=r
    . LOCK
    QUIT r
IM(x,y)    NEW (x,y)
    SET u=x,v=y,c=1,a=0
    FOR  SET q=v\u,t=c,c=a-(q*c),a=t,t=u,u=v-(q*u),v=t QUIT:u=0
    SET a=a#y
    SET:a<0 a=a+y
    QUIT a
PM(a,b,m)    NEW (a,b,m)
    SET r=1,aa=a
    FOR  SET:b#2 r=r*aa#m SET b=b\2 QUIT:b=0  SET aa=aa*aa#m
    QUIT r
 

0 0
Vitaliy Serdtsev · Jul 16, 2019

Intel i5-2400

10000 digits ~ 58 sec.

calcPI(n) public { 
  s $lb(len,nines,predigit,r)=$lb(10*n\3,0,0,"")
  
  f i=1:1:len s a(i)=2
  
  f j=1:1:n {
    s q=0
    f i=len:-1:1 s x=10*a(i)+(q*i), a(i)=x#(2*i-1), q=x\(2*i-1)
    s a(1)=q#10, q=q\10
    i q=9 {
      s nines=nines+1
    }elseif q=10 {
      s r=r_(predigit+1)_$$repeat^%qarfunc(0,nines), predigit=0, nines=0
    }else{
      s r=r_predigit, predigit=q
      s:nines r=r_$$repeat^%qarfunc(9,nines), nines=0
    }
  }
  q r_predigit
}
0 0