fix docs & improve test

docs/src/api.md

@@ -84,7 +84,6 @@ wigner_init_float
 PFRational
 gcd(::PFRational, ::PFRational)
 lcm(::PFRational, ::PFRational)
-sgcd(::PFRational, ::PFRational)
 wigner_init_pf
 pf_binomial
 eCG

src/PFRational/PFRational.jl

@@ -6,7 +6,7 @@ A positive rational number in prime factorization form.
 ```math
 q = \prod_{i=1}^{n} p_i^{e_i}
 ```
-where `p_i` is the `i`-th prime number and `e_i` is the exponent of `p_i`, `e_i` can be negative.
+where ``p_i`` is the ``i``-th prime number and ``e_i`` can be negative.
 """
 struct PFRational{T<:Integer} <: Real
     e::Vector{T}
@@ -24,14 +24,6 @@ pf_alloc!(x::PFRational{T}, n::Integer) where {T<:Integer} = begin
     resize!(x.e, n)
     fill!(x.e, T(0))
 end
-function pf_fit!(x::PFRational)
-    n = length(x.e)
-    while n > 0 && x.e[n] == 0
-        n -= 1
-    end
-    resize!(x.e, n)
-    return x
-end
 
 Base.inv(x::PFRational) = PFRational(-x.e)
 
@@ -56,7 +48,6 @@ end
 end
 
 ^(x::PFRational, n::Integer) = PFRational(x.e .* n)
-square(x::PFRational) = PFRational(x.e .* 2)
 square!(x::PFRational) = begin
     for i in eachindex(x.e)
         x.e[i] *= 2
@@ -91,30 +82,6 @@ Base.lcm(x::PFRational, y::PFRational) = begin
     return PFRational{T}(e)
 end
 
-"""
-    scgd(x::PFRational, y::PFRational)
-here we define `scgd` of rational, which means gcd of both numerator and denominator
-"""
-scgd(x::PFRational, y::PFRational) = begin
-    T = promote_type(eltype(x.e), eltype(y.e))
-    e = zeros(T, max(length(x.e), length(y.e)))
-    copyto!(e, x.e)
-    for i in eachindex(y.e)
-        e_min, e_max = minmax(e[i], y.e[i])
-        if e_min > 0
-            e[i] = e_min
-        elseif e_max < 0
-            e[i] = e_max
-        else
-            e[i] = 0
-        end
-    end
-    for i in length(y.e)+1:length(x.e)
-        e[i] = 0
-    end
-    return PFRational{T}(e)
-end
-
 # assume length(y) <= length(x)
 _copy!(x::PFRational{T}, y::PFRational{T}) where {T<:Integer} = begin
     for i in eachindex(y.e)
@@ -159,24 +126,6 @@ _lcm!(x::PFRational{T}, y::PFRational{T}) where {T<:Integer} = begin
     end
 end
 
-# assume length(y) <= length(x)
-_sgcd!(x::PFRational{T}, y::PFRational{T}) where {T<:Integer} = begin
-    for i in eachindex(y.e)
-        e_min, e_max = minmax(x.e[i], y.e[i])
-        if e_min > 0
-            x.e[i] = e_min
-        elseif e_max < 0
-            x.e[i] = e_max
-        else
-            x.e[i] = 0
-        end
-    end
-    for i in length(y.e)+1:length(x.e)
-        x.e[i] = 0
-    end
-end
-
-
 # positive integers in prime factorization form
 _pf_integers = PFRational{Int16}[one(PFRational{Int16})]
 # binomial coefficients in prime factorization form

src/PFRational/primetable.jl

@@ -7,7 +7,6 @@ mutable struct PrimePowItem
         new(size, offset, ulongdata, bigdata)
     end
 end
-PrimePowItem() = PrimePowItem(0, 0, Culong[], BigInt[])
 
 function make_pow_item(p::Culong, k::Int)
     offset = floor(Int, log(big(p), big(typemax(Culong))))
@@ -23,11 +22,9 @@ function make_pow_item(p::Culong, k::Int)
     else
         bigdata = BigInt[]
     end
-    PrimePowItem(k, offset, ulongdata, bigdata)
+    PrimePowItem(max(k, offset), offset, ulongdata, bigdata)
 end
 
-make_pow_item(p::Culong) = make_pow_item(p, 0)
-
 # first n primes
 _prime_table = Culong[2, 3, 5, 7]
 # first n primes and their powers

test/runtests.jl

@@ -36,6 +36,7 @@ include("test_pf.jl")
 @testset "test eCG0" begin test_eCG0(1:3) end
 @testset "test e3j" begin test_e3j(1:3) end
 @testset "test e6j" begin test_e6j(1:3) end
+@testset "test big binomial" begin test_big_binomial() end
 
 try
     gsl3j(1, 1, 1, 0, 0, 0)

test/test_pf.jl

@@ -84,3 +84,17 @@ function test_e6j(test_range::AbstractArray)
         @test convert(Float64, float(d)) == ef
     end end end end end end
 end
+
+
+function test_big_binomial()
+    wigner_init_pf(100, "nmax", 0)
+    wigner_init_float(100, "nmax", 0)
+    n = 100
+    for k in 0:n
+        x = pf_binomial(n, k)
+        y = binomial(big(n), big(k))
+        @test numerator(x) == y
+        f = fbinomial(n, k)
+        @test convert(Float64, numerator(x)) ≈ f
+    end
+end