Computing the value of sine function accurately

[Theia]

Hi all
Simple question: How I can compute the value of \(a = \sin \left( 2017 \sqrt[5]{2} \right) \) under following assumptions:

• No use of advanced numerical libraries is allowed.
• Only accepted operations are: comparisons, absolute value, addition, subtraction, multiplication and division.
• More advanced operations must be returned to the use of allowed operations.
• End result is needed by 25 digit accuracy.

Ps. I know, \( a \approx -1 \), so equivalent way is to find the number \(b = 1 + a. \) But what is the easiest way to do that? :unsure:

Thank you!

 

[I like Serena]

How about:

1. Evaluate $\sqrt[5] 2$ by Newton Raphson with $f(x)=x^5-2$ and $x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)},\,x_0=2$.
2. Evaluate $\tau=2\pi$ by one of its many algorithms. For instance $\tau = 8\tan^{-1}(1)=8\Big(1-\frac 13+\frac 15-\ldots\Big)$.
3. Evaluate $x\overset{def}=\Big(2017\sqrt[5] 2 \bmod{\tau}\Big)-\frac 34\tau$.
4. Evaluate $a=\sin\Big(x+\frac 34\tau\Big)=-\cos x=-\Big(1-\frac{x^2}{2!}+\frac{x^4}{4!}-\ldots\Big)$.

Intermediate results will need more than 25 digit accuracy.
A ballpark figure of 30 digit accuracy should suffice in this case.

 

[Theia]

Interesting! At least worth to try and see how it works! :giggle: Thank you!

 

[Theia]

My approach was:

I typed a very robust library to handle n-digit arithmetic. And then followed the scheme Klaas van Aarsen wrote.

# -*- coding: utf-8 -*-

#procedures for N-digit computation

#number of decimals
N = 30

def tofloat(string):
  #to convert number from string to N-float
  minus = string.find('-')
  if(minus > -1): #number is negative
    minus = -1
    #remove it from the string
    string = string.replace('-', '', 1)
  else:
    minus = 1
  dot = string.find('.')
  if(dot == -1): #if decimal point does not exist
    #return integer part
    ip = int(string)
    dp = 0
    number = [ip, dp, minus]
    return number
  #if decimal point exists
  ip = int(string[:dot]) #integer part
  dps = string[dot+1:] #decimal part as string
  ldps = len(dps) #length of decimal part
  if(ldps == N): #no extra digits
    dp = int(dps)
  elif(ldps > N): #extra digits
    dp = int(dps[:N]) #first N digits
    e = int(dps[N]) #next digit
    if(e > 4):
      dp += 1
    if(dp > ONE1):
      dp -= ONE
      ip += 1
  else: #too few digits
    a = N - ldps
    zeros = '0' *a
    dp = int(dps + zeros)
  number = [ip, dp, minus]
  return number
    
def show(number,out=1):
  #routines to print N-digit number
  ip = number[0]
  dp = number[1]
  minus = number[2]
  string = ''
  if(minus == -1):
    string = '-'
  string = string + str(ip)
  if(dp > 0):
    #how many digits in dp?
    dps = str(dp)
    a = N - len(dps) # zeros missing
    if(a > 0):
      zeros = '0' *a
    else:
      zeros = ''
    string = string + '.' + zeros + dps
  if(out == 1):
    print(string)
    return
  return string

def sum(a,b):
  #summation of numbers a and b
  na = a[0]*ONE + a[1]
  if(a[2] < 0):
    na = -na
  nb = b[0]*ONE + b[1]
  if(b[2] < 0):
    nb = -nb
  n = na + nb
  sig = 1
  #analyze the result
  if(n < 0):
    sig = -1
    n = -n
  number = [n // ONE,n % ONE,sig]
  return number

def subs(a,b):
  #substraction a - b
  number = sum(a,[b[0],b[1],-b[2]])
  return number

def prod(a,b):
  #product a*b
  #determine the sign of the result
  sig = a[2]*b[2]
  ip = a[0]*b[0]
  dp1 = a[0]*b[1]
  dp2 = a[1]*b[0]
  dp3 = a[1]*b[1]
  dp31 = dp3 // ONE
  #need to find the first dropped digit
  f = dp3 - ONE*dp31
  #one1d = ONE // 10
  f = f // ONE1D
  if(f > 4):
    dp31 += 1
  dp = dp1 + dp2 + dp31
  if(dp > ONE1):
    k = dp // ONE
    dp = dp % ONE
    ip += k
  number = [ip, dp, sig]
  return number

def div(a,b):
  #division a/b
  #sign of the result
  sig = a[2]*b[2]
  na = a[0]*ONE + a[1]
  na *= NSIFT1
  nb = b[0]*ONE + b[1]
  n = na // nb
  #last digit is for rounding
  r = n % 10
  n -= r
  n = n // 10
  #build number
  ip = n // ONE
  dp = n % ONE
  if(r > 4):
    dp += 1
    if(dp > ONE1):
      dp -= ONE
      ip += 1
  number = [ip,dp,sig]
  return number

NZEROS = '0' *N
NSIFT1 = int('1' + NZEROS + '0')
ONE = int('1' + NZEROS)
ONE1 = ONE - 1
ONE1D = ONE // 10

# -*- coding: utf-8 -*-

import ndigit as nd

#compute the value of a = sin(2017*(2)^(1/5))

def newite(x):
  #computes 5th root of x
  o = [1,0,1]
  r = nd.prod(x,nd.tofloat('0.3'))
  p4 = nd.tofloat('4')
  while((o[0] > 0) or ((o[0] == 0) and (o[1] > 20))):
    r2 = nd.prod(r,r)
    r4 = nd.prod(r2,r2)
    r5 = nd.prod(r,r4)
    o = nd.subs(r5,x)
    n = nd.prod(p4,r4)
    d = nd.div(o,n)
    r = nd.subs(r,d)
  return r

def macatan(x,err = 20):
  #computes atan(x) with error less than err/N units
  x2 = nd.prod(x,x)
  term = [5,0,1]
  su = x
  n = [1,0,1]
  sig = -1
  while((term[0] > 0) or ((term[0] == 0) and (term[1] > err))):
    x = nd.prod(x,x2)
    n = nd.sum(n,[2,0,1])
    term = nd.div(x,n)
    if(sig < 0):
      su = nd.subs(su,term)
      sig = 1
    else:
      su = nd.sum(su,term)
      sig = -1
  return su

#compute 2^(1/5) by Newton iteration
fth2 = newite(nd.tofloat('2'))

#compute pi using Machin's serie:
#pi = 4*(4*atan(1/5) - atan(1/239))
#   = 16*atan(1/5) - 4*atan(1/239)

x1 = nd.tofloat('0.2')
x2 = [239,0,1]
x2 = nd.div([1,0,1],x2) # 1/239

pi1 = macatan(x1,1)
pi1 = nd.prod(pi1,[16,0,1])
pi2 = macatan(x2)
pi2 = nd.prod(pi2,[4,0,1])

pi = nd.subs(pi1,pi2) # pi
pi2 = nd.sum(pi,pi) # 2*pi

#compute the number 2017*2^(1/5)...
x = nd.prod([2017,0,1],fth2)
#and remove 2pi-multiples
n2pi = nd.div(x,pi2)
xremove = nd.prod(pi2,[n2pi[0],0,1])
x = nd.subs(x,xremove)
#now x ~ 3pi/2. Difference is
pi32 = nd.sum(pi,pi2) # 3*pi
pi32 = nd.prod(pi32,nd.tofloat('0.5')) # 3*pi/2
x = nd.subs(x,pi32)

#compute the a = -cos(x) = -1 + x^2/2 - ...
#difference from -1

x2 = nd.prod(x,x)
x = [1,0,1]
sig = 1
v = 0
nim = [1,0,1]
term = [0,45677,1]
su = [0,0,1]
while(term[1] > 20):
  x = nd.prod(x,x2)
  v += 1
  nim = nd.prod(nim,[v,0,1])
  v += 1
  nim = nd.prod(nim,[v,0,1])
  term = nd.div(x,nim)
  if(sig > 0):
    su = nd.sum(su,term)
    sig = -1
  else:
    su = nd.subs(su,term)
    sig = 1

print('difference from -1:', nd.show(su,0))
a = nd.sum([1,0,-1],su)
print('value:', nd.show(a,0))

#correct value:
#-0.999999999999999978567771261060983