Is my reasoning correct?
I can see at least five major problems in your reasoning.
Do most genes have two alleles?
The claim that most genes have two alleles
is wrong. Actually, the entire top answer in the linked quora post is wrong. The person answering confounds simplified cases used to teach intro genetics class with what nature really looks like.
As about one every 100 nucleotides is a SNP (bionumbers) and given that a gene is on average about 30k sites (bionumbers), the average number of alleles is probably closer to 300 although that would be counting synonymous mutations which is a bit unfair. About 20% of substitutions are synonymous (bionumbers), which would lead to an estimate of 240 alleles per gene on average. Even if you consider that 90% of these mutations have no effect on the phenotype (which is maybe an overestimation, bionumbers), you're left with 24 alleles per gene. Of course, this is a very vague estimate but it is, I think, much better than 2.
Non-coding sequences
Also, the majority of phenotypic distinction is caused by non-coding (non - gene) sequences such as regulatory sequences. You cannot ignore them in the calculation.
According to this pie chart (not sure if we can trust this chart though), 5% of the genome is made of regulatory sequences (against 1.5% of the genome that is coding).
One genotype ≠ one phenotype
If you have in mind that one genotype = one phenotype, then it would be very wrong too. There is environmental variance, epigenetic variance and other sources of variance (e.g. developmental noise) that underly phenotypic variance (see this post to understand the sources of phenotypic variance in populations).
indels and CNV
You are only considering substitutions. There are plenty of indels and ohter CNV too.
The existence of this variation does not only complicate the question, it actually renders the question a bit undefined because if you allow the genome to take any size, then number of possibilities is necessarily infinite.
Diploidy
Your calculations assume individuals are haploids. But humans are diplontic.
Gene are not atomic
Genes are not units that cannot be split. Recombination do happen in genes and it would be wrong to take the number of variants at genes (or at any other type of sequences) as basis for your calculation
Attempt at a calculation
So, ignoring many of the problems, here is my attempt at a calculation.
There are about $1.5 \cdot 10^7$ SNPs in the human genome (bionumbers). Based on this pie chart (not sure if we can trust this chart though), I will assume that about 6.5% of them are in coding or regulatory regions and I will assume (based on that paper) that 90% of them have no fitness effect and by extension I will assume it also mean no phenotypic effect (although the actually percentage should been increased). This will result in $1.5 \cdot 10^7 \cdot 0.065 \cdot 0.1 = 97500$ SNPs of interest. It results into $2^{97500} ≈ 10^{32500}$ possible haplotypes of interest. The number of possible genotypes of interest is therefore $\left(10^{32500}\right)^2 = 10^{65000}$. I also assumed absence of epistatic interactions otherwise the number would become so fantastically closer to infinity!
Interest in the calculation
I fail to find much interest in making such calculation because it entirely depends upon the model we want to consider. There is no offense against the OP's question here (and in fact the OP shows that he understands the importance of the considered model), I just want to make sure readers understand that these calculations are based on a very arbitrary model and one cannot make any sense of the given result without understanding the model.
The calculations from the Quora post that OP linked assumes that every SNP could be polymorphic. As long as the intention behind the calculation is not clearly stated, this is a perfectly reasonable model. Here, in my attempt, I assume that the number of SNPs is fixed, but any combination of it is possible. This can or cannot be a fine assumptions. Also, the estimates of the number of SNPs are of course, estimates of the number of SNPs that are at a frequency that is not too low to be detected! So, if we were, for some reason, to consider all the SNPs (by assuming some drift-mutation equilibrium), then the result would differ quite a bit. But again, it all depends upon the specific model one wants to consider.