/* ;; Copyright (C) 2021 Yasuaki Honda ;; modular-j.mac ;; This file is part of ModularJ package. ;; ;; This program is free software; you can redistribute it and/or ;; modify it under the terms of the GNU General Public License as ;; published by the Free Software Foundation; either version 2 of ;; the License, or (at your option) any later version. ;; ;; This program is distributed in the hope that it will be ;; useful, but WITHOUT ANY WARRANTY; without even the implied ;; warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR ;; PURPOSE. See the GNU General Public License for more details. ;; ;; You should have received a copy of the GNU General Public ;; License along with this program; if not, see ;; . */ /* General Description Klein's elliptic modular j-invariant, or j function in short, has a notable property that the special value for the imaginary quadratic irrational numbers can be computed exactly. The following is the algorithm implemented in Maxima in this notebook. Let's assume 𝑡 is an imaginary quadratic irrational number, the argument for the j function. Compute the determinant D of t. List up all the reduced imaginary quadratic irrational numbers x0,,,xn whose determinant is D. Note that the number of such irrationals n is called the class number of the determinant D. Compute the aproximate class equation (x-j(x0))...(x-j(xn)) where j(xi) is computed using q-expansion of j function. Find the exact class equation by rounding up the coefficients of the aproximate class equation and solve it using solve() algebraically. Find the correspondance between aproximate j(xi) and the exact j(xi) to know the exact value of j(xi). API Manual InFD(p) p is a complex number. InFD(p) checks if p is in the fundamental domain of j(t). FindInFD(p) p is a complex number. FindInFD(p) returns a complex number in the fundamental domain and equivalent to p. It shows the necessary transformations, too. float_modular_j(t) t is a complex float number. float_modular_j(t) numerically computes the value of j(t). FindDetIrr(t) t is an imaginary quadratic irrational number. FindDetIrr(t) computes the determinant of t. FindReducedQuadIrrD(D) D is an integer. FindReducedQuadIrrD(D) computes the list of all the reduced imaginary quadratic irrational numbers whose determinant is D. FindClassNumberD(D) D is an integer. FindClassNumberD(D) computes the class number of the determinant D. FindClassPolyD(var,D) var is a variable, D is an integer. FindClassPolyD(D) computes the Hilbert class polynomial of variable var and the determinant equal to D. exact_modular_j(t) t is an imaginary quadratic irrational number of the form (a1+a2*sqrt(a3))/a4 where a1,,,a4 are integers, a2,a4 positive and a4 negative. exact_modular_j(t) computes the exact value (which is an algebraic number) of j(t) using the above algorithm. */ debug_modular_j:false; InFD(p):=block([absp,realp],absp:abs(p),realp:realpart(p), if imagpart(p) <= 0 then return (false) else return ((absp>1 and realp>=-1/2 and realp<1/2) or (absp=1 and -1/2<=realp and realp<=0)))$ Ta(p):=block([realp,ans],realp:realpart(p),ans:realp-floor(realp+1/2), if ans=1/2 then ans:ans-1, return(realp-ans))$ FindInFD(p):= ev(block([a],a:Ta(p), if a#0 then print("Ta:",a), p:p-a, if abs(p)<1 then (print("S"),return(FindInFD(-1/p))), if abs(p)>1 then return(p), if abs(p)=1 then if realpart(p)<=0 then return(p) else return(-1/p)),ratsimp,algebraic:true)$ j_coeff:[196884, 21493760, 864299970, 20245856256, 333202640600, 4252023300096, 44656994071935, 401490886656000, 3176440229784420, 22567393309593600, 146211911499519294, 874313719685775360, 4872010111798142520, 25497827389410525184, 126142916465781843075]$ float_modular_j(t):=ev(1/q+744+sum(j_coeff[i]*q^i,i,1,length(j_coeff)),q:rectform(exp(2*%pi*%i*t)),eval,numer)$ bfloat_modular_j(t):=ev(1/q+744+sum(j_coeff[i]*q^i,i,1,length(j_coeff)),q:rectform(exp(2*%pi*%i*t)),fpprec:32,eval,bfloat)$ Checkabc(a,b,c):=(((-a=3. if e is mulformed, return false. Otherwise, [exp1^(exp2-1/2),exp1] is returned. */ if (not(atom(e)) and op(e)="^" and integerp(args(e)[1]) and not(atom(args(e)[2])) and op(args(e)[2])="/" and args(args(e)[2])[2]=2 ) then [args(e)[1]^(args(e)[2]-1/2),args(e)[1]] else false$ irr_number(t):=block([a1,a2,a3,a4,tmp], if imagpart(t)=0 then error("irr_number: argument is not an imaginary quadratic irrational number ",t), a4:ratdenom(t), t:t*a4, a1:realpart(t), tmp:imagpart(t), /* tmp pattern: int, sqrt(int), int^(int/2), int*sqrt(int),int*int^(int/2) */ if integerp(tmp) then [a2,a3]:[tmp,1] elseif op(tmp)='sqrt then [a2,a3]:[1,args(tmp)[1]] elseif op(tmp)="^" then block([ret:irr_number_parse_exp(tmp)], if ret=false then error("irr_number: argument is not well formated as an imaginary quadratic irrational number ",t) else [a2,a3]:ret) elseif op(tmp)="*" and not(atom(args(tmp)[2])) and op(args(tmp)[2])='sqrt then [a2,a3]:[args(tmp)[1],args(args(tmp)[2])[1]] elseif op(tmp)="*" and not(atom(args(tmp)[2])) and op(args(tmp)[2])="^" then block([ret:irr_number_parse_exp(args(tmp)[2])], if ret=false then error("irr_number: argument is not well formated as an imaginary quadratic irrational number ",t) else [a2,a3]:[args(tmp)[1]*ret[1],ret[2]]) else error("irr_number: argument is not well formated as an imaginary quadratic irrational number ",t), return([a1,a2,-a3,a4]))$ FindDetIrr(t):=block([tparts], tparts:irr_number(t), D:ev(4*a2^2*a3*a4^2/gcd(gcd(a4^2,2*a1*a4),a1^2-a2^2*a3)^2,a1:tparts[1],a2:tparts[2],a3:tparts[3],a4:tparts[4]), return(D))$ exact_modular_j(t):=block([D,fpjt,cpoly,sol,numsol,x], D:FindDetIrr(t), if debug_modular_j=true then print("D=",D), cpoly:FindClassPolyD(x,D), if cpoly=0 then error("exact_modular_j(): Cannot compute class polynomial for", t), if debug_modular_j=true then print("cpoly=",cpoly), if hipow(cpoly,x)>=5 then error("exact_modular_j(): Degree of class polynomial is higher than 5, ", cpoly), sol:solve(cpoly), numsol:map(lambda([xx],rhs(float(xx))),sol), fpjt:float_modular_j(t), block([mindex,minval],mindex:1,minval:abs(numsol[1]-fpjt), for i:1 thru length(numsol) do if abs(numsol[i]-fpjt)