TABLELESS CALCULATION OF CIRCULAR FUNCTIONSON DYADIC RATIONALS
a continued radical approachkourzanov@acm.org[1]1 . 1

INTRODUCTION• Trigonometry use is widespread

INTRODUCTION• Trigonometry use is widespread• Circular functions and their inverses (as well as square roots) are of great practical value

INTRODUCTION• Trigonometry use is widespread• Circular functions and their inverses (as well as square roots) are of great practical value• Complex numbers, frequency estimation, FFT, DFT and other DSP processing algorithms etc.

QUESTIONS
1. Do we need special implementations of sqrt, sin, cos etc.?
2. Can we not reduce one implementation to others?
3. Algorithms for inverse functions derivable from forward ones?
4. Avoid Taylor/Maclaurin and pre-computed tables a-la CORDIC?
5. Can we enumerate all circular triginometric values?

CONTEXTn• Dyadic rationals x = m, m = 2k , n ∈ Z, k ∈ Nnπ• Focus on cos 2m, by symmetry, restrict x ∈ [0, π2 ]CONSTRAINTS• maintaining explicit tables containing ﬁxed-point values for e.g., FFT or a CORDIC can be too
expensive for resource-constrained implementations.
• For the software, the eﬃciency is limited by the memory hierarchy bottleneck• For the hardware, excessive (de)multiplexing that is associated to large (frequently > 4Kbit)
tables can outweigh the advantages of avoiding explicit multipliers• There is a trade-oﬀ between multiplier-based approach with series and multiplexer-based
approach with tabling that a resource-constrained application may not want to make.

MATHEMATICAL BACKGROUND• Basic identity cos π4 = sin π4 = 1

=
2

1
1
2 = 2

x+π

1+cos x2

• Half-angle formulæ cos x2 = (−1)⌊ 2π ⌋2MATHEMATICAL BACKGROUND• Basic identity cos π4 = sin π4 = 1

=
2

1
1
2 = 2

x+π

1+cos x2

• Half-angle formulæ cos x2 = (−1)⌊ 2π ⌋

2

⌊

x

• Continued radical expansion where ski (x) = (−1) π2k−ixcos k = sk0 (x)2+ 12 ⌋x
1 + cos 2k−1
2MATHEMATICAL BACKGROUND• Basic identity cos π4 = sin π4 = 1

=
2

1
1
2 = 2

x+π

1+cos x2

• Half-angle formulæ cos x2 = (−1)⌊ 2π ⌋

2

⌊

x

• Continued radical expansion where ski (x) = (−1) π2k−ixcos k = sk0 (x)2= sk0 (x)+ 12 ⌋x
1 + cos 2k−1
21
+ sk1 (x)2x
1 + cos 2k−2
22 ∗ 2

MATHEMATICAL BACKGROUND• Basic identity cos π4 = sin π4 = 1

=
2

1
1
2 = 2

x+π

1+cos x2

• Half-angle formulæ cos x2 = (−1)⌊ 2π ⌋

2

⌊

x

• Continued radical expansion where ski (x) = (−1) π2k−ixcos k = sk0 (x)2= sk0 (x)= sk0 (x)x
1 + cos 2k−1
21
+ sk1 (x)
2

1
k (x)
+
s
1
20 ∗ 2+ 12 ⌋x
1 + cos 2k−2
22 ∗ 2

1
k (x) …
+s2
22 ∗ 2

MATHEMATICAL BACKGROUND• Basic identity cos π4 = sin π4 = 1

=
2

1
1
2 = 2

x+π

1+cos x2

• Half-angle formulæ cos x2 = (−1)⌊ 2π ⌋

2

⌊

x

• Continued radical expansion where ski (x) = (−1) π2k−ix
1 + cos 2k−1
2x
cos k = sk0 (x)2= sk0 (x)+ 12 ⌋

1+ sk1 (x)2

= sk0 (x)

1
k (x)
+
s
1
20 ∗ 2= sk0 (x)

2
k (x)
+
s
1
1
22

x
1 + cos 2k−2
22 ∗ 2

1
k (x) …
+s2
22 ∗ 2

2
k (x)
+
s
2
2
22

2
+…
223

ENCODING COS VALUES IN HEXADECAIC INTERVAL
ski (n)nm(+,+,+)

1
16

(+,+)

1
8

(+,+,-)

3
16

(+)

1
4

(+,-,-)

5
16

(+,-)

3
8

(+,-,+)

7
16

()

1
2

(-,-,+)

9
16

(-,-)

5
8

(-,-,-)

11
16

(-)

3
4

(-,+,-)

13
16

(-,+)

7
8

(-,+,+)

1516nπ
cos 2m

2−1 +

2−3 +

2−1 +
2−1 +

2−1 +

2−7 +

2−3 +

2−7

2−3 +

2−7 −

2−1 +

2−3

2−3 −

2−7 −

2−1 +
2−1 +

sig sig' tree

2−3 −

2−3 −

2−7

2−7 +

2−1
2−1 −

2−3 −

2−1 −
2−1 −

2−1 −

2−3 −

2−7

2−3 −

2−7 −

2−1 −

2−3

2−3 +

2−7 −

2−1 −
2−1 −

2−7 +

2−3 +

2−3 +

2−7

2−7 +

2−15 016

082−15 117

0

42−15 319

192−15 218022−15 622

3112−15 7231

52−15 521

2102−15 420

PRELIMINARIESFRACTAL SEQUENCES FOR sig
function S(x::Rational,init=0,sign=0)::Vector{Int}ss=Vector{Int}([init]); m=denominator(x)>>1
while numerator(x)!=m
if numerator(x)<mss=[mod(sum(ss),2)==0 ? sign : 1-sign; ss]
x=numerator(x)//m
else
ss=[mod(sum(ss),2)==0 ? 1-sign : sign; ss]
x=(numerator(x)-m)//mend
m=denominator(x)>>1
end
return ss
endPRELIMINARIESFRACTAL SEQUENCES FOR sig
function S(x::Rational,init=0,sign=0)::Vector{Int}ss=Vector{Int}([init]); m=denominator(x)>>1
while numerator(x)!=m
if numerator(x)<mss=[mod(sum(ss),2)==0 ? sign : 1-sign; ss]
x=numerator(x)//m
else
ss=[mod(sum(ss),2)==0 ? 1-sign : sign; ss]
x=(numerator(x)-m)//mend
m=denominator(x)>>1
end
return ss
end

H(k)= if k==1 return [0]elsehk=H(k-1)
r=zeros(Int,2^k-1)
r[2:2:end]=hk
r[1:2:end>>1]=hk[1:2:end]r[1 .+ (end+1)>>1:2:end]=reverse(hk[1:2:end]).+2^(k-2)
return r
endPRELIMINARIESFRACTAL SEQUENCES FOR sig
function S(x::Rational,init=0,sign=0)::Vector{Int}ss=Vector{Int}([init]); m=denominator(x)>>1
while numerator(x)!=m
if numerator(x)<mss=[mod(sum(ss),2)==0 ? sign : 1-sign; ss]
x=numerator(x)//m
else
ss=[mod(sum(ss),2)==0 ? 1-sign : sign; ss]
x=(numerator(x)-m)//mend
m=denominator(x)>>1
end
return ss
end

H(k)= if k==1 return [0]elsehk=H(k-1)
r=zeros(Int,2^k-1)
r[2:2:end]=hk

H1(k,seed=0)=
if k==1
return [seed]else
h=H1(k-1,seed)[h;flat(zip(h[1 .+ end>>1:end],
reverse(h[1:2:end]).+2^(k-2)

r[1:2:end>>1]=hk[1:2:end]r[1 .+ (end+1)>>1:2:end]=reverse(hk[1:2:end]).+2^(k-2)
return r
end

))]
end

ANOTHER OEIS -BASED FORMULATION
G(k,seed=0)=if k==0 return [seed]

K1(k,seed=1)=if k==1; return [seed]else
g=G(k-1,seed)else
hk=K1(k-1,seed)[g;reverse(g).+2^(k-1)]
end[flat(zip(1:2^(k-1)-1,hk));2^(k-1)]
end

Gray sequence
A003188(n)=if n==1 return 0else k=trunc(Int,log2(n-1))
g=A003188.(1:2^k)[g;reverse(g).+2^k][n]endKimberling paraphrases
A003602(n)=if n==1; return 1
else k=trunc(Int,log2(n-1))
ak=A003602.(1:2^k)[flat(zip(1:2^k,ak));2^k][n]endsig = A003188 ( A003602 ), i.e., modulation of the Grey sequence by Kimberling’s paraphrases[2]verify(k) = G(k)[K1(k)]==H1(k)SQRT IN SCHEME(define (sqrt-loop n start)(let loop ([num n][res 0][bit start])(cond ([zerofx? bit] res)([>=fx num (+fx res bit)](loop (-fx num (+fx res bit))(+fx bit (bit-rsh res 1))(bit-rsh bit 2)))(else
(loop number(bit-rsh res 1)(bit-rsh bit 2))))))(def-syntax *maxbit* 60)(define (sqrt-start n)(do ([bit (bit-lsh 1 (-fx *maxbit* 2))(bit-rsh bit 2)])([<=fx bit n] bit)))(Define (sqrt fix: n)(sqrt-loop n (sqrt-start n)))SQRT IN SCHEME(SIMPLIFIED)(define (sqrt-loop n start)(let loop ([num n](Define (sqrt-loop fix: n start)(Let loop (fix:[num n][res 0][bit start])(cond ([zerofx? bit] res)([>=fx num (+fx res bit)](loop (-fx num (+fx res bit))(+fx bit (bit-rsh res 1))(bit-rsh bit 2)))[res 0][bit start])(cond ([zero? bit] res)([>= num (+ res bit)](loop (- num (+ res bit))⟶(+ bit (rsh res 1))(rsh bit 2)))

(else

(else

(loop number(bit-rsh res 1)(loop num(rsh res 1)(bit-rsh bit 2))))))(def-syntax *maxbit* 60)(define (sqrt-start n)(do ([bit (bit-lsh 1 (-fx *maxbit* 2))(bit-rsh bit 2)])([<=fx bit n] bit)))(rsh bit 2))))))

(Define (sqrt fix: n)(sqrt-loop n (sqrt-start n)))

DEVELOPMENTCOS USING CONTINUATIONS(Define (dyadic-cos-cont n :fix m :fix sig :fix f :flo cont :proc);; exact(cond ([= n 0] (cont 1))([= n m] (cont 0));; results(else (Let* (fix:[m (rsh m 1)][sign (if [= (bit-and sig 1) 1] -1 1)])(cond ([= n m](cont (* sign (sqrt (* 2 f)))))([< n m](dyadic-cos-cont n m(bit-or (bit-lsh sig 1)(if (even? (bit-count sig)) 0 1))(sqr f)(λ (x) (cont (* sign (sqrt (+ (* 2 f) x)))))))([> n m](dyadic-cos-cont (- n m) m
(bit-or (bit-lsh sig 1)(if (even? (bit-count sig)) 1 0))(sqr f)(λ (x) (cont (* sign (sqrt (+ (* 2 f) x)))))
))
)))))

KEY INSIGHT #1nπ• sig signature encodes the value of cos 2m∈ R as a N valuen= 2−k , k ≥ 1• sig has a value of 0 for all m• because of performance impact of closures, depth needs to be tracked explicitly(or with lf as on the next slide)• with depth, sig is a surjective decoding of a N value as a forward circular R value

so enumerating sig at increasing depth does enumerate all values

H1.(1:6)
1: [0]
2: [0, 0, 1]
3: [0, 0, 1, 0, 3, 1, 2]
4: [0, 0, 1, 0, 3, 1, 2, 0, 6, 3, 7, 1, 5, 2, 4]
5: [0, 0, 1, 0, 3, 1, 2, 0, 6, 3, 7, 1, 5, 2, 4, 0, 12, 6, 13, 3, 15, 7, 14, 1, 10, 5, 11, 2, 9, 4, 8]
6: [0, 0, 1, 0, 3, 1, 2, 0, 6, 3, 7, 1, 5, 2, 4, 0, 12, 6, 13, 3, 15, 7, 14, 1, 10, 5, 11, 2, 9, 4, 8, 0, 24, 12, 25, 6, 27, 13, 26, 3, 30, 1

DEFUNCTIONALIZATION(Define (dyadic-cos-defun fix: n m sig f lf)(cond ([= n 0] 1)([= n m] 0)(else (Let* ([m :fix (rsh m 1)])(cond ([= n m] (dyadic-apply sig f lf))([< n m] (dyadic-cos-defun n m;; exact;; results(bit-or (bit-lsh sig 1) (if (even? (bit-count sig)) 0 1))(lsh f lf)(lsh lf 1)))([> n m] (dyadic-cos-defun (- n m) m
(bit-or (bit-lsh sig 1) (if (even? (bit-count sig)) 1 0))(lsh f lf)(lsh lf 1))))))))

(Define (dyadic-apply fix: sig f lf)(Let loop (fix:[s sig][f f]
[lf (rsh lf 1)][a :flo 0])(if (zero? lf) a(loop (rsh s 1)
(rsh f lf)(rsh lf 1)(* (if [= (and s 1) 1] -1 1)(sqrt (+ (/ 2 (toflo f)) a)))))))

CONTROL OPTIMIZATIONS #1
Strength reduction to avoid branching(if [= (bit-and s 1) 1]
-1
1)

⟶

{1 ↦ −1, 0 ↦ 1} ≡ 1 − 2(s mod 2)⟶(>-> s :fix
(and _ 1)(lsh _ 1)(- 1 _))

CONTROL OPTIMIZATIONS #2
Syntactic rewriting of code to calculate ski (n)(bit-or(bit-lsh s 1)
(if (even? (bit-count s)) 1 0))⟶(>-> s :fix(lsh _ 1)(or _ (if (even? (bit-count s)) 1 0)))⟶{1 ↦ 0, 0 ↦ 1} ≡ 1 − s mod 2⟶(>-> s(lsh _ 1)(bit-count s);; [Nth 5];; [Nth 4];; [Nth 3](and _ 1)(- 1 _);; [Nth 2];; [Nth 1](or _ [Nth 4])) ;; refer to the above sub-expressions 1 and 4

REPRESENTATION OPTIMIZATIONS
Recall that cos π4 = 12

2. This hints at a possible simpliﬁcation of sk0• f0 = 2−0 = 1• f1 = 2−1• f2 = 2−3

2fi−1−2i +1for i > 0, using ski for ski (n)• thus fi = 2 = 2

f1 + sk1

f2 + sk2 …

REPRESENTATION OPTIMIZATIONS
Recall that cos π4 = 12

2. This hints at a possible simpliﬁcation of sk0

f1 + sk1

f2 + sk2 …

2k
2 + s2
2
2

2
3 + …
2
2

• f0 = 2−0 = 1• f1 = 2−1• f2 = 2−3

2fi−1−2i +1for i > 0, using ski for ski (n)• thus fi = 2 = 2

sk0

2
k
1 + s12
22
k
2 + s22
2

2
sk0
3 + … =
22
2

22
2 21 + 22 sk1
2

REPRESENTATION OPTIMIZATIONS
Recall that cos π4 = 12

2. This hints at a possible simpliﬁcation of sk0

f1 + sk1

f2 + sk2 …

2k
2 + s2
2
2

2
3 + …
2
2

• f0 = 2−0 = 1• f1 = 2−1• f2 = 2−3

2fi−1−2i +1for i > 0, using ski for ski (n)• thus fi = 2 = 2

sk0
sk0
=
2

2
k
1 + s12
22
k
2 + s22
2

2 + sk124
2 22 + 24 sk2
2

2
sk0
3 + … =
22
2
2
sk0
3 + … =
22
2

22
2 21 + 22 sk1
2
2 + sk1

2 + sk228
2 23 + 28 sk3 …
2

REPRESENTATION OPTIMIZATIONS
Recall that cos π4 = 12

2. This hints at a possible simpliﬁcation of sk0

f1 + sk1

f2 + sk2 …

2k
2 + s2
2
2

2
3 + …
2
2

• f0 = 2−0 = 1• f1 = 2−1• f2 = 2−3

2fi−1−2i +1for i > 0, using ski for ski (n)• thus fi = 2 = 2

sk0
sk0
=
2

2
k
1 + s12
22
k
2 + s22
2

2 + sk124
2 22 + 24 sk2
2
sk0
=
2

2
sk0
3 + … =
22
2
2
sk0
3 + … =
22
2
2 + sk1

2 + sk222
2 21 + 22 sk1
2
2 + sk1

2 + sk2

2 + sk3 2 + …

28
2 23 + 28 sk3 …
2

KEY INSIGHT #2• sig' signature just preﬁxes a 1-bit before each value of sig• this makes it an injective mapping from forward circular trigonometric value to Nnπ
• all codes (elements of sig') obtained this way can be converted back to values of cos 2m• depth is tracked by sig'[n] being > 3

N.B. because of the extra preﬁx, bit-counting logic is reversed below[2][4, 2, 5][8, 4, 9, 2, 11, 5, 10][16, 8, 17, 4, 19, 9, 18, 2, 22, 11, 23, 5, 21, 10, 20][32, 16, 33, 8, 35, 17, 34, 4, 38, 19, 39, 9, 37, 18, 36, 2, 44, 22, 45, 11, 47, 23, 46, 5, 42, 21, 43, 10, 41, 20, 40][64, 32, 65, 16, 67, 33, 66, 8, 70, 35, 71, 17, 69, 34, 68, 4, 76, 38, 77, 19, 79, 39, 78, 9, 74, 37, 75, 18, 73, 36, 72, 2, 88, 44, 89, 22,

Enumeration becomes simple: at each depth, take all odd elements:[2, 4, 5, 8, 9, 11, 10, 16, 17, 19, 18, 22, 23, 21, 20, 32, 33, 35, 34, 38, 39, 37, 36, 44, 45, 47, 46, 42, 43, 41, 40, 64, 65, 67, 66]

FINAL ENCODE/DECODE LOOP
(Define (dyadic-trig x :rat sig' :fix)(cond ([= x 0] (* 4 (- 1 (and sig' 1)))) ;; exact([= x 1] (* 4 (and sig' 1)));; results(else (Let encode (fix:[n (numerator x)][m (denominator x)][s sig'])(Let* ([m :fix (rsh m 1)])(cond ([= n m] (dyadic-decode s))([< n m] (encode n m (>-> s (lsh _ 1)
(bit-count s)(and _ 1)(- 1 _);; bit-counting;; logic reversed(or _ [Nth 4]))))([> n m] (encode (- n m) m (>-> s (lsh _ 1)
(bit-count s)(and _ 1);; bit-counting(or _ [Nth 3]);; logic reversed)))))))))(Define (dyadic-decode fix: sig')(Let loop (fix:[s sig'] [a :flo 0])(if (<= s 3) (+ 2 a)(loop (rsh s 1)(* (>-> s (and _ 1) (lsh _ 1) (- 1 _))(sqrt (+ 2 a)))))))

FORWARD AND BACKWARD
FUNCTIONS
Because of the following identity, calculating sin amounts to the same
iteration as for cos, modulo a diﬀerent starting value of sig'

sin x =

1 − cos2 x =1
1 − 2 (2 ±
2

1
…) =
2

4
4 − 2 (2 ±2(Define (dyadic-cos x) ((∘ half sqrt) (dyadic-trig x #b10)))(Define (dyadic-sin x) ((∘ half sqrt) (dyadic-trig x #b11)))(Define (dyadic-tan x) ((∘ sqrt sub1 reciprocal [κ / _ 4]) (dyadic-trig x #b10)))(Define (dyadic-cot x) ((∘ sqrt sub1 reciprocal [κ / _ 4]) (dyadic-trig x #b11)))(Define (dyadic-cos² x) ([κ / _ 4] (dyadic-trig x #b10)))(Define (dyadic-sin² x) ([κ / _ 4] (dyadic-trig x #b11)))(Define (dyadic-tan² x) ((∘ sub1 reciprocal [κ / _ 4]) (dyadic-trig x #b10)))(Define (dyadic-cot² x) ((∘ sub1 reciprocal [κ / _ 4]) (dyadic-trig x #b11)))

1
…) =
2

2∓

…

FURTHER GENERALIZATION(Define (quadrature x :rat inv :fix mir :fix)(cond ((or [> x #r:1/4] [< x #r:-1/4])(quadrature (- x (* (sign x) #r:1/2))(neg inv)(neg mir)))([< x 0](quadrature (neg x) inv (neg mir)))(else(values x inv mir))))(Define (Cos x :rat) (let-values ([(y inv mir) (quadrature x 1 1)])(* inv (dyadic-cos (* y 4)))))
(Define (Sin x :rat) (let-values ([(y inv mir) (quadrature x 1 1)])(* mir (dyadic-sin (* y 4)))))
(Define (Tan x :rat) (let-values ([(y inv mir) (quadrature x 1 1)])(* mir inv (dyadic-tan (* y 4)))))
(Define (Cot x :rat) (let-values ([(y inv mir) (quadrature x 1 1)])(* mir inv (dyadic-cot (* y 4)))))FIRST ATTEMPT AT INVERSE FUNCTIONS(Define (dyadic-arccos x :flo ss :fix *eps* :flo)(cond ([= x 0.] (- 1 ss))([= x 1.] ss)(else(Let loop ([y :rat #r:1/2]
fix:[s ss][d 0])(Let ([m :fix (>-> (denominator y) :fix (lsh _ 1))][xe :flo ((∘ half sqrt) (dyadic-apply' s d))][sss (>-> ss (lsh _ 1) (- 1 _) (toflo _))])(cond ([< (abs (- x xe)) *eps*] y)([< (* sss x) (* sss xe)](loop (+ y #r:1/m)(>-> s (lsh _ 1)
(bit-count s)(and _ 1)(- 1 _)(or _ [Nth 4]))(+ d 1)))([> (* sss x) (* sss xe)](loop (- y #r:1/m)(>-> s (lsh _ 1)
(bit-count s)(and _ 1)(or _ [Nth 3]))(+ d 1)))))))));; exact;; results

KEY INSIGHT #3• the overhead of binary search depends on the argument• target value used to drive the decision process• Let’s fuse the binary search together with the evaluation of the target value used to drive thedecision process, and get rid of the overhead

sk0
x=
2

2 + sk1

2 + sk2

2 + sk3 2 + …

1. view the argument as signature sig' value “decoded” into x ∈ R
2. undo all steps of the forward algorithm in reverse:• initial remainder is set to the input value, doubled and squared:(∘ half sqrt)−1 = (∘ sqr double)• at every un-nesting, determine:▪ if the remainder is approximately equal to 2, terminate the search▪ if less than 2, the sign ski = −1
▪ if greater than 2, the sign is ski = 1▪ In either case, 2 is subtracted and the remainder is squared

IMPROVED LOOP
(Define (dyadic-find x :flo ss :fix *eps* :flo)(cond ([= x 0.] (- 1 ss))([= x 4.] ss)(else (Let loop ([y :rat #r:1/2];; exact;; results[x :flo x]
fix:[s ss])
(Let ([m :fix (>-> (denominator y) :fix (lsh _ 1))])(cond ([< (abs (- x 2)) *eps*] y)([< x 2.](loop (+ y (* (>-> s (bit-count _)(and _ 1)(lsh _ 1)(- 1 _))#r:1/m))(>-> x (- _ 2) (sqr _))(>-> s (lsh _ 1) (or _ 1))))([> x 2.](loop (- y (* (>-> s (bit-count _)(and _ 1)(lsh _ 1)(- 1 _))#r:1/m))(>-> x (- _ 2) (sqr _))(>-> s (lsh _ 1) (or _ 0))))))))))

… OTHER INVERSE FUNCTIONS
Because (∘ sqr double sqrt reciprocal add1 sqr)=(∘ [κ / 4 _] add1 sqr)
(def-syntax *eps* 1e-10)(Define [dyadic-acos x :flo] (dyadic-find ((∘ sqr double) x) 0 *eps*))(Define [dyadic-asin x :flo] (dyadic-find ((∘ sqr double) x) 1 *eps*))(Define [dyadic-atan x :flo] (dyadic-find ((∘ [κ / 4 _] add1 sqr) x) 0 *eps*))(Define [dyadic-acot x :flo] (dyadic-find ((∘ [κ / 4 _] add1 sqr) x) 1 *eps*))(Define (ACos x :flo) (if [< x 0.](- #r:1/2 (/ (dyadic-acos (neg x)) 4))(/ (dyadic-acos x 4)))
(Define (ASin x :flo) (if [< x 0.](neg (/ (dyadic-asin (neg x)) 4))(/ (dyadic-asin x 4)))
(Define (ATan x :flo) (if [< x 0.](neg (/ (dyadic-atan (neg x)) 4))(/ (dyadic-atan x 4)))
(Define (ACot x :flo) (if [< x 0.](- #r:1/2 (/ (dyadic-acot (neg x)) 4))(/ (dyadic-acot x 4)))

EVALUATIONANSWERS
1. Do we need special implementations of sqrt, sin, cos etc.?
Need sqrt, forward funcs ≡ continued radicals
2. Algorithms for inverse trigonometric functions derivable from forward ones?
Direct inverse of forward steps, in reverse
3. Avoid Taylor/Maclaurin series and pre-computed tables a-la CORDIC?
Tables or expensive series not required
4. Can we enumerate all circular values in the interval [0, π2 ]?
Yes, similar to Cantor’s enumeration of Q

CONCLUSIONSThis application of functional programming is beneﬁcial in showing
what can materialise when thinking algorithmically, in addition tothinking mathematically. In particular, a derivation of an algorithm forthe inverse can be straightforward given the algebraic representation.

• Computing cispi (aka sincospi) now avoids duplicate calculation• forward circular trigonometry is not competitive due to overhead of software-based sqrt
▪ this will change once sqrt is implemented in hardware• backward circular trigonometry is competitive due to only sqr being used
▪ only uses bit-shifts, so mapping to hardware is trivial▪ advantageous when #iterations is < 9

APPENDIX
1. although my regular job is still NXP, this particular work has been performed on a private basis↩
2. we thank AZ for this observation ↩

8 . 1